diff options
Diffstat (limited to 'winsup/mingw/mingwex/math/cephes_emath.c')
-rw-r--r-- | winsup/mingw/mingwex/math/cephes_emath.c | 1318 |
1 files changed, 1318 insertions, 0 deletions
diff --git a/winsup/mingw/mingwex/math/cephes_emath.c b/winsup/mingw/mingwex/math/cephes_emath.c new file mode 100644 index 00000000000..ab798a2d255 --- /dev/null +++ b/winsup/mingw/mingwex/math/cephes_emath.c @@ -0,0 +1,1318 @@ +/* This file is extracted from S L Moshier's ioldoubl.c, + * modified for use in MinGW + * + * Extended precision arithmetic functions for long double I/O. + * This program has been placed in the public domain. + */ + + + +/* + * Revision history: + * + * 5 Jan 84 PDP-11 assembly language version + * 6 Dec 86 C language version + * 30 Aug 88 100 digit version, improved rounding + * 15 May 92 80-bit long double support + * + * Author: S. L. Moshier. + * + * 6 Oct 02 Modified for MinGW by inlining utility routines, + * removing global variables and splitting out strtold + * from _IO_ldtoa and _IO_ldtostr. + * + * Danny Smith <dannysmith@users.sourceforge.net> + * + */ + + +#include "cephes_emath.h" + +/* + * The constants are for 64 bit precision. + */ + + +/* Move in external format number, + * converting it to internal format. + */ +void __emovi(const short unsigned int * __restrict__ a, + short unsigned int * __restrict__ b) +{ +register const unsigned short *p; +register unsigned short *q; +int i; + +q = b; +p = a + (NE-1); /* point to last word of external number */ +/* get the sign bit */ +if( *p & 0x8000 ) + *q++ = 0xffff; +else + *q++ = 0; +/* get the exponent */ +*q = *p--; +*q++ &= 0x7fff; /* delete the sign bit */ +#ifdef INFINITY +if( (*(q-1) & 0x7fff) == 0x7fff ) + { +#ifdef NANS + if( __eisnan(a) ) + { + *q++ = 0; + for( i=3; i<NI; i++ ) + *q++ = *p--; + return; + } +#endif + for( i=2; i<NI; i++ ) + *q++ = 0; + return; + } +#endif +/* clear high guard word */ +*q++ = 0; +/* move in the significand */ +for( i=0; i<NE-1; i++ ) + *q++ = *p--; +/* clear low guard word */ +*q = 0; +} + + +/* +; Add significands +; x + y replaces y +*/ + +void __eaddm(const short unsigned int * __restrict__ x, + short unsigned int * __restrict__ y) +{ +register unsigned long a; +int i; +unsigned int carry; + +x += NI-1; +y += NI-1; +carry = 0; +for( i=M; i<NI; i++ ) + { + a = (unsigned long )(*x) + (unsigned long )(*y) + carry; + if( a & 0x10000 ) + carry = 1; + else + carry = 0; + *y = (unsigned short )a; + --x; + --y; + } +} + +/* +; Subtract significands +; y - x replaces y +*/ + +void __esubm(const short unsigned int * __restrict__ x, + short unsigned int * __restrict__ y) +{ +unsigned long a; +int i; +unsigned int carry; + +x += NI-1; +y += NI-1; +carry = 0; +for( i=M; i<NI; i++ ) + { + a = (unsigned long )(*y) - (unsigned long )(*x) - carry; + if( a & 0x10000 ) + carry = 1; + else + carry = 0; + *y = (unsigned short )a; + --x; + --y; + } +} + + +/* Multiply significand of e-type number b +by 16-bit quantity a, e-type result to c. */ + +static void __m16m(short unsigned int a, + short unsigned int * __restrict__ b, + short unsigned int * __restrict__ c) +{ +register unsigned short *pp; +register unsigned long carry; +unsigned short *ps; +unsigned short p[NI]; +unsigned long aa, m; +int i; + +aa = a; +pp = &p[NI-2]; +*pp++ = 0; +*pp = 0; +ps = &b[NI-1]; + +for( i=M+1; i<NI; i++ ) + { + if( *ps == 0 ) + { + --ps; + --pp; + *(pp-1) = 0; + } + else + { + m = (unsigned long) aa * *ps--; + carry = (m & 0xffff) + *pp; + *pp-- = (unsigned short )carry; + carry = (carry >> 16) + (m >> 16) + *pp; + *pp = (unsigned short )carry; + *(pp-1) = carry >> 16; + } + } +for( i=M; i<NI; i++ ) + c[i] = p[i]; +} + + +/* Divide significands. Neither the numerator nor the denominator +is permitted to have its high guard word nonzero. */ + + +int __edivm(short unsigned int * __restrict__ den, + short unsigned int * __restrict__ num) +{ +int i; +register unsigned short *p; +unsigned long tnum; +unsigned short j, tdenm, tquot; +unsigned short tprod[NI+1]; +unsigned short equot[NI]; + +p = &equot[0]; +*p++ = num[0]; +*p++ = num[1]; + +for( i=M; i<NI; i++ ) + { + *p++ = 0; + } +__eshdn1( num ); +tdenm = den[M+1]; +for( i=M; i<NI; i++ ) + { + /* Find trial quotient digit (the radix is 65536). */ + tnum = (((unsigned long) num[M]) << 16) + num[M+1]; + + /* Do not execute the divide instruction if it will overflow. */ + if( (tdenm * 0xffffUL) < tnum ) + tquot = 0xffff; + else + tquot = tnum / tdenm; + + /* Prove that the divide worked. */ +/* + tcheck = (unsigned long )tquot * tdenm; + if( tnum - tcheck > tdenm ) + tquot = 0xffff; +*/ + /* Multiply denominator by trial quotient digit. */ + __m16m( tquot, den, tprod ); + /* The quotient digit may have been overestimated. */ + if( __ecmpm( tprod, num ) > 0 ) + { + tquot -= 1; + __esubm( den, tprod ); + if( __ecmpm( tprod, num ) > 0 ) + { + tquot -= 1; + __esubm( den, tprod ); + } + } + __esubm( tprod, num ); + equot[i] = tquot; + __eshup6(num); + } +/* test for nonzero remainder after roundoff bit */ +p = &num[M]; +j = 0; +for( i=M; i<NI; i++ ) + { + j |= *p++; + } +if( j ) + j = 1; + +for( i=0; i<NI; i++ ) + num[i] = equot[i]; + +return( (int )j ); +} + + + +/* Multiply significands */ +int __emulm(const short unsigned int * __restrict__ a, + short unsigned int * __restrict__ b) +{ +const unsigned short *p; +unsigned short *q; +unsigned short pprod[NI]; +unsigned short equot[NI]; +unsigned short j; +int i; + +equot[0] = b[0]; +equot[1] = b[1]; +for( i=M; i<NI; i++ ) + equot[i] = 0; + +j = 0; +p = &a[NI-1]; +q = &equot[NI-1]; +for( i=M+1; i<NI; i++ ) + { + if( *p == 0 ) + { + --p; + } + else + { + __m16m( *p--, b, pprod ); + __eaddm(pprod, equot); + } + j |= *q; + __eshdn6(equot); + } + +for( i=0; i<NI; i++ ) + b[i] = equot[i]; + +/* return flag for lost nonzero bits */ +return( (int)j ); +} + + + +/* + * Normalize and round off. + * + * The internal format number to be rounded is "s". + * Input "lost" indicates whether the number is exact. + * This is the so-called sticky bit. + * + * Input "subflg" indicates whether the number was obtained + * by a subtraction operation. In that case if lost is nonzero + * then the number is slightly smaller than indicated. + * + * Input "exp" is the biased exponent, which may be negative. + * the exponent field of "s" is ignored but is replaced by + * "exp" as adjusted by normalization and rounding. + * + * Input "rcntrl" is the rounding control. + * + * Input "rnprc" is precison control (64 or NBITS). + */ + +void __emdnorm(short unsigned int *s, int lost, int subflg, long int exp, int rcntrl, int rndprc) +{ +int i, j; +unsigned short r; +int rw = NI-1; /* low guard word */ +int re = NI-2; +const unsigned short rmsk = 0xffff; +const unsigned short rmbit = 0x8000; +#if NE == 6 +unsigned short rbit[NI] = {0,0,0,0,0,0,0,1,0}; +#else +unsigned short rbit[NI] = {0,0,0,0,0,0,0,0,0,0,0,1,0}; +#endif + +/* Normalize */ +j = __enormlz( s ); + +/* a blank significand could mean either zero or infinity. */ +#ifndef INFINITY +if( j > NBITS ) + { + __ecleazs( s ); + return; + } +#endif +exp -= j; +#ifndef INFINITY +if( exp >= 32767L ) + goto overf; +#else +if( (j > NBITS) && (exp < 32767L) ) + { + __ecleazs( s ); + return; + } +#endif +if( exp < 0L ) + { + if( exp > (long )(-NBITS-1) ) + { + j = (int )exp; + i = __eshift( s, j ); + if( i ) + lost = 1; + } + else + { + __ecleazs( s ); + return; + } + } +/* Round off, unless told not to by rcntrl. */ +if( rcntrl == 0 ) + goto mdfin; +if (rndprc == 64) + { + rw = 7; + re = 6; + rbit[NI-2] = 0; + rbit[6] = 1; + } + +/* Shift down 1 temporarily if the data structure has an implied + * most significant bit and the number is denormal. + * For rndprc = 64 or NBITS, there is no implied bit. + * But Intel long double denormals lose one bit of significance even so. + */ +#if IBMPC +if( (exp <= 0) && (rndprc != NBITS) ) +#else +if( (exp <= 0) && (rndprc != 64) && (rndprc != NBITS) ) +#endif + { + lost |= s[NI-1] & 1; + __eshdn1(s); + } +/* Clear out all bits below the rounding bit, + * remembering in r if any were nonzero. + */ +r = s[rw] & rmsk; +if( rndprc < NBITS ) + { + i = rw + 1; + while( i < NI ) + { + if( s[i] ) + r |= 1; + s[i] = 0; + ++i; + } + } +s[rw] &= ~rmsk; +if( (r & rmbit) != 0 ) + { + if( r == rmbit ) + { + if( lost == 0 ) + { /* round to even */ + if( (s[re] & 1) == 0 ) + goto mddone; + } + else + { + if( subflg != 0 ) + goto mddone; + } + } + __eaddm( rbit, s ); + } +mddone: +#if IBMPC +if( (exp <= 0) && (rndprc != NBITS) ) +#else +if( (exp <= 0) && (rndprc != 64) && (rndprc != NBITS) ) +#endif + { + __eshup1(s); + } +if( s[2] != 0 ) + { /* overflow on roundoff */ + __eshdn1(s); + exp += 1; + } +mdfin: +s[NI-1] = 0; +if( exp >= 32767L ) + { +#ifndef INFINITY +overf: +#endif +#ifdef INFINITY + s[1] = 32767; + for( i=2; i<NI-1; i++ ) + s[i] = 0; +#else + s[1] = 32766; + s[2] = 0; + for( i=M+1; i<NI-1; i++ ) + s[i] = 0xffff; + s[NI-1] = 0; + if( (rndprc < 64) || (rndprc == 113) ) + s[rw] &= ~rmsk; +#endif + return; + } +if( exp < 0 ) + s[1] = 0; +else + s[1] = (unsigned short )exp; +} + + +/* +; Multiply. +; +; unsigned short a[NE], b[NE], c[NE]; +; emul( a, b, c ); c = b * a +*/ +void __emul(const short unsigned int *a, + const short unsigned int *b, + short unsigned int *c) +{ +unsigned short ai[NI], bi[NI]; +int i, j; +long lt, lta, ltb; + +#ifdef NANS +/* NaN times anything is the same NaN. */ +if( __eisnan(a) ) + { + __emov(a,c); + return; + } +if( __eisnan(b) ) + { + __emov(b,c); + return; + } +/* Zero times infinity is a NaN. */ +if( (__eisinf(a) && __eiiszero(b)) + || (__eisinf(b) && __eiiszero(a)) ) + { + mtherr( "emul", DOMAIN ); + __enan_NBITS( c ); + return; + } +#endif +/* Infinity times anything else is infinity. */ +#ifdef INFINITY +if( __eisinf(a) || __eisinf(b) ) + { + if( __eisneg(a) ^ __eisneg(b) ) + *(c+(NE-1)) = 0x8000; + else + *(c+(NE-1)) = 0; + __einfin(c); + return; + } +#endif +__emovi( a, ai ); +__emovi( b, bi ); +lta = ai[E]; +ltb = bi[E]; +if( ai[E] == 0 ) + { + for( i=1; i<NI-1; i++ ) + { + if( ai[i] != 0 ) + { + lta -= __enormlz( ai ); + goto mnzer1; + } + } + __eclear(c); + return; + } +mnzer1: + +if( bi[E] == 0 ) + { + for( i=1; i<NI-1; i++ ) + { + if( bi[i] != 0 ) + { + ltb -= __enormlz( bi ); + goto mnzer2; + } + } + __eclear(c); + return; + } +mnzer2: + +/* Multiply significands */ +j = __emulm( ai, bi ); +/* calculate exponent */ +lt = lta + ltb - (EXONE - 1); +__emdnorm( bi, j, 0, lt, 64, NBITS ); +/* calculate sign of product */ +if( ai[0] == bi[0] ) + bi[0] = 0; +else + bi[0] = 0xffff; +__emovo( bi, c ); +} + + +/* move out internal format to ieee long double */ +void __toe64(short unsigned int *a, short unsigned int *b) +{ +register unsigned short *p, *q; +unsigned short i; + +#ifdef NANS +if( __eiisnan(a) ) + { + __enan_64( b ); + return; + } +#endif +#ifdef IBMPC +/* Shift Intel denormal significand down 1. */ +if( a[E] == 0 ) + __eshdn1(a); +#endif +p = a; +#ifdef MIEEE +q = b; +#else +q = b + 4; /* point to output exponent */ +#if 1 +/* NOTE: if data type is 96 bits wide, clear the last word here. */ +*(q+1)= 0; +#endif +#endif + +/* combine sign and exponent */ +i = *p++; +#ifdef MIEEE +if( i ) + *q++ = *p++ | 0x8000; +else + *q++ = *p++; +*q++ = 0; +#else +if( i ) + *q-- = *p++ | 0x8000; +else + *q-- = *p++; +#endif +/* skip over guard word */ +++p; +/* move the significand */ +#ifdef MIEEE +for( i=0; i<4; i++ ) + *q++ = *p++; +#else +#ifdef INFINITY +if (__eiisinf (a)) + { + /* Intel long double infinity. */ + *q-- = 0x8000; + *q-- = 0; + *q-- = 0; + *q = 0; + return; + } +#endif +for( i=0; i<4; i++ ) + *q-- = *p++; +#endif +} + + +/* Compare two e type numbers. + * + * unsigned short a[NE], b[NE]; + * ecmp( a, b ); + * + * returns +1 if a > b + * 0 if a == b + * -1 if a < b + * -2 if either a or b is a NaN. + */ +int __ecmp(const short unsigned int * __restrict__ a, + const short unsigned int * __restrict__ b) +{ +unsigned short ai[NI], bi[NI]; +register unsigned short *p, *q; +register int i; +int msign; + +#ifdef NANS +if (__eisnan (a) || __eisnan (b)) + return( -2 ); +#endif +__emovi( a, ai ); +p = ai; +__emovi( b, bi ); +q = bi; + +if( *p != *q ) + { /* the signs are different */ +/* -0 equals + 0 */ + for( i=1; i<NI-1; i++ ) + { + if( ai[i] != 0 ) + goto nzro; + if( bi[i] != 0 ) + goto nzro; + } + return(0); +nzro: + if( *p == 0 ) + return( 1 ); + else + return( -1 ); + } +/* both are the same sign */ +if( *p == 0 ) + msign = 1; +else + msign = -1; +i = NI-1; +do + { + if( *p++ != *q++ ) + { + goto diff; + } + } +while( --i > 0 ); + +return(0); /* equality */ + + + +diff: + +if( *(--p) > *(--q) ) + return( msign ); /* p is bigger */ +else + return( -msign ); /* p is littler */ +} + +/* +; Shift significand +; +; Shifts significand area up or down by the number of bits +; given by the variable sc. +*/ +int __eshift(short unsigned int *x, int sc) +{ +unsigned short lost; +unsigned short *p; + +if( sc == 0 ) + return( 0 ); + +lost = 0; +p = x + NI-1; + +if( sc < 0 ) + { + sc = -sc; + while( sc >= 16 ) + { + lost |= *p; /* remember lost bits */ + __eshdn6(x); + sc -= 16; + } + + while( sc >= 8 ) + { + lost |= *p & 0xff; + __eshdn8(x); + sc -= 8; + } + + while( sc > 0 ) + { + lost |= *p & 1; + __eshdn1(x); + sc -= 1; + } + } +else + { + while( sc >= 16 ) + { + __eshup6(x); + sc -= 16; + } + + while( sc >= 8 ) + { + __eshup8(x); + sc -= 8; + } + + while( sc > 0 ) + { + __eshup1(x); + sc -= 1; + } + } +if( lost ) + lost = 1; +return( (int )lost ); +} + + + +/* +; normalize +; +; Shift normalizes the significand area pointed to by argument +; shift count (up = positive) is returned. +*/ +int __enormlz(short unsigned int *x) +{ +register unsigned short *p; +int sc; + +sc = 0; +p = &x[M]; +if( *p != 0 ) + goto normdn; +++p; +if( *p & 0x8000 ) + return( 0 ); /* already normalized */ +while( *p == 0 ) + { + __eshup6(x); + sc += 16; +/* With guard word, there are NBITS+16 bits available. + * return true if all are zero. + */ + if( sc > NBITS ) + return( sc ); + } +/* see if high byte is zero */ +while( (*p & 0xff00) == 0 ) + { + __eshup8(x); + sc += 8; + } +/* now shift 1 bit at a time */ +while( (*p & 0x8000) == 0) + { + __eshup1(x); + sc += 1; + if( sc > (NBITS+16) ) + { + mtherr( "enormlz", UNDERFLOW ); + return( sc ); + } + } +return( sc ); + +/* Normalize by shifting down out of the high guard word + of the significand */ +normdn: + +if( *p & 0xff00 ) + { + __eshdn8(x); + sc -= 8; + } +while( *p != 0 ) + { + __eshdn1(x); + sc -= 1; + + if( sc < -NBITS ) + { + mtherr( "enormlz", OVERFLOW ); + return( sc ); + } + } +return( sc ); +} + + +/* Move internal format number out, + * converting it to external format. + */ +void __emovo(const short unsigned int * __restrict__ a, + short unsigned int * __restrict__ b) +{ +register const unsigned short *p; +register unsigned short *q; +unsigned short i; + +p = a; +q = b + (NE-1); /* point to output exponent */ +/* combine sign and exponent */ +i = *p++; +if( i ) + *q-- = *p++ | 0x8000; +else + *q-- = *p++; +#ifdef INFINITY +if( *(p-1) == 0x7fff ) + { +#ifdef NANS + if( __eiisnan(a) ) + { + __enan_NBITS( b ); + return; + } +#endif + __einfin(b); + return; + } +#endif +/* skip over guard word */ +++p; +/* move the significand */ +for( i=0; i<NE-1; i++ ) + *q-- = *p++; +} + + +#if USE_LDTOA + + +void __eiremain(short unsigned int *den, short unsigned int *num, + short unsigned int *equot ) +{ +long ld, ln; +unsigned short j; + +ld = den[E]; +ld -= __enormlz( den ); +ln = num[E]; +ln -= __enormlz( num ); +__ecleaz( equot ); +while( ln >= ld ) + { + if( __ecmpm(den,num) <= 0 ) + { + __esubm(den, num); + j = 1; + } + else + { + j = 0; + } + __eshup1(equot); + equot[NI-1] |= j; + __eshup1(num); + ln -= 1; + } +__emdnorm( num, 0, 0, ln, 0, NBITS ); +} + + +void __eadd1(const short unsigned int * __restrict__ a, + const short unsigned int * __restrict__ b, + short unsigned int * __restrict__ c, + int subflg) +{ +unsigned short ai[NI], bi[NI], ci[NI]; +int i, lost, j, k; +long lt, lta, ltb; + +#ifdef INFINITY +if( __eisinf(a) ) + { + __emov(a,c); + if( subflg ) + __eneg(c); + return; + } +if( __eisinf(b) ) + { + __emov(b,c); + return; + } +#endif +__emovi( a, ai ); +__emovi( b, bi ); +if( sub ) + ai[0] = ~ai[0]; + +/* compare exponents */ +lta = ai[E]; +ltb = bi[E]; +lt = lta - ltb; +if( lt > 0L ) + { /* put the larger number in bi */ + __emovz( bi, ci ); + __emovz( ai, bi ); + __emovz( ci, ai ); + ltb = bi[E]; + lt = -lt; + } +lost = 0; +if( lt != 0L ) + { + if( lt < (long )(-NBITS-1) ) + goto done; /* answer same as larger addend */ + k = (int )lt; + lost = __eshift( ai, k ); /* shift the smaller number down */ + } +else + { +/* exponents were the same, so must compare significands */ + i = __ecmpm( ai, bi ); + if( i == 0 ) + { /* the numbers are identical in magnitude */ + /* if different signs, result is zero */ + if( ai[0] != bi[0] ) + { + __eclear(c); + return; + } + /* if same sign, result is double */ + /* double denomalized tiny number */ + if( (bi[E] == 0) && ((bi[3] & 0x8000) == 0) ) + { + __eshup1( bi ); + goto done; + } + /* add 1 to exponent unless both are zero! */ + for( j=1; j<NI-1; j++ ) + { + if( bi[j] != 0 ) + { +/* This could overflow, but let emovo take care of that. */ + ltb += 1; + break; + } + } + bi[E] = (unsigned short )ltb; + goto done; + } + if( i > 0 ) + { /* put the larger number in bi */ + __emovz( bi, ci ); + __emovz( ai, bi ); + __emovz( ci, ai ); + } + } +if( ai[0] == bi[0] ) + { + __eaddm( ai, bi ); + subflg = 0; + } +else + { + __esubm( ai, bi ); + subflg = 1; + } +__emdnorm( bi, lost, subflg, ltb, 64, NBITS); + +done: +__emovo( bi, c ); +} + + +/* y = largest integer not greater than x + * (truncated toward minus infinity) + * + * unsigned short x[NE], y[NE] + * + * efloor( x, y ); + */ + + +void __efloor(short unsigned int *x, short unsigned int *y) +{ +register unsigned short *p; +int e, expon, i; +unsigned short f[NE]; +const unsigned short bmask[] = { +0xffff, +0xfffe, +0xfffc, +0xfff8, +0xfff0, +0xffe0, +0xffc0, +0xff80, +0xff00, +0xfe00, +0xfc00, +0xf800, +0xf000, +0xe000, +0xc000, +0x8000, +0x0000, +}; + +__emov( x, f ); /* leave in external format */ +expon = (int )f[NE-1]; +e = (expon & 0x7fff) - (EXONE - 1); +if( e <= 0 ) + { + __eclear(y); + goto isitneg; + } +/* number of bits to clear out */ +e = NBITS - e; +__emov( f, y ); +if( e <= 0 ) + return; + +p = &y[0]; +while( e >= 16 ) + { + *p++ = 0; + e -= 16; + } +/* clear the remaining bits */ +*p &= bmask[e]; +/* truncate negatives toward minus infinity */ +isitneg: + +if( (unsigned short )expon & (unsigned short )0x8000 ) + { + for( i=0; i<NE-1; i++ ) + { + if( f[i] != y[i] ) + { + __esub( __eone, y, y ); + break; + } + } + } +} + +/* +; Subtract external format numbers. +; +; unsigned short a[NE], b[NE], c[NE]; +; esub( a, b, c ); c = b - a +*/ + + +void __esub(const short unsigned int * a, + const short unsigned int * b, + short unsigned int * c) +{ + +#ifdef NANS +if( __eisnan(a) ) + { + __emov (a, c); + return; + } +if( __eisnan(b) ) + { + __emov(b,c); + return; + } +/* Infinity minus infinity is a NaN. + * Test for subtracting infinities of the same sign. + */ +if( __eisinf(a) && __eisinf(b) && ((__eisneg (a) ^ __eisneg (b)) == 0)) + { + mtherr( "esub", DOMAIN ); + __enan_NBITS( c ); + return; + } +#endif +__eadd1( a, b, c, 1 ); +} + + + +/* +; Divide. +; +; unsigned short a[NI], b[NI], c[NI]; +; ediv( a, b, c ); c = b / a +*/ + +void __ediv(const short unsigned int *a, + const short unsigned int *b, + short unsigned int *c) +{ +unsigned short ai[NI], bi[NI]; +int i; +long lt, lta, ltb; + +#ifdef NANS +/* Return any NaN input. */ +if( __eisnan(a) ) + { + __emov(a,c); + return; + } +if( __eisnan(b) ) + { + __emov(b,c); + return; + } +/* Zero over zero, or infinity over infinity, is a NaN. */ +if( (__eiszero(a) && __eiszero(b)) + || (__eisinf (a) && __eisinf (b)) ) + { + mtherr( "ediv", DOMAIN ); + __enan_NBITS( c ); + return; + } +#endif +/* Infinity over anything else is infinity. */ +#ifdef INFINITY +if( __eisinf(b) ) + { + if( __eisneg(a) ^ __eisneg(b) ) + *(c+(NE-1)) = 0x8000; + else + *(c+(NE-1)) = 0; + __einfin(c); + return; + } +if( __eisinf(a) ) + { + __eclear(c); + return; + } +#endif +__emovi( a, ai ); +__emovi( b, bi ); +lta = ai[E]; +ltb = bi[E]; +if( bi[E] == 0 ) + { /* See if numerator is zero. */ + for( i=1; i<NI-1; i++ ) + { + if( bi[i] != 0 ) + { + ltb -= __enormlz( bi ); + goto dnzro1; + } + } + __eclear(c); + return; + } +dnzro1: + +if( ai[E] == 0 ) + { /* possible divide by zero */ + for( i=1; i<NI-1; i++ ) + { + if( ai[i] != 0 ) + { + lta -= __enormlz( ai ); + goto dnzro2; + } + } + if( ai[0] == bi[0] ) + *(c+(NE-1)) = 0; + else + *(c+(NE-1)) = 0x8000; + __einfin(c); + mtherr( "ediv", SING ); + return; + } +dnzro2: + +i = __edivm( ai, bi ); +/* calculate exponent */ +lt = ltb - lta + EXONE; +__emdnorm( bi, i, 0, lt, 64, NBITS ); +/* set the sign */ +if( ai[0] == bi[0] ) + bi[0] = 0; +else + bi[0] = 0Xffff; +__emovo( bi, c ); +} + +void __e64toe(short unsigned int *pe, short unsigned int *y) +{ +unsigned short yy[NI]; +unsigned short *p, *q, *e; +int i; + +e = pe; +p = yy; +for( i=0; i<NE-5; i++ ) + *p++ = 0; +#ifdef IBMPC +for( i=0; i<5; i++ ) + *p++ = *e++; +#endif +#ifdef DEC +for( i=0; i<5; i++ ) + *p++ = *e++; +#endif +#ifdef MIEEE +p = &yy[0] + (NE-1); +*p-- = *e++; +++e; +for( i=0; i<4; i++ ) + *p-- = *e++; +#endif + +#ifdef IBMPC +/* For Intel long double, shift denormal significand up 1 + -- but only if the top significand bit is zero. */ +if((yy[NE-1] & 0x7fff) == 0 && (yy[NE-2] & 0x8000) == 0) + { + unsigned short temp[NI+1]; + __emovi(yy, temp); + __eshup1(temp); + __emovo(temp,y); + return; + } +#endif +#ifdef INFINITY +/* Point to the exponent field. */ +p = &yy[NE-1]; +if( *p == 0x7fff ) + { +#ifdef NANS +#ifdef IBMPC + for( i=0; i<4; i++ ) + { + if((i != 3 && pe[i] != 0) + /* Check for Intel long double infinity pattern. */ + || (i == 3 && pe[i] != 0x8000)) + { + __enan_NBITS( y ); + return; + } + } +#else + for( i=1; i<=4; i++ ) + { + if( pe[i] != 0 ) + { + __enan_NBITS( y ); + return; + } + } +#endif +#endif /* NANS */ + __eclear( y ); + __einfin( y ); + if( *p & 0x8000 ) + __eneg(y); + return; + } +#endif +p = yy; +q = y; +for( i=0; i<NE; i++ ) + *q++ = *p++; +} + +#endif /* USE_LDTOA */ |