diff options
Diffstat (limited to 'winsup/mingw/mingwex/ldtoa.c')
-rw-r--r-- | winsup/mingw/mingwex/ldtoa.c | 614 |
1 files changed, 614 insertions, 0 deletions
diff --git a/winsup/mingw/mingwex/ldtoa.c b/winsup/mingw/mingwex/ldtoa.c new file mode 100644 index 000000000..09dc7911b --- /dev/null +++ b/winsup/mingw/mingwex/ldtoa.c @@ -0,0 +1,614 @@ +/* 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> + * + */ + + +#ifdef USE_LDTOA + +#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 + +#if NE == 10 +static const unsigned short __emtens[NTEN + 1][NE] = +{ + {0x2030, 0xcffc, 0xa1c3, 0x8123, + 0x2de3, 0x9fde, 0xd2ce, 0x04c8, 0xa6dd, 0x0ad8,}, /* 10**-4096 */ + {0x8264, 0xd2cb, 0xf2ea, 0x12d4, + 0x4925, 0x2de4, 0x3436, 0x534f, 0xceae, 0x256b,}, /* 10**-2048 */ + {0xf53f, 0xf698, 0x6bd3, 0x0158, + 0x87a6, 0xc0bd, 0xda57, 0x82a5, 0xa2a6, 0x32b5,}, + {0xe731, 0x04d4, 0xe3f2, 0xd332, + 0x7132, 0xd21c, 0xdb23, 0xee32, 0x9049, 0x395a,}, + {0xa23e, 0x5308, 0xfefb, 0x1155, + 0xfa91, 0x1939, 0x637a, 0x4325, 0xc031, 0x3cac,}, + {0xe26d, 0xdbde, 0xd05d, 0xb3f6, + 0xac7c, 0xe4a0, 0x64bc, 0x467c, 0xddd0, 0x3e55,}, + {0x2a20, 0x6224, 0x47b3, 0x98d7, + 0x3f23, 0xe9a5, 0xa539, 0xea27, 0xa87f, 0x3f2a,}, + {0x0b5b, 0x4af2, 0xa581, 0x18ed, + 0x67de, 0x94ba, 0x4539, 0x1ead, 0xcfb1, 0x3f94,}, + {0xbf71, 0xa9b3, 0x7989, 0xbe68, + 0x4c2e, 0xe15b, 0xc44d, 0x94be, 0xe695, 0x3fc9,}, + {0x3d4d, 0x7c3d, 0x36ba, 0x0d2b, + 0xfdc2, 0xcefc, 0x8461, 0x7711, 0xabcc, 0x3fe4,}, + {0xc155, 0xa4a8, 0x404e, 0x6113, + 0xd3c3, 0x652b, 0xe219, 0x1758, 0xd1b7, 0x3ff1,}, + {0xd70a, 0x70a3, 0x0a3d, 0xa3d7, + 0x3d70, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,}, + {0xcccd, 0xcccc, 0xcccc, 0xcccc, + 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,}, /* 10**-1 */ +}; + +#else +static const unsigned short __emtens[NTEN+1][NE] = { +{0x2de4,0x9fde,0xd2ce,0x04c8,0xa6dd,0x0ad8,}, /* 10**-4096 */ +{0x4925,0x2de4,0x3436,0x534f,0xceae,0x256b,}, /* 10**-2048 */ +{0x87a6,0xc0bd,0xda57,0x82a5,0xa2a6,0x32b5,}, +{0x7133,0xd21c,0xdb23,0xee32,0x9049,0x395a,}, +{0xfa91,0x1939,0x637a,0x4325,0xc031,0x3cac,}, +{0xac7d,0xe4a0,0x64bc,0x467c,0xddd0,0x3e55,}, +{0x3f24,0xe9a5,0xa539,0xea27,0xa87f,0x3f2a,}, +{0x67de,0x94ba,0x4539,0x1ead,0xcfb1,0x3f94,}, +{0x4c2f,0xe15b,0xc44d,0x94be,0xe695,0x3fc9,}, +{0xfdc2,0xcefc,0x8461,0x7711,0xabcc,0x3fe4,}, +{0xd3c3,0x652b,0xe219,0x1758,0xd1b7,0x3ff1,}, +{0x3d71,0xd70a,0x70a3,0x0a3d,0xa3d7,0x3ff8,}, +{0xcccd,0xcccc,0xcccc,0xcccc,0xcccc,0x3ffb,}, /* 10**-1 */ +}; +#endif + +/* This routine will not return more than NDEC+1 digits. */ +void __etoasc(short unsigned int * __restrict__ x, + char * __restrict__ string, + const int ndigits, const int outformat, + int* outexp) +{ +long digit; +unsigned short y[NI], t[NI], u[NI], w[NI], equot[NI]; +const unsigned short *r, *p; +const unsigned short *ten; +unsigned short sign; +int i, j, k, expon, ndigs; +char *s, *ss; +unsigned short m; + +ndigs = ndigits; +#ifdef NANS +if( __eisnan(x) ) + { + sprintf( string, " NaN " ); + expon = 9999; + goto bxit; + } +#endif +__emov( x, y ); /* retain external format */ +if( y[NE-1] & 0x8000 ) + { + sign = 0xffff; + y[NE-1] &= 0x7fff; + } +else + { + sign = 0; + } +expon = 0; +ten = &__etens[NTEN][0]; +__emov( __eone, t ); +/* Test for zero exponent */ +if( y[NE-1] == 0 ) + { + for( k=0; k<NE-1; k++ ) + { + if( y[k] != 0 ) + goto tnzro; /* denormalized number */ + } + goto isone; /* legal all zeros */ + } +tnzro: + +/* Test for infinity. + */ +if( y[NE-1] == 0x7fff ) + { + if( sign ) + sprintf( string, " -Infinity " ); + else + sprintf( string, " Infinity " ); + expon = 9999; + goto bxit; + } + +/* Test for exponent nonzero but significand denormalized. + * This is an error condition. + */ +if( (y[NE-1] != 0) && ((y[NE-2] & 0x8000) == 0) ) + { + mtherr( "etoasc", DOMAIN ); + sprintf( string, "NaN" ); + expon = 9999; + goto bxit; + } + +/* Compare to 1.0 */ +i = __ecmp( __eone, y ); +if( i == 0 ) + goto isone; + +if( i < 0 ) + { /* Number is greater than 1 */ +/* Convert significand to an integer and strip trailing decimal zeros. */ + __emov( y, u ); + u[NE-1] = EXONE + NBITS - 1; + + p = &__etens[NTEN-4][0]; + m = 16; +do + { + __ediv( p, u, t); + __efloor( t, w ); + for( j=0; j<NE-1; j++ ) + { + if( t[j] != w[j] ) + goto noint; + } + __emov( t, u ); + expon += (int )m; +noint: + p += NE; + m >>= 1; + } +while( m != 0 ); + +/* Rescale from integer significand */ + u[NE-1] += y[NE-1] - (unsigned int )(EXONE + NBITS - 1); + __emov( u, y ); +/* Find power of 10 */ + __emov( __eone, t ); + m = MAXP; + p = &__etens[0][0]; + while( __ecmp( ten, u ) <= 0 ) + { + if( __ecmp( p, u ) <= 0 ) + { + __ediv( p, u, u ); + __emul( p, t, t ); + expon += (int )m; + } + m >>= 1; + if( m == 0 ) + break; + p += NE; + } + } +else + { /* Number is less than 1.0 */ +/* Pad significand with trailing decimal zeros. */ + if( y[NE-1] == 0 ) + { + while( (y[NE-2] & 0x8000) == 0 ) + { + __emul( ten, y, y ); + expon -= 1; + } + } + else + { + __emovi( y, w ); + for( i=0; i<NDEC+1; i++ ) + { + if( (w[NI-1] & 0x7) != 0 ) + break; +/* multiply by 10 */ + __emovz( w, u ); + __eshdn1( u ); + __eshdn1( u ); + __eaddm( w, u ); + u[1] += 3; + while( u[2] != 0 ) + { + __eshdn1(u); + u[1] += 1; + } + if( u[NI-1] != 0 ) + break; + if( __eone[NE-1] <= u[1] ) + break; + __emovz( u, w ); + expon -= 1; + } + __emovo( w, y ); + } + k = -MAXP; + p = &__emtens[0][0]; + r = &__etens[0][0]; + __emov( y, w ); + __emov( __eone, t ); + while( __ecmp( __eone, w ) > 0 ) + { + if( __ecmp( p, w ) >= 0 ) + { + __emul( r, w, w ); + __emul( r, t, t ); + expon += k; + } + k /= 2; + if( k == 0 ) + break; + p += NE; + r += NE; + } + __ediv( t, __eone, t ); + } +isone: +/* Find the first (leading) digit. */ +__emovi( t, w ); +__emovz( w, t ); +__emovi( y, w ); +__emovz( w, y ); +__eiremain( t, y, equot); +digit = equot[NI-1]; +while( (digit == 0) && (__eiszero(y) == 0) ) + { + __eshup1( y ); + __emovz( y, u ); + __eshup1( u ); + __eshup1( u ); + __eaddm( u, y ); + __eiremain( t, y, equot); + digit = equot[NI-1]; + expon -= 1; + } +s = string; +if( sign ) + *s++ = '-'; +else + *s++ = ' '; +/* Examine number of digits requested by caller. */ +if( outformat == 3 ) + ndigs += expon; +/* +else if( ndigs < 0 ) + ndigs = 0; +*/ +if( ndigs > NDEC ) + ndigs = NDEC; +if( digit == 10 ) + { + *s++ = '1'; + *s++ = '.'; + if( ndigs > 0 ) + { + *s++ = '0'; + ndigs -= 1; + } + expon += 1; + if( ndigs < 0 ) + { + ss = s; + goto doexp; + } + } +else + { + *s++ = (char )digit + '0'; + *s++ = '.'; + } +/* Generate digits after the decimal point. */ +for( k=0; k<=ndigs; k++ ) + { +/* multiply current number by 10, without normalizing */ + __eshup1( y ); + __emovz( y, u ); + __eshup1( u ); + __eshup1( u ); + __eaddm( u, y ); + __eiremain( t, y, equot); + *s++ = (char )equot[NI-1] + '0'; + } +digit = equot[NI-1]; +--s; +ss = s; +/* round off the ASCII string */ +if( digit > 4 ) + { +/* Test for critical rounding case in ASCII output. */ + if( digit == 5 ) + { + if( __eiiszero(y) == 0 ) + goto roun; /* round to nearest */ + if( (*(s-1) & 1) == 0 ) + goto doexp; /* round to even */ + } +/* Round up and propagate carry-outs */ +roun: + --s; + k = *s & 0x7f; +/* Carry out to most significant digit? */ + if( ndigs < 0 ) + { + /* This will print like "1E-6". */ + *s = '1'; + + expon += 1; + goto doexp; + } + else if( k == '.' ) + { + --s; + k = *s; + k += 1; + *s = (char )k; +/* Most significant digit carries to 10? */ + if( k > '9' ) + { + expon += 1; + *s = '1'; + } + goto doexp; + } +/* Round up and carry out from less significant digits */ + k += 1; + *s = (char )k; + if( k > '9' ) + { + *s = '0'; + goto roun; + } + } +doexp: +#if defined (__GO32__) || defined (__MINGW32__) +if( expon >= 0 ) + sprintf( ss, "e+%02d", expon ); +else + sprintf( ss, "e-%02d", -expon ); +#else + sprintf( ss, "E%d", expon ); +#endif +bxit: + +if (outexp) + *outexp = expon; +} + + +/* FIXME: Not thread safe */ +static char outstr[128]; + + char * +_IO_ldtoa(long double d, int mode, int ndigits, int *decpt, + int *sign, char **rve) +{ +unsigned short e[NI]; +char *s, *p; +int k; +int outexpon = 0; + +union + { + unsigned short int us[6]; + long double ld; + } xx; +xx.ld = d; +__e64toe(xx.us, e ); +if( __eisneg(e) ) + *sign = 1; +else + *sign = 0; +/* Mode 3 is "f" format. */ +if( mode != 3 ) + ndigits -= 1; +/* Mode 0 is for %.999 format, which is supposed to give a + minimum length string that will convert back to the same binary value. + For now, just ask for 20 digits which is enough but sometimes too many. */ +if( mode == 0 ) + ndigits = 20; +/* This sanity limit must agree with the corresponding one in etoasc, to + keep straight the returned value of outexpon. */ +if( ndigits > NDEC ) + ndigits = NDEC; + +__etoasc( e, outstr, ndigits, mode, &outexpon ); +s = outstr; +if( __eisinf(e) || __eisnan(e) ) + { + *decpt = 9999; + goto stripspaces; + } +*decpt = outexpon + 1; + +/* Transform the string returned by etoasc into what the caller wants. */ + +/* Look for decimal point and delete it from the string. */ +s = outstr; +while( *s != '\0' ) + { + if( *s == '.' ) + goto yesdecpt; + ++s; + } +goto nodecpt; + +yesdecpt: + +/* Delete the decimal point. */ +while( *s != '\0' ) + { + *s = *(s+1); + ++s; + } + +nodecpt: + +/* Back up over the exponent field. */ +while( *s != 'E' && *s != 'e' && s > outstr) + --s; +*s = '\0'; + +stripspaces: + +/* Strip leading spaces and sign. */ +p = outstr; +while( *p == ' ' || *p == '-') + ++p; + +/* Find new end of string. */ +s = outstr; +while( (*s++ = *p++) != '\0' ) + ; +--s; + +/* Strip trailing zeros. */ +if( mode == 2 ) + k = 1; +else if( ndigits > outexpon ) + k = ndigits; +else + k = outexpon; + +while( *(s-1) == '0' && ((s - outstr) > k)) + *(--s) = '\0'; + +/* In f format, flush small off-scale values to zero. + Rounding has been taken care of by etoasc. */ +if( mode == 3 && ((ndigits + outexpon) < 0)) + { + s = outstr; + *s = '\0'; + *decpt = 0; + } + +if( rve ) + *rve = s; +return outstr; +} + +void +_IO_ldtostr(long double *x, char *string, int ndigs, int flags, char fmt) +{ +unsigned short w[NI]; +char *t, *u; +int outexpon = 0; +int outformat = -1; +char dec_sym = *(localeconv()->decimal_point); + +__e64toe( (unsigned short *)x, w ); +__etoasc( w, string, ndigs, outformat, &outexpon ); + +if( ndigs == 0 && flags == 0 ) + { + /* Delete the decimal point unless alternate format. */ + t = string; + while( *t != '.' ) + ++t; + u = t + 1; + while( *t != '\0' ) + *t++ = *u++; + } +if (*string == ' ') + { + t = string; + u = t + 1; + while( *t != '\0' ) + *t++ = *u++; + } +if (fmt == 'E') + { + t = string; + while( *t != 'e' ) + ++t; + *t = 'E'; + } +if (dec_sym != '.') + { + t = string; + while (*t != '.') + ++t; + *t = dec_sym; + } +} + +#endif /* USE_LDTOA */ |