#include "MyStdLib.h"
#include "MyString.h"
#include "math.h"
#include "fp_private.h"
#include "fenv.h"
#include "bcd.h"
#define INTPRECISION ((roundpre)(3))
#define SIGSIZEDIV 5 // deliver 160 good bits
#define SIGDIVQOUT SIGSIZE+SIGSIZEDIV
struct parts { int32_t highsig, lowsig;
};
typedef struct parts parts;
typedef union {
double extval;
struct parts intval;
} eparts;
typedef union {
uint32_t along;
float flt;
} fparts;
struct prod { uint32_t high;
uint32_t low;
};
typedef struct prod prod;
static void excessDigits (decimal *d, char rndChar, int excess, int32_t *expadjust);
int32_t apml ( double x );
double tenpower (const int32_t n);
double dec2d ( const decimal *d );
double big2d ( int sgn, big *b );
float big2float ( int sgn, big *b );
long big2long ( int sgn, big *b );
short big2short ( int sgn, big *b );
void dec2big ( const decimal *d, big *b );
void getnan ( double x, decimal *d );
void bigb2d ( double x, int32_t scale, decimal *d );
void ext2big ( double x, big *b );
uint32_t bigshifter ( big *a, big *d, int n);
void addbigs ( big *a, big *b, big *c );
void subbigs ( big *a, big *b, big *c );
void getsig ( double x, decimal *d );
int divuL (uint32_t dvsr, uint32_t *dvdnd);
struct prod muluL (uint32_t a, uint32_t b);
int subDL (prod a, uint32_t *b);
void subL (uint32_t *a, uint32_t *b, int length); int addL (uint32_t *a, uint32_t *b, int length); int firstdiff (uint32_t *a, uint32_t *b, int length);
void num2dec( const decform *f, double x, decimal *d ) { decform df;
int32_t logx, len, scale;
int32_t rnd;
fenv_t env;
double y, tens;
int dummy;
d->sgn = (char) __signbitd (x);
switch ( __fpclassifyd ( x ) ) {
#if 0
case FP_SNAN:
feraiseexcept (FE_INVALID);
case FP_QNAN:
getnan (x, d);
break;
#else
case FP_NAN: {
hexdouble z;
z.d = x;
if (!( z.i.hi & dQuietNan ))
feraiseexcept (FE_INVALID);
getnan (x, d);
break;
}
#endif
case FP_INFINITE:
d->sig.length = 1; d->sig.text [0] = 'I';
d->sig.text [1] = 0; break;
case FP_ZERO:
d->sig.length = 1; d->sig.text [0] = '0';
d->sig.text [1] = 0; break;
case FP_NORMAL:
case FP_SUBNORMAL:
rnd = fegetround (); dummy = feholdexcept (&env); df = *f;
logx = apml (x);
if ((df.style == FLOATDECIMAL) && (df.digits < 1)) df.digits = 1;
do {
if (df.style == FLOATDECIMAL) len = df.digits;
else len = df.digits + logx + 1;
if (len > MAXDIGITS) len = MAXDIGITS;
scale = len - logx - 1;
if (len < 1) len = 1;
feclearexcept (FE_INEXACT);
if (labs (scale) > MAXTEN) bigb2d (x, scale, d); else { tens = tenpower (scale);
dummy = fesetround (FE_TOWARDZERO); if (scale < 0) y = x / tens;
else if (scale > 0) y = x * tens;
else y = x;
dummy = fesetround (rnd); if (fetestexcept (FE_INEXACT) && (logb (y) > 51)) {
feclearexcept (FE_INEXACT); bigb2d (x, scale, d); }
else getsig (y, d); }
logx++;
} while (d->sig.length > len);
feupdateenv (&env);
d->exp = (short) -scale;
break;
}
}
double dec2num (const decimal *d) {
eparts s;
int i, i2, hex;
if ((d->sig.text [0] == '0') || (d->sig.length <= 0)) { s.intval.highsig = (d->sgn) ? 0x80000000 : 0 ;
s.intval.lowsig = 0;
return s.extval;
}
else if (d->sig.text [0] == 'I') { s.intval.highsig = (d->sgn) ? 0xfff00000 : 0x7ff00000 ;
s.intval.lowsig = 0;
return s.extval;
}
else if (d->sig.text [0] == 'N') { i2 = (d->sig.length < 5) ? d->sig.length + 3 : 8;
s.intval.highsig = 0;
for (i = 1; i <= i2; i++) {
s.intval.highsig <<= 4;
hex = (i < d->sig.length) ? d->sig.text [i] - '0' : 0 ;
if (hex > 9) hex += 9; hex &= 0x000f; s.intval.highsig |= hex;
}
if (!(s.intval.highsig &= 0x3fffffff)) s.intval.highsig = 0x00150000; s.intval.highsig |= 0x40000000;
s.intval.lowsig = s.intval.highsig << 23;
s.intval.highsig = 0x7ff00000 | (s.intval.highsig >> 11);
if (d->sgn) s.intval.highsig |= 0x80000000;
return s.extval;
}
else {
return dec2d ( d ); }
}
float dec2f ( const decimal *d ) {
int i, i2, hex;
fparts y;
big b;
if ((d->sig.text [0] == '0') || (d->sig.length <= 0)) { y.along = (d->sgn) ? 0x80000000 : 0 ;
return y.flt;
}
else if (d->sig.text [0] == 'I') { y.along = (d->sgn) ? 0xff800000 : 0x7f800000 ;
return y.flt;
}
else if (d->sig.text [0] == 'N') { i2 = (d->sig.length < 3) ? d->sig.length + 3 : 6;
y.along = 0;
for (i = 1; i <= i2; i++) {
y.along <<= 4;
hex = (i < d->sig.length) ? d->sig.text [i] - '0' : 0 ;
if (hex > 9) hex += 9; hex &= 0x000f; y.along |= hex;
}
if (!(y.along &= 0x003fffff)) y.along = 0x00001500;
y.along |= 0x7fc00000; if (d->sgn) y.along |= 0x80000000;
return y.flt;
}
else {
dec2big (d, &b);
return big2float ( (int) d->sgn, &b );
}
}
void dec2big ( const decimal *d, big *b ) {
big a;
int32_t i, i2, j, j2;
int sticky;
j2 = 0; a.sig.lng [j2] = 0; i2 = d->sig.length - 1; for (i = 0; i <= i2; i++) { a.sig.lng [0] += (d->sig.text [i] - '0'); if (i < i2) for (j = 0; j <= j2; j++) a.sig.lng [j] *= 10;
for (j = 0; j <= j2; j++) {
if (a.sig.lng [j] & 0xffff0000) { if (j == j2) { j2++;
a.sig.lng [j2] = 0;
}
a.sig.lng [j + 1] += ((a.sig.lng [j] >> 16) & 0x0000ffff); a.sig.lng [j] &= 0x0000ffff; }
}
}
for (i = 0; i < SIGSIZE; i++) b->sig.lng [i] = 0;
for (j = 0; j <= j2; j++) b->sig.shrt [j] = (short) a.sig.lng [j2 - j];
b->exp = 16*j2 + 15;
j2 /= 2;
while ((b->sig.lng [0] & 0x80000000) == 0) { b->exp--;
for (j = 0; j < j2; j++) {
b->sig.lng [j] <<= 1;
if ((b->sig.lng [j + 1] & 0x80000000) != 0)
b->sig.lng [j]++;
}
b->sig.lng [j2] <<= 1;
}
sticky = (b->sig.lng [1] & 0x000007ff) != 0;
i = 2;
while ((!sticky) && (i < SIGSIZE)) sticky = b->sig.lng [i++] != 0;
if (d->exp < 0) { bigtenpower (d->exp, &a);
adivb2c (b, &a, b);
}
else if (d->exp > 0) {
bigtenpower (d->exp, &a);
axb2c ( b, &a, b, FALSE );
}
}
short dec2s ( const decimal *d ) { big b;
if ((d->sig.text [0] == '0') || (d->sig.length <= 0)) { return 0;
}
else if ((d->sig.text [0] == 'I') || (d->sig.text [0] == 'N')) {
feraiseexcept (FE_INVALID);
return 0x8000;
}
else {
dec2big (d, &b);
return big2short ( (int) d->sgn, &b );
}
}
long dec2l ( const decimal *d ) { big b;
if ((d->sig.text [0] == '0') || (d->sig.length <= 0)) { return 0;
}
else if ((d->sig.text [0] == 'I') || (d->sig.text [0] == 'N')) {
feraiseexcept (FE_INVALID);
return 0x80000000;
}
else {
dec2big (d, &b);
return big2long ( (int) d->sgn, &b );
}
}
void dec2str(const decform *f,const decimal *d,char *s) {
unsigned char ExpDigs [9];
int i, i0, n = 0, exp, k;
ldiv_t qr;
n = strlen ((char*) d->sig.text);
if (n > SIGDIGLEN) {
s = strcpy (s, "NAN(017)"); return;
}
s = strncpy (s, (char*) d->sig.text, n);
s [n] = '\0';
exp = d->exp;
if (s [0] == '?') s = strcpy (s, "?");
else {
if (s [0] == '0') { s = strcpy (s, "0");
exp = 0; }
if (s [0] == 'I') s = strcpy (s, "INF");
else {
if (s [0] == 'N') {
i = strlen (s) - 1;
if (i > 4) i = 4;
i0 = i - 1;
if ( i0 < 1) i0 = 1;
qr.quot = 0;
while (i >= i0) {
qr.quot *= 16;
qr.quot += (s [i0] - '0');
if (s [i0] > '9') qr.quot += ('0' - 'A' + 10);
i0++;
};
s = strcpy (s, "NAN(XXX)");
for (i = 6; i >= 4; i--) {
qr = ldiv (qr.quot, 10);
s [i] = (unsigned char) ('0' + qr.rem);
}
}
else {
if (f->style == FLOATDECIMAL) { if (f->digits > (DECSTROUTLEN - 8)) n = DECSTROUTLEN - 8;
else {
if (f->digits < d->sig.length) n = d->sig.length;
else n = f->digits;
}
for (i = d->sig.length + 1; i <= n; i++) s = strcat (s, "0");
if (n > 1) {
k = strlen (s);
memmove (&s[2], &s[1], k);
s [1] = '.';
}
i = abs(exp + d->sig.length - 1); ExpDigs [0] = 0; do {
qr = ldiv (i, 10);
k = strlen ((char*) ExpDigs) + 1;
memmove ((char*) &ExpDigs [1], (char*) &ExpDigs [0], k);
ExpDigs [0] = (unsigned char) ('0' + qr.rem);
i = qr.quot;
} while (i > 0);
if ((exp + d->sig.length) < 1) s = strcat (s, "e-");
else s = strcat (s, "e+");
s = strcat (s, (char*) ExpDigs);
}
else { n = (f->digits < (-exp)) ? - exp : f->digits;
if (n < 1) {
if ((d->sig.length + exp + 1) <= DECSTROUTLEN) {
for (i = 1; i <= exp; i++) s = strcat(s, "0");
}
else s = strcpy (s, "?");
}
else {
if ((d->sig.length + exp) > 0) {
if ((d->sig.length + exp + n + 2) <= DECSTROUTLEN) {
for (i = 1; i <= (exp + n); i++) s = strcat(s, "0");
k = strlen (s);
i = k - n;
memmove (&s[i + 1], &s[i], n + 1);
s [i] = '.';
}
else s = strcpy (s, "?");
}
else {
if ((n + 3) <= DECSTROUTLEN) {
for (i = d->sig.length + 1; i <= (-exp); i++) {
k = strlen (s) + 1;
memmove (&s[1], &s[0], k);
s [0] = '0';
}
for (i = - exp + 1; i <= n; i++) s = strcat(s, "0");
k = strlen (s) + 1;
memmove (&s[2], &s[0], k);
s [0] = '0';
s [1] = '.';
}
else s = strcpy (s, "?");
}
}
} }
}
}
if (s [0] != '?') {
if ((d->sgn == 1) || (f->style == FLOATDECIMAL)) {
k = strlen (s) + 1;
memmove (&s[1], &s[0], k);
s [0] = (d->sgn) ? '-' : ' ';
}
}
}
void str2dec(const char *s,short *ix,decimal *d,short *vp) {
const char white[] = {' ', '\xca', '\t', '\0'}; const char zero[] = "0";
const char digits[] = "0123456789";
const char digitsdot[] = "0123456789.";
const char hex[] = "0123456789ABCDEF";
char rndChar = 0;
int i; int iok; int j;
int expsign, i1, i2, n1, n2;
int excess = 0; uint32_t nancode;
int32_t expadjust, expdigits;
*vp = d->sgn = d->exp = 0; d->sig.length = strlen (strcpy ((char*) d->sig.text, "N0011"));
iok = i = *ix; i += strspn (&s [i], white); if (s [i] == '+') i++; else if (s [i] == '-') { i++;
d->sgn = 1;
}
if ((j = strspn (&s [i], digitsdot))) { if ((j == 1) && (s [i] == '.')) { *vp = (short) (s [++i] == '\0'); return;
}
if (strspn (&s [i], ".") > 1) return; expadjust = i1 = i2 = n2 = 0;
i += strspn (&s [i], zero);
if ((n1 = strspn (&s [i], digits))) { i1 = i;
i += n1;
if (s [i] == '.') { i++;
if ((n2 = strspn (&s [i], digits))) {
i2 = i;
expadjust = -n2;
i += n2;
}
}
}
else { if (s [i] == '.') { i++;
j = strspn (&s [i], zero); i += j;
if ((n2 = strspn (&s [i], digits))) {
i2 = i; expadjust = -(j+n2);
i += n2;
}
}
}
if (n2)
while (s [i2 + n2 - 1] == '0') { n2--;
expadjust++;
}
if ((!n2) && n1) while (s [i1 + n1 - 1] == '0') { n1--;
expadjust++;
}
if ((j = n1 + n2) > SIGDIGLEN) {
rndChar = (n1 > SIGDIGLEN) ? s [i1 + SIGDIGLEN]: s [i2 + SIGDIGLEN - n1]; excess = j - SIGDIGLEN; n2 -= excess;
if (n2 < 0) {
n1 += n2;
n2 = 0;
}
expadjust += excess;
}
if (!(j = n1 + n2)) d->sig.length = strlen (strcpy ((char*) d->sig.text, "0"));
else if (j <= SIGDIGLEN) {
if (n1) strncpy ((char*) d->sig.text, &s [i1], n1);
if (n2) strncpy ((char*) &d->sig.text [n1], &s [i2], n2);
d->sig.text [j] = '\0';
if (excess) excessDigits (d, rndChar, excess, &expadjust);
if (expadjust > 32767) d->exp = 32767; else {
if (expadjust < -32767) d->exp = 0x8000; else d->exp = expadjust; }
}
else return; iok = i;
if ((s [i] == 'E') || (s [i] == 'e')){
i++;
expsign = 0;
if (s [i] == '+') i++; else if (s [i] == '-') {
i++;
expsign = 1;
}
expdigits = 0;
if ((j = strspn (&s [i], digits))) { while (j) {
if (expdigits > 0x0c000000L) expdigits = 0x0c000000L; expdigits *= 10;
expdigits += (s [i++] - '0');
j--;
}
if (expsign) expdigits = -expdigits; iok = i;
}
if (d->sig.text [0] != '0') {
expdigits += expadjust;
if (expdigits > 32767) d->exp = 32767; else {
if (expdigits < -32767) d->exp = 0x8000; else d->exp = expdigits; }
}
}
d->sig.length = strlen ((char*) d->sig.text);
}
else if ((s [i] == 'N') || (s [i] == 'n')) { i++;
if ((s [i] == 'A') || (s [i] == 'a')) {
i++;
if ((s [i] == 'N') || (s [i] == 'n')) {
iok = ++i;
d->sig.length = strlen (strcpy ((char*) d->sig.text, "N4000")); if (s [i] == '(') {
i++;
nancode = 0;
if ((j = strspn (&s [i], digits))) { while (j) {
nancode *= 10;
nancode += (s [i++] - '0');
j--;
}
}
if (s [i] == ')') {
iok = ++i;
for (j = 4; j >= 3; j--) {
d->sig.text [j] = hex [0x0000000f & nancode];
nancode >>= 4;
}
}
}
}
}
}
else if ((s [i] == 'I') || (s [i] == 'i')) { i++;
if ((s [i] == 'N') || (s [i] == 'n')) {
i++;
if ((s [i] == 'F') || (s [i] == 'f')) {
iok = ++i;
d->sig.length = strlen (strcpy ((char*) d->sig.text, "I")); }
}
}
else if (s [i] == '°') {
iok = ++i;
d->sig.length = strlen (strcpy ((char*) d->sig.text, "I")); }
*ix = (short) iok; *vp = (short) (s [i] == '\0'); }
void excessDigits (decimal *d, char rndChar, int excess, int32_t *expadjust) {
int i;
feraiseexcept (FE_INEXACT);
switch (fegetround ()) {
case FE_TONEAREST :
if ((excess > 1) || (rndChar != '5')) { if (rndChar < '5') return;
}
else { if (strchr ("13579", d->sig.text [SIGDIGLEN - 1]) == 0)
return; }
break;
case FE_UPWARD :
if (d->sgn) return;
break;
case FE_DOWNWARD :
if (!d->sgn) return;
break;
case FE_TOWARDZERO :
return;
break;
}
i = SIGDIGLEN - 1;
while (((++d->sig.text [i]) > '9') && (i >= 0)) { d->sig.text [i--] = '0'; }
if (i < 0) { d->sig.text [0] = '1'; (*expadjust)++; }
}
int32_t apml ( double x )
{
eparts s;
register int32_t e, g = 0x4D10;
s.extval = x;
if ((e = ((s.intval.highsig & 0x7FF00000) >> 20))) { s.intval.highsig = (((s.intval.highsig & 0x000FFFFF) | 0x00100000) << 11) |
((s.intval.lowsig >> 21) & 0x000007ff);
}
else {
e++;
s.intval.highsig = ((s.intval.highsig & 0x000FFFFF) << 11) |
((s.intval.lowsig >> 21) & 0x000007ff);
}
e -= 0x3FF; while (s.intval.highsig >= 0) { e--;
s.intval.highsig = s.intval.highsig << 1;
if (s.intval.lowsig < 0) s.intval.highsig++; s.intval.lowsig = s.intval.lowsig << 1;
}
if (e < 0) g++; return (g * e + ((g * ((s.intval.highsig & 0x7FFFFFFF) >> 15)) >> 16)) >> 16;
}
double tenpower (const int32_t n)
{
ldiv_t qr;
double x, y;
const double tens [] = {
1, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,
1e20, 1e21, 1e22};
uint32_t i = labs (n);
if (i <= MAXTEN) return tens [i];
qr = ldiv (i, MAXTEN);
i = qr.quot;
y = tens [qr.rem];
x = tens [MAXTEN];
do {
if (i & 1) y *= x;
if (i >>= 1) x *= x;
} while (i);
return y;
}
void bigtenpower (const int32_t n, big *y )
{
big xtable [] = {
{ 0x00000049, {{ 0x87867832, 0x6EAC9000, 0, 0, 0, 0, 0, 0 }} },
{ 0x00000092, {{ 0x8F7E32CE, 0x7BEA5C6F, 0xE4820023, 0xA2000000, 0, 0, 0, 0 }} },
{ 0x00000124, {{ 0xA0DC75F1, 0x778E39D6, 0x696361AE, 0x3DB1C721, 0x373B6B34, 0x318119EB,
0x65080000, 0 }} },
{ 0x00000248, {{ 0xCA28A291, 0x859BBF93, 0x7D7B8F75, 0x03CFDCFE, 0xD11F91FF, 0x10629770,
0x56291E60, 0x0B4D00DB }} },
{ 0x00000491, {{ 0x9FA42700, 0xDB900AD2, 0x5EBF18B6, 0xD27795FF, 0x9DF3E0BD, 0x5F019366,
0xF477189F, 0xA875F113 }} },
{ 0x00000922, {{ 0xC71AA36A, 0x1F8F01CB, 0x9DAD43F2, 0x30E1226E, 0x83689C3C, 0xBD362290,
0x50C00F72, 0x12E04C7D }} },
{ 0x00001245, {{ 0x9ADA6CD4, 0x96EF0E05, 0x2F1A208F, 0xDEDFF747, 0x57E155F4, 0xAE05D035,
0xE00E35AB, 0xDA0C5952 }} },
{ 0x0000248A, {{ 0xBB570A9A, 0x9BD977CC, 0x4C808753, 0xBB22FEF8, 0x6FC5802C, 0xDE0B3272,
0x10E980A1, 0xC0CCFD83 }} } };
int j;
ldiv_t qr;
uint32_t i;
fenv_t e;
int dummy;
dummy = feholdexcept (&e); dummy = fesetround (FE_TONEAREST);
i = labs (n);
if (i > 5631) i = 5631;
qr = ldiv (i, MAXTEN);
i = qr.quot;
ext2big (tenpower (qr.rem), y);
if (i == 0) {
feupdateenv (&e);
return;
}
j = 0;
do {
if (i & 1) axb2c ( y, &xtable [j], y, TRUE );
j++;
i >>= 1;
} while (i);
feupdateenv (&e);
}
void getsig ( double x, decimal *d )
{
eparts s;
register int e, i, carry, round, sticky;
unsigned char length, c[128];
s.extval = x;
c [length = 0] = 0;
if ((e = ((s.intval.highsig & 0x7FF00000) >> 20))) { s.intval.highsig = (((s.intval.highsig & 0x000FFFFF) | 0x00100000) << 11) |
((s.intval.lowsig >> 21) & 0x000007ff);
}
else {
e++;
s.intval.highsig = ((s.intval.highsig & 0x000FFFFF) << 11) |
((s.intval.lowsig >> 21) & 0x000007ff);
}
e -= 0x3FF; s.intval.lowsig = s.intval.lowsig << 11;
while ((e < 0) && ((s.intval.highsig & 0xC0000000) != 0)) {
if (s.intval.highsig & 1) s.intval.lowsig = 1;
s.intval.highsig = (s.intval.highsig >> 1) & 0x7FFFFFFF;
e++;
}
do {
if ((s.intval.highsig & 0x80000000) != 0) carry = 1;
else carry = 0;
if (--e < 0) { round = (s.intval.highsig & 0x40000000) != 0;
sticky = ((s.intval.highsig & 0x3FFFFFFF) != 0) ||
(s.intval.lowsig != 0) || fetestexcept (FE_INEXACT);
switch ( fegetround () ) {
case FE_TONEAREST:
if (round && (carry || sticky))
carry++;
break;
case FE_UPWARD:
if ((d->sgn == 0) && ((round != 0) || (sticky != 0)))
carry++;
break;
case FE_DOWNWARD:
if ((d->sgn == 1) && ((round != 0) || (sticky != 0)))
carry++;
break;
case FE_TOWARDZERO:
break;
}
if ((round != 0) || (sticky != 0)) feraiseexcept (FE_INEXACT);
}
for ( i = length; i >= 0; i-- ) { c [i] = c [i] + c [i] + carry;
carry = 0;
while (c [i] > 9) {
c [i] -= 10;
carry++;
}
}
if (carry != 0) { for ( i = ++length; i; i-- ) c [i] = c [i - 1];
c [0] = carry;
}
s.intval.highsig = s.intval.highsig << 1;
if ((s.intval.lowsig & 0x80000000) != 0) s.intval.highsig = s.intval.highsig + 1; s.intval.lowsig = s.intval.lowsig << 1;
} while ((e >= 0) && (length < MAXDIGITS)); for ( i = 0; i <= length; i++ ) c [i] += '0';
d->sig.length = length + 1;
for (i = 0; i <= length; i++) d->sig.text [i] = c [i];
d->sig.text [length + 1] = 0; return;
}
void biggetsig ( big *s, decimal *d )
{
register int i, carry, round, sticky;
unsigned char length, c[128];
c [length = 0] = 0;
while ((s->exp < 0) && ((s->sig.lng [0] & 0xC0000000) != 0)) { if (s->sig.lng [0] & 1) s->sig.lng [1] = 1; s->sig.lng [0] = (s->sig.lng [0] >> 1) & 0x7FFFFFFF;
s->exp++;
}
do { if ((s->sig.lng [0] & 0x80000000) != 0) carry = 1;
else carry = 0;
if (--s->exp < 0) { round = (s->sig.lng [0] & 0x40000000) != 0;
sticky = ((s->sig.lng [0] & 0x3FFFFFFF) != 0) || fetestexcept (FE_INEXACT);
i = 1;
while ((!sticky) && (i < SIGSIZE)) { sticky = s->sig.lng [i] != 0;
i++;
}
switch ( fegetround () ) {
case FE_TONEAREST:
if ((round != 0) && ((carry != 0) || (sticky != 0))) carry++;
break;
case FE_UPWARD:
if ((d->sgn == 0) && ((round != 0) || (sticky != 0)))
carry++;
break;
case FE_DOWNWARD:
if ((d->sgn == 1) && ((round != 0) || (sticky != 0)))
carry++;
break;
case FE_TOWARDZERO:
break;
}
if ((round != 0) || (sticky != 0)) feraiseexcept (FE_INEXACT);
}
for ( i = length; i >= 0; i-- ) { c [i] = c [i] + c [i] + carry;
carry = 0;
while (c [i] > 9) {
c [i] -= 10;
carry++;
}
}
if (carry != 0) { for ( i = ++length; i; i-- ) c [i] = c [i - 1];
c [0] = carry;
}
for (i = 0; i < SIGSIZEM; i++)
s->sig.lng [i] = (s->sig.lng [i] << 1) + ((s->sig.lng [i + 1] & 0x80000000) != 0);
s->sig.lng [i] = s->sig.lng [i] << 1;
} while ((s->exp >= 0) && (length <= MAXDIGITS)); for ( i = 0; i <= length; i++ ) c [i] += '0';
d->sig.length = length + 1;
for (i = 0; i <= length; i++) d->sig.text [i] = c [i];
d->sig.text [length + 1] = 0; return;
}
double dec2d ( const decimal *d) {
big a, b;
int32_t i, i2, j, j2;
eparts s;
int sticky;
int32_t bias = 1023;
j2 = 0; a.sig.lng [j2] = 0; i2 = d->sig.length - 1; for (i = 0; i <= i2; i++) { a.sig.lng [0] += (d->sig.text [i] - '0'); if (i < i2) for (j = 0; j <= j2; j++) a.sig.lng [j] *= 10;
for (j = 0; j <= j2; j++) {
if (a.sig.lng [j] & 0xffff0000) { if (j == j2) { j2++;
a.sig.lng [j2] = 0;
}
a.sig.lng [j + 1] += ((a.sig.lng [j] >> 16) & 0x0000ffff); a.sig.lng [j] &= 0x0000ffff; }
}
}
for (i = 0; i < SIGSIZE; i++) b.sig.lng [i] = 0;
for (j = 0; j <= j2; j++) b.sig.shrt [j] = (short) a.sig.lng [j2 - j];
b.exp = 16*j2 + 15;
j2 /= 2;
while ((b.sig.lng [0] & 0x80000000) == 0) { b.exp--;
for (j = 0; j < j2; j++) {
b.sig.lng [j] <<= 1;
if ((b.sig.lng [j + 1] & 0x80000000) != 0)
b.sig.lng [j]++;
}
b.sig.lng [j2] <<= 1;
}
sticky = (b.sig.lng [1] & 0x000007ff) != 0;
i = 2;
while ((!sticky) && (i < SIGSIZE)) sticky = b.sig.lng [i++] != 0;
if ((!sticky) && (abs(d->exp) <= MAXTEN)) { s.intval.highsig = ((b.exp + bias) << 20) | ((b.sig.lng [0] >> 11) & 0x000fffff);
if (d->sgn) s.intval.highsig |= 0x80000000; s.intval.lowsig = (b.sig.lng [0] << 21) | ((b.sig.lng [1] >> 11) & 0x001fffff);
if (d->exp < 0) s.extval /= tenpower (d->exp); if (d->exp > 0) s.extval *= tenpower (d->exp);
return s.extval;
}
else {
if (d->exp < 0) { bigtenpower (d->exp, &a);
adivb2c (&b, &a, &b);
}
else if (d->exp > 0) {
bigtenpower (d->exp, &a);
axb2c ( &b, &a, &b, FALSE );
}
return big2d ( (int) d->sgn, &b );
}
}
double big2d ( int sgn, big *b )
{
int32_t i, j;
eparts s;
int carry, round, sticky, bump;
int32_t maxexp = 1023;
int32_t bias = maxexp;
int32_t minexp = -1022;
int numbits = 54;
b->sig.lng [2] = (b->sig.lng [1] << 21) | (((b->sig.lng [2] >> 11) | b->sig.lng [2]) & 0x001fffff);
b->sig.lng [1] = (b->sig.lng [0] << 21) | ((b->sig.lng [1] >> 11) & 0x001fffff);
b->sig.lng [0] = (b->sig.lng [0] >> 11) & 0x001fffff;
if (b->exp < minexp) { j = numbits;
while ((j) && (b->exp < minexp)) { j--;
b->exp++;
bump = 0;
for (i = 0; i < SIGSIZE; i++) {
carry = (b->sig.lng [i] & 0x00000001) != 0;
b->sig.lng [i] = 0x7fffffff & (b->sig.lng [i] >> 1);
if (bump) b->sig.lng [i] |= 0x80000000;
bump = carry;
}
if (carry) b->sig.lng [SIGSIZEM] |= 0x00000001;
}
b->exp = minexp;
}
carry = (b->sig.lng [1] & 0x00000001) != 0;
round = (b->sig.lng [2] & 0x80000000) != 0;
sticky = (b->sig.lng [2] & 0x7fffffff) != 0;
i = 3;
while ((!sticky) && (i < SIGSIZE)) sticky = b->sig.lng [i++] != 0; bump = 0;
switch ( fegetround () ) {
case FE_TONEAREST:
if ((round != 0) && ((carry != 0) || (sticky != 0))) bump++;
break;
case FE_UPWARD:
if ((sgn == 0) && ((round != 0) || (sticky != 0)))
bump++;
break;
case FE_DOWNWARD:
if ((sgn == 1) && ((round != 0) || (sticky != 0)))
bump++;
break;
case FE_TOWARDZERO:
break;
}
if ((round != 0) || (sticky != 0)) feraiseexcept (FE_INEXACT);
if (bump) { if (!(++b->sig.lng [1]))
if ((++b->sig.lng [0]) == 0x00200000) {
b->sig.lng [0] = 0x00100000;
b->exp++;
}
}
if ((b->sig.lng [0] & 0xfff00000) == 0) {
b->exp--;
if ((round != 0) || (sticky != 0)) feraiseexcept (FE_UNDERFLOW);
}
if (b->exp <= maxexp) { s.intval.highsig = ((b->exp + bias) << 20) | (b->sig.lng [0] & 0x000fffff);
if (sgn) s.intval.highsig |= 0x80000000; s.intval.lowsig = b->sig.lng [1];
}
else { s.intval.lowsig = 0;
s.intval.highsig = (sgn) ? 0xfff00000 : 0x7ff00000 ;
feraiseexcept (FE_INEXACT + FE_OVERFLOW); if (((fegetround () == FE_UPWARD) && (sgn == 1)) ||
((fegetround () == FE_DOWNWARD) && (sgn == 0)) ||
(fegetround () == FE_TOWARDZERO)) {
s.intval.lowsig = 0xffffffff;
s.intval.highsig = (sgn) ? 0xffefffff : 0x7fefffff ;
}
}
return s.extval;
}
float big2float ( int sgn, big *b ) {
int32_t i, j;
fparts y;
int carry, round, sticky, bump;
int32_t maxexp = 127;
int32_t bias = maxexp;
int32_t minexp = -126;
int numbits = 25;
b->sig.lng [1] = (b->sig.lng [0] << 24) | (((b->sig.lng [1] >> 8) | b->sig.lng [1]) & 0x00ffffff);
b->sig.lng [0] = (b->sig.lng [0] >> 8) & 0x00ffffff;
if (b->exp < minexp) { j = numbits;
while ((j) && (b->exp < minexp)) { j--;
b->exp++;
bump = 0;
for (i = 0; i < SIGSIZE; i++) {
carry = (b->sig.lng [i] & 0x00000001) != 0;
b->sig.lng [i] = 0x7fffffff & (b->sig.lng [i] >> 1);
if (bump) b->sig.lng [i] |= 0x80000000;
bump = carry;
}
if (carry) b->sig.lng [SIGSIZEM] |= 0x00000001;
}
b->exp = minexp;
}
carry = (b->sig.lng [0] & 0x00000001) != 0;
round = (b->sig.lng [1] & 0x80000000) != 0;
sticky = (b->sig.lng [1] & 0x7fffffff) != 0;
i = 2;
while ((!sticky) && (i < SIGSIZE)) sticky = b->sig.lng [i++] != 0; bump = 0;
switch ( fegetround () ) {
case FE_TONEAREST:
if ((round != 0) && ((carry != 0) || (sticky != 0))) bump++;
break;
case FE_UPWARD:
if ((sgn == 0) && ((round != 0) || (sticky != 0)))
bump++;
break;
case FE_DOWNWARD:
if ((sgn == 1) && ((round != 0) || (sticky != 0)))
bump++;
break;
case FE_TOWARDZERO:
break;
}
if ((round != 0) || (sticky != 0)) feraiseexcept (FE_INEXACT);
if (bump) { if ((++b->sig.lng [0]) == 0x01000000) {
b->sig.lng [0] = 0x00800000;
b->exp++;
}
}
if ((b->sig.lng [0] & 0xff800000) == 0) {
b->exp--;
if ((round != 0) || (sticky != 0)) feraiseexcept (FE_UNDERFLOW);
}
if (b->exp <= maxexp) { y.along = ((b->exp + bias) << 23) | (b->sig.lng [0] & 0x007fffff);
if (sgn) y.along |= 0x80000000; }
else { y.along = (sgn) ? 0xfff00000 : 0x7ff00000 ;
feraiseexcept (FE_INEXACT + FE_OVERFLOW); if (((fegetround () == FE_UPWARD) && (sgn == 1)) ||
((fegetround () == FE_DOWNWARD) && (sgn == 0)) ||
(fegetround () == FE_TOWARDZERO)) {
y.along = (sgn) ? 0xffefffff : 0x7fefffff ;
}
}
return y.flt;
}
long big2long ( int sgn, big *b ) {
int32_t i, j;
int carry, round, sticky, bump;
int32_t maxexp = 31;
int32_t minexp = 31;
int numbits = 33;
if (b->exp < minexp) { j = numbits;
while ((j) && (b->exp < minexp)) { j--;
b->exp++;
bump = 0;
for (i = 0; i < SIGSIZE; i++) {
carry = (b->sig.lng [i] & 0x00000001) != 0;
b->sig.lng [i] = 0x7fffffff & (b->sig.lng [i] >> 1);
if (bump) b->sig.lng [i] |= 0x80000000;
bump = carry;
}
if (carry) b->sig.lng [SIGSIZEM] |= 0x00000001;
}
b->exp = minexp;
}
carry = (b->sig.lng [0] & 0x00000001) != 0;
round = (b->sig.lng [1] & 0x80000000) != 0;
sticky = (b->sig.lng [1] & 0x7fffffff) != 0;
i = 2;
while ((!sticky) && (i < SIGSIZE)) sticky = b->sig.lng [i++] != 0; bump = 0;
switch ( fegetround () ) {
case FE_TONEAREST:
if ((round != 0) && ((carry != 0) || (sticky != 0))) bump++;
break;
case FE_UPWARD:
if ((sgn == 0) && ((round != 0) || (sticky != 0)))
bump++;
break;
case FE_DOWNWARD:
if ((sgn == 1) && ((round != 0) || (sticky != 0)))
bump++;
break;
case FE_TOWARDZERO:
break;
}
if (bump) { if ((++b->sig.lng [0]) == 0x01000000) {
b->sig.lng [0] = 0x00800000;
b->exp++;
}
}
if (sgn) {
if ((b->exp > maxexp) ||
(((b->sig.lng [0] & 0x80000000) != 0) && (b->sig.lng [0] != 0x80000000))) {
feraiseexcept (FE_INVALID);
return 0x80000000;
}
else {
if ((round != 0) || (sticky != 0)) feraiseexcept (FE_INEXACT);
return -b->sig.lng [0];
}
}
else {
if ((b->exp > maxexp) || ((b->sig.lng [0] & 0x80000000) != 0)) {
feraiseexcept (FE_INVALID);
return 0x80000000;
}
else {
if ((round != 0) || (sticky != 0)) feraiseexcept (FE_INEXACT);
return b->sig.lng [0];
}
}
}
short big2short ( int sgn, big *b ) {
int32_t i, j;
int carry, round, sticky, bump;
int32_t maxexp = 31;
int32_t minexp = 31;
int numbits = 33;
if (b->exp < minexp) { j = numbits;
while ((j) && (b->exp < minexp)) { j--;
b->exp++;
bump = 0;
for (i = 0; i < SIGSIZE; i++) {
carry = (b->sig.lng [i] & 0x00000001) != 0;
b->sig.lng [i] = 0x7fffffff & (b->sig.lng [i] >> 1);
if (bump) b->sig.lng [i] |= 0x80000000;
bump = carry;
}
if (carry) b->sig.lng [SIGSIZEM] |= 0x00000001;
}
b->exp = minexp;
}
carry = (b->sig.lng [0] & 0x00000001) != 0;
round = (b->sig.lng [1] & 0x80000000) != 0;
sticky = (b->sig.lng [1] & 0x7fffffff) != 0;
i = 2;
while ((!sticky) && (i < SIGSIZE)) sticky = b->sig.lng [i++] != 0; bump = 0;
switch ( fegetround () ) {
case FE_TONEAREST:
if ((round != 0) && ((carry != 0) || (sticky != 0))) bump++;
break;
case FE_UPWARD:
if ((sgn == 0) && ((round != 0) || (sticky != 0)))
bump++;
break;
case FE_DOWNWARD:
if ((sgn == 1) && ((round != 0) || (sticky != 0)))
bump++;
break;
case FE_TOWARDZERO:
break;
}
if (bump) { if ((++b->sig.lng [0]) == 0x01000000) {
b->sig.lng [0] = 0x00800000;
b->exp++;
}
}
if (sgn) {
if ((b->exp > maxexp) || (((b->sig.lng [0] & 0xffff8000) != 0) && (b->sig.lng [0] != 0x00008000))) {
feraiseexcept (FE_INVALID);
return 0x8000;
}
else {
if ((round != 0) || (sticky != 0)) feraiseexcept (FE_INEXACT);
return -b->sig.shrt [1];
}
}
else {
if ((b->exp > maxexp) || ((b->sig.lng [0] & 0xffff8000) != 0)) {
feraiseexcept (FE_INVALID);
return 0x8000;
}
else {
if ((round != 0) || (sticky != 0)) feraiseexcept (FE_INEXACT);
return b->sig.shrt [1];
}
}
}
int big2twobigs ( big *b, big *c ) {
int i, j, j1, j2, kleft, kright, sgn, carry;
int32_t mask0;
c->exp = b->exp;
if ((sgn = (b->sig.lng [1] & 0x00000400) != 0)) { carry = 1;
for (i = SIGSIZEM; i > 1; i--)
carry = ((c->sig.lng [i-1] = ~b->sig.lng [i] + carry) == 0) && carry;
c->sig.lng [0] = ~b->sig.lng [1] + carry;
b->sig.lng [1] &= 0xfffff800;
if ((b->sig.lng [1] += 0x00000800) == 0)
if ((++b->sig.lng [0]) == 0) {
b->sig.lng [0] = 0x80000000;
b->exp++;
}
}
else { for (i = SIGSIZEM; i > 1; i--) c->sig.lng [i-1] = b->sig.lng [i];
c->sig.lng [0] = b->sig.lng [1];
b->sig.lng [1] &= 0xfffff800;
}
for (i = 2; i < SIGSIZE; i++) b->sig.lng [i] = 0;
c->sig.lng [SIGSIZEM] = 0;
if (c->sig.lng [0] &= 0x000007ff) {
mask0 = 0x00000400;
kleft = 21;
j1 = 0;
while ((mask0 & c->sig.lng [j1]) == 0) {
kleft++;
mask0 = mask0 >> 1;
}
}
else {
j1 = 1;
while ((j1 < SIGSIZE) && (c->sig.lng [j1] == 0)) j1++;
if (j1 == SIGSIZE) { c->exp = 0;
c->sig.lng [0] = 0;
return 0; }
else {
mask0 = 0x80000000;
kleft = 0;
while ((mask0 & c->sig.lng [j1]) == 0) {
kleft++;
mask0 = mask0 >> 1;
}
}
}
c->exp -= (32*(j1 + 1) + kleft);
kright = 32 - kleft;
j2 = SIGSIZEM - j1 - 1;
for (j = 0; j <= j2; j++) {
c->sig.lng [j] = (c->sig.lng [j + j1] << kleft) | (c->sig.lng [j + j1 + 1] >> kright);
}
for (j = j2; j < SIGSIZEM; j++) c->sig.lng [j] = 0;
return sgn;
}
void getnan ( double x, decimal *d )
{
eparts s;
register int32_t n;
s.extval = x;
d->sig.length = 5; d->sig.text [0] = 'N';
d->sig.text [1] = (unsigned char) ('4' + ((s.intval.highsig & 0x00020000) >> 17));
n = (s.intval.highsig & 0x0001E000) >> 13;
if (n > 9) n = n + 7;
d->sig.text [2] = (unsigned char) ('0' + n);
n = (s.intval.highsig & 0x00001E00) >> 9;
if (n > 9) n = n + 7;
d->sig.text [3] = (unsigned char) ('0' + n);
n = (s.intval.highsig & 0x000001E0) >> 5;
if (n > 9) n = n + 7;
d->sig.text [4] = (unsigned char) ('0' + n);
d->sig.text [5] = 0; }
void bigb2d ( double x, int32_t scale, decimal *d ) {
big y, z;
ext2big (x, &y);
if (scale < 0) {
bigtenpower (scale, &z);
adivb2c (&y, &z, &y);
}
else if (scale > 0) {
bigtenpower (scale, &z);
axb2c ( &y, &z, &y, FALSE );
}
biggetsig (&y, d);
}
void ext2big ( double x, big *b )
{
eparts s;
register int i, j;
uint32_t carry;
s.extval = x;
if ((b->exp = ((s.intval.highsig & 0x7FF00000) >> 20))) { b->sig.lng [0] = (((s.intval.highsig & 0x000FFFFF) | 0x00100000) << 11) |
((s.intval.lowsig >> 21) & 0x000007ff);
}
else {
b->exp++;
b->sig.lng [0] = ((s.intval.highsig & 0x000FFFFF) << 11) |
((s.intval.lowsig >> 21) & 0x000007ff);
}
b->exp -= 0x3FF; b->sig.lng [1] = s.intval.lowsig << 11;
for (i = 2; i < SIGSIZE; i++) b->sig.lng [i] = 0;
j = 32*SIGSIZE;
while ((j-- > 0) && ((b->sig.lng [0] & 0x80000000) == 0)) {
for (i = 0; i < SIGSIZE; i++) {
carry = (i < SIGSIZEM) ? (b->sig.lng [i + 1] & 0x80000000) != 0 : 0;
b->sig.lng [i] = (b->sig.lng [i] << 1) + carry;
}
b->exp--;
}
}
void axb2c ( big *a, big *b, big *c, int finishRounding )
{
int i, j, k, i2, j2, k2, carry, round, sticky;
uint32_t temp, product [SIGSIZE4M];
c->exp = a->exp + b->exp + 1;
i2 = SIGSIZE2M;
while (a->sig.shrt [i2] == 0) i2--; j2 = SIGSIZE2M;
while (b->sig.shrt [j2] == 0) j2--; k2 = i2 + j2;
for (k = 0; k < SIGSIZE4M; k++) product [k] = 0; for (i = 0; i <= i2; i++) { for (j = 0; j <= j2; j++) {
k = i + j;
temp = a->sig.shrt [i] * b->sig.shrt [j];
if (k != 0) {
product [k - 1] += ((temp >> 16) & 0x0000FFFFu);
product [k] += (temp & 0x0000FFFFu);
}
else product [k] = temp; }
}
for (k = k2; k > 0; k--) { product [k - 1] += (product [k] >> 16);
product [k] &= 0x0000FFFFu;
}
while ((product [0] & 0x80000000u) == 0) { c->exp--;
for (k = 0; k < k2; k++) {
product [k] <<= 1;
if ((product [k + 1] & 0x00008000u) != 0)
product [k]++;
}
product [k2] <<= 1;
}
c->sig.lng [0] = product [0];
carry = (product [SIGSIZE2M2] & 0x00000001u) != 0;
round = (product [SIGSIZE2M ] & 0x00008000u) != 0;
sticky = (product [SIGSIZE2M] & 0x00007FFFu) != 0;
k = SIGSIZE2;
while ((!sticky) && (k <= k2)) { sticky = (product [k] & 0x00007FFFu) != 0;
k++;
}
if (finishRounding) {
if ((round != 0) || (sticky != 0)) feraiseexcept (FE_INEXACT);
if ((round != 0) && ((carry != 0) || (sticky != 0))) { k = SIGSIZE2M2;
product [k]++;
while (((product [k] & 0x0000FFFFu) == 0) && (k > 0)) { k--;
product [k]++;
}
if ((k == 0) && (product [0] == 0)) { product [0] = 0x80000000u; c->exp++; }
}
}
else if ((round != 0) || (sticky != 0)) { feraiseexcept (FE_INEXACT);
product [k] |= 1; } c->sig.lng [0] = product [0]; for (k = 2; k <= SIGSIZE2M; k++)
c->sig.shrt [k] = (short) product [k - 1] & 0x0000FFFFu;
}
void adivb2c ( big *a, big *b, big *c )
{
int i, j, k, carry, sticky, excess = FALSE;
int asize = SIGSIZEM, bsize = SIGSIZEM, csize, dqsize;
uint32_t dq [SIGDIVQOUT];
struct prod prodresult;
c->exp = a->exp - b->exp - 1;
while (a->sig.lng [asize] == 0) asize--; while (b->sig.lng [bsize] == 0) bsize--; for (i = 0; i <= asize; i++) dq [i] = a->sig.lng [i];
dq [dqsize = asize + 1] = 0; while (dqsize < (SIGSIZEDIV + bsize)) dq [++dqsize] = 0;
csize = (asize > bsize) ? asize : bsize; i = firstdiff (&(a->sig.lng [0]), &(b->sig.lng [0]), csize); if (a->sig.lng [i] >= b->sig.lng [i]) { c->exp++;
carry = FALSE;
for (i = dqsize; i > 0; i--) {
dq [i] >>= 1;
if ((dq [i - 1] & 0x00000001u) != 0) dq [i] |= 0x80000000u;
else dq [i] &= 0x7fffffffu;
}
dq [0] >>= 1;
dq [0] &= 0x7fffffffu;
}
for (i = 0; i < SIGSIZEDIV; i++) {
if (divuL (b->sig.lng [0], &dq [i])) {
dq [i] = 0;
subL ( &(b->sig.lng [1]), &dq [i + 1], bsize-1);
do {
dq [i]--;
k = addL ( &(b->sig.lng [0]), &dq [i + 1], bsize);
k = firstdiff (&(b->sig.lng [0]), &dq [i + 1], bsize);
} while (b->sig.lng [k] <= dq [i + 1 + k]);
}
else {
excess = FALSE;
for (j = bsize; j > 0; j--) {
prodresult = muluL (b->sig.lng [j], dq [i]);
k = i + j;
if (subDL (prodresult, &dq [k])) {
do {
if (--k > i) dq [k]--;
else excess = TRUE;
}
while ((k > i) && (dq [k] == 0xffffffffu));
}
}
if (excess) {
do dq [i]--;
while (!addL ( &(b->sig.lng [0]), &dq [i + 1], bsize));
}
}
}
i = SIGSIZEDIV;
sticky = dq [i] != 0;
while ((!sticky) && (i < dqsize)) sticky = dq [++i] != 0;
for (i = 0; i < SIGSIZEDIV; i++) c->sig.lng [i] = dq [i];
c->sig.lng [SIGSIZEDIV] = sticky;
for (i = SIGSIZEDIV + 1; i < SIGSIZE; i++) c->sig.lng [i] = 0;
}
uint32_t bigshifter ( big *a, big *d, int n) {
ldiv_t qr;
uint32_t buf [SIGSIZE + 2];
int j, j1, j2, kleft;
for (j = 0; j < SIGSIZEP2; j++) buf [j] = 0; qr = ldiv (n, 32);
j2 = j1 = (qr.quot < SIGSIZEP2) ? qr.quot : SIGSIZEP2; for (j = 0; j < SIGSIZE; j++) {
if (j1 < SIGSIZEP2) buf [j1++] = a->sig.lng [j];
else buf [SIGSIZEP2 - 1] |= a->sig.lng [j];
}
if (qr.rem != 0) {
kleft = 32 - qr.rem;
buf [SIGSIZE + 1] |= buf [SIGSIZE] << kleft;
for (j = SIGSIZE; j > j2; j--)
buf [j] = (buf [j - 1] << kleft) | (buf [j] >> qr.rem);
if (j2 <= SIGSIZE) buf [j2] = buf [j2] >> qr.rem;
}
for (j = 0; j < SIGSIZE; j++) d->sig.lng [j] = buf [j];
return buf [SIGSIZE] | (buf [SIGSIZE + 1] != 0);
}
void addbigs ( big *a, big *b, big *c ) {
big tbig;
uint32_t temp, ovrflw, round, carry, sticky;
int i;
if (a->exp >= b->exp) {
c->exp = a->exp;
temp = bigshifter (b, &tbig, a->exp - b->exp);
ovrflw = addL ( &(a->sig.lng [0]), &(tbig.sig.lng [0]), SIGSIZEM);
}
else {
c->exp = b->exp;
temp = bigshifter (a, &tbig, b->exp - a->exp);
ovrflw = addL ( &(b->sig.lng [0]), &(tbig.sig.lng [0]), SIGSIZEM);
}
if (ovrflw != 0) { c->exp++;
sticky = temp != 0;
temp = bigshifter (&tbig, &tbig, 1) | sticky;
tbig.sig.lng [0] |= 0x80000000;
}
carry = (tbig.sig.lng [SIGSIZEM] & 1) != 0;
round = (temp & 0x80000000) != 0;
sticky = (temp & 0x7FFFFFFF) != 0;
if (round || sticky) feraiseexcept (FE_INEXACT);
if (round && (carry || sticky)) { i = SIGSIZEM;
while ((++tbig.sig.lng [i] == 0) && (i > 0)) i--;
if (i == 0)
if (++tbig.sig.lng [i] == 0) {
tbig.sig.lng [i] = 0x80000000;
c->exp++;
}
}
for (i = 0; i < SIGSIZE; i++) c->sig.lng [i] = tbig.sig.lng [i];
}
void subbigs ( big *a, big *b, big *c ) {
big tbig;
uint32_t temp, round, carry, sticky;
int i, j;
c->exp = a->exp;
temp = bigshifter (b, &tbig, a->exp - b->exp);
carry = (tbig.sig.lng [SIGSIZEM] & 1) != 0;
round = (temp & 0x80000000) != 0;
sticky = (temp & 0x7FFFFFFF) != 0;
if (round || sticky) feraiseexcept (FE_INEXACT);
if (round && (carry || sticky)) { i = SIGSIZEM;
while ((++tbig.sig.lng [i] == 0) && (i > 0)) i--;
if (i == 0) ++tbig.sig.lng [i];
}
for (i = 0; i < SIGSIZE; i++) c->sig.lng [i] = a->sig.lng [i];
subL ( &(tbig.sig.lng [0]), &(c->sig.lng [0]), SIGSIZEM);
j = 32*SIGSIZE;
while ((j-- > 0) && ((c->sig.lng [0] & 0x80000000) == 0)) {
for (i = 0; i < SIGSIZE; i++) {
carry = (i < SIGSIZEM) ? (c->sig.lng [i + 1] & 0x80000000) != 0 : 0;
c->sig.lng [i] = (c->sig.lng [i] << 1) + carry;
}
c->exp--;
}
}
int divuL (uint32_t dvsr, uint32_t *dvdnd)
{
register uint32_t quo = 0;
register int i;
if (dvsr <= dvdnd [0]) return 1; else {
for (i = 1; i <= 32; i++) {
quo <<= 1;
if ((dvdnd [0] & 0x80000000u) != 0) {
quo += 1;
dvdnd [0] <<= 1;
if ((dvdnd [1] & 0x80000000u) != 0) dvdnd [0]++;
dvdnd [0] -= dvsr;
}
else {
dvdnd [0] <<= 1;
if ((dvdnd [1] & 0x80000000u) != 0) dvdnd [0]++;
if (dvsr <= dvdnd [0]) {
quo += 1;
dvdnd [0] -= dvsr;
}
}
dvdnd [1] <<= 1;
}
dvdnd [1] = dvdnd [0];
dvdnd [0] = quo;
return 0; }
}
struct prod muluL (uint32_t a, uint32_t b)
{
union {
uint32_t along;
uint16_t ashort [2];
} u1, u2, p00, p01, p10, p11, r;
struct prod temp;
u1.along = a;
u2.along = b;
p00.along = (uint32_t) u1.ashort [0] * (uint32_t) u2.ashort [0];
p01.along = (uint32_t) u1.ashort [0] * (uint32_t) u2.ashort [1];
p10.along = (uint32_t) u1.ashort [1] * (uint32_t) u2.ashort [0];
p11.along = (uint32_t) u1.ashort [1] * (uint32_t) u2.ashort [1];
r.ashort [1] = p11.ashort [1];
u1.along = (uint32_t) p11.ashort [0]
+ (uint32_t) p01.ashort [1]
+ (uint32_t) p10.ashort [1];
r.ashort [0] = u1.ashort [1];
temp.low = r.along;
u2.along = (uint32_t) p00.ashort [1]
+ (uint32_t) p01.ashort [0]
+ (uint32_t) p10.ashort [0]
+ (uint32_t) u1.ashort [0];
r.ashort [1] = u2.ashort [1];
r.ashort [0] = p00.ashort [0] + u2.ashort [0];
temp.high = r.along;
return temp;
}
int subDL (prod a, uint32_t *b)
{
register int borrow, carry;
borrow = a.low > b [1];
b [1] -= a.low;
carry = (a.high > b [0]) || ((a.high == b [0]) && borrow);
b [0] -= (a.high + (uint32_t) borrow);
return carry;
}
void subL (uint32_t *a, uint32_t *b, int length) {
register int k, borrow = FALSE, carry;
for (k = length; k >= 0; k--) {
carry = (a [k] > b [k]) || ((a [k] == b [k]) && borrow);
b [k] -= (a [k] + (uint32_t) borrow);
borrow = carry;
}
}
int addL (uint32_t *a, uint32_t *b, int length) {
register int k, carry = FALSE;
register uint32_t c;
for (k = length; k >= 0; k--) {
c = b [k] + a [k] + (uint32_t) carry;
carry = (((0x80000000u & a [k]) != 0) && ((0x80000000u & b [k]) != 0))
|| (((0x80000000u & a [k]) != 0) && ((0x80000000u & c) == 0))
|| (((0x80000000u & b [k]) != 0) && ((0x80000000u & c) == 0));
b [k] = c;
}
return carry;
}
int firstdiff (uint32_t *a, uint32_t *b, int length)
{ register int i = 0;
while ((a [i] == b [i]) && (i < length)) i++;
return i;
}