#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))
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];
u_long b[4];
u_long c[4];
int neg = 0;
if (a_i < 0)
{
neg = 1;
M_NEG(a_i, a_f);
}
if (b_i < 0)
{
neg = !neg;
M_NEG(b_i, b_f);
}
a[0] = a_f & LOW_MASK;
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;
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] = 0;
for (i = 0; i < 4; i++)
for (j = 0; j < 4; j++)
{
u_long result_low, result_high;
result_low = (u_long)a[i] * (u_long)b[j];
if ((i+j) & 1)
{
result_high = result_low >> (FRACTION_PREC/2);
result_low <<= FRACTION_PREC/2;
}
else
{
result_high = 0;
}
if (((c[(i+j)/2] >> 1) + (result_low >> 1)) & (u_int32)((unsigned)1<<(FRACTION_PREC - 1)))
result_high++;
c[(i+j)/2] += result_low;
if (((c[(i+j+1)/2] >> 1) + (result_high >> 1)) & (u_int32)((unsigned)1<<(FRACTION_PREC - 1)))
c[1+(i+j)/2]++;
c[(i+j+1)/2] += result_high;
}
#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])
{
i = ((unsigned)1 << (FRACTION_PREC-1)) - 1;
f = ~(unsigned)0;
}
else
{
i = c[2];
f = c[1];
}
if (neg)
{
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
}