diff options
Diffstat (limited to 'winsup/mingw/mingwex/strtold.c')
| -rw-r--r-- | winsup/mingw/mingwex/strtold.c | 384 |
1 files changed, 384 insertions, 0 deletions
diff --git a/winsup/mingw/mingwex/strtold.c b/winsup/mingw/mingwex/strtold.c new file mode 100644 index 000000000..81db8a6f2 --- /dev/null +++ b/winsup/mingw/mingwex/strtold.c @@ -0,0 +1,384 @@ +/* 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 "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, sgnflg, nexp, exp, prec, lost; +int k, trail, c; +long lexp; +unsigned short nsign; +const unsigned short *p; +char *sp, *lstr; +char *s; + +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; + +lost = 0; +nsign = 0; +decflg = 0; +sgnflg = 0; +nexp = 0; +exp = 0; +prec = 0; +__ecleaz( yy ); +trail = 0; + +nxtcom: +k = *s - '0'; +if( (k >= 0) && (k <= 9) ) + { +/* Ignore leading zeros */ + if( (prec == 0) && (decflg == 0) && (k == 0) ) + goto donchr; +/* 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-- = 'z'; + trail = 1; + if( *s == 'z' ) + goto donchr; + } +/* 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; + } + prec += 1; + goto donchr; + } +if (*s == dec_sym) + { + if( decflg ) + goto daldone; + ++decflg; + } +else + switch( *s ) + { + case 'z': + break; + case 'E': + case 'e': + goto expnt; + case '-': + nsign = 0xffff; + if( sgnflg ) + goto daldone; + ++sgnflg; + break; + case '+': + if( sgnflg ) + goto daldone; + ++sgnflg; + break; + case 'i': + case 'I': + { + s++; + if (*s != 'n' && *s != 'N') + goto zero; + s++; + if (*s != 'f' && *s != 'F') + goto zero; + s++; + if ((*s == 'i' || *s == 'I') && (s[1] == 'n' || s[1] == 'N') + && (s[2] == 'i' || s[2] == 'I') + && (s[3] == 't' || s[3] == 'T') + && (s[4] == 'y' || s[4] == 'Y')) + s += 5; + goto infinite; + } + case 'n': + case 'N': + { + s++; + if (*s != 'a' && *s != 'A') + goto zero; + s++; + if (*s != 'n' && *s != 'N') + goto zero; + s++; + __enan_NI16( yy ); + goto aexit; + } + default: + goto daldone; + } +donchr: +++s; +goto nxtcom; + +/* Exponent interpretation */ +expnt: + +esign = 1; +exp = 0; +++s; +/* check for + or - */ +if( *s == '-' ) + { + esign = -1; + ++s; + } +if( *s == '+' ) + ++s; +while( (*s >= '0') && (*s <= '9') && exp < 4978) + { + exp *= 10; + exp += *s++ - '0'; + } +if( esign < 0 ) + exp = -exp; +if( exp > 4932 ) + { + errno = ERANGE; +infinite: + __ecleaz(yy); + yy[E] = 0x7fff; /* infinity */ + goto aexit; + } +if( exp < -4977 ) + { + errno = ERANGE; +zero: + __ecleaz(yy); + goto aexit; + } + +daldone: +nexp = exp - 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 ); +exp = 1; +do + { + if( exp & nexp ) + __emul( p, xt, xt ); + p -= NE; + exp = exp + exp; + } +while( exp <= 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 ); +return (lenldstr + s - lstr); +} + + +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; +} |