Annotation of embedaddon/ntp/libparse/mfp_mul.c, revision 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