contrib/ntp/libntp/mfp_mul.c

   1 /*
   2  * /src/NTP/ntp-4/libntp/mfp_mul.c,v 4.3 1999/02/21 12:17:37 kardel RELEASE_19990228_A
   3  *
   4  * $Created: Sat Aug 16 20:35:08 1997 $
   5  *
   6  * Copyright (C) 1997, 1998 by Frank Kardel
   7  */
   8 #include <stdio.h>
   9 #include "ntp_stdlib.h"
  10 #include "ntp_types.h"
  11 #include "ntp_fp.h"
  12
  13 #define LOW_MASK  (u_int32)((1<<(FRACTION_PREC/2))-1)
  14 #define HIGH_MASK (u_int32)(LOW_MASK << (FRACTION_PREC/2))
  15
  16 void
  17 mfp_mul(
  18         int32   *o_i,
  19         u_int32 *o_f,
  20         int32    a_i,
  21         u_int32  a_f,
  22         int32    b_i,
  23         u_int32  b_f
  24         )
  25 {
  26   int32 i, j;
  27   u_int32  f;
  28   u_long a[4];                  /* operand a */
  29   u_long b[4];                  /* operand b */
  30   u_long  c[4];                 /* result c */
  31
  32   int neg = 0;
  33
  34   if (a_i < 0)                  /* examine sign situation */
  35     {
  36       neg = 1;
  37       M_NEG(a_i, a_f);
  38     }
  39
  40   if (b_i < 0)                  /* examine sign situation */
  41     {
  42       neg = !neg;
  43       M_NEG(b_i, b_f);
  44     }
  45
  46   a[0] = a_f & LOW_MASK;        /* prepare a operand */
  47   a[1] = (a_f & HIGH_MASK) >> (FRACTION_PREC/2);
  48   a[2] = a_i & LOW_MASK;
  49   a[3] = (a_i & HIGH_MASK) >> (FRACTION_PREC/2);
  50
  51   b[0] = b_f & LOW_MASK;        /* prepare b operand */
  52   b[1] = (b_f & HIGH_MASK) >> (FRACTION_PREC/2);
  53   b[2] = b_i & LOW_MASK;
  54   b[3] = (b_i & HIGH_MASK) >> (FRACTION_PREC/2);
  55
  56   c[0] = c[1] = c[2] = c[3] = 0;
  57
  58   for (i = 0; i < 4; i++)       /* we do assume 32 * 32 = 64 bit multiplication */
  59     for (j = 0; j < 4; j++)
  60       {
  61         u_long result_low, result_high;
  62
  63         result_low = (u_long)a[i] * (u_long)b[j];       /* partial product */
  64
  65         if ((i+j) & 1)          /* splits across two result registers */
  66           {
  67             result_high   = result_low >> (FRACTION_PREC/2);
  68             result_low  <<= FRACTION_PREC/2;
  69           }
  70         else
  71           {                     /* stays in a result register - except for overflows */
  72             result_high = 0;
  73           }
  74
  75         if (((c[(i+j)/2] >> 1) + (result_low >> 1)) & (u_int32)((unsigned)1<<(FRACTION_PREC - 1)))
  76           result_high++;        /* propagate overflows */
  77
  78         c[(i+j)/2]   += result_low; /* add up partial products */
  79
  80         if (((c[(i+j+1)/2] >> 1) + (result_high >> 1)) & (u_int32)((unsigned)1<<(FRACTION_PREC - 1)))
  81           c[1+(i+j)/2]++;       /* propagate overflows */
  82
  83         c[(i+j+1)/2] += result_high;
  84       }
  85
  86 #ifdef DEBUG
  87   if (debug > 6)
  88     printf("mfp_mul: 0x%04lx%04lx%04lx%04lx * 0x%04lx%04lx%04lx%04lx = 0x%08lx%08lx%08lx%08lx\n",
  89          a[3], a[2], a[1], a[0], b[3], b[2], b[1], b[0], c[3], c[2], c[1], c[0]);
  90 #endif
  91
  92   if (c[3])                     /* overflow */
  93     {
  94       i = ((unsigned)1 << (FRACTION_PREC-1)) - 1;
  95       f = ~(unsigned)0;
  96     }
  97   else
  98     {                           /* take produkt - discarding extra precision */
  99       i = c[2];
 100       f = c[1];
 101     }
 102
 103   if (neg)                      /* recover sign */
 104     {
 105       M_NEG(i, f);
 106     }
 107
 108   *o_i = i;
 109   *o_f = f;
 110
 111 #ifdef DEBUG
 112   if (debug > 6)
 113     printf("mfp_mul: %s * %s => %s\n",
 114            mfptoa((u_long)a_i, a_f, 6),
 115            mfptoa((u_long)b_i, b_f, 6),
 116            mfptoa((u_long)i, f, 6));
 117 #endif
 118 }
 119
 120 /*
 121  * mfp_mul.c,v
 122  * Revision 4.3  1999/02/21 12:17:37  kardel
 123  * 4.91f reconcilation
 124  *
 125  * Revision 4.2  1998/12/20 23:45:28  kardel
 126  * fix types and warnings
 127  *
 128  * Revision 4.1  1998/05/24 07:59:57  kardel
 129  * conditional debug support
 130  *
 131  * Revision 4.0  1998/04/10 19:46:38  kardel
 132  * Start 4.0 release version numbering
 133  *
 134  * Revision 1.1  1998/04/10 19:27:47  kardel
 135  * initial NTP VERSION 4 integration of PARSE with GPS166 binary support
 136  *
 137  * Revision 1.1  1997/10/06 21:05:46  kardel
 138  * new parse structure
 139  *
 140  */