contrib/ntp/libparse/mfp_mul.c

   1 /*
   2  * /src/NTP/ntp4-dev/libparse/mfp_mul.c,v 4.9 2005/07/17 20:34:40 kardel RELEASE_20050717_A
   3  *
   4  * mfp_mul.c,v 4.9 2005/07/17 20:34:40 kardel RELEASE_20050717_A
   5  *
   6  * $Created: Sat Aug 16 20:35:08 1997 $
   7  *
   8  * Copyright (c) 1997-2005 by Frank Kardel <kardel <AT> ntp.org>
   9  *
  10  * Redistribution and use in source and binary forms, with or without
  11  * modification, are permitted provided that the following conditions
  12  * are met:
  13  * 1. Redistributions of source code must retain the above copyright
  14  *    notice, this list of conditions and the following disclaimer.
  15  * 2. Redistributions in binary form must reproduce the above copyright
  16  *    notice, this list of conditions and the following disclaimer in the
  17  *    documentation and/or other materials provided with the distribution.
  18  * 3. Neither the name of the author nor the names of its contributors
  19  *    may be used to endorse or promote products derived from this software
  20  *    without specific prior written permission.
  21  *
  22  * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
  23  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  24  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  25  * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
  26  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  27  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  28  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  29  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  30  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  31  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  32  * SUCH DAMAGE.
  33  *
  34  */
  35 #include <stdio.h>
  36 #include "ntp_stdlib.h"
  37 #include "ntp_types.h"
  38 #include "ntp_fp.h"
  39
  40 #define LOW_MASK  (u_int32)((1<<(FRACTION_PREC/2))-1)
  41 #define HIGH_MASK (u_int32)(LOW_MASK << (FRACTION_PREC/2))
  42
  43 /*
  44  * for those who worry about overflows (possibly triggered by static analysis tools):
  45  *
  46  * Largest value of a 2^n bit number is 2^n-1.
  47  * Thus the result is: (2^n-1)*(2^n-1) = 2^2n - 2^n - 2^n + 1 < 2^2n
  48  * Here overflow can not happen for 2 reasons:
  49  * 1) the code actually multiplies the absolute values of two signed
  50  *    64bit quantities.thus effectively multiplying 2 63bit quantities.
  51  * 2) Carry propagation is from low to high, building principle is
  52  *    addition, so no storage for the 2^2n term from above is needed.
  53  */
  54
  55 void
  56 mfp_mul(
  57         int32   *o_i,
  58         u_int32 *o_f,
  59         int32    a_i,
  60         u_int32  a_f,
  61         int32    b_i,
  62         u_int32  b_f
  63         )
  64 {
  65   int32 i, j;
  66   u_int32  f;
  67   u_long a[4];                  /* operand a */
  68   u_long b[4];                  /* operand b */
  69   u_long c[5];                  /* result c - 5 items for performance - see below */
  70   u_long carry;
  71
  72   int neg = 0;
  73
  74   if (a_i < 0)                  /* examine sign situation */
  75     {
  76       neg = 1;
  77       M_NEG(a_i, a_f);
  78     }
  79
  80   if (b_i < 0)                  /* examine sign situation */
  81     {
  82       neg = !neg;
  83       M_NEG(b_i, b_f);
  84     }
  85
  86   a[0] = a_f & LOW_MASK;        /* prepare a operand */
  87   a[1] = (a_f & HIGH_MASK) >> (FRACTION_PREC/2);
  88   a[2] = a_i & LOW_MASK;
  89   a[3] = (a_i & HIGH_MASK) >> (FRACTION_PREC/2);
  90
  91   b[0] = b_f & LOW_MASK;        /* prepare b operand */
  92   b[1] = (b_f & HIGH_MASK) >> (FRACTION_PREC/2);
  93   b[2] = b_i & LOW_MASK;
  94   b[3] = (b_i & HIGH_MASK) >> (FRACTION_PREC/2);
  95
  96   c[0] = c[1] = c[2] = c[3] = c[4] = 0;
  97
  98   for (i = 0; i < 4; i++)       /* we do assume 32 * 32 = 64 bit multiplication */
  99     for (j = 0; j < 4; j++)
 100       {
 101         u_long result_low, result_high;
 102         int low_index = (i+j)/2;      /* formal [0..3]  - index for low long word */
 103         int mid_index = 1+low_index;  /* formal [1..4]! - index for high long word
 104                                          will generate unecessary add of 0 to c[4]
 105                                          but save 15 'if (result_high) expressions' */
 106         int high_index = 1+mid_index; /* formal [2..5]! - index for high word overflow
 107                                          - only assigned on overflow (limits range to 2..3) */
 108
 109         result_low = (u_long)a[i] * (u_long)b[j];       /* partial product */
 110
 111         if ((i+j) & 1)          /* splits across two result registers */
 112           {
 113             result_high   = result_low >> (FRACTION_PREC/2);
 114             result_low  <<= FRACTION_PREC/2;
 115             carry         = (unsigned)1<<(FRACTION_PREC/2);
 116           }
 117         else
 118           {                     /* stays in a result register - except for overflows */
 119             result_high = 0;
 120             carry       = 1;
 121           }
 122
 123         if (((c[low_index] >> 1) + (result_low >> 1) + ((c[low_index] & result_low & carry) != 0)) &
 124             (u_int32)((unsigned)1<<(FRACTION_PREC - 1))) {
 125           result_high++;        /* propagate overflows */
 126         }
 127
 128         c[low_index]   += result_low; /* add up partial products */
 129
 130         if (((c[mid_index] >> 1) + (result_high >> 1) + ((c[mid_index] & result_high & 1) != 0)) &
 131             (u_int32)((unsigned)1<<(FRACTION_PREC - 1))) {
 132           c[high_index]++;              /* propagate overflows of high word sum */
 133         }
 134
 135         c[mid_index] += result_high;  /* will add a 0 to c[4] once but saves 15 if conditions */
 136       }
 137
 138 #ifdef DEBUG
 139   if (debug > 6)
 140     printf("mfp_mul: 0x%04lx%04lx%04lx%04lx * 0x%04lx%04lx%04lx%04lx = 0x%08lx%08lx%08lx%08lx\n",
 141          a[3], a[2], a[1], a[0], b[3], b[2], b[1], b[0], c[3], c[2], c[1], c[0]);
 142 #endif
 143
 144   if (c[3])                     /* overflow */
 145     {
 146       i = ((unsigned)1 << (FRACTION_PREC-1)) - 1;
 147       f = ~(unsigned)0;
 148     }
 149   else
 150     {                           /* take produkt - discarding extra precision */
 151       i = c[2];
 152       f = c[1];
 153     }
 154
 155   if (neg)                      /* recover sign */
 156     {
 157       M_NEG(i, f);
 158     }
 159
 160   *o_i = i;
 161   *o_f = f;
 162
 163 #ifdef DEBUG
 164   if (debug > 6)
 165     printf("mfp_mul: %s * %s => %s\n",
 166            mfptoa((u_long)a_i, a_f, 6),
 167            mfptoa((u_long)b_i, b_f, 6),
 168            mfptoa((u_long)i, f, 6));
 169 #endif
 170 }
 171
 172 /*
 173  * History:
 174  *
 175  * mfp_mul.c,v
 176  * Revision 4.9  2005/07/17 20:34:40  kardel
 177  * correct carry propagation implementation
 178  *
 179  * Revision 4.8  2005/07/12 16:17:26  kardel
 180  * add explanation why we do not write into c[4]
 181  *
 182  * Revision 4.7  2005/04/16 17:32:10  kardel
 183  * update copyright
 184  *
 185  * Revision 4.6  2004/11/14 15:29:41  kardel
 186  * support PPSAPI, upgrade Copyright to Berkeley style
 187  *
 188  * Revision 4.3  1999/02/21 12:17:37  kardel
 189  * 4.91f reconcilation
 190  *
 191  * Revision 4.2  1998/12/20 23:45:28  kardel
 192  * fix types and warnings
 193  *
 194  * Revision 4.1  1998/05/24 07:59:57  kardel
 195  * conditional debug support
 196  *
 197  * Revision 4.0  1998/04/10 19:46:38  kardel
 198  * Start 4.0 release version numbering
 199  *
 200  * Revision 1.1  1998/04/10 19:27:47  kardel
 201  * initial NTP VERSION 4 integration of PARSE with GPS166 binary support
 202  *
 203  * Revision 1.1  1997/10/06 21:05:46  kardel
 204  * new parse structure
 205  *
 206  */