2 * /src/NTP/ntp-4/libntp/mfp_mul.c,v 4.3 1999/02/21 12:17:37 kardel RELEASE_19990228_A
4 * $Created: Sat Aug 16 20:35:08 1997 $
6 * Copyright (C) 1997, 1998 by Frank Kardel
9 #include "ntp_stdlib.h"
10 #include "ntp_types.h"
13 #define LOW_MASK (u_int32)((1<<(FRACTION_PREC/2))-1)
14 #define HIGH_MASK (u_int32)(LOW_MASK << (FRACTION_PREC/2))
28 u_long a[4]; /* operand a */
29 u_long b[4]; /* operand b */
30 u_long c[4]; /* result c */
34 if (a_i < 0) /* examine sign situation */
40 if (b_i < 0) /* examine sign situation */
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);
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);
56 c[0] = c[1] = c[2] = c[3] = 0;
58 for (i = 0; i < 4; i++) /* we do assume 32 * 32 = 64 bit multiplication */
59 for (j = 0; j < 4; j++)
61 u_long result_low, result_high;
63 result_low = (u_long)a[i] * (u_long)b[j]; /* partial product */
65 if ((i+j) & 1) /* splits across two result registers */
67 result_high = result_low >> (FRACTION_PREC/2);
68 result_low <<= FRACTION_PREC/2;
71 { /* stays in a result register - except for overflows */
75 if (((c[(i+j)/2] >> 1) + (result_low >> 1)) & (u_int32)((unsigned)1<<(FRACTION_PREC - 1)))
76 result_high++; /* propagate overflows */
78 c[(i+j)/2] += result_low; /* add up partial products */
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 */
83 c[(i+j+1)/2] += result_high;
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]);
92 if (c[3]) /* overflow */
94 i = ((unsigned)1 << (FRACTION_PREC-1)) - 1;
98 { /* take produkt - discarding extra precision */
103 if (neg) /* recover sign */
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));
122 * Revision 4.3 1999/02/21 12:17:37 kardel
123 * 4.91f reconcilation
125 * Revision 4.2 1998/12/20 23:45:28 kardel
126 * fix types and warnings
128 * Revision 4.1 1998/05/24 07:59:57 kardel
129 * conditional debug support
131 * Revision 4.0 1998/04/10 19:46:38 kardel
132 * Start 4.0 release version numbering
134 * Revision 1.1 1998/04/10 19:27:47 kardel
135 * initial NTP VERSION 4 integration of PARSE with GPS166 binary support
137 * Revision 1.1 1997/10/06 21:05:46 kardel
138 * new parse structure