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