fis-gtm/sr_port/eb_muldiv.c

317 lines
7.3 KiB
C

/****************************************************************
* *
* Copyright 2001 Sanchez Computer Associates, Inc. *
* *
* This source code contains the intellectual property *
* of its copyright holder(s), and is made available *
* under a license. If you do not know the terms of *
* the license, please stop and do not read further. *
* *
****************************************************************/
/* eb_muldiv - emulate extended precision (18-digit) multiplication and division */
#include "mdef.h"
#include "arit.h"
#include "eb_muldiv.h"
#define FLO_HI 1e9
#define FLO_LO 1e8
#define FLO_BIAS 1e3
#define DFLOAT2MINT(X,DF) (X[1] = DF, X[0] = (DF - (double)X[1])*FLO_HI)
#define RADIX 10000
/* eb_int_mul - multiply two GT.M INT's
*
* input
* v1, u1 - GT.M INT's
*
* output
* if the product will fit into a GT.M INT format:
* function result = FALSE => no promotion necessary
* p[] = INT product of (v1*m1)
* else (implies overflow out of INT format):
* function result = TRUE => promotion to extended precision necessary
* p[] = undefined
*/
bool eb_int_mul (int4 v1, int4 u1, int4 p[])
{
double pf;
int4 tp[2], promote;
promote = TRUE; /* promote if overflow or too many significant fractional digits */
pf = (double)u1*(double)v1/FLO_BIAS;
if ((pf < FLO_HI) && (pf > -FLO_HI))
{
DFLOAT2MINT(tp, pf);
if (tp[0] == 0) /* don't need extra precision */
{
promote = FALSE;
p[0] = tp[0]; p[1] = tp[1];
}
}
return promote;
}
/* eb_mul - multiply two GT.M extended precision numeric values
*
* input
* v[], u[] - GT.M extended precision numeric value mantissas
*
* output
* function result = scale factor of result
* p[] = GT.M extended precision mantissa of (u*v)
*/
int4 eb_mul (int4 v[], int4 u[], int4 p[]) /* p = u*v */
{
short i, j;
int4 acc, carry, m1[5], m2[5], prod[9], scale;
/* Throughout, larger index => more significance. */
for (i = 0 ; i < 9 ; i++)
prod[i] = 0;
/* Break up 2-4-(3/1)-4-4 */
m1[0] = v[0] % RADIX;
m2[0] = u[0] % RADIX;
m1[1] = (v[0]/RADIX) % RADIX;
m2[1] = (u[0]/RADIX) % RADIX;
m1[2] = (v[1] % (RADIX/10))*10 + v[0]/(RADIX*RADIX);
m2[2] = (u[1] % (RADIX/10))*10 + u[0]/(RADIX*RADIX);
m1[3] = (v[1]/(RADIX/10)) % RADIX;
m2[3] = (u[1]/(RADIX/10)) % RADIX;
m1[4] = v[1]/((RADIX/10)*RADIX);
m2[4] = u[1]/((RADIX/10)*RADIX);
for (j = 0 ; j <= 4 ; j++)
{
if (m2[j] != 0)
{
for (i = 0, carry = 0 ; i <= 4 ; i++)
{
acc = m1[i]*m2[j] + prod[i+j] + carry;
prod[i+j] = acc % RADIX;
carry = acc / RADIX;
}
if ( 9 > i+j)
prod[i+j] = carry;
else
if (0 != carry)
assert(FALSE);
}
}
if (prod[8] >= RADIX/10)
{
/* Assemble back 4-4-1/3-4-2 */
scale = 0; /* no scaling needed */
p[0] = ((prod[6]%1000)*RADIX + prod[5])*(RADIX/ 100) + (prod[4]/ 100);
p[1] = ( prod[8] *RADIX + prod[7])*(RADIX/1000) + (prod[6]/1000);
}
else /* prod[8] < RADIX/10 [means not normalized] */
{
/* Assemble back 3-4-2/2-4-3 */
scale = -1; /* to compensate for normalization */
p[0] = ((prod[6]%100)*RADIX + prod[5])*(RADIX/ 10) + (prod[4]/ 10);
p[1] = ( prod[8] *RADIX + prod[7])*(RADIX/100) + (prod[6]/100);
}
return scale;
}
/* eb_mvint_div - divide to GT.M INT's
*
* input
* v, u - INT's to be divided
*
* output
* if the quotient will fit into a GT.M INT:
* function value = FALSE => no promotion necessary
* q[] = INT quotient of (v/u)
* else (implies overflow out of GT.M INT forat):
* function value = TRUE => promotion to extended precision necessary
* q[] = undefined
*/
bool eb_mvint_div (int4 v, int4 u, int4 q[])
{
double qf;
int4 tq[2], promote;
promote = TRUE; /* promote if overflow or too many significant fractional digits */
qf = (double)v*FLO_BIAS/(double)u;
if ((qf < FLO_HI) && (qf > -FLO_HI))
{
DFLOAT2MINT(tq, qf);
if (tq[0] == 0) /* don't need extra word of precision */
{
promote = FALSE;
q[0] = tq[0]; q[1] = tq[1];
}
}
return promote;
}
/* eb_int_div - integer division of two GT.M INT's
*
* input
* v1, u1 - GT.M INT's to be divided
*
* output
* if result fits into a GT.M INT:
* function value = FALSE => no promotion necessary
* q[] = INT result of (v1\u1)
* else (implies some sort of overflow):
* function result = TRUE => promotion to extended precision necessary
* q[] = undefined
*/
bool eb_int_div (int4 v1, int4 u1, int4 q[])
{
double qf;
qf= (double)v1*FLO_BIAS/(double)u1;
if (qf < FLO_HI && qf > -FLO_HI)
{
DFLOAT2MINT(q,qf);
q[1]= (q[1]/MV_BIAS)*MV_BIAS;
return FALSE;
}
else
{
return TRUE;
}
}
/* eb_div - divide two GT.M extended precision numeric values
*
* input
* x[], y[] - GT.M extended precision numeric value mantissas
*
* output
* function result = scale factor of result
* q[] = GT.M extended precision mantissa of (y/x)
*/
int4 eb_div (int4 x[], int4 y[], int4 q[]) /* q = y/x */
{
int4 borrow, carry, i, j, scale, prod, qx[5], xx[5], yx[10];
for (i = 0 ; i < 5 ; i++)
yx[i] = 0;
if (x[1] < y[1] || (x[1] == y[1] && x[0] <= y[0])) /* i.e., if x <= y */
{
/* Break y apart 3-4-2/2-4-3 */
scale = 1;
yx[5] = (y[0]%(RADIX/10))*10;
yx[6] = (y[0]/(RADIX/10))%RADIX;
yx[7] = (y[1]%(RADIX/100))*(RADIX/100) + y[0]/((RADIX/10)*RADIX);
yx[8] = (y[1]/(RADIX/100))%RADIX;
yx[9] = y[1]/((RADIX/100)*RADIX);
}
else
{
/* Break y apart 4-4-1/3-4-2 */
scale = 0;
yx[5] = (y[0]%(RADIX/100))*100;
yx[6] = (y[0]/(RADIX/100))%RADIX;
yx[7] = (y[1]%(RADIX/1000))*(RADIX/10) + y[0]/((RADIX/100)*RADIX);
yx[8] = (y[1]/(RADIX/1000))%RADIX;
yx[9] = y[1]/((RADIX/1000)*RADIX);
}
/* Break x apart 4-4-1/3-4-2 */
xx[0] = (x[0]%(RADIX/100))*100;
xx[1] = (x[0]/(RADIX/100))%RADIX;
xx[2] = (x[1]%(RADIX/1000))*(RADIX/10) + x[0]/((RADIX/100)*RADIX);
xx[3] = (x[1]/(RADIX/1000))%RADIX;
xx[4] = x[1]/((RADIX/1000)*RADIX);
assert (yx[9] <= xx[4]);
for (i = 4 ; i >= 0 ; i--)
{
qx[i] = (yx[i+5]*RADIX + yx[i+4]) / xx[4];
if (qx[i] != 0)
{
/* Multiply x by qx[i] and subtract from remainder. */
for (j = 0, borrow = 0 ; j <= 4 ; j++)
{
prod = qx[i]*xx[j] + borrow;
borrow = prod/RADIX;
yx[i+j] -= (prod%RADIX);
if (yx[i+j] < 0)
{
yx[i+j] += RADIX;
borrow ++;
}
}
yx[i+5] -= borrow;
while (yx[i+5] < 0)
{
qx[i] --; /* estimate too high */
for (j = 0, carry = 0 ; j <= 4 ; j++)
{
yx[i+j] += (xx[j] + carry);
carry = yx[i+j]/RADIX;
yx[i+j] %= RADIX;
}
yx[i+5] += carry;
}
}
assert (0 <= qx[i] && qx[i] < RADIX); /* make sure in range */
assert (yx[i+5] == 0); /* check that remainder doesn't overflow */
}
/* Assemble q 4-4-1/3-4-2 */
q[0] = ((qx[2]%1000)*RADIX + qx[1])*100 + (qx[0]/ 100);
q[1] = ( qx[4] *RADIX + qx[3])* 10 + (qx[2]/1000);
assert ( (FLO_LO <= q[1] && q[1] < FLO_HI)
|| (q[1] == 0 && q[0] == 0 && y[1] == 0 && y[0] == 0) );
return scale;
}
/* eb_int_mod - INT modulus of two GT.M INT's
*
* input
* v1, u1 - GT.M INT's
*
* output
* p[] = INT value of (v1 mod u1) == (v1 - (u1*floor(v1/u1)))
*/
void eb_int_mod (int4 v1, int4 u1, int4 p[])
{
int4 quo, rat, neg;
if (u1 == 0 || v1 == 0)
{
p[1]= 0;
}
else
{
quo = v1/u1;
rat = v1 != quo*u1;
neg = (v1 < 0 && u1 > 0) || (v1 > 0 && u1 < 0);
p[1] = v1 - u1*(quo - (neg && rat));
}
return;
}