iSharkFly-Open
/
fis-gtm
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
fis-gtm/sr_port/op_exp.c

308 lines
8.4 KiB
C
Raw Permalink Blame History

/****************************************************************
* *
* Copyright 2001, 2012 Fidelity Information Services, 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. *
* *
****************************************************************/
#include "mdef.h"
#include "gtm_string.h"
#include <math.h>
#include "arit.h"
#include "stringpool.h"
#include "op.h"
#include "mvalconv.h"
#define ACCUR_PERCENT 0.00000000000000055
#define MAX_M_INT 999999999
LITREF int4 ten_pwr[];
LITREF mval literal_one;
LITREF mval literal_zero;
STATICFNDEF void op_exp_flgchk(mval *mv);
error_def(ERR_NEGFRACPWR);
error_def(ERR_NUMOFLOW);
void op_exp(mval *u, mval* v, mval *p)
{
mval u1, *u1_p;
double accuracy, exponent;
double x, x1, y, z, z2, z3, z4, z5, id, il;
int im0, im1, ie, i, j, j1;
boolean_t fraction = FALSE, in = FALSE;
boolean_t neg = FALSE, even = TRUE;
mval w, zmv;
int4 n, n1;
int4 z1_rnd, z2_rnd, pten;
u1_p = &u1;
memcpy(u1_p, u, SIZEOF(mval));
MV_FORCE_NUM(u1_p);
MV_FORCE_NUM(v);
if ((0 == v->m[1]) && (0 == v->m[0]))
{ /* 0**0 = 1 */
*p = literal_one;
return;
}
if (0 != (v->mvtype & MV_INT))
{ /* Integer-ish exponent (could have up to 3 digits to right of decimal pt) */
n = v->m[1];
if (0 == n)
{ /* anything**0 = 1 where anything != 0 */
*p = literal_one;
return;
}
if (0 != (u1_p->mvtype & MV_INT))
{ /* Integer-ish base */
if (0 == u1_p->m[1])
{ /* 0**anything = 0 */
*p = literal_zero;
return;
}
} else if ((0 == u1_p->m[1]) && (0 == u1_p->m[0]))
{ /* 0**anything = 0 */
*p = literal_zero;
return;
}
n1 = n / MV_BIAS;
if ((n1 * MV_BIAS) == n)
{ /* True non-fractional exponent */
if (0 == v->m[1])
{ /* Duplicate of check on line 58? */
*p = literal_one;
return;
}
if (0 > n1)
{ /* Create inverse due to negative exponent */
op_div((mval *)&literal_one, u1_p, &w);
n1 = -n1;
} else
w = *u1_p;
zmv = literal_one;
for ( ; ; )
{ /* Compute integer exponent */
if (n1 & 1)
op_mul(&zmv, &w, &zmv);
n1 /= 2;
if (!n1)
break;
op_mul(&w, &w, &w);
}
*p = zmv;
return;
} else
{ /* Have non-integer exponent (has fractional component) */
if (0 != (u1_p->mvtype & MV_INT))
{ /* Base is integer-ish */
if (0 > u1_p->m[1])
{ /* Base is negative, invalid exponent expression */
rts_error(VARLSTCNT(1) ERR_NEGFRACPWR);
return;
}
} else
{ /* Base is NOT integer-sh */
if (u1_p->sgn)
{ /* Base is negative, invalid exponent expression */
rts_error(VARLSTCNT(1) ERR_NEGFRACPWR);
return;
}
}
}
} else
{ /* Exponent NOT integer-ish */
if (0 != (u1_p->mvtype & MV_INT))
{ /* Base is integer-ish */
if (0 > u1_p->m[1])
{ /* Base is negative - make positive but record was negative */
u1_p->m[1] = -u1_p->m[1];
neg = TRUE;
}
if (0 == u1_p->m[1])
{ /* 0**anything is 0 */
*p = literal_zero;
return;
}
} else if (u1_p->sgn)
{ /* Base is NOT integer-ish and is negative - clear sign and record with flag */
u1_p->sgn = 0;
neg = TRUE;
if ((0 == u1_p->m[1]) && (0 == u1_p->m[0]))
{ /* 0**anything is zero */
*p = literal_zero;
return;
}
}
if (NUM_DEC_DG_2L > (ie = (v->e - MV_XBIAS))) /* Note assignment */
{ /* Need to determine 2 things:
* 1. Whether this exponent has a fractional part (vs just being large)
* 2. If the base value is negative, result is negative based on whether exponent is even/odd
* Since all exponent digits are visible, determine these items.
*/
if (0 < ie)
{ /* Decimal point shifting to right - see if fully integer*/
if (NUM_DEC_DG_1L > ie)
{ /* Exponent fits in int4 */
for (i = 1, j = 10; (NUM_DEC_DG_1L - ie) > i; j *= 10, i++)
; /* Determine exponential scale */
im1 = v->m[1];
if (0 != (i = im1 % j)) /* Note assignment */
fraction = TRUE;
else
if (0 != (i = im1 % (2 * j)))
even = FALSE;
} else
{
im0 = v->m[0];
if (NUM_DEC_DG_1L == ie)
{ /* Var ie is at max for no loss in single int */
if (0 == im0)
{
im1 = v->m[1];
if (0 != (i = im1 % 2)) /* Note assignment */
even = FALSE;
} else
fraction = TRUE;
} else
{
for (i=1, j=10; (NUM_DEC_DG_2L - ie) > i; j *= 10, i++)
;
if (0 == (i = im0 % j)) /* Note assignment */
{
if (0 != (i = im0 % (2 * j))) /* Note assignment */
even = FALSE;
} else
fraction = TRUE;
}
}
} else /* Var ie is negative (fraction) or 0. This is a non-integerish number so ie
* must be a fraction. ie must be > 0 to be non-integer.
*/
fraction = TRUE;
} else
{
if (NUM_DEC_DG_2L == ie)
{ /* Var ie is max to have no significant digit loss */
im0 = v->m[0];
if (0 != (i = im0 % 2)) /* Note assignment */
even = FALSE;
}
}
if (fraction && neg)
{ /* Fractional exponent and negative base not valid */
rts_error(VARLSTCNT(1) ERR_NEGFRACPWR);
return;
}
}
x = mval2double(u1_p); /* Convert base and exponent to double for feeding to pow*() function */
y = mval2double(v);
z = pow(x, y);
# ifdef UNIX
if (HUGE_VAL == z) /* Infinity return value check */
{
rts_error(VARLSTCNT(1) ERR_NUMOFLOW);
return;
}
# endif
assert(!neg); /* Should be taken care of in one of the op_mul() using sections dealing with whole exponents */
p->sgn = 0; /* Positive numbers only from here on out */
if (0 == z)
{
*p = literal_zero;
return;
}
/* Remaining code's main purpose is to convert the double float value to our internal format taking care that
* actual integer values (e.g. 4**.5 is 2) are correctly rendered. This can be tricky not only because floating
* point values are often imprecise (e.g. 3.0 represented as 2.999999999999999) but because the double-float
* type has fewer digits of accuracy than GTM carries.
*/
accuracy = z * ACCUR_PERCENT;
/* Scale accuracy up in case it turned to zero above */
for (z2 = ACCUR_PERCENT; (0.0 == accuracy); z2 *= 10, accuracy = z * z2)
;
if ((fabs(z) <= accuracy) || (0 >= z))
{ /* Zero equivalency test - some small chance pow() could return a very small negative number. Any
* possible negative numbers must have whole number exponent which is handled earlier via op_mul().
*/
*p = literal_zero;
return;
}
/* Integer check - GT.M pseudo int MV_INT check - must be decimal form 999999.999 */
z2 = floor((z + .0005) * MV_BIAS);
if (MANT_HI > z2)
{ /* Have proper range, check if accuracy warrants conversion */
n1 = (int)z2;
z2 = (double)n1 / MV_BIAS;
if (fabs(z - z2) < accuracy)
{ /* We can treat this as a GT.M int */
((mval_b *)p)->sgne = 0;
p->mvtype = (MV_NM | MV_INT);
p->m[0] = 0;
p->m[1] = (p->sgn) ? -n1 : n1;
return;
}
}
/* Store as a type-2 non-'integer', i.e. put z in the form z1_rnd * 10^n + z2_rnd * 10^(n-9).
* Could add checks for zero/infinity here to avoid lengthy (300ish iterations) while loops below.
*/
n = 0;
while (1e16 <= z)
{
n += 5;
z *= 1e-5;
}
while (1e1 > z)
{
n -= 5;
z *= 1e5;
}
while (1e9 <= z)
{
n++;
z *= .1;
}
while (1e8 > z)
{
n--;
z *= 10.0;
}
z1_rnd = (int)z; /* casts should keep z1_rnd and z2_rnd in proper range, hence the assert below */
z -= (double)z1_rnd;
z *= 1e9;
z2_rnd = (int)z;
assert(((MANT_LO <= z1_rnd) && (MANT_HI > z1_rnd)) && ((0 <= z2_rnd) && (MANT_HI > z2_rnd)));
if ((0 == z2_rnd) && (-11 <= n) && (-3 >= n) && (0 == (z1_rnd % ten_pwr[-3 - n])))
{ /* This is a second-chance at detection an integer in case something slipped through the earlier
* check. Not expecting it to ever be invoked but is here as a safety net.
*/
z1_rnd /= ten_pwr[-3-n];
((mval_b *)p)->sgne = 0;
p->mvtype = (MV_NM | MV_INT);
p->m[1] = z1_rnd;
return;
}
exponent = MV_XBIAS + n + 9;
if (exponent >= EXPHI)
{
rts_error(VARLSTCNT(1) ERR_NUMOFLOW);
return;
}
if (exponent < EXPLO)
{
*p = literal_zero;
return;
}
p->mvtype = MV_NM;
p->e = exponent;
p->m[1] = z1_rnd;
p->m[0] = z2_rnd;
return;
}
Reference in New Issue View Git Blame Copy Permalink