1181834Sroberto/* 2181834Sroberto * /src/NTP/ntp4-dev/libparse/mfp_mul.c,v 4.9 2005/07/17 20:34:40 kardel RELEASE_20050717_A 3181834Sroberto * 4181834Sroberto * mfp_mul.c,v 4.9 2005/07/17 20:34:40 kardel RELEASE_20050717_A 5181834Sroberto * 6181834Sroberto * $Created: Sat Aug 16 20:35:08 1997 $ 7181834Sroberto * 8181834Sroberto * Copyright (c) 1997-2005 by Frank Kardel <kardel <AT> ntp.org> 9181834Sroberto * 10181834Sroberto * Redistribution and use in source and binary forms, with or without 11181834Sroberto * modification, are permitted provided that the following conditions 12181834Sroberto * are met: 13181834Sroberto * 1. Redistributions of source code must retain the above copyright 14181834Sroberto * notice, this list of conditions and the following disclaimer. 15181834Sroberto * 2. Redistributions in binary form must reproduce the above copyright 16181834Sroberto * notice, this list of conditions and the following disclaimer in the 17181834Sroberto * documentation and/or other materials provided with the distribution. 18181834Sroberto * 3. Neither the name of the author nor the names of its contributors 19181834Sroberto * may be used to endorse or promote products derived from this software 20181834Sroberto * without specific prior written permission. 21181834Sroberto * 22181834Sroberto * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND 23181834Sroberto * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 24181834Sroberto * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 25181834Sroberto * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE 26181834Sroberto * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 27181834Sroberto * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 28181834Sroberto * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 29181834Sroberto * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 30181834Sroberto * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 31181834Sroberto * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 32181834Sroberto * SUCH DAMAGE. 33181834Sroberto * 34181834Sroberto */ 35290001Sglebius#include <config.h> 36181834Sroberto#include <stdio.h> 37181834Sroberto#include "ntp_stdlib.h" 38181834Sroberto#include "ntp_types.h" 39181834Sroberto#include "ntp_fp.h" 40181834Sroberto 41181834Sroberto#define LOW_MASK (u_int32)((1<<(FRACTION_PREC/2))-1) 42181834Sroberto#define HIGH_MASK (u_int32)(LOW_MASK << (FRACTION_PREC/2)) 43181834Sroberto 44181834Sroberto/* 45181834Sroberto * for those who worry about overflows (possibly triggered by static analysis tools): 46181834Sroberto * 47181834Sroberto * Largest value of a 2^n bit number is 2^n-1. 48181834Sroberto * Thus the result is: (2^n-1)*(2^n-1) = 2^2n - 2^n - 2^n + 1 < 2^2n 49181834Sroberto * Here overflow can not happen for 2 reasons: 50181834Sroberto * 1) the code actually multiplies the absolute values of two signed 51181834Sroberto * 64bit quantities.thus effectively multiplying 2 63bit quantities. 52181834Sroberto * 2) Carry propagation is from low to high, building principle is 53181834Sroberto * addition, so no storage for the 2^2n term from above is needed. 54181834Sroberto */ 55181834Sroberto 56181834Srobertovoid 57181834Srobertomfp_mul( 58181834Sroberto int32 *o_i, 59181834Sroberto u_int32 *o_f, 60181834Sroberto int32 a_i, 61181834Sroberto u_int32 a_f, 62181834Sroberto int32 b_i, 63181834Sroberto u_int32 b_f 64181834Sroberto ) 65181834Sroberto{ 66181834Sroberto int32 i, j; 67181834Sroberto u_int32 f; 68181834Sroberto u_long a[4]; /* operand a */ 69181834Sroberto u_long b[4]; /* operand b */ 70181834Sroberto u_long c[5]; /* result c - 5 items for performance - see below */ 71181834Sroberto u_long carry; 72181834Sroberto 73181834Sroberto int neg = 0; 74181834Sroberto 75181834Sroberto if (a_i < 0) /* examine sign situation */ 76181834Sroberto { 77181834Sroberto neg = 1; 78181834Sroberto M_NEG(a_i, a_f); 79181834Sroberto } 80181834Sroberto 81181834Sroberto if (b_i < 0) /* examine sign situation */ 82181834Sroberto { 83181834Sroberto neg = !neg; 84181834Sroberto M_NEG(b_i, b_f); 85181834Sroberto } 86181834Sroberto 87181834Sroberto a[0] = a_f & LOW_MASK; /* prepare a operand */ 88181834Sroberto a[1] = (a_f & HIGH_MASK) >> (FRACTION_PREC/2); 89181834Sroberto a[2] = a_i & LOW_MASK; 90181834Sroberto a[3] = (a_i & HIGH_MASK) >> (FRACTION_PREC/2); 91181834Sroberto 92181834Sroberto b[0] = b_f & LOW_MASK; /* prepare b operand */ 93181834Sroberto b[1] = (b_f & HIGH_MASK) >> (FRACTION_PREC/2); 94181834Sroberto b[2] = b_i & LOW_MASK; 95181834Sroberto b[3] = (b_i & HIGH_MASK) >> (FRACTION_PREC/2); 96181834Sroberto 97181834Sroberto c[0] = c[1] = c[2] = c[3] = c[4] = 0; 98181834Sroberto 99181834Sroberto for (i = 0; i < 4; i++) /* we do assume 32 * 32 = 64 bit multiplication */ 100181834Sroberto for (j = 0; j < 4; j++) 101181834Sroberto { 102181834Sroberto u_long result_low, result_high; 103181834Sroberto int low_index = (i+j)/2; /* formal [0..3] - index for low long word */ 104181834Sroberto int mid_index = 1+low_index; /* formal [1..4]! - index for high long word 105181834Sroberto will generate unecessary add of 0 to c[4] 106181834Sroberto but save 15 'if (result_high) expressions' */ 107181834Sroberto int high_index = 1+mid_index; /* formal [2..5]! - index for high word overflow 108181834Sroberto - only assigned on overflow (limits range to 2..3) */ 109181834Sroberto 110181834Sroberto result_low = (u_long)a[i] * (u_long)b[j]; /* partial product */ 111181834Sroberto 112181834Sroberto if ((i+j) & 1) /* splits across two result registers */ 113181834Sroberto { 114181834Sroberto result_high = result_low >> (FRACTION_PREC/2); 115181834Sroberto result_low <<= FRACTION_PREC/2; 116181834Sroberto carry = (unsigned)1<<(FRACTION_PREC/2); 117181834Sroberto } 118181834Sroberto else 119181834Sroberto { /* stays in a result register - except for overflows */ 120181834Sroberto result_high = 0; 121181834Sroberto carry = 1; 122181834Sroberto } 123181834Sroberto 124181834Sroberto if (((c[low_index] >> 1) + (result_low >> 1) + ((c[low_index] & result_low & carry) != 0)) & 125181834Sroberto (u_int32)((unsigned)1<<(FRACTION_PREC - 1))) { 126181834Sroberto result_high++; /* propagate overflows */ 127181834Sroberto } 128181834Sroberto 129181834Sroberto c[low_index] += result_low; /* add up partial products */ 130181834Sroberto 131181834Sroberto if (((c[mid_index] >> 1) + (result_high >> 1) + ((c[mid_index] & result_high & 1) != 0)) & 132181834Sroberto (u_int32)((unsigned)1<<(FRACTION_PREC - 1))) { 133181834Sroberto c[high_index]++; /* propagate overflows of high word sum */ 134181834Sroberto } 135181834Sroberto 136181834Sroberto c[mid_index] += result_high; /* will add a 0 to c[4] once but saves 15 if conditions */ 137181834Sroberto } 138181834Sroberto 139181834Sroberto#ifdef DEBUG 140181834Sroberto if (debug > 6) 141181834Sroberto printf("mfp_mul: 0x%04lx%04lx%04lx%04lx * 0x%04lx%04lx%04lx%04lx = 0x%08lx%08lx%08lx%08lx\n", 142181834Sroberto a[3], a[2], a[1], a[0], b[3], b[2], b[1], b[0], c[3], c[2], c[1], c[0]); 143181834Sroberto#endif 144181834Sroberto 145181834Sroberto if (c[3]) /* overflow */ 146181834Sroberto { 147181834Sroberto i = ((unsigned)1 << (FRACTION_PREC-1)) - 1; 148181834Sroberto f = ~(unsigned)0; 149181834Sroberto } 150181834Sroberto else 151181834Sroberto { /* take produkt - discarding extra precision */ 152181834Sroberto i = c[2]; 153181834Sroberto f = c[1]; 154181834Sroberto } 155181834Sroberto 156181834Sroberto if (neg) /* recover sign */ 157181834Sroberto { 158181834Sroberto M_NEG(i, f); 159181834Sroberto } 160181834Sroberto 161181834Sroberto *o_i = i; 162181834Sroberto *o_f = f; 163181834Sroberto 164181834Sroberto#ifdef DEBUG 165181834Sroberto if (debug > 6) 166181834Sroberto printf("mfp_mul: %s * %s => %s\n", 167181834Sroberto mfptoa((u_long)a_i, a_f, 6), 168181834Sroberto mfptoa((u_long)b_i, b_f, 6), 169181834Sroberto mfptoa((u_long)i, f, 6)); 170181834Sroberto#endif 171181834Sroberto} 172181834Sroberto 173181834Sroberto/* 174181834Sroberto * History: 175181834Sroberto * 176181834Sroberto * mfp_mul.c,v 177181834Sroberto * Revision 4.9 2005/07/17 20:34:40 kardel 178181834Sroberto * correct carry propagation implementation 179181834Sroberto * 180181834Sroberto * Revision 4.8 2005/07/12 16:17:26 kardel 181181834Sroberto * add explanation why we do not write into c[4] 182181834Sroberto * 183181834Sroberto * Revision 4.7 2005/04/16 17:32:10 kardel 184181834Sroberto * update copyright 185181834Sroberto * 186181834Sroberto * Revision 4.6 2004/11/14 15:29:41 kardel 187181834Sroberto * support PPSAPI, upgrade Copyright to Berkeley style 188181834Sroberto * 189181834Sroberto * Revision 4.3 1999/02/21 12:17:37 kardel 190181834Sroberto * 4.91f reconcilation 191181834Sroberto * 192181834Sroberto * Revision 4.2 1998/12/20 23:45:28 kardel 193181834Sroberto * fix types and warnings 194181834Sroberto * 195181834Sroberto * Revision 4.1 1998/05/24 07:59:57 kardel 196181834Sroberto * conditional debug support 197181834Sroberto * 198181834Sroberto * Revision 4.0 1998/04/10 19:46:38 kardel 199181834Sroberto * Start 4.0 release version numbering 200181834Sroberto * 201181834Sroberto * Revision 1.1 1998/04/10 19:27:47 kardel 202181834Sroberto * initial NTP VERSION 4 integration of PARSE with GPS166 binary support 203181834Sroberto * 204181834Sroberto * Revision 1.1 1997/10/06 21:05:46 kardel 205181834Sroberto * new parse structure 206181834Sroberto * 207181834Sroberto */ 208