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 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