diff options
author | Danny Smith <dannysmith@users.sourceforge.net> | 2002-11-26 03:11:06 +0300 |
---|---|---|
committer | Danny Smith <dannysmith@users.sourceforge.net> | 2002-11-26 03:11:06 +0300 |
commit | 5362be5926f9b07cb8c61fc10e0ed4485fc8bf86 (patch) | |
tree | 6759ce6b8c7e1d2112b9545e7f26a3ea1155b5df /winsup/mingw/mingwex/math | |
parent | 4e85569d1146c519fdb51758df7d8ca5c92e4311 (diff) |
Add strtold and wcstold to libmingwex.a
* mingwex/strtold.c: New file.
* mingwex/wcstold.c: New file.
* mingwex/ldtoa.c: New file.
* mingwex/math/cephes_emath.h: New file.
* mingwex/math/cephes_emath.c: New file.
* mingwex/Makefile.in (DISTFILES): Add new files.
(MATH_DISTFILES): Ditto.
(STDLIB_OBJS): New. Define as strtold.c wcstold.c.
(MATH_OBJS): Add cephes_emath.o.
(LIB_OBJS): Add $(STDLIB_OBJS).
* include/stdlib.h (strtold, wcstold): Add prototypes.
* include/wchar.h (wcstold): Add prototype.
Add missing ChangeLog entry for 2002-11-09.
Diffstat (limited to 'winsup/mingw/mingwex/math')
-rw-r--r-- | winsup/mingw/mingwex/math/cephes_emath.c | 1318 | ||||
-rw-r--r-- | winsup/mingw/mingwex/math/cephes_emath.h | 713 |
2 files changed, 2031 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 000000000..ab798a2d2 --- /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 */ diff --git a/winsup/mingw/mingwex/math/cephes_emath.h b/winsup/mingw/mingwex/math/cephes_emath.h new file mode 100644 index 000000000..133b0e03f --- /dev/null +++ b/winsup/mingw/mingwex/math/cephes_emath.h @@ -0,0 +1,713 @@ +#ifndef _CEPHES_EMATH_H +#define _CEPHES_EMATH_H + +/* 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> + * + */ + + + +/* ieee.c + * + * Extended precision IEEE binary floating point arithmetic routines + * + * Numbers are stored in C language as arrays of 16-bit unsigned + * short integers. The arguments of the routines are pointers to + * the arrays. + * + * + * External e type data structure, simulates Intel 8087 chip + * temporary real format but possibly with a larger significand: + * + * NE-1 significand words (least significant word first, + * most significant bit is normally set) + * exponent (value = EXONE for 1.0, + * top bit is the sign) + * + * + * Internal data structure of a number (a "word" is 16 bits): + * + * ei[0] sign word (0 for positive, 0xffff for negative) + * ei[1] biased __exponent (value = EXONE for the number 1.0) + * ei[2] high guard word (always zero after normalization) + * ei[3] + * to ei[NI-2] significand (NI-4 significand words, + * most significant word first, + * most significant bit is set) + * ei[NI-1] low guard word (0x8000 bit is rounding place) + * + * + * + * Routines for external format numbers + * + * __asctoe64( string, &d ) ASCII string to long double + * __asctoeg( string, e, prec ) ASCII string to specified precision + * __e64toe( &d, e ) IEEE long double precision to e type + * __eadd( a, b, c ) c = b + a + * __eclear(e) e = 0 + * __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. + * __ediv( a, b, c ) c = b / a + * __efloor( a, b ) truncate to integer, toward -infinity + * __efrexp( a, exp, s ) extract exponent and significand + * __eifrac( e, &l, frac ) e to long integer and e type fraction + * __euifrac( e, &l, frac ) e to unsigned long integer and e type fraction + * __einfin( e ) set e to infinity, leaving its sign alone + * __eldexp( a, n, b ) multiply by 2**n + * __emov( a, b ) b = a + * __emul( a, b, c ) c = b * a + * __eneg(e) e = -e + * __eround( a, b ) b = nearest integer value to a + * __esub( a, b, c ) c = b - a + * __e24toasc( &f, str, n ) single to ASCII string, n digits after decimal + * __e53toasc( &d, str, n ) double to ASCII string, n digits after decimal + * __e64toasc( &d, str, n ) long double to ASCII string + * __etoasc( e, str, n ) e to ASCII string, n digits after decimal + * __etoe24( e, &f ) convert e type to IEEE single precision + * __etoe53( e, &d ) convert e type to IEEE double precision + * __etoe64( e, &d ) convert e type to IEEE long double precision + * __eisneg( e ) 1 if sign bit of e != 0, else 0 + * __eisinf( e ) 1 if e has maximum exponent (non-IEEE) + * or is infinite (IEEE) + * __eisnan( e ) 1 if e is a NaN + * __esqrt( a, b ) b = square root of a + * + * + * Routines for internal format numbers + * + * __eaddm( ai, bi ) add significands, bi = bi + ai + * __ecleaz(ei) ei = 0 + * __ecleazs(ei) set ei = 0 but leave its sign alone + * __ecmpm( ai, bi ) compare significands, return 1, 0, or -1 + * __edivm( ai, bi ) divide significands, bi = bi / ai + * __emdnorm(ai,l,s,exp) normalize and round off + * __emovi( a, ai ) convert external a to internal ai + * __emovo( ai, a ) convert internal ai to external a + * __emovz( ai, bi ) bi = ai, low guard word of bi = 0 + * __emulm( ai, bi ) multiply significands, bi = bi * ai + * __enormlz(ei) left-justify the significand + * __eshdn1( ai ) shift significand and guards down 1 bit + * __eshdn8( ai ) shift down 8 bits + * __eshdn6( ai ) shift down 16 bits + * __eshift( ai, n ) shift ai n bits up (or down if n < 0) + * __eshup1( ai ) shift significand and guards up 1 bit + * __eshup8( ai ) shift up 8 bits + * __eshup6( ai ) shift up 16 bits + * __esubm( ai, bi ) subtract significands, bi = bi - ai + * + * + * The result is always normalized and rounded to NI-4 word precision + * after each arithmetic operation. + * + * Exception flags are NOT fully supported. + * + * Define INFINITY in mconf.h for support of infinity; otherwise a + * saturation arithmetic is implemented. + * + * Define NANS for support of Not-a-Number items; otherwise the + * arithmetic will never produce a NaN output, and might be confused + * by a NaN input. + * If NaN's are supported, the output of ecmp(a,b) is -2 if + * either a or b is a NaN. This means asking if(ecmp(a,b) < 0) + * may not be legitimate. Use if(ecmp(a,b) == -1) for less-than + * if in doubt. + * Signaling NaN's are NOT supported; they are treated the same + * as quiet NaN's. + * + * Denormals are always supported here where appropriate (e.g., not + * for conversion to DEC numbers). + */ + +#include <stdio.h> +#include <stdlib.h> +#include <string.h> +#include <errno.h> +#include <math.h> +#include <locale.h> +#include <ctype.h> + +#define alloca __builtin_alloca + +/* Don't build non-ANSI _IO_ldtoa. It is not thread safe. */ +#ifndef USE_LDTOA +#define USE_LDTOA 0 +#endif + + + /* Number of 16 bit words in external x type format */ +#define NE 6 + + /* Number of 16 bit words in internal format */ +#define NI (NE+3) + + /* Array offset to exponent */ +#define E 1 + + /* Array offset to high guard word */ +#define M 2 + + /* Number of bits of precision */ +#define NBITS ((NI-4)*16) + + /* Maximum number of decimal digits in ASCII conversion + * = NBITS*log10(2) + */ +#define NDEC (NBITS*8/27) + + /* The exponent of 1.0 */ +#define EXONE (0x3fff) + + +#define mtherr(x,y) + + +extern long double strtold (const char * __restrict__ s, char ** __restrict__ se); +extern int __asctoe64(const char * __restrict__ ss, + short unsigned int * __restrict__ y); +extern void __emul(const short unsigned int * a, + const short unsigned int * b, + short unsigned int * c); +extern int __ecmp(const short unsigned int * __restrict__ a, + const short unsigned int * __restrict__ b); +extern int __enormlz(short unsigned int *x); +extern int __eshift(short unsigned int *x, int sc); +extern void __eaddm(const short unsigned int * __restrict__ x, + short unsigned int * __restrict__ y); +extern void __esubm(const short unsigned int * __restrict__ x, + short unsigned int * __restrict__ y); +extern void __emdnorm(short unsigned int *s, int lost, int subflg, + long int exp, int rcntrl, const int rndprc); +extern void __toe64(short unsigned int * __restrict__ a, + short unsigned int * __restrict__ b); +extern int __edivm(short unsigned int * __restrict__ den, + short unsigned int * __restrict__ num); +extern int __emulm(const short unsigned int * __restrict__ a, + short unsigned int * __restrict__ b); +extern void __emovi(const short unsigned int * __restrict__ a, + short unsigned int * __restrict__ b); +extern void __emovo(const short unsigned int * __restrict__ a, + short unsigned int * __restrict__ b); + +#if USE_LDTOA + +extern char * _IO_ldtoa(long double, int, int, int *, int *, char **); +extern void _IO_ldtostr(long double *x, char *string, int ndigs, + int flags, char fmt); + +extern void __eiremain(short unsigned int * __restrict__ den, + short unsigned int *__restrict__ num, + short unsigned int *__restrict__ equot); +extern void __efloor(short unsigned int *x, short unsigned int *y); +extern void __eadd1(const short unsigned int * __restrict__ a, + const short unsigned int * __restrict__ b, + short unsigned int * __restrict__ c, + int subflg); +extern void __esub(const short unsigned int *a, const short unsigned int *b, + short unsigned int *c); +extern void __ediv(const short unsigned int *a, const short unsigned int *b, + short unsigned int *c); +extern void __e64toe(short unsigned int *pe, short unsigned int *y); + + +#endif + +static __inline__ int __eisneg(const short unsigned int *x); +static __inline__ int __eisinf(const short unsigned int *x); +static __inline__ int __eisnan(const short unsigned int *x); +static __inline__ int __eiszero(const short unsigned int *a); +static __inline__ void __emovz(register const short unsigned int * __restrict__ a, + register short unsigned int * __restrict__ b); +static __inline__ void __eclear(register short unsigned int *x); +static __inline__ void __ecleaz(register short unsigned int *xi); +static __inline__ void __ecleazs(register short unsigned int *xi); +static __inline__ int __eiisinf(const short unsigned int *x); +static __inline__ int __eiisnan(const short unsigned int *x); +static __inline__ int __eiiszero(const short unsigned int *x); +static __inline__ void __enan_64(short unsigned int *nan); +static __inline__ void __enan_NBITS (short unsigned int *nan); +static __inline__ void __enan_NI16 (short unsigned int *nan); +static __inline__ void __einfin(register short unsigned int *x); +static __inline__ void __eneg(short unsigned int *x); +static __inline__ void __eshup1(register short unsigned int *x); +static __inline__ void __eshup8(register short unsigned int *x); +static __inline__ void __eshup6(register short unsigned int *x); +static __inline__ void __eshdn1(register short unsigned int *x); +static __inline__ void __eshdn8(register short unsigned int *x); +static __inline__ void __eshdn6(register short unsigned int *x); + + + +/* Intel IEEE, low order words come first: + */ +#define IBMPC 1 + +/* Define 1 for ANSI C atan2() function + * See atan.c and clog.c. + */ +#define ANSIC 1 + +/*define VOLATILE volatile*/ +#define VOLATILE + +/* For 12-byte long doubles on an i386, pad a 16-bit short 0 + * to the end of real constants initialized by integer arrays. + * + * #define XPD 0, + * + * Otherwise, the type is 10 bytes long and XPD should be + * defined blank. + * + * #define XPD + */ +#define XPD 0, +/* #define XPD */ +#define NANS +#define INFINITY + +/* NaN's require infinity support. */ +#ifdef NANS +#ifndef INFINITY +#define INFINITY +#endif +#endif + +/* This handles 64-bit long ints. */ +#define LONGBITS (8 * sizeof(long)) + + +#define NTEN 12 +#define MAXP 4096 + +extern const unsigned short __etens[NTEN + 1][NE]; + +/* +; Clear out entire external format number. +; +; unsigned short x[]; +; eclear( x ); +*/ + +static __inline__ void __eclear(register short unsigned int *x) +{ + memset(x, 0, NE * sizeof(unsigned short)); +} + + +/* Move external format number from a to b. + * + * emov( a, b ); + */ + +static __inline__ void __emov(register const short unsigned int * __restrict__ a, + register short unsigned int * __restrict__ b) +{ + memcpy(b, a, NE * sizeof(unsigned short)); +} + + +/* +; Negate external format number +; +; unsigned short x[NE]; +; eneg( x ); +*/ + +static __inline__ void __eneg(short unsigned int *x) +{ + +#ifdef NANS +if( __eisnan(x) ) + return; +#endif +x[NE-1] ^= 0x8000; /* Toggle the sign bit */ +} + + +/* Return 1 if external format number is negative, + * else return zero. + */ +static __inline__ int __eisneg(const short unsigned int *x) +{ + +#ifdef NANS +if( __eisnan(x) ) + return( 0 ); +#endif +if( x[NE-1] & 0x8000 ) + return( 1 ); +else + return( 0 ); +} + + +/* Return 1 if external format number has maximum possible exponent, + * else return zero. + */ +static __inline__ int __eisinf(const short unsigned int *x) +{ + +if( (x[NE-1] & 0x7fff) == 0x7fff ) + { +#ifdef NANS + if( __eisnan(x) ) + return( 0 ); +#endif + return( 1 ); + } +else + return( 0 ); +} + +/* Check if e-type number is not a number. + */ +static __inline__ int __eisnan(const short unsigned int *x) +{ +#ifdef NANS +int i; +/* NaN has maximum __exponent */ +if( (x[NE-1] & 0x7fff) == 0x7fff ) +/* ... and non-zero significand field. */ + for( i=0; i<NE-1; i++ ) + { + if( *x++ != 0 ) + return (1); + } +#endif +return (0); +} + +/* +; Fill __entire number, including __exponent and significand, with +; largest possible number. These programs implement a saturation +; value that is an ordinary, legal number. A special value +; "infinity" may also be implemented; this would require tests +; for that value and implementation of special rules for arithmetic +; operations involving inifinity. +*/ + +static __inline__ void __einfin(register short unsigned int *x) +{ +register int i; + +#ifdef INFINITY +for( i=0; i<NE-1; i++ ) + *x++ = 0; +*x |= 32767; +#else +for( i=0; i<NE-1; i++ ) + *x++ = 0xffff; +*x |= 32766; +*(x-5) = 0; +#endif +} + +/* Clear out internal format number. + */ + +static __inline__ void __ecleaz(register short unsigned int *xi) +{ + memset(xi, 0, NI * sizeof(unsigned short)); +} + +/* same, but don't touch the sign. */ + +static __inline__ void __ecleazs(register short unsigned int *xi) +{ + ++xi; + memset(xi, 0, (NI-1) * sizeof(unsigned short)); +} + + + +/* Move internal format number from a to b. + */ +static __inline__ void __emovz(register const short unsigned int * __restrict__ a, + register short unsigned int * __restrict__ b) +{ + memcpy(b, a, (NI-1) * sizeof(unsigned short)); + b[NI-1]=0; +} + +/* Return nonzero if internal format number is a NaN. + */ + +static __inline__ int __eiisnan (const short unsigned int *x) +{ +int i; + +if( (x[E] & 0x7fff) == 0x7fff ) + { + for( i=M+1; i<NI; i++ ) + { + if( x[i] != 0 ) + return(1); + } + } +return(0); +} + +/* Return nonzero if external format number is zero. */ + +static __inline__ int +__eiszero(const short unsigned int * a) +{ +if (*((long double*) a) == 0) + return (1); +return (0); +} + +/* Return nonzero if internal format number is zero. */ + +static __inline__ int +__eiiszero(const short unsigned int * ai) +{ + int i; + /* skip the sign word */ + for( i=1; i<NI-1; i++ ) + { + if( ai[i] != 0 ) + return (0); + } + return (1); +} + + +/* Return nonzero if internal format number is infinite. */ + +static __inline__ int +__eiisinf (const unsigned short *x) +{ + +#ifdef NANS + if (__eiisnan (x)) + return (0); +#endif + if ((x[E] & 0x7fff) == 0x7fff) + return (1); + return (0); +} + +/* +; Compare significands of numbers in internal format. +; Guard words are included in the comparison. +; +; unsigned short a[NI], b[NI]; +; cmpm( a, b ); +; +; for the significands: +; returns +1 if a > b +; 0 if a == b +; -1 if a < b +*/ +static __inline__ int __ecmpm(register const short unsigned int * __restrict__ a, + register const short unsigned int * __restrict__ b) +{ +int i; + +a += M; /* skip up to significand area */ +b += M; +for( i=M; i<NI; i++ ) + { + if( *a++ != *b++ ) + goto difrnt; + } +return(0); + +difrnt: +if( *(--a) > *(--b) ) + return(1); +else + return(-1); +} + + +/* +; Shift significand down by 1 bit +*/ + +static __inline__ void __eshdn1(register short unsigned int *x) +{ +register unsigned short bits; +int i; + +x += M; /* point to significand area */ + +bits = 0; +for( i=M; i<NI; i++ ) + { + if( *x & 1 ) + bits |= 1; + *x >>= 1; + if( bits & 2 ) + *x |= 0x8000; + bits <<= 1; + ++x; + } +} + +/* +; Shift significand up by 1 bit +*/ + +static __inline__ void __eshup1(register short unsigned int *x) +{ +register unsigned short bits; +int i; + +x += NI-1; +bits = 0; + +for( i=M; i<NI; i++ ) + { + if( *x & 0x8000 ) + bits |= 1; + *x <<= 1; + if( bits & 2 ) + *x |= 1; + bits <<= 1; + --x; + } +} + + + +/* +; Shift significand down by 8 bits +*/ + +static __inline__ void __eshdn8(register short unsigned int *x) +{ +register unsigned short newbyt, oldbyt; +int i; + +x += M; +oldbyt = 0; +for( i=M; i<NI; i++ ) + { + newbyt = *x << 8; + *x >>= 8; + *x |= oldbyt; + oldbyt = newbyt; + ++x; + } +} + +/* +; Shift significand up by 8 bits +*/ + +static __inline__ void __eshup8(register short unsigned int *x) +{ +int i; +register unsigned short newbyt, oldbyt; + +x += NI-1; +oldbyt = 0; + +for( i=M; i<NI; i++ ) + { + newbyt = *x >> 8; + *x <<= 8; + *x |= oldbyt; + oldbyt = newbyt; + --x; + } +} + +/* +; Shift significand up by 16 bits +*/ + +static __inline__ void __eshup6(register short unsigned int *x) +{ +int i; +register unsigned short *p; + +p = x + M; +x += M + 1; + +for( i=M; i<NI-1; i++ ) + *p++ = *x++; + +*p = 0; +} + +/* +; Shift significand down by 16 bits +*/ + +static __inline__ void __eshdn6(register short unsigned int *x) +{ +int i; +register unsigned short *p; + +x += NI-1; +p = x + 1; + +for( i=M; i<NI-1; i++ ) + *(--p) = *(--x); + +*(--p) = 0; +} + +/* +; Add significands +; x + y replaces y +*/ + +static __inline__ void __enan_64(unsigned short* nan) +{ + static const unsigned short nan64[6] + = {0, 0, 0, 0xc000, 0xffff, 0}; + nan = (unsigned short*) nan64; + return; +} + +static __inline__ void __enan_NBITS(unsigned short* nan) +{ + int i; + for( i=0; i<NE-2; i++ ) + *nan++ = 0; + *nan++ = 0xc000; + *nan++ = 0x7fff; + return; +} + +static __inline__ void __enan_NI16(unsigned short* nan) +{ + int i; + *nan++ = 0; + *nan = 0x7fff; + *nan = 0; + *nan = 0xc000; + for( i=4; i<NI; i++ ) + *nan++ = 0; + return; +} + + +#endif /* _CEPHES_EMATH_H */ + |