diff options
Diffstat (limited to 'libparse/mfp_mul.c')
| -rw-r--r-- | libparse/mfp_mul.c | 207 |
1 files changed, 207 insertions, 0 deletions
diff --git a/libparse/mfp_mul.c b/libparse/mfp_mul.c new file mode 100644 index 0000000..43c1e39 --- /dev/null +++ b/libparse/mfp_mul.c @@ -0,0 +1,207 @@ +/* + * /src/NTP/ntp4-dev/libparse/mfp_mul.c,v 4.9 2005/07/17 20:34:40 kardel RELEASE_20050717_A + * + * mfp_mul.c,v 4.9 2005/07/17 20:34:40 kardel RELEASE_20050717_A + * + * $Created: Sat Aug 16 20:35:08 1997 $ + * + * Copyright (c) 1997-2005 by Frank Kardel <kardel <AT> ntp.org> + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions + * are met: + * 1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the distribution. + * 3. Neither the name of the author nor the names of its contributors + * may be used to endorse or promote products derived from this software + * without specific prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND + * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE + * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL + * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS + * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) + * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT + * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY + * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF + * SUCH DAMAGE. + * + */ +#include <config.h> +#include <stdio.h> +#include "ntp_stdlib.h" +#include "ntp_types.h" +#include "ntp_fp.h" + +#define LOW_MASK (u_int32)((1<<(FRACTION_PREC/2))-1) +#define HIGH_MASK (u_int32)(LOW_MASK << (FRACTION_PREC/2)) + +/* + * for those who worry about overflows (possibly triggered by static analysis tools): + * + * Largest value of a 2^n bit number is 2^n-1. + * Thus the result is: (2^n-1)*(2^n-1) = 2^2n - 2^n - 2^n + 1 < 2^2n + * Here overflow can not happen for 2 reasons: + * 1) the code actually multiplies the absolute values of two signed + * 64bit quantities.thus effectively multiplying 2 63bit quantities. + * 2) Carry propagation is from low to high, building principle is + * addition, so no storage for the 2^2n term from above is needed. + */ + +void +mfp_mul( + int32 *o_i, + u_int32 *o_f, + int32 a_i, + u_int32 a_f, + int32 b_i, + u_int32 b_f + ) +{ + int32 i, j; + u_int32 f; + u_long a[4]; /* operand a */ + u_long b[4]; /* operand b */ + u_long c[5]; /* result c - 5 items for performance - see below */ + u_long carry; + + int neg = 0; + + if (a_i < 0) /* examine sign situation */ + { + neg = 1; + M_NEG(a_i, a_f); + } + + if (b_i < 0) /* examine sign situation */ + { + neg = !neg; + M_NEG(b_i, b_f); + } + + a[0] = a_f & LOW_MASK; /* prepare a operand */ + a[1] = (a_f & HIGH_MASK) >> (FRACTION_PREC/2); + a[2] = a_i & LOW_MASK; + a[3] = (a_i & HIGH_MASK) >> (FRACTION_PREC/2); + + b[0] = b_f & LOW_MASK; /* prepare b operand */ + b[1] = (b_f & HIGH_MASK) >> (FRACTION_PREC/2); + b[2] = b_i & LOW_MASK; + b[3] = (b_i & HIGH_MASK) >> (FRACTION_PREC/2); + + c[0] = c[1] = c[2] = c[3] = c[4] = 0; + + for (i = 0; i < 4; i++) /* we do assume 32 * 32 = 64 bit multiplication */ + for (j = 0; j < 4; j++) + { + u_long result_low, result_high; + int low_index = (i+j)/2; /* formal [0..3] - index for low long word */ + int mid_index = 1+low_index; /* formal [1..4]! - index for high long word + will generate unecessary add of 0 to c[4] + but save 15 'if (result_high) expressions' */ + int high_index = 1+mid_index; /* formal [2..5]! - index for high word overflow + - only assigned on overflow (limits range to 2..3) */ + + result_low = (u_long)a[i] * (u_long)b[j]; /* partial product */ + + if ((i+j) & 1) /* splits across two result registers */ + { + result_high = result_low >> (FRACTION_PREC/2); + result_low <<= FRACTION_PREC/2; + carry = (unsigned)1<<(FRACTION_PREC/2); + } + else + { /* stays in a result register - except for overflows */ + result_high = 0; + carry = 1; + } + + if (((c[low_index] >> 1) + (result_low >> 1) + ((c[low_index] & result_low & carry) != 0)) & + (u_int32)((unsigned)1<<(FRACTION_PREC - 1))) { + result_high++; /* propagate overflows */ + } + + c[low_index] += result_low; /* add up partial products */ + + if (((c[mid_index] >> 1) + (result_high >> 1) + ((c[mid_index] & result_high & 1) != 0)) & + (u_int32)((unsigned)1<<(FRACTION_PREC - 1))) { + c[high_index]++; /* propagate overflows of high word sum */ + } + + c[mid_index] += result_high; /* will add a 0 to c[4] once but saves 15 if conditions */ + } + +#ifdef DEBUG + if (debug > 6) + printf("mfp_mul: 0x%04lx%04lx%04lx%04lx * 0x%04lx%04lx%04lx%04lx = 0x%08lx%08lx%08lx%08lx\n", + a[3], a[2], a[1], a[0], b[3], b[2], b[1], b[0], c[3], c[2], c[1], c[0]); +#endif + + if (c[3]) /* overflow */ + { + i = ((unsigned)1 << (FRACTION_PREC-1)) - 1; + f = ~(unsigned)0; + } + else + { /* take produkt - discarding extra precision */ + i = c[2]; + f = c[1]; + } + + if (neg) /* recover sign */ + { + M_NEG(i, f); + } + + *o_i = i; + *o_f = f; + +#ifdef DEBUG + if (debug > 6) + printf("mfp_mul: %s * %s => %s\n", + mfptoa((u_long)a_i, a_f, 6), + mfptoa((u_long)b_i, b_f, 6), + mfptoa((u_long)i, f, 6)); +#endif +} + +/* + * History: + * + * mfp_mul.c,v + * Revision 4.9 2005/07/17 20:34:40 kardel + * correct carry propagation implementation + * + * Revision 4.8 2005/07/12 16:17:26 kardel + * add explanation why we do not write into c[4] + * + * Revision 4.7 2005/04/16 17:32:10 kardel + * update copyright + * + * Revision 4.6 2004/11/14 15:29:41 kardel + * support PPSAPI, upgrade Copyright to Berkeley style + * + * Revision 4.3 1999/02/21 12:17:37 kardel + * 4.91f reconcilation + * + * Revision 4.2 1998/12/20 23:45:28 kardel + * fix types and warnings + * + * Revision 4.1 1998/05/24 07:59:57 kardel + * conditional debug support + * + * Revision 4.0 1998/04/10 19:46:38 kardel + * Start 4.0 release version numbering + * + * Revision 1.1 1998/04/10 19:27:47 kardel + * initial NTP VERSION 4 integration of PARSE with GPS166 binary support + * + * Revision 1.1 1997/10/06 21:05:46 kardel + * new parse structure + * + */ |