-
/* Extended precision arithmetic functions for long double I/O.
* This program has been placed in the public domain.
*/
/* The exponent of 1.0 */
#define EXONE (0x3fff)
-/* Control structure for long doublue conversion including rounding precision values.
+ /* Maximum exponent digits - base 10 */
+ #define MAX_EXP_DIGITS 5
+
+/* Control structure for long double conversion including rounding precision values.
* rndprc can be set to 80 (if NE=6), 64, 56, 53, or 24 bits.
*/
typedef struct
static void efloor(short unsigned int *x, short unsigned int *y, LDPARMS *ldp);
static void etoasc(short unsigned int *x, char *string, int ndigs, int outformat, LDPARMS *ldp);
+union uconv
+{
+ unsigned short pe;
+ long double d;
+};
+
#if LDBL_MANT_DIG == 24
static void e24toe(short unsigned int *pe, short unsigned int *y, LDPARMS *ldp);
#elif LDBL_MANT_DIG == 53
*
* Exception flags are NOT fully supported.
*
- * Define INFINITY in mconf.h for support of infinity; otherwise a
+ * Define USE_INFINITY in mconf.h for support of infinity; otherwise a
* saturation arithmetic is implemented.
*
* Define NANS for support of Not-a-Number items; otherwise the
#define VOLATILE
#define NANS
-#define INFINITY
+#define USE_INFINITY
/* NaN's require infinity support. */
#ifdef NANS
#ifndef INFINITY
-#define INFINITY
+#define USE_INFINITY
#endif
#endif
{
register int i;
-#ifdef INFINITY
+#ifdef USE_INFINITY
for( i=0; i<NE-1; i++ )
*x++ = 0;
*x |= 32767;
/* get the exponent */
*q = *p--;
*q++ &= 0x7fff; /* delete the sign bit */
-#ifdef INFINITY
+#ifdef USE_INFINITY
if( (*(q-1) & 0x7fff) == 0x7fff )
{
#ifdef NANS
*q-- = *p++ | 0x8000;
else
*q-- = *p++;
-#ifdef INFINITY
+#ifdef USE_INFINITY
if( *(p-1) == 0x7fff )
{
#ifdef NANS
/* Return nonzero if internal format number is infinite. */
static int
-eiisinf (x)
- unsigned short x[];
+eiisinf (unsigned short x[])
{
#ifdef NANS
j = enormlz( s );
/* a blank significand could mean either zero or infinity. */
-#ifndef INFINITY
+#ifndef USE_INFINITY
if( j > NBITS )
{
ecleazs( s );
}
#endif
exp -= j;
-#ifndef INFINITY
+#ifndef USE_INFINITY
if( exp >= 32767L )
goto overf;
#else
s[NI-1] = 0;
if( exp >= 32767L )
{
-#ifndef INFINITY
+#ifndef USE_INFINITY
overf:
#endif
-#ifdef INFINITY
+#ifdef USE_INFINITY
s[1] = 32767;
for( i=2; i<NI-1; i++ )
s[i] = 0;
int i, lost, j, k;
long lt, lta, ltb;
-#ifdef INFINITY
+#ifdef USE_INFINITY
if( eisinf(a) )
{
emov(a,c);
}
#endif
/* Infinity over anything else is infinity. */
-#ifdef INFINITY
+#ifdef USE_INFINITY
if( eisinf(b) )
{
if( eisneg(a) ^ eisneg(b) )
}
#endif
/* Infinity times anything else is infinity. */
-#ifdef INFINITY
+#ifdef USE_INFINITY
if( eisinf(a) || eisinf(b) )
{
if( eisneg(a) ^ eisneg(b) )
if( r & 0x8000 )
yy[0] = 0xffff;
r &= 0x7fff;
-#ifdef INFINITY
+#ifdef USE_INFINITY
if( r == 0x7fff )
{
#ifdef NANS
return;
}
#endif
-#ifdef INFINITY
+#ifdef USE_INFINITY
/* Point to the exponent field. */
p = &yy[NE-1];
if( (*p & 0x7fff) == 0x7fff )
eneg(y);
return;
}
-#endif /* INFINITY */
+#endif /* USE_INFINITY */
p = yy;
q = y;
for( i=0; i<NE; i++ )
for( i=0; i<4; i++ )
*q++ = *p++;
#else
-#ifdef INFINITY
+#ifdef USE_INFINITY
#ifdef IBMPC
if (eiisinf (a))
{
return;
}
#endif /* IBMPC */
-#endif /* INFINITY */
+#endif /* USE_INFINITY */
for( i=0; i<4; i++ )
*q-- = *p++;
#endif
yy[0] = 0xffff;
yy[M] = (r & 0x0f) | 0x10;
r &= ~0x800f; /* strip sign and 4 significand bits */
-#ifdef INFINITY
+#ifdef USE_INFINITY
if( r == 0x7ff0 )
{
#ifdef NANS
i = *p++;
if( i >= (unsigned int )2047 )
{ /* Saturate at largest number less than infinity. */
-#ifdef INFINITY
+#ifdef USE_INFINITY
*y |= 0x7ff0;
#ifdef IBMPC
*(--y) = 0;
*y++ = 0;
*y++ = 0;
#endif /* IBMPC */
-#else /* !INFINITY */
+#else /* !USE_INFINITY */
*y |= (unsigned short )0x7fef;
#ifdef IBMPC
*(--y) = 0xffff;
*y++ = 0xffff;
*y++ = 0xffff;
#endif
-#endif /* !INFINITY */
+#endif /* !USE_INFINITY */
return;
}
if( i == 0 )
yy[0] = 0xffff;
yy[M] = (r & 0x7f) | 0200;
r &= ~0x807f; /* strip sign and 7 significand bits */
-#ifdef INFINITY
+#ifdef USE_INFINITY
if( r == 0x7f80 )
{
#ifdef NANS
i = *p++;
if( i >= 255 )
{ /* Saturate at largest number less than infinity. */
-#ifdef INFINITY
+#ifdef USE_INFINITY
*y |= (unsigned short )0x7f80;
#ifdef IBMPC
*(--y) = 0;
++y;
*y = 0;
#endif
-#else /* !INFINITY */
+#else /* !USE_INFINITY */
*y |= (unsigned short )0x7f7f;
#ifdef IBMPC
*(--y) = 0xffff;
++y;
*y = 0xffff;
#endif
-#endif /* !INFINITY */
+#endif /* !USE_INFINITY */
return;
}
if( i == 0 )
{
unsigned short e[NI];
char *s, *p;
-int k;
+int i, j, k;
+int orig_ndigits;
LDPARMS rnd;
LDPARMS *ldp = &rnd;
char *outstr;
+char outbuf[NDEC + MAX_EXP_DIGITS + 10];
+union uconv du;
+du.d = d;
+orig_ndigits = ndigits;
rnd.rlast = -1;
rnd.rndprc = NBITS;
}
#if LDBL_MANT_DIG == 24
-e24toe( (unsigned short *)&d, e, ldp );
+e24toe( &du.pe, e, ldp );
#elif LDBL_MANT_DIG == 53
-e53toe( (unsigned short *)&d, e, ldp );
+e53toe( &du.pe, e, ldp );
#elif LDBL_MANT_DIG == 64
-e64toe( (unsigned short *)&d, e, ldp );
+e64toe( &du.pe, e, ldp );
#else
-e113toe( (unsigned short *)&d, e, ldp );
+e113toe( &du.pe, e, ldp );
#endif
if( eisneg(e) )
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;
-/* reentrancy addition to use mprec storage pool */
-_REENT_MP_RESULT(ptr) = Balloc (ptr, 3);
-_REENT_MP_RESULT_K(ptr) = 3;
-outstr = (char *)_REENT_MP_RESULT(ptr);
-
-etoasc( e, outstr, ndigits, mode, ldp );
-s = outstr;
+etoasc( e, outbuf, ndigits, mode, ldp );
+s = outbuf;
if( eisinf(e) || eisnan(e) )
{
*decpt = 9999;
/* Transform the string returned by etoasc into what the caller wants. */
/* Look for decimal point and delete it from the string. */
-s = outstr;
+s = outbuf;
while( *s != '\0' )
{
if( *s == '.' )
nodecpt:
/* Back up over the exponent field. */
-while( *s != 'E' && s > outstr)
+while( *s != 'E' && s > outbuf)
--s;
*s = '\0';
stripspaces:
/* Strip leading spaces and sign. */
-p = outstr;
+p = outbuf;
while( *p == ' ' || *p == '-')
++p;
/* Find new end of string. */
-s = outstr;
+s = outbuf;
while( (*s++ = *p++) != '\0' )
;
--s;
else
k = ldp->outexpon;
-while( *(s-1) == '0' && ((s - outstr) > k))
+while( *(s-1) == '0' && ((s - outbuf) > 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 + ldp->outexpon) < 0))
{
- s = outstr;
+ s = outbuf;
*s = '\0';
*decpt = 0;
}
+/* reentrancy addition to use mprec storage pool */
+/* we want to have enough space to hold the formatted result */
+
+if (mode == 3) /* f format, account for sign + dec digits + decpt + frac */
+ i = *decpt + orig_ndigits + 3;
+else /* account for sign + max precision digs + E + exp sign + exponent */
+ i = orig_ndigits + MAX_EXP_DIGITS + 4;
+
+j = sizeof (__ULong);
+for (_REENT_MP_RESULT_K(ptr) = 0; sizeof (_Bigint) - sizeof (__ULong) + j <= i; j <<= 1)
+ _REENT_MP_RESULT_K(ptr)++;
+_REENT_MP_RESULT(ptr) = Balloc (ptr, _REENT_MP_RESULT_K(ptr));
+
+/* Copy from internal temporary buffer to permanent buffer. */
+outstr = (char *)_REENT_MP_RESULT(ptr);
+strcpy (outstr, outbuf);
+
if( rve )
- *rve = s;
+ *rve = outstr + (s - outbuf);
+
return outstr;
}
LDPARMS rnd;
LDPARMS *ldp = &rnd;
+union uconv du;
+
rnd.rlast = -1;
rnd.rndprc = NBITS;
-
+du.d = *d;
#if LDBL_MANT_DIG == 24
-e24toe( (unsigned short *)d, e, ldp );
+e24toe( &du.pe, e, ldp );
#elif LDBL_MANT_DIG == 53
-e53toe( (unsigned short *)d, e, ldp );
+e53toe( &du.pe, e, ldp );
#elif LDBL_MANT_DIG == 64
-e64toe( (unsigned short *)d, e, ldp );
+e64toe( &du.pe, e, ldp );
#else
-e113toe( (unsigned short *)d, e, ldp );
+e113toe( &du.pe, e, ldp );
#endif
if( (e[NE-1] & 0x7fff) == 0x7fff )
emovo( y, t, ldp );
if( ecmp(t,ezero) != 0 )
goto roun; /* round to nearest */
- if( (*(s-1) & 1) == 0 )
+ if( ndigs < 0 || (*(s-1-(*(s-1)=='.')) & 1) == 0 )
goto doexp; /* round to even */
}
/* Round up and propagate carry-outs */
long double _strtold (char *s, char **se)
{
- long double x;
+ union uconv x;
LDPARMS rnd;
LDPARMS *ldp = &rnd;
int lenldstr;
rnd.rlast = -1;
rnd.rndprc = NBITS;
- lenldstr = asctoeg( s, (unsigned short *)&x, LDBL_MANT_DIG, ldp );
+ lenldstr = asctoeg( s, &x.pe, LDBL_MANT_DIG, ldp );
if (se)
*se = s + lenldstr;
- return x;
+ return x.d;
}
#define REASONABLE_LEN 200
for (i=0; i < n; i++)
*nan++ = *p++;
}
-
-
-
-
OSDN Git Service
