Annotation of embedaddon/ntp/libparse/mfp_mul.c, revision 1.1.1.1
1.1 misho 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: */
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to mfp_mul.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [ELWIX - Embedded LightWeight unIX -] / embedaddon / ntp / libparse