newlib-cygwin/winsup/mingw/mingwex/ldtoa.c

615 lines
13 KiB
C

/* 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 */