385 lines
7.8 KiB
C
385 lines
7.8 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>
|
|
*/
|
|
|
|
|
|
#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;
|
|
}
|