dist/ntp/libparse/mfp_mul.c

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