mfp_mul.c revision 285612
1228753Smm/* 2228753Smm * /src/NTP/ntp4-dev/libparse/mfp_mul.c,v 4.9 2005/07/17 20:34:40 kardel RELEASE_20050717_A 3228753Smm * 4228753Smm * mfp_mul.c,v 4.9 2005/07/17 20:34:40 kardel RELEASE_20050717_A 5228753Smm * 6228753Smm * $Created: Sat Aug 16 20:35:08 1997 $ 7228753Smm * 8228753Smm * Copyright (c) 1997-2005 by Frank Kardel <kardel <AT> ntp.org> 9228753Smm * 10228753Smm * Redistribution and use in source and binary forms, with or without 11228753Smm * modification, are permitted provided that the following conditions 12228753Smm * are met: 13228753Smm * 1. Redistributions of source code must retain the above copyright 14228753Smm * notice, this list of conditions and the following disclaimer. 15228753Smm * 2. Redistributions in binary form must reproduce the above copyright 16228753Smm * notice, this list of conditions and the following disclaimer in the 17228753Smm * documentation and/or other materials provided with the distribution. 18228753Smm * 3. Neither the name of the author nor the names of its contributors 19228753Smm * may be used to endorse or promote products derived from this software 20228753Smm * without specific prior written permission. 21228753Smm * 22228753Smm * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND 23228753Smm * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 24228753Smm * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 25228753Smm * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE 26228753Smm * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 27228753Smm * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 28228753Smm * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 29228763Smm * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 30228753Smm * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 31228753Smm * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 32228753Smm * SUCH DAMAGE. 33228753Smm * 34228753Smm */ 35228753Smm#include <config.h> 36228753Smm#include <stdio.h> 37228753Smm#include "ntp_stdlib.h" 38228753Smm#include "ntp_types.h" 39228753Smm#include "ntp_fp.h" 40228753Smm 41228753Smm#define LOW_MASK (u_int32)((1<<(FRACTION_PREC/2))-1) 42228753Smm#define HIGH_MASK (u_int32)(LOW_MASK << (FRACTION_PREC/2)) 43228753Smm 44228753Smm/* 45228753Smm * for those who worry about overflows (possibly triggered by static analysis tools): 46228753Smm * 47228753Smm * Largest value of a 2^n bit number is 2^n-1. 48228753Smm * Thus the result is: (2^n-1)*(2^n-1) = 2^2n - 2^n - 2^n + 1 < 2^2n 49228753Smm * Here overflow can not happen for 2 reasons: 50228753Smm * 1) the code actually multiplies the absolute values of two signed 51228753Smm * 64bit quantities.thus effectively multiplying 2 63bit quantities. 52228753Smm * 2) Carry propagation is from low to high, building principle is 53228753Smm * addition, so no storage for the 2^2n term from above is needed. 54232153Smm */ 55228753Smm 56228753Smmvoid 57228753Smmmfp_mul( 58228753Smm int32 *o_i, 59228753Smm u_int32 *o_f, 60228753Smm int32 a_i, 61228753Smm u_int32 a_f, 62228753Smm int32 b_i, 63228753Smm u_int32 b_f 64248616Smm ) 65228753Smm{ 66305192Smm int32 i, j; 67232153Smm u_int32 f; 68228753Smm u_long a[4]; /* operand a */ 69228753Smm u_long b[4]; /* operand b */ 70228753Smm u_long c[5]; /* result c - 5 items for performance - see below */ 71248616Smm u_long carry; 72228753Smm 73228753Smm int neg = 0; 74228753Smm 75228753Smm if (a_i < 0) /* examine sign situation */ 76248616Smm { 77302001Smm neg = 1; 78228753Smm M_NEG(a_i, a_f); 79248616Smm } 80228753Smm 81228753Smm if (b_i < 0) /* examine sign situation */ 82228753Smm { 83228753Smm neg = !neg; 84228753Smm M_NEG(b_i, b_f); 85302001Smm } 86228753Smm 87228753Smm a[0] = a_f & LOW_MASK; /* prepare a operand */ 88228753Smm a[1] = (a_f & HIGH_MASK) >> (FRACTION_PREC/2); 89228753Smm a[2] = a_i & LOW_MASK; 90228753Smm a[3] = (a_i & HIGH_MASK) >> (FRACTION_PREC/2); 91248616Smm 92228753Smm b[0] = b_f & LOW_MASK; /* prepare b operand */ 93228753Smm b[1] = (b_f & HIGH_MASK) >> (FRACTION_PREC/2); 94228753Smm b[2] = b_i & LOW_MASK; 95324418Smm b[3] = (b_i & HIGH_MASK) >> (FRACTION_PREC/2); 96228753Smm 97228753Smm c[0] = c[1] = c[2] = c[3] = c[4] = 0; 98228753Smm 99228753Smm for (i = 0; i < 4; i++) /* we do assume 32 * 32 = 64 bit multiplication */ 100228753Smm for (j = 0; j < 4; j++) 101228753Smm { 102228753Smm u_long result_low, result_high; 103228753Smm int low_index = (i+j)/2; /* formal [0..3] - index for low long word */ 104228753Smm int mid_index = 1+low_index; /* formal [1..4]! - index for high long word 105228753Smm will generate unecessary add of 0 to c[4] 106228753Smm but save 15 'if (result_high) expressions' */ 107228753Smm int high_index = 1+mid_index; /* formal [2..5]! - index for high word overflow 108228753Smm - only assigned on overflow (limits range to 2..3) */ 109228753Smm 110228753Smm result_low = (u_long)a[i] * (u_long)b[j]; /* partial product */ 111228753Smm 112228753Smm if ((i+j) & 1) /* splits across two result registers */ 113228753Smm { 114228753Smm result_high = result_low >> (FRACTION_PREC/2); 115228753Smm result_low <<= FRACTION_PREC/2; 116228753Smm carry = (unsigned)1<<(FRACTION_PREC/2); 117228753Smm } 118228753Smm else 119228753Smm { /* stays in a result register - except for overflows */ 120228753Smm result_high = 0; 121228753Smm carry = 1; 122232153Smm } 123228753Smm 124228753Smm if (((c[low_index] >> 1) + (result_low >> 1) + ((c[low_index] & result_low & carry) != 0)) & 125228753Smm (u_int32)((unsigned)1<<(FRACTION_PREC - 1))) { 126228753Smm result_high++; /* propagate overflows */ 127228753Smm } 128228753Smm 129228753Smm c[low_index] += result_low; /* add up partial products */ 130228753Smm 131228753Smm if (((c[mid_index] >> 1) + (result_high >> 1) + ((c[mid_index] & result_high & 1) != 0)) & 132228753Smm (u_int32)((unsigned)1<<(FRACTION_PREC - 1))) { 133228753Smm c[high_index]++; /* propagate overflows of high word sum */ 134228753Smm } 135228753Smm 136228753Smm c[mid_index] += result_high; /* will add a 0 to c[4] once but saves 15 if conditions */ 137228753Smm } 138228753Smm 139228753Smm#ifdef DEBUG 140228753Smm if (debug > 6) 141228753Smm printf("mfp_mul: 0x%04lx%04lx%04lx%04lx * 0x%04lx%04lx%04lx%04lx = 0x%08lx%08lx%08lx%08lx\n", 142228753Smm a[3], a[2], a[1], a[0], b[3], b[2], b[1], b[0], c[3], c[2], c[1], c[0]); 143228753Smm#endif 144228753Smm 145228753Smm if (c[3]) /* overflow */ 146228753Smm { 147228753Smm i = ((unsigned)1 << (FRACTION_PREC-1)) - 1; 148228753Smm f = ~(unsigned)0; 149228753Smm } 150228753Smm else 151228753Smm { /* take produkt - discarding extra precision */ 152228753Smm i = c[2]; 153228753Smm f = c[1]; 154228753Smm } 155228753Smm 156228753Smm if (neg) /* recover sign */ 157228753Smm { 158228753Smm M_NEG(i, f); 159228753Smm } 160228753Smm 161228753Smm *o_i = i; 162228753Smm *o_f = f; 163228753Smm 164228753Smm#ifdef DEBUG 165228753Smm if (debug > 6) 166228753Smm printf("mfp_mul: %s * %s => %s\n", 167228753Smm mfptoa((u_long)a_i, a_f, 6), 168228753Smm mfptoa((u_long)b_i, b_f, 6), 169228753Smm mfptoa((u_long)i, f, 6)); 170228753Smm#endif 171228753Smm} 172228753Smm 173228753Smm/* 174228753Smm * History: 175228753Smm * 176228753Smm * mfp_mul.c,v 177228753Smm * Revision 4.9 2005/07/17 20:34:40 kardel 178228753Smm * correct carry propagation implementation 179228753Smm * 180228753Smm * Revision 4.8 2005/07/12 16:17:26 kardel 181228753Smm * add explanation why we do not write into c[4] 182228753Smm * 183228753Smm * Revision 4.7 2005/04/16 17:32:10 kardel 184228753Smm * update copyright 185228753Smm * 186228753Smm * Revision 4.6 2004/11/14 15:29:41 kardel 187228753Smm * support PPSAPI, upgrade Copyright to Berkeley style 188228753Smm * 189228753Smm * Revision 4.3 1999/02/21 12:17:37 kardel 190228753Smm * 4.91f reconcilation 191228753Smm * 192228753Smm * Revision 4.2 1998/12/20 23:45:28 kardel 193228753Smm * fix types and warnings 194228753Smm * 195228753Smm * Revision 4.1 1998/05/24 07:59:57 kardel 196228753Smm * conditional debug support 197228753Smm * 198228753Smm * Revision 4.0 1998/04/10 19:46:38 kardel 199228753Smm * Start 4.0 release version numbering 200228753Smm * 201232153Smm * Revision 1.1 1998/04/10 19:27:47 kardel 202228753Smm * initial NTP VERSION 4 integration of PARSE with GPS166 binary support 203228753Smm * 204228753Smm * Revision 1.1 1997/10/06 21:05:46 kardel 205228753Smm * new parse structure 206228753Smm * 207228753Smm */ 208228753Smm