crt: Use the __mingw_strtold implementation for strtold, too The second implementation, which was under the 'strtold' name, was lacking (it didn't support e.g. strings like "0x125p-1"), while the '__mingw_strtold' one seemed correct and complete. Signed-off-by: Martin Storsjö <martin@martin.st>
diff --git a/mingw-w64-crt/Makefile.am b/mingw-w64-crt/Makefile.am index 8f3b4b8..edfbfc5 100644 --- a/mingw-w64-crt/Makefile.am +++ b/mingw-w64-crt/Makefile.am
@@ -435,7 +435,7 @@ gdtoa/strtopx.c gdtoa/sum.c gdtoa/ulp.c \ \ math/cbrt.c \ - math/cbrtf.c math/cbrtl.c math/cephes_emath.c math/copysign.c math/copysignf.c \ + math/cbrtf.c math/cbrtl.c math/copysign.c math/copysignf.c \ math/coshf.c math/coshl.c math/erfl.c \ math/expf.c \ math/fabs.c math/fabsf.c math/fabsl.c math/fdim.c math/fdimf.c math/fdiml.c \ @@ -457,7 +457,7 @@ math/sqrtf.c math/sqrtl.c math/tanhf.c math/tanhl.c math/tgamma.c \ math/tgammaf.c math/tgammal.c math/truncl.c \ math/powi.def.h math/sqrt.def.h \ - math/cephes_emath.h math/cephes_mconf.h math/fp_consts.h math/abs64.c \ + math/cephes_mconf.h math/fp_consts.h math/abs64.c \ \ misc/mb_wc_common.h \ misc/mingw_longjmp.S \ @@ -474,7 +474,7 @@ misc/mingw_matherr.c misc/mingw_mbwc_convert.c misc/mingw_usleep.c misc/mingw_wcstod.c misc/mingw_wcstof.c \ misc/mkstemp.c misc/seterrno.c misc/sleep.c \ misc/strnlen.c misc/strsafe.c \ - misc/strtoimax.c misc/strtold.c misc/strtoumax.c misc/tdelete.c misc/tfind.c \ + misc/strtoimax.c misc/strtoumax.c misc/tdelete.c misc/tfind.c \ misc/tsearch.c misc/twalk.c \ misc/wcsnlen.c misc/wcstof.c \ misc/wcstoimax.c misc/wcstold.c misc/wcstoumax.c misc/wctob.c misc/wctrans.c \
diff --git a/mingw-w64-crt/gdtoa/strtodnrp.c b/mingw-w64-crt/gdtoa/strtodnrp.c index a7b9972..97d74fc 100644 --- a/mingw-w64-crt/gdtoa/strtodnrp.c +++ b/mingw-w64-crt/gdtoa/strtodnrp.c
@@ -93,4 +93,8 @@ long double __cdecl __mingw_strtold (const char * __restrict__ src, char ** __restrict__ endptr) __attribute__((alias("__strtod"))); + +long double __cdecl +strtold (const char * __restrict__ src, char ** __restrict__ endptr) + __attribute__((alias("__strtod"))); #endif
diff --git a/mingw-w64-crt/gdtoa/strtopx.c b/mingw-w64-crt/gdtoa/strtopx.c index 64db719..eb0d63b 100644 --- a/mingw-w64-crt/gdtoa/strtopx.c +++ b/mingw-w64-crt/gdtoa/strtopx.c
@@ -127,4 +127,8 @@ __mingw_strtold (const char * __restrict__ src, char ** __restrict__ endptr) __attribute__((alias("__strtold"))); +long double __cdecl +strtold (const char * __restrict__ src, char ** __restrict__ endptr) + __attribute__((alias("__strtold"))); + #endif
diff --git a/mingw-w64-crt/math/cephes_emath.c b/mingw-w64-crt/math/cephes_emath.c deleted file mode 100644 index 41e297c..0000000 --- a/mingw-w64-crt/math/cephes_emath.c +++ /dev/null
@@ -1,1283 +0,0 @@ -/** - * This file has no copyright assigned and is placed in the Public Domain. - * This file is part of the mingw-w64 runtime package. - * No warranty is given; refer to the file DISCLAIMER.PD within this package. - */ -#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 "expo" is the biased exponent, which may be negative. - * the exponent field of "s" is ignored but is replaced by - * "expo" 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, int expo, 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 - expo -= j; -#ifndef INFINITY - if (expo >= 32767L) - goto overf; -#else - if ((j > NBITS) && (expo < 32767L)) - { - __ecleazs(s); - return; - } -#endif - if (expo < 0L) - { - if (expo > (long)(-NBITS - 1)) - { - j = (int)expo; - 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 ((expo <= 0) && (rndprc != NBITS)) -#else - if ((expo <= 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 ^ 0xffff); - 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 ((expo <= 0) && (rndprc != NBITS)) -#else - if ((expo <= 0) && (rndprc != 64) && (rndprc != NBITS)) -#endif - { - __eshup1(s); - } - if (s[2] != 0) - { /* overflow on roundoff */ - __eshdn1(s); - expo += 1; - } -mdfin: - s[NI - 1] = 0; - if (expo >= 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 ^ 0xffff); -#endif - return; - } - if (expo < 0) - s[1] = 0; - else - s[1] = (unsigned short)expo; -} - - -/* -; 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/mingw-w64-crt/math/cephes_emath.h b/mingw-w64-crt/math/cephes_emath.h deleted file mode 100644 index 58a8e13..0000000 --- a/mingw-w64-crt/math/cephes_emath.h +++ /dev/null
@@ -1,719 +0,0 @@ -/** - * This file has no copyright assigned and is placed in the Public Domain. - * This file is part of the mingw-w64 runtime package. - * No warranty is given; refer to the file DISCLAIMER.PD within this package. - */ -#ifndef _CEPHES_EMATH_H -#define _CEPHES_EMATH_H - -/** - * This is a workaround for a gcc bug - */ -#define __restrict__ - -/* 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> - -#undef alloca -#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(fname, code) - - -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, - 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 *nanptr); -static __inline__ void __enan_NBITS (short unsigned int *nanptr); -static __inline__ void __enan_NI16 (short unsigned int *nanptr); -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 1 - -/* 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 - -/* -; 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) -{ - union { - long double ld; - unsigned short sh[8]; - } av; - av.ld = 0.0; - memcpy (av.sh, a, 12); - if (av.ld == 0.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* nanptr) -{ - int i; - for (i = 0; i < 3; i++) - *nanptr++ = 0; - *nanptr++ = 0xc000; - *nanptr++ = 0x7fff; - *nanptr = 0; - return; -} - -static __inline__ void __enan_NBITS(unsigned short* nanptr) -{ - int i; - for (i = 0; i < NE - 2; i++) - *nanptr++ = 0; - *nanptr++ = 0xc000; - *nanptr = 0x7fff; - return; -} - -static __inline__ void __enan_NI16(unsigned short* nanptr) -{ - int i; - *nanptr++ = 0; - *nanptr++ = 0x7fff; - *nanptr++ = 0; - *nanptr++ = 0xc000; - for (i = 4; i < NI; i++) - *nanptr++ = 0; - return; -} - -#endif /* _CEPHES_EMATH_H */ -
diff --git a/mingw-w64-crt/misc/strtold.c b/mingw-w64-crt/misc/strtold.c deleted file mode 100644 index 2b99eba..0000000 --- a/mingw-w64-crt/misc/strtold.c +++ /dev/null
@@ -1,398 +0,0 @@ -/** - * This file has no copyright assigned and is placed in the Public Domain. - * This file is part of the mingw-w64 runtime package. - * No warranty is given; refer to the file DISCLAIMER.PD within this package. - */ - -#include "math/cephes_emath.h" - -#if NE == 10 -/* 1.0E0 */ -static const unsigned short __eone[NE] = { - 0x0000, 0x0000, 0x0000, 0x0000, - 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3fff, -}; -#else -static const unsigned short __eone[NE] = { - 0, 0000000,0000000,0000000,0100000,0x3fff, -}; -#endif - -#if NE == 10 -static const unsigned short __etens[NTEN + 1][NE] = -{ - {0x6576, 0x4a92, 0x804a, 0x153f, - 0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,}, /* 10**4096 */ - {0x6a32, 0xce52, 0x329a, 0x28ce, - 0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,}, /* 10**2048 */ - {0x526c, 0x50ce, 0xf18b, 0x3d28, - 0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,}, - {0x9c66, 0x58f8, 0xbc50, 0x5c54, - 0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,}, - {0x851e, 0xeab7, 0x98fe, 0x901b, - 0xddbb, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,}, - {0x0235, 0x0137, 0x36b1, 0x336c, - 0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,}, - {0x50f8, 0x25fb, 0xc76b, 0x6b71, - 0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,}, - {0x0000, 0x0000, 0x0000, 0x0000, - 0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,}, - {0x0000, 0x0000, 0x0000, 0x0000, - 0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,}, - {0x0000, 0x0000, 0x0000, 0x0000, - 0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,}, - {0x0000, 0x0000, 0x0000, 0x0000, - 0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,}, - {0x0000, 0x0000, 0x0000, 0x0000, - 0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,}, - {0x0000, 0x0000, 0x0000, 0x0000, - 0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,}, /* 10**1 */ -}; -#else -static const unsigned short __etens[NTEN+1][NE] = { - {0xc94c,0x979a,0x8a20,0x5202,0xc460,0x7525,},/* 10**4096 */ - {0xa74d,0x5de4,0xc53d,0x3b5d,0x9e8b,0x5a92,},/* 10**2048 */ - {0x650d,0x0c17,0x8175,0x7586,0xc976,0x4d48,}, - {0xcc65,0x91c6,0xa60e,0xa0ae,0xe319,0x46a3,}, - {0xddbc,0xde8d,0x9df9,0xebfb,0xaa7e,0x4351,}, - {0xc66f,0x8cdf,0x80e9,0x47c9,0x93ba,0x41a8,}, - {0x3cbf,0xa6d5,0xffcf,0x1f49,0xc278,0x40d3,}, - {0xf020,0xb59d,0x2b70,0xada8,0x9dc5,0x4069,}, - {0x0000,0x0000,0x0400,0xc9bf,0x8e1b,0x4034,}, - {0x0000,0x0000,0x0000,0x2000,0xbebc,0x4019,}, - {0x0000,0x0000,0x0000,0x0000,0x9c40,0x400c,}, - {0x0000,0x0000,0x0000,0x0000,0xc800,0x4005,}, - {0x0000,0x0000,0x0000,0x0000,0xa000,0x4002,}, /* 10**1 */ -}; -#endif - -int __asctoe64(const char * __restrict__ ss, short unsigned int * __restrict__ y) -{ - unsigned short yy[NI], xt[NI], tt[NI]; - int esign, decflg, nexp, expo, lost; - int k, c; - int valid_lead_string = 0; - int have_non_zero_mant = 0; - int prec = 0; - /* int trail = 0; */ - int lexp; - unsigned short nsign = 0; - const unsigned short *p; - char *sp, *lstr; - char *s; - - const char dec_sym = *(localeconv ()->decimal_point); - - int lenldstr = 0; - - /* Copy the input string. */ - c = strlen (ss) + 2; - lstr = (char *) alloca (c); - s = (char *) ss; - while( isspace ((int)(unsigned char)*s)) /* skip leading spaces */ - { - ++s; - ++lenldstr; - } - sp = lstr; - for (k = 0; k < c; k++) - { - if ((*sp++ = *s++) == '\0') - break; - } - *sp = '\0'; - s = lstr; - - if (*s == '-') - { - nsign = 0xffff; - ++s; - } - else if (*s == '+') - { - ++s; - } - - if (_strnicmp("INF", s , 3) == 0) - { - valid_lead_string = 1; - s += 3; - if ( _strnicmp ("INITY", s, 5) == 0) - s += 5; - __ecleaz(yy); - yy[E] = 0x7fff; /* infinity */ - goto aexit; - } - else if(_strnicmp ("NAN", s, 3) == 0) - { - valid_lead_string = 1; - s += 3; - __enan_NI16( yy ); - goto aexit; - } - - /* FIXME: Handle case of strtold ("NAN(n_char_seq)",endptr) */ - - /* Now get some digits. */ - lost = 0; - decflg = 0; - nexp = 0; - expo = 0; - __ecleaz( yy ); - - /* Ignore leading zeros */ - while (*s == '0') - { - valid_lead_string = 1; - s++; - } - -nxtcom: - - k = *s - '0'; - if ((k >= 0) && (k <= 9)) - { -#if 0 -/* The use of a special char as a flag for trailing zeroes causes problems when input - actually contains the char */ -/* Identify and strip trailing zeros after the decimal point. */ - if ((trail == 0) && (decflg != 0)) - { - sp = s; - while ((*sp >= '0') && (*sp <= '9')) - ++sp; - --sp; - while (*sp == '0') - { - *sp-- = (char)-1; - trail++; - } - if( *s == (char)-1 ) - goto donchr; - } -#endif - -/* If enough digits were given to more than fill up the yy register, - * continuing until overflow into the high guard word yy[2] - * guarantees that there will be a roundoff bit at the top - * of the low guard word after normalization. - */ - if (yy[2] == 0) - { - if( decflg ) - nexp += 1; /* count digits after decimal point */ - __eshup1( yy ); /* multiply current number by 10 */ - __emovz( yy, xt ); - __eshup1( xt ); - __eshup1( xt ); - __eaddm( xt, yy ); - __ecleaz( xt ); - xt[NI-2] = (unsigned short )k; - __eaddm( xt, yy ); - } - else - { - /* Mark any lost non-zero digit. */ - lost |= k; - /* Count lost digits before the decimal point. */ - if (decflg == 0) - nexp -= 1; - } - have_non_zero_mant |= k; - prec ++; - /* goto donchr; */ - } - else if (*s == dec_sym) - { - if( decflg ) - goto daldone; - ++decflg; - } - else if ((*s == 'E') || (*s == 'e') ) - { - if (prec || valid_lead_string) - goto expnt; - else - goto daldone; - } -#if 0 - else if (*s == (char)-1) - goto donchr; -#endif - else /* an invalid char */ - goto daldone; - - /* donchr: */ - ++s; - goto nxtcom; - -/* Exponent interpretation */ -expnt: - - esign = 1; - expo = 0; - /* Save position in case we need to fall back. */ - sp = s; - ++s; - /* check for + or - */ - if (*s == '-') - { - esign = -1; - ++s; - } - if (*s == '+') - ++s; - - /* Check for valid exponent. */ - if (!(*s >= '0' && *s <= '9')) - { - s = sp; - goto daldone; - } - - while ((*s >= '0') && (*s <= '9')) - { - /* Stop modifying exp if we are going to overflow anyway, - but keep parsing the string. */ - if (expo < 4978) - { - expo *= 10; - expo += *s - '0'; - } - s++; - } - - if (esign < 0) - expo = -expo; - - if (expo > 4977) /* maybe overflow */ - { - __ecleaz(yy); - if (have_non_zero_mant) - yy[E] = 0x7fff; - goto aexit; - } - else if (expo < -4977) /* underflow */ - { - __ecleaz(yy); - goto aexit; - } - -daldone: - - nexp = expo - nexp; - - /* Pad trailing zeros to minimize power of 10, per IEEE spec. */ - while ((nexp > 0) && (yy[2] == 0)) - { - __emovz( yy, xt ); - __eshup1( xt ); - __eshup1( xt ); - __eaddm( yy, xt ); - __eshup1( xt ); - if (xt[2] != 0) - break; - nexp -= 1; - __emovz( xt, yy ); - } - if ((k = __enormlz(yy)) > NBITS) - { - __ecleaz(yy); - goto aexit; - } - lexp = (EXONE - 1 + NBITS) - k; - __emdnorm( yy, lost, 0, lexp, 64, NBITS ); - /* convert to external format */ - - /* Multiply by 10**nexp. If precision is 64 bits, - * the maximum relative error incurred in forming 10**n - * for 0 <= n <= 324 is 8.2e-20, at 10**180. - * For 0 <= n <= 999, the peak relative error is 1.4e-19 at 10**947. - * For 0 >= n >= -999, it is -1.55e-19 at 10**-435. - */ - lexp = yy[E]; - if (nexp == 0) - { - k = 0; - goto expdon; - } - esign = 1; - if (nexp < 0) - { - nexp = -nexp; - esign = -1; - if (nexp > 4096) - { /* Punt. Can't handle this without 2 divides. */ - __emovi( __etens[0], tt ); - lexp -= tt[E]; - k = __edivm( tt, yy ); - lexp += EXONE; - nexp -= 4096; - } - } - p = &__etens[NTEN][0]; - __emov( __eone, xt ); - expo = 1; - do - { - if (expo & nexp) - __emul( p, xt, xt ); - p -= NE; - expo = expo + expo; - } - while (expo <= MAXP); - - __emovi( xt, tt ); - if (esign < 0) - { - lexp -= tt[E]; - k = __edivm( tt, yy ); - lexp += EXONE; - } - else - { - lexp += tt[E]; - k = __emulm( tt, yy ); - lexp -= EXONE - 1; - } - -expdon: - - /* Round and convert directly to the destination type */ - - __emdnorm( yy, k, 0, lexp, 64, 64 ); - -aexit: - - yy[0] = nsign; - - __toe64( yy, y ); - - /* Check for overflow, undeflow */ - if (have_non_zero_mant && - (*((long double*) y) == 0.0L || isinf (*((long double*) y)))) - errno = ERANGE; - - if (prec || valid_lead_string) - return (lenldstr + (s - lstr)); - - return 0; -} - - -long double strtold (const char * __restrict__ s, char ** __restrict__ se) -{ - int lenldstr; - union - { - unsigned short int us[6]; - long double ld; - } xx = {{0}}; - - lenldstr = __asctoe64( s, xx.us); - if (se) - *se = (char*)s + lenldstr; - - return xx.ld; -} -