summaryrefslogtreecommitdiffstats
path: root/winsup/mingw/mingwex/ldtoa.c
diff options
context:
space:
mode:
Diffstat (limited to 'winsup/mingw/mingwex/ldtoa.c')
-rw-r--r--winsup/mingw/mingwex/ldtoa.c614
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 */