207 lines
6.4 KiB
C
207 lines
6.4 KiB
C
|
/*
|
||
|
|
* /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 <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
|
||
|
|
*
|
||
|
|
*/
|