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 <config.h>
  36 #include <stdio.h>
  37 #include "ntp_stdlib.h"
  38 #include "ntp_types.h"
  39 #include "ntp_fp.h"
  40
  41 #define LOW_MASK  (u_int32)((1<<(FRACTION_PREC/2))-1)
  42 #define HIGH_MASK (u_int32)(LOW_MASK << (FRACTION_PREC/2))
  43
  44 /*
  45  * for those who worry about overflows (possibly triggered by static analysis tools):
  46  *
  47  * Largest value of a 2^n bit number is 2^n-1.
  48  * Thus the result is: (2^n-1)*(2^n-1) = 2^2n - 2^n - 2^n + 1 < 2^2n
  49  * Here overflow can not happen for 2 reasons:
  50  * 1) the code actually multiplies the absolute values of two signed
  51  *    64bit quantities.thus effectively multiplying 2 63bit quantities.
  52  * 2) Carry propagation is from low to high, building principle is
  53  *    addition, so no storage for the 2^2n term from above is needed.
  54  */
  55
  56 void
  57 mfp_mul(
  58         int32   *o_i,
  59         u_int32 *o_f,
  60         int32    a_i,
  61         u_int32  a_f,
  62         int32    b_i,
  63         u_int32  b_f
  64         )
  65 {
  66   int32 i, j;
  67   u_int32  f;
  68   u_long a[4];                  /* operand a */
  69   u_long b[4];                  /* operand b */
  70   u_long c[5];                  /* result c - 5 items for performance - see below */
  71   u_long carry;
  72
  73   int neg = 0;
  74
  75   if (a_i < 0)                  /* examine sign situation */
  76     {
  77       neg = 1;
  78       M_NEG(a_i, a_f);
  79     }
  80
  81   if (b_i < 0)                  /* examine sign situation */
  82     {
  83       neg = !neg;
  84       M_NEG(b_i, b_f);
  85     }
  86
  87   a[0] = a_f & LOW_MASK;        /* prepare a operand */
  88   a[1] = (a_f & HIGH_MASK) >> (FRACTION_PREC/2);
  89   a[2] = a_i & LOW_MASK;
  90   a[3] = (a_i & HIGH_MASK) >> (FRACTION_PREC/2);
  91
  92   b[0] = b_f & LOW_MASK;        /* prepare b operand */
  93   b[1] = (b_f & HIGH_MASK) >> (FRACTION_PREC/2);
  94   b[2] = b_i & LOW_MASK;
  95   b[3] = (b_i & HIGH_MASK) >> (FRACTION_PREC/2);
  96
  97   c[0] = c[1] = c[2] = c[3] = c[4] = 0;
  98
  99   for (i = 0; i < 4; i++)       /* we do assume 32 * 32 = 64 bit multiplication */
 100     for (j = 0; j < 4; j++)
 101       {
 102         u_long result_low, result_high;
 103         int low_index = (i+j)/2;      /* formal [0..3]  - index for low long word */
 104         int mid_index = 1+low_index;  /* formal [1..4]! - index for high long word
 105                                          will generate unecessary add of 0 to c[4]
 106                                          but save 15 'if (result_high) expressions' */
 107         int high_index = 1+mid_index; /* formal [2..5]! - index for high word overflow
 108                                          - only assigned on overflow (limits range to 2..3) */
 109
 110         result_low = (u_long)a[i] * (u_long)b[j];       /* partial product */
 111
 112         if ((i+j) & 1)          /* splits across two result registers */
 113           {
 114             result_high   = result_low >> (FRACTION_PREC/2);
 115             result_low  <<= FRACTION_PREC/2;
 116             carry         = (unsigned)1<<(FRACTION_PREC/2);
 117           }
 118         else
 119           {                     /* stays in a result register - except for overflows */
 120             result_high = 0;
 121             carry       = 1;
 122           }
 123
 124         if (((c[low_index] >> 1) + (result_low >> 1) + ((c[low_index] & result_low & carry) != 0)) &
 125             (u_int32)((unsigned)1<<(FRACTION_PREC - 1))) {
 126           result_high++;        /* propagate overflows */
 127         }
 128
 129         c[low_index]   += result_low; /* add up partial products */
 130
 131         if (((c[mid_index] >> 1) + (result_high >> 1) + ((c[mid_index] & result_high & 1) != 0)) &
 132             (u_int32)((unsigned)1<<(FRACTION_PREC - 1))) {
 133           c[high_index]++;              /* propagate overflows of high word sum */
 134         }
 135
 136         c[mid_index] += result_high;  /* will add a 0 to c[4] once but saves 15 if conditions */
 137       }
 138
 139 #ifdef DEBUG
 140   if (debug > 6)
 141     printf("mfp_mul: 0x%04lx%04lx%04lx%04lx * 0x%04lx%04lx%04lx%04lx = 0x%08lx%08lx%08lx%08lx\n",
 142          a[3], a[2], a[1], a[0], b[3], b[2], b[1], b[0], c[3], c[2], c[1], c[0]);
 143 #endif
 144
 145   if (c[3])                     /* overflow */
 146     {
 147       i = ((unsigned)1 << (FRACTION_PREC-1)) - 1;
 148       f = ~(unsigned)0;
 149     }
 150   else
 151     {                           /* take produkt - discarding extra precision */
 152       i = c[2];
 153       f = c[1];
 154     }
 155
 156   if (neg)                      /* recover sign */
 157     {
 158       M_NEG(i, f);
 159     }
 160
 161   *o_i = i;
 162   *o_f = f;
 163
 164 #ifdef DEBUG
 165   if (debug > 6)
 166     printf("mfp_mul: %s * %s => %s\n",
 167            mfptoa((u_long)a_i, a_f, 6),
 168            mfptoa((u_long)b_i, b_f, 6),
 169            mfptoa((u_long)i, f, 6));
 170 #endif
 171 }
 172
 173 /*
 174  * History:
 175  *
 176  * mfp_mul.c,v
 177  * Revision 4.9  2005/07/17 20:34:40  kardel
 178  * correct carry propagation implementation
 179  *
 180  * Revision 4.8  2005/07/12 16:17:26  kardel
 181  * add explanation why we do not write into c[4]
 182  *
 183  * Revision 4.7  2005/04/16 17:32:10  kardel
 184  * update copyright
 185  *
 186  * Revision 4.6  2004/11/14 15:29:41  kardel
 187  * support PPSAPI, upgrade Copyright to Berkeley style
 188  *
 189  * Revision 4.3  1999/02/21 12:17:37  kardel
 190  * 4.91f reconcilation
 191  *
 192  * Revision 4.2  1998/12/20 23:45:28  kardel
 193  * fix types and warnings
 194  *
 195  * Revision 4.1  1998/05/24 07:59:57  kardel
 196  * conditional debug support
 197  *
 198  * Revision 4.0  1998/04/10 19:46:38  kardel
 199  * Start 4.0 release version numbering
 200  *
 201  * Revision 1.1  1998/04/10 19:27:47  kardel
 202  * initial NTP VERSION 4 integration of PARSE with GPS166 binary support
 203  *
 204  * Revision 1.1  1997/10/06 21:05:46  kardel
 205  * new parse structure
 206  *
 207  */