xref: /freebsd/contrib/ntp/libparse/mfp_mul.c (revision 313376588638950ba1e93c403dd8c97bc52fd3a2)
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  */
207