libs/newlib-cygwin
libs/newlib-cygwin
4
0
Fork 0
mirror of git://sourceware.org/git/newlib-cygwin.git synced 2025-01-19 04:49:25 +08:00
Code Issues Projects Releases Wiki Activity
newlib-cygwin/newlib/libc/stdlib/ldtoa.c
Corinna Vinschen 4d90e53359 ldtoa: fix dropping too many digits from output 
ldtoa cuts the number of digits it returns based on a computation of
number of supported bits (144) divide by log10(2).  Not only is the
integer approximation of log10(2) ~= 8/27 missing a digit here, it
also fails to take really small double and long double values into
account.

Allow for the full potential precision of long double values.  At the
same time, change the local string array allocation to request only as
much bytes as necessary to support the caller-requested number of
digits, to keep the stack size low on small targets.

In the long run a better fix would be to switch to gdtoa, as the BSD
variants, as well as Mingw64 do.

Signed-off-by: Corinna Vinschen <corinna@vinschen.de>
2021-11-04 13:14:17 +01:00

3879 lines
77 KiB
C
Raw Blame History

This file contains invisible Unicode characters

This file contains invisible Unicode characters that are indistinguishable to humans but may be processed differently by a computer. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

/* Extended precision arithmetic functions for long double I/O.
* This program has been placed in the public domain.
*/
#include <_ansi.h>
#include <reent.h>
#include <string.h>
#include <stdlib.h>
#include "mprec.h"
/* These are the externally visible entries. */
/* linux name: long double _IO_strtold (char *, char **); */
long double _strtold (char *, char **);
char *_ldtoa_r (struct _reent *, long double, int, int, int *, int *,
char **);
int _ldcheck (long double *);
#if 0
void _IO_ldtostr (long double *, char *, int, int, char);
#endif
/* Number of 16 bit words in external x type format */
#define NE 10
/* Number of 16 bit words in internal format */
#define NI (NE+3)
/* Array offset to exponent */
#define E 1
/* Array offset to high guard word */
#define M 2
/* Number of bits of precision */
#define NBITS ((NI-4)*16)
/* Maximum number of decimal digits in ASCII conversion
* Take full possible size of output into account
*/
#define NDEC 1023
/* The exponent of 1.0 */
#define EXONE (0x3fff)
/* 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
{
int rlast;
int rndprc;
int rw;
int re;
int outexpon;
unsigned short rmsk;
unsigned short rmbit;
unsigned short rebit;
unsigned short rbit[NI];
unsigned short equot[NI];
} LDPARMS;
static void esub (const short unsigned int *a, const short unsigned int *b,
short unsigned int *c, LDPARMS * ldp);
static void emul (const short unsigned int *a, const short unsigned int *b,
short unsigned int *c, LDPARMS * ldp);
static void ediv (const short unsigned int *a, const short unsigned int *b,
short unsigned int *c, LDPARMS * ldp);
static int ecmp (const short unsigned int *a, const short unsigned int *b);
static int enormlz (short unsigned int *x);
static int eshift (short unsigned int *x, int sc);
static void eshup1 (register short unsigned int *x);
static void eshup8 (register short unsigned int *x);
static void eshup6 (register short unsigned int *x);
static void eshdn1 (register short unsigned int *x);
static void eshdn8 (register short unsigned int *x);
static void eshdn6 (register short unsigned int *x);
static void eneg (short unsigned int *x);
static void emov (register const short unsigned int *a,
register short unsigned int *b);
static void eclear (register short unsigned int *x);
static void einfin (register short unsigned int *x, register LDPARMS * ldp);
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
static void e53toe (short unsigned int *pe, short unsigned int *y,
LDPARMS * ldp);
#elif LDBL_MANT_DIG == 64
static void e64toe (short unsigned int *pe, short unsigned int *y,
LDPARMS * ldp);
#else
static void e113toe (short unsigned int *pe, short unsigned int *y,
LDPARMS * ldp);
#endif
/* econst.c */
/* e type constants used by high precision check routines */
#if NE == 10
/* 0.0 */
static const unsigned short ezero[NE] = { 0x0000, 0x0000, 0x0000, 0x0000,
0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,
};
/* 1.0E0 */
static const unsigned short eone[NE] = { 0x0000, 0x0000, 0x0000, 0x0000,
0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3fff,
};
#else
/* 0.0 */
static const unsigned short ezero[NE] = {
0, 0000000, 0000000, 0000000, 0000000, 0000000,
};
/* 1.0E0 */
static const unsigned short eone[NE] = {
0, 0000000, 0000000, 0000000, 0100000, 0x3fff,
};
#endif
/* Debugging routine for displaying errors */
#ifdef DEBUG
/* Notice: the order of appearance of the following
* messages is bound to the error codes defined
* in mconf.h.
*/
static const char *const ermsg[7] = {
"unknown", /* error code 0 */
"domain", /* error code 1 */
"singularity", /* et seq. */
"overflow",
"underflow",
"total loss of precision",
"partial loss of precision"
};
#define mtherr(name, code) printf( "\n%s %s error\n", name, ermsg[code] );
#else
#define mtherr(name, code)
#endif
/* ieee.c
*
* Extended precision IEEE binary floating point arithmetic routines
*
* Numbers are stored in C language as arrays of 16-bit unsigned
* short integers. The arguments of the routines are pointers to
* the arrays.
*
*
* External e type data structure, simulates Intel 8087 chip
* temporary real format but possibly with a larger significand:
*
* NE-1 significand words (least significant word first,
* most significant bit is normally set)
* exponent (value = EXONE for 1.0,
* top bit is the sign)
*
*
* Internal data structure of a number (a "word" is 16 bits):
*
* ei[0] sign word (0 for positive, 0xffff for negative)
* ei[1] biased exponent (value = EXONE for the number 1.0)
* ei[2] high guard word (always zero after normalization)
* ei[3]
* to ei[NI-2] significand (NI-4 significand words,
* most significant word first,
* most significant bit is set)
* ei[NI-1] low guard word (0x8000 bit is rounding place)
*
*
*
* Routines for external format numbers
*
* asctoe( string, e ) ASCII string to extended double e type
* asctoe64( string, &d ) ASCII string to long double
* asctoe53( string, &d ) ASCII string to double
* asctoe24( string, &f ) ASCII string to single
* asctoeg( string, e, prec, ldp ) ASCII string to specified precision
* e24toe( &f, e, ldp ) IEEE single precision to e type
* e53toe( &d, e, ldp ) IEEE double precision to e type
* e64toe( &d, e, ldp ) IEEE long double precision to e type
* e113toe( &d, e, ldp ) IEEE long double precision to e type
* eabs(e) absolute value
* eadd( a, b, c ) c = b + a
* eclear(e) e = 0
* ecmp (a, b) Returns 1 if a > b, 0 if a == b,
* -1 if a < b, -2 if either a or b is a NaN.
* ediv( a, b, c, ldp ) c = b / a
* efloor( a, b, ldp ) truncate to integer, toward -infinity
* efrexp( a, exp, s ) extract exponent and significand
* eifrac( e, &l, frac ) e to long integer and e type fraction
* euifrac( e, &l, frac ) e to unsigned long integer and e type fraction
* einfin( e, ldp ) set e to infinity, leaving its sign alone
* eldexp( a, n, b ) multiply by 2**n
* emov( a, b ) b = a
* emul( a, b, c, ldp ) c = b * a
* eneg(e) e = -e
* eround( a, b ) b = nearest integer value to a
* esub( a, b, c, ldp ) c = b - a
* e24toasc( &f, str, n ) single to ASCII string, n digits after decimal
* e53toasc( &d, str, n ) double to ASCII string, n digits after decimal
* e64toasc( &d, str, n ) long double to ASCII string
* etoasc(e,str,n,fmt,ldp)e to ASCII string, n digits after decimal
* etoe24( e, &f ) convert e type to IEEE single precision
* etoe53( e, &d ) convert e type to IEEE double precision
* etoe64( e, &d ) convert e type to IEEE long double precision
* ltoe( &l, e ) long (32 bit) integer to e type
* ultoe( &l, e ) unsigned long (32 bit) integer to e type
* eisneg( e ) 1 if sign bit of e != 0, else 0
* eisinf( e ) 1 if e has maximum exponent (non-IEEE)
* or is infinite (IEEE)
* eisnan( e ) 1 if e is a NaN
* esqrt( a, b ) b = square root of a
*
*
* Routines for internal format numbers
*
* eaddm( ai, bi ) add significands, bi = bi + ai
* ecleaz(ei) ei = 0
* ecleazs(ei) set ei = 0 but leave its sign alone
* ecmpm( ai, bi ) compare significands, return 1, 0, or -1
* edivm( ai, bi, ldp ) divide significands, bi = bi / ai
* emdnorm(ai,l,s,exp,ldp) normalize and round off
* emovi( a, ai ) convert external a to internal ai
* emovo( ai, a, ldp ) convert internal ai to external a
* emovz( ai, bi ) bi = ai, low guard word of bi = 0
* emulm( ai, bi, ldp ) multiply significands, bi = bi * ai
* enormlz(ei) left-justify the significand
* eshdn1( ai ) shift significand and guards down 1 bit
* eshdn8( ai ) shift down 8 bits
* eshdn6( ai ) shift down 16 bits
* eshift( ai, n ) shift ai n bits up (or down if n < 0)
* eshup1( ai ) shift significand and guards up 1 bit
* eshup8( ai ) shift up 8 bits
* eshup6( ai ) shift up 16 bits
* esubm( ai, bi ) subtract significands, bi = bi - ai
*
*
* The result is always normalized and rounded to NI-4 word precision
* after each arithmetic operation.
*
* Exception flags are NOT fully supported.
*
* 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
* arithmetic will never produce a NaN output, and might be confused
* by a NaN input.
* If NaN's are supported, the output of ecmp(a,b) is -2 if
* either a or b is a NaN. This means asking if(ecmp(a,b) < 0)
* may not be legitimate. Use if(ecmp(a,b) == -1) for less-than
* if in doubt.
* Signaling NaN's are NOT supported; they are treated the same
* as quiet NaN's.
*
* Denormals are always supported here where appropriate (e.g., not
* for conversion to DEC numbers).
*/
/*
* 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
* 22 Nov 00 Revised to fit into newlib by Jeff Johnston <jjohnstn@redhat.com>
*
* Author: S. L. Moshier.
*
* Copyright (c) 1984,2000 S.L. Moshier
*
* Permission to use, copy, modify, and distribute this software for any
* purpose without fee is hereby granted, provided that this entire notice
* is included in all copies of any software which is or includes a copy
* or modification of this software and in all copies of the supporting
* documentation for such software.
*
* THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
* WARRANTY. IN PARTICULAR, THE AUTHOR MAKES NO REPRESENTATION
* OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY OF THIS
* SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
*
*/
#include <stdio.h>
/* #include "\usr\include\stdio.h" */
/*#include "ehead.h"*/
/*#include "mconf.h"*/
/* mconf.h
*
* Common include file for math routines
*
*
*
* SYNOPSIS:
*
* #include "mconf.h"
*
*
*
* DESCRIPTION:
*
* This file contains definitions for error codes that are
* passed to the common error handling routine mtherr()
* (which see).
*
* The file also includes a conditional assembly definition
* for the type of computer arithmetic (IEEE, DEC, Motorola
* IEEE, or UNKnown).
*
* For Digital Equipment PDP-11 and VAX computers, certain
* IBM systems, and others that use numbers with a 56-bit
* significand, the symbol DEC should be defined. In this
* mode, most floating point constants are given as arrays
* of octal integers to eliminate decimal to binary conversion
* errors that might be introduced by the compiler.
*
* For computers, such as IBM PC, that follow the IEEE
* Standard for Binary Floating Point Arithmetic (ANSI/IEEE
* Std 754-1985), the symbol IBMPC should be defined. These
* numbers have 53-bit significands. In this mode, constants
* are provided as arrays of hexadecimal 16 bit integers.
*
* To accommodate other types of computer arithmetic, all
* constants are also provided in a normal decimal radix
* which one can hope are correctly converted to a suitable
* format by the available C language compiler. To invoke
* this mode, the symbol UNK is defined.
*
* An important difference among these modes is a predefined
* set of machine arithmetic constants for each. The numbers
* MACHEP (the machine roundoff error), MAXNUM (largest number
* represented), and several other parameters are preset by
* the configuration symbol. Check the file const.c to
* ensure that these values are correct for your computer.
*
* For ANSI C compatibility, define ANSIC equal to 1. Currently
* this affects only the atan2() function and others that use it.
*/
/* Constant definitions for math error conditions
*/
#define DOMAIN 1 /* argument domain error */
#define SING 2 /* argument singularity */
#define OVERFLOW 3 /* overflow range error */
#define UNDERFLOW 4 /* underflow range error */
#define TLOSS 5 /* total loss of precision */
#define PLOSS 6 /* partial loss of precision */
#define EDOM 33
#define ERANGE 34
typedef struct
{
double r;
double i;
} cmplx;
/* Type of computer arithmetic */
#ifndef DEC
#ifdef __IEEE_LITTLE_ENDIAN
#define IBMPC 1
#else /* !__IEEE_LITTLE_ENDIAN */
#define MIEEE 1
#endif /* !__IEEE_LITTLE_ENDIAN */
#endif /* !DEC */
/* Define 1 for ANSI C atan2() function
* See atan.c and clog.c.
*/
#define ANSIC 1
/*define VOLATILE volatile*/
#define VOLATILE
#define NANS
#define USE_INFINITY
/* NaN's require infinity support. */
#ifdef NANS
#ifndef INFINITY
#define USE_INFINITY
#endif
#endif
/* This handles 64-bit long ints. */
#define LONGBITS (8 * sizeof(long))
static void eaddm (short unsigned int *x, short unsigned int *y);
static void esubm (short unsigned int *x, short unsigned int *y);
static void emdnorm (short unsigned int *s, int lost, int subflg,
long int exp, int rcntrl, LDPARMS * ldp);
#if 0 /* Broken, unusable implementation of strtold */
static int asctoeg (char *ss, short unsigned int *y, int oprec,
LDPARMS * ldp);
#endif
static void enan (short unsigned int *nan, int size);
#if LDBL_MANT_DIG == 24
static void toe24 (short unsigned int *x, short unsigned int *y);
#elif LDBL_MANT_DIG == 53
static void toe53 (short unsigned int *x, short unsigned int *y);
#elif LDBL_MANT_DIG == 64
static void toe64 (short unsigned int *a, short unsigned int *b);
#else
static void toe113 (short unsigned int *a, short unsigned int *b);
#endif
static void eiremain (short unsigned int *den, short unsigned int *num,
LDPARMS * ldp);
static int ecmpm (register short unsigned int *a,
register short unsigned int *b);
static int edivm (short unsigned int *den, short unsigned int *num,
LDPARMS * ldp);
static int emulm (short unsigned int *a, short unsigned int *b,
LDPARMS * ldp);
static int eisneg (const short unsigned int *x);
static int eisinf (const short unsigned int *x);
static void emovi (const short unsigned int *a, short unsigned int *b);
static void emovo (short unsigned int *a, short unsigned int *b,
LDPARMS * ldp);
static void emovz (register short unsigned int *a,
register short unsigned int *b);
static void ecleaz (register short unsigned int *xi);
static void eadd1 (const short unsigned int *a, const short unsigned int *b,
short unsigned int *c, int subflg, LDPARMS * ldp);
static int eisnan (const short unsigned int *x);
static int eiisnan (short unsigned int *x);
#ifdef DEC
static void etodec (), todec (), dectoe ();
#endif
/*
; Clear out entire external format number.
;
; unsigned short x[];
; eclear( x );
*/
static void
eclear (register short unsigned int *x)
{
register int i;
for (i = 0; i < NE; i++)
*x++ = 0;
}
/* Move external format number from a to b.
*
* emov( a, b );
*/
static void
emov (register const short unsigned int *a, register short unsigned int *b)
{
register int i;
for (i = 0; i < NE; i++)
*b++ = *a++;
}
/*
; Negate external format number
;
; unsigned short x[NE];
; eneg( x );
*/
static void
eneg (short unsigned int *x)
{
#ifdef NANS
if (eisnan (x))
return;
#endif
x[NE - 1] ^= 0x8000; /* Toggle the sign bit */
}
/* Return 1 if external format number is negative,
* else return zero.
*/
static int
eisneg (const short unsigned int *x)
{
#ifdef NANS
if (eisnan (x))
return (0);
#endif
if (x[NE - 1] & 0x8000)
return (1);
else
return (0);
}
/* Return 1 if external format number has maximum possible exponent,
* else return zero.
*/
static int
eisinf (const short unsigned int *x)
{
if ((x[NE - 1] & 0x7fff) == 0x7fff)
{
#ifdef NANS
if (eisnan (x))
return (0);
#endif
return (1);
}
else
return (0);
}
/* Check if e-type number is not a number.
*/
static int
eisnan (const short unsigned int *x)
{
#ifdef NANS
int i;
/* NaN has maximum exponent */
if ((x[NE - 1] & 0x7fff) != 0x7fff)
return (0);
/* ... and non-zero significand field. */
for (i = 0; i < NE - 1; i++)
{
if (*x++ != 0)
return (1);
}
#endif
return (0);
}
/*
; Fill entire number, including exponent and significand, with
; largest possible number. These programs implement a saturation
; value that is an ordinary, legal number. A special value
; "infinity" may also be implemented; this would require tests
; for that value and implementation of special rules for arithmetic
; operations involving inifinity.
*/
static void
einfin (register short unsigned int *x, register LDPARMS * ldp)
{
register int i;
#ifdef USE_INFINITY
for (i = 0; i < NE - 1; i++)
*x++ = 0;
*x |= 32767;
ldp = ldp;
#else
for (i = 0; i < NE - 1; i++)
*x++ = 0xffff;
*x |= 32766;
if (ldp->rndprc < NBITS)
{
if (ldp->rndprc == 113)
{
*(x - 9) = 0;
*(x - 8) = 0;
}
if (ldp->rndprc == 64)
{
*(x - 5) = 0;
}
if (ldp->rndprc == 53)
{
*(x - 4) = 0xf800;
}
else
{
*(x - 4) = 0;
*(x - 3) = 0;
*(x - 2) = 0xff00;
}
}
#endif
}
/* Move in external format number,
* converting it to internal format.
*/
static void
emovi (const short unsigned int *a, short unsigned int *b)
{
register const unsigned short *p;
register unsigned short *q;
int i;
q = b;
p = a + (NE - 1); /* point to last word of external number */
/* get the sign bit */
if (*p & 0x8000)
*q++ = 0xffff;
else
*q++ = 0;
/* get the exponent */
*q = *p--;
*q++ &= 0x7fff; /* delete the sign bit */
#ifdef USE_INFINITY
if ((*(q - 1) & 0x7fff) == 0x7fff)
{
#ifdef NANS
if (eisnan (a))
{
*q++ = 0;
for (i = 3; i < NI; i++)
*q++ = *p--;
return;
}
#endif
for (i = 2; i < NI; i++)
*q++ = 0;
return;
}
#endif
/* clear high guard word */
*q++ = 0;
/* move in the significand */
for (i = 0; i < NE - 1; i++)
*q++ = *p--;
/* clear low guard word */
*q = 0;
}
/* Move internal format number out,
* converting it to external format.
*/
static void
emovo (short unsigned int *a, short unsigned int *b, LDPARMS * ldp)
{
register unsigned short *p, *q;
unsigned short i;
p = a;
q = b + (NE - 1); /* point to output exponent */
/* combine sign and exponent */
i = *p++;
if (i)
*q-- = *p++ | 0x8000;
else
*q-- = *p++;
#ifdef USE_INFINITY
if (*(p - 1) == 0x7fff)
{
#ifdef NANS
if (eiisnan (a))
{
enan (b, NBITS);
return;
}
#endif
einfin (b, ldp);
return;
}
#endif
/* skip over guard word */
++p;
/* move the significand */
for (i = 0; i < NE - 1; i++)
*q-- = *p++;
}
/* Clear out internal format number.
*/
static void
ecleaz (register short unsigned int *xi)
{
register int i;
for (i = 0; i < NI; i++)
*xi++ = 0;
}
/* same, but don't touch the sign. */
static void
ecleazs (register short unsigned int *xi)
{
register int i;
++xi;
for (i = 0; i < NI - 1; i++)
*xi++ = 0;
}
/* Move internal format number from a to b.
*/
static void
emovz (register short unsigned int *a, register short unsigned int *b)
{
register int i;
for (i = 0; i < NI - 1; i++)
*b++ = *a++;
/* clear low guard word */
*b = 0;
}
/* Return nonzero if internal format number is a NaN.
*/
static int
eiisnan (short unsigned int *x)
{
int i;
if ((x[E] & 0x7fff) == 0x7fff)
{
for (i = M + 1; i < NI; i++)
{
if (x[i] != 0)
return (1);
}
}
return (0);
}
#if LDBL_MANT_DIG == 64
/* Return nonzero if internal format number is infinite. */
static int
eiisinf (unsigned short x[])
{
#ifdef NANS
if (eiisnan (x))
return (0);
#endif
if ((x[E] & 0x7fff) == 0x7fff)
return (1);
return (0);
}
#endif /* LDBL_MANT_DIG == 64 */
/*
; Compare significands of numbers in internal format.
; Guard words are included in the comparison.
;
; unsigned short a[NI], b[NI];
; cmpm( a, b );
;
; for the significands:
; returns +1 if a > b
; 0 if a == b
; -1 if a < b
*/
static int
ecmpm (register short unsigned int *a, register short unsigned int *b)
{
int i;
a += M; /* skip up to significand area */
b += M;
for (i = M; i < NI; i++)
{
if (*a++ != *b++)
goto difrnt;
}
return (0);
difrnt:
if (*(--a) > *(--b))
return (1);
else
return (-1);
}
/*
; Shift significand down by 1 bit
*/
static void
eshdn1 (register short unsigned int *x)
{
register unsigned short bits;
int i;
x += M; /* point to significand area */
bits = 0;
for (i = M; i < NI; i++)
{
if (*x & 1)
bits |= 1;
*x >>= 1;
if (bits & 2)
*x |= 0x8000;
bits <<= 1;
++x;
}
}
/*
; Shift significand up by 1 bit
*/
static void
eshup1 (register short unsigned int *x)
{
register unsigned short bits;
int i;
x += NI - 1;
bits = 0;
for (i = M; i < NI; i++)
{
if (*x & 0x8000)
bits |= 1;
*x <<= 1;
if (bits & 2)
*x |= 1;
bits <<= 1;
--x;
}
}
/*
; Shift significand down by 8 bits
*/
static void
eshdn8 (register short unsigned int *x)
{
register unsigned short newbyt, oldbyt;
int i;
x += M;
oldbyt = 0;
for (i = M; i < NI; i++)
{
newbyt = *x << 8;
*x >>= 8;
*x |= oldbyt;
oldbyt = newbyt;
++x;
}
}
/*
; Shift significand up by 8 bits
*/
static void
eshup8 (register short unsigned int *x)
{
int i;
register unsigned short newbyt, oldbyt;
x += NI - 1;
oldbyt = 0;
for (i = M; i < NI; i++)
{
newbyt = *x >> 8;
*x <<= 8;
*x |= oldbyt;
oldbyt = newbyt;
--x;
}
}
/*
; Shift significand up by 16 bits
*/
static void
eshup6 (register short unsigned int *x)
{
int i;
register unsigned short *p;
p = x + M;
x += M + 1;
for (i = M; i < NI - 1; i++)
*p++ = *x++;
*p = 0;
}
/*
; Shift significand down by 16 bits
*/
static void
eshdn6 (register short unsigned int *x)
{
int i;
register unsigned short *p;
x += NI - 1;
p = x + 1;
for (i = M; i < NI - 1; i++)
*(--p) = *(--x);
*(--p) = 0;
}
/*
; Add significands
; x + y replaces y
*/
static void
eaddm (short unsigned int *x, short unsigned int *y)
{
register unsigned long a;
int i;
unsigned int carry;
x += NI - 1;
y += NI - 1;
carry = 0;
for (i = M; i < NI; i++)
{
a = (unsigned long) (*x) + (unsigned long) (*y) + carry;
if (a & 0x10000)
carry = 1;
else
carry = 0;
*y = (unsigned short) a;
--x;
--y;
}
}
/*
; Subtract significands
; y - x replaces y
*/
static void
esubm (short unsigned int *x, short unsigned int *y)
{
unsigned long a;
int i;
unsigned int carry;
x += NI - 1;
y += NI - 1;
carry = 0;
for (i = M; i < NI; i++)
{
a = (unsigned long) (*y) - (unsigned long) (*x) - carry;
if (a & 0x10000)
carry = 1;
else
carry = 0;
*y = (unsigned short) a;
--x;
--y;
}
}
/* Divide significands */
/* Multiply significand of e-type number b
by 16-bit quantity a, e-type result to c. */
static void
m16m (short unsigned int a, short unsigned int *b, short unsigned int *c)
{
register unsigned short *pp;
register unsigned long carry;
unsigned short *ps;
unsigned short p[NI];
unsigned long aa, m;
int i;
aa = a;
pp = &p[NI - 2];
*pp++ = 0;
*pp = 0;
ps = &b[NI - 1];
for (i = M + 1; i < NI; i++)
{
if (*ps == 0)
{
--ps;
--pp;
*(pp - 1) = 0;
}
else
{
m = (unsigned long) aa **ps--;
carry = (m & 0xffff) + *pp;
*pp-- = (unsigned short) carry;
carry = (carry >> 16) + (m >> 16) + *pp;
*pp = (unsigned short) carry;
*(pp - 1) = carry >> 16;
}
}
for (i = M; i < NI; i++)
c[i] = p[i];
}
/* Divide significands. Neither the numerator nor the denominator
is permitted to have its high guard word nonzero. */
static int
edivm (short unsigned int *den, short unsigned int *num, LDPARMS * ldp)
{
int i;
register unsigned short *p;
unsigned long tnum;
unsigned short j, tdenm, tquot;
unsigned short tprod[NI + 1];
unsigned short *equot = ldp->equot;
p = &equot[0];
*p++ = num[0];
*p++ = num[1];
for (i = M; i < NI; i++)
{
*p++ = 0;
}
eshdn1 (num);
tdenm = den[M + 1];
for (i = M; i < NI; i++)
{
/* Find trial quotient digit (the radix is 65536). */
tnum = (((unsigned long) num[M]) << 16) + num[M + 1];
/* Do not execute the divide instruction if it will overflow. */
if ((tdenm * 0xffffUL) < tnum)
tquot = 0xffff;
else
tquot = tnum / tdenm;
/* Prove that the divide worked. */
/*
tcheck = (unsigned long )tquot * tdenm;
if( tnum - tcheck > tdenm )
tquot = 0xffff;
*/
/* Multiply denominator by trial quotient digit. */
m16m (tquot, den, tprod);
/* The quotient digit may have been overestimated. */
if (ecmpm (tprod, num) > 0)
{
tquot -= 1;
esubm (den, tprod);
if (ecmpm (tprod, num) > 0)
{
tquot -= 1;
esubm (den, tprod);
}
}
/*
if( ecmpm( tprod, num ) > 0 )
{
eshow( "tprod", tprod );
eshow( "num ", num );
printf( "tnum = %08lx, tden = %04x, tquot = %04x\n",
tnum, den[M+1], tquot );
}
*/
esubm (tprod, num);
/*
if( ecmpm( num, den ) >= 0 )
{
eshow( "num ", num );
eshow( "den ", den );
printf( "tnum = %08lx, tden = %04x, tquot = %04x\n",
tnum, den[M+1], tquot );
}
*/
equot[i] = tquot;
eshup6 (num);
}
/* test for nonzero remainder after roundoff bit */
p = &num[M];
j = 0;
for (i = M; i < NI; i++)
{
j |= *p++;
}
if (j)
j = 1;
for (i = 0; i < NI; i++)
num[i] = equot[i];
return ((int) j);
}
/* Multiply significands */
static int
emulm (short unsigned int *a, short unsigned int *b, LDPARMS * ldp)
{
unsigned short *p, *q;
unsigned short pprod[NI];
unsigned short j;
int i;
unsigned short *equot = ldp->equot;
equot[0] = b[0];
equot[1] = b[1];
for (i = M; i < NI; i++)
equot[i] = 0;
j = 0;
p = &a[NI - 1];
q = &equot[NI - 1];
for (i = M + 1; i < NI; i++)
{
if (*p == 0)
{
--p;
}
else
{
m16m (*p--, b, pprod);
eaddm (pprod, equot);
}
j |= *q;
eshdn6 (equot);
}
for (i = 0; i < NI; i++)
b[i] = equot[i];
/* return flag for lost nonzero bits */
return ((int) j);
}
/*
static void eshow(str, x)
char *str;
unsigned short *x;
{
int i;
printf( "%s ", str );
for( i=0; i<NI; i++ )
printf( "%04x ", *x++ );
printf( "\n" );
}
*/
/*
* Normalize and round off.
*
* The internal format number to be rounded is "s".
* Input "lost" indicates whether the number is exact.
* This is the so-called sticky bit.
*
* Input "subflg" indicates whether the number was obtained
* by a subtraction operation. In that case if lost is nonzero
* then the number is slightly smaller than indicated.
*
* Input "exp" is the biased exponent, which may be negative.
* the exponent field of "s" is ignored but is replaced by
* "exp" as adjusted by normalization and rounding.
*
* Input "rcntrl" is the rounding control.
*/
static void
emdnorm (short unsigned int *s, int lost, int subflg, long int exp,
int rcntrl, LDPARMS * ldp)
{
int i, j;
unsigned short r;
/* Normalize */
j = enormlz (s);
/* a blank significand could mean either zero or infinity. */
#ifndef USE_INFINITY
if (j > NBITS)
{
ecleazs (s);
return;
}
#endif
exp -= j;
#ifndef USE_INFINITY
if (exp >= 32767L)
goto overf;
#else
if ((j > NBITS) && (exp < 32767L))
{
ecleazs (s);
return;
}
#endif
if (exp < 0L)
{
if (exp > (long) (-NBITS - 1))
{
j = (int) exp;
i = eshift (s, j);
if (i)
lost = 1;
}
else
{
ecleazs (s);
return;
}
}
/* Round off, unless told not to by rcntrl. */
if (rcntrl == 0)
goto mdfin;
/* Set up rounding parameters if the control register changed. */
if (ldp->rndprc != ldp->rlast)
{
ecleaz (ldp->rbit);
switch (ldp->rndprc)
{
default:
case NBITS:
ldp->rw = NI - 1; /* low guard word */
ldp->rmsk = 0xffff;
ldp->rmbit = 0x8000;
ldp->rebit = 1;
ldp->re = ldp->rw - 1;
break;
case 113:
ldp->rw = 10;
ldp->rmsk = 0x7fff;
ldp->rmbit = 0x4000;
ldp->rebit = 0x8000;
ldp->re = ldp->rw;
break;
case 64:
ldp->rw = 7;
ldp->rmsk = 0xffff;
ldp->rmbit = 0x8000;
ldp->rebit = 1;
ldp->re = ldp->rw - 1;
break;
/* For DEC arithmetic */
case 56:
ldp->rw = 6;
ldp->rmsk = 0xff;
ldp->rmbit = 0x80;
ldp->rebit = 0x100;
ldp->re = ldp->rw;
break;
case 53:
ldp->rw = 6;
ldp->rmsk = 0x7ff;
ldp->rmbit = 0x0400;
ldp->rebit = 0x800;
ldp->re = ldp->rw;
break;
case 24:
ldp->rw = 4;
ldp->rmsk = 0xff;
ldp->rmbit = 0x80;
ldp->rebit = 0x100;
ldp->re = ldp->rw;
break;
}
ldp->rbit[ldp->re] = ldp->rebit;
ldp->rlast = ldp->rndprc;
}
/* Shift down 1 temporarily if the data structure has an implied
* most significant bit and the number is denormal.
* For rndprc = 64 or NBITS, there is no implied bit.
* But Intel long double denormals lose one bit of significance even so.
*/
#if IBMPC
if ((exp <= 0) && (ldp->rndprc != NBITS))
#else
if ((exp <= 0) && (ldp->rndprc != 64) && (ldp->rndprc != NBITS))
#endif
{
lost |= s[NI - 1] & 1;
eshdn1 (s);
}
/* Clear out all bits below the rounding bit,
* remembering in r if any were nonzero.
*/
r = s[ldp->rw] & ldp->rmsk;
if (ldp->rndprc < NBITS)
{
i = ldp->rw + 1;
while (i < NI)
{
if (s[i])
r |= 1;
s[i] = 0;
++i;
}
}
s[ldp->rw] &= ~ldp->rmsk;
if ((r & ldp->rmbit) != 0)
{
if (r == ldp->rmbit)
{
if (lost == 0)
{ /* round to even */
if ((s[ldp->re] & ldp->rebit) == 0)
goto mddone;
}
else
{
if (subflg != 0)
goto mddone;
}
}
eaddm (ldp->rbit, s);
}
mddone:
#if IBMPC
if ((exp <= 0) && (ldp->rndprc != NBITS))
#else
if ((exp <= 0) && (ldp->rndprc != 64) && (ldp->rndprc != NBITS))
#endif
{
eshup1 (s);
}
if (s[2] != 0)
{ /* overflow on roundoff */
eshdn1 (s);
exp += 1;
}
mdfin:
s[NI - 1] = 0;
if (exp >= 32767L)
{
#ifndef USE_INFINITY
overf:
#endif
#ifdef USE_INFINITY
s[1] = 32767;
for (i = 2; i < NI - 1; i++)
s[i] = 0;
#else
s[1] = 32766;
s[2] = 0;
for (i = M + 1; i < NI - 1; i++)
s[i] = 0xffff;
s[NI - 1] = 0;
if ((ldp->rndprc < 64) || (ldp->rndprc == 113))
{
s[ldp->rw] &= ~ldp->rmsk;
if (ldp->rndprc == 24)
{
s[5] = 0;
s[6] = 0;
}
}
#endif
return;
}
if (exp < 0)
s[1] = 0;
else
s[1] = (unsigned short) exp;
}
/*
; Subtract external format numbers.
;
; unsigned short a[NE], b[NE], c[NE];
; LDPARMS *ldp;
; esub( a, b, c, ldp ); c = b - a
*/
static void
esub (const short unsigned int *a, const short unsigned int *b,
short unsigned int *c, LDPARMS * ldp)
{
#ifdef NANS
if (eisnan (a))
{
emov (a, c);
return;
}
if (eisnan (b))
{
emov (b, c);
return;
}
/* Infinity minus infinity is a NaN.
* Test for subtracting infinities of the same sign.
*/
if (eisinf (a) && eisinf (b) && ((eisneg (a) ^ eisneg (b)) == 0))
{
mtherr ("esub", DOMAIN);
enan (c, NBITS);
return;
}
#endif
eadd1 (a, b, c, 1, ldp);
}
static void
eadd1 (const short unsigned int *a, const short unsigned int *b,
short unsigned int *c, int subflg, LDPARMS * ldp)
{
unsigned short ai[NI], bi[NI], ci[NI];
int i, lost, j, k;
long lt, lta, ltb;
#ifdef USE_INFINITY
if (eisinf (a))
{
emov (a, c);
if (subflg)
eneg (c);
return;
}
if (eisinf (b))
{
emov (b, c);
return;
}
#endif
emovi (a, ai);
emovi (b, bi);
if (subflg)
ai[0] = ~ai[0];
/* compare exponents */
lta = ai[E];
ltb = bi[E];
lt = lta - ltb;
if (lt > 0L)
{ /* put the larger number in bi */
emovz (bi, ci);
emovz (ai, bi);
emovz (ci, ai);
ltb = bi[E];
lt = -lt;
}
lost = 0;
if (lt != 0L)
{
if (lt < (long) (-NBITS - 1))
goto done; /* answer same as larger addend */
k = (int) lt;
lost = eshift (ai, k); /* shift the smaller number down */
}
else
{
/* exponents were the same, so must compare significands */
i = ecmpm (ai, bi);
if (i == 0)
{ /* the numbers are identical in magnitude */
/* if different signs, result is zero */
if (ai[0] != bi[0])
{
eclear (c);
return;
}
/* if same sign, result is double */
/* double denomalized tiny number */
if ((bi[E] == 0) && ((bi[3] & 0x8000) == 0))
{
eshup1 (bi);
goto done;
}
/* add 1 to exponent unless both are zero! */
for (j = 1; j < NI - 1; j++)
{
if (bi[j] != 0)
{
/* This could overflow, but let emovo take care of that. */
ltb += 1;
break;
}
}
bi[E] = (unsigned short) ltb;
goto done;
}
if (i > 0)
{ /* put the larger number in bi */
emovz (bi, ci);
emovz (ai, bi);
emovz (ci, ai);
}
}
if (ai[0] == bi[0])
{
eaddm (ai, bi);
subflg = 0;
}
else
{
esubm (ai, bi);
subflg = 1;
}
emdnorm (bi, lost, subflg, ltb, 64, ldp);
done:
emovo (bi, c, ldp);
}
/*
; Divide.
;
; unsigned short a[NE], b[NE], c[NE];
; LDPARMS *ldp;
; ediv( a, b, c, ldp ); c = b / a
*/
static void
ediv (const short unsigned int *a, const short unsigned int *b,
short unsigned int *c, LDPARMS * ldp)
{
unsigned short ai[NI], bi[NI];
int i;
long lt, lta, ltb;
#ifdef NANS
/* Return any NaN input. */
if (eisnan (a))
{
emov (a, c);
return;
}
if (eisnan (b))
{
emov (b, c);
return;
}
/* Zero over zero, or infinity over infinity, is a NaN. */
if (((ecmp (a, ezero) == 0) && (ecmp (b, ezero) == 0))
|| (eisinf (a) && eisinf (b)))
{
mtherr ("ediv", DOMAIN);
enan (c, NBITS);
return;
}
#endif
/* Infinity over anything else is infinity. */
#ifdef USE_INFINITY
if (eisinf (b))
{
if (eisneg (a) ^ eisneg (b))
*(c + (NE - 1)) = 0x8000;
else
*(c + (NE - 1)) = 0;
einfin (c, ldp);
return;
}
if (eisinf (a))
{
eclear (c);
return;
}
#endif
emovi (a, ai);
emovi (b, bi);
lta = ai[E];
ltb = bi[E];
if (bi[E] == 0)
{ /* See if numerator is zero. */
for (i = 1; i < NI - 1; i++)
{
if (bi[i] != 0)
{
ltb -= enormlz (bi);
goto dnzro1;
}
}
eclear (c);
return;
}
dnzro1:
if (ai[E] == 0)
{ /* possible divide by zero */
for (i = 1; i < NI - 1; i++)
{
if (ai[i] != 0)
{
lta -= enormlz (ai);
goto dnzro2;
}
}
if (ai[0] == bi[0])
*(c + (NE - 1)) = 0;
else
*(c + (NE - 1)) = 0x8000;
einfin (c, ldp);
mtherr ("ediv", SING);
return;
}
dnzro2:
i = edivm (ai, bi, ldp);
/* calculate exponent */
lt = ltb - lta + EXONE;
emdnorm (bi, i, 0, lt, 64, ldp);
/* set the sign */
if (ai[0] == bi[0])
bi[0] = 0;
else
bi[0] = 0Xffff;
emovo (bi, c, ldp);
}
/*
; Multiply.
;
; unsigned short a[NE], b[NE], c[NE];
; LDPARMS *ldp
; emul( a, b, c, ldp ); c = b * a
*/
static void
emul (const short unsigned int *a, const short unsigned int *b,
short unsigned int *c, LDPARMS * ldp)
{
unsigned short ai[NI], bi[NI];
int i, j;
long lt, lta, ltb;
#ifdef NANS
/* NaN times anything is the same NaN. */
if (eisnan (a))
{
emov (a, c);
return;
}
if (eisnan (b))
{
emov (b, c);
return;
}
/* Zero times infinity is a NaN. */
if ((eisinf (a) && (ecmp (b, ezero) == 0))
|| (eisinf (b) && (ecmp (a, ezero) == 0)))
{
mtherr ("emul", DOMAIN);
enan (c, NBITS);
return;
}
#endif
/* Infinity times anything else is infinity. */
#ifdef USE_INFINITY
if (eisinf (a) || eisinf (b))
{
if (eisneg (a) ^ eisneg (b))
*(c + (NE - 1)) = 0x8000;
else
*(c + (NE - 1)) = 0;
einfin (c, ldp);
return;
}
#endif
emovi (a, ai);
emovi (b, bi);
lta = ai[E];
ltb = bi[E];
if (ai[E] == 0)
{
for (i = 1; i < NI - 1; i++)
{
if (ai[i] != 0)
{
lta -= enormlz (ai);
goto mnzer1;
}
}
eclear (c);
return;
}
mnzer1:
if (bi[E] == 0)
{
for (i = 1; i < NI - 1; i++)
{
if (bi[i] != 0)
{
ltb -= enormlz (bi);
goto mnzer2;
}
}
eclear (c);
return;
}
mnzer2:
/* Multiply significands */
j = emulm (ai, bi, ldp);
/* calculate exponent */
lt = lta + ltb - (EXONE - 1);
emdnorm (bi, j, 0, lt, 64, ldp);
/* calculate sign of product */
if (ai[0] == bi[0])
bi[0] = 0;
else
bi[0] = 0xffff;
emovo (bi, c, ldp);
}
#if LDBL_MANT_DIG > 64
static void
e113toe (short unsigned int *pe, short unsigned int *y, LDPARMS * ldp)
{
register unsigned short r;
unsigned short *e, *p;
unsigned short yy[NI];
int denorm, i;
e = pe;
denorm = 0;
ecleaz (yy);
#ifdef IBMPC
e += 7;
#endif
r = *e;
yy[0] = 0;
if (r & 0x8000)
yy[0] = 0xffff;
r &= 0x7fff;
#ifdef USE_INFINITY
if (r == 0x7fff)
{
#ifdef NANS
#ifdef IBMPC
for (i = 0; i < 7; i++)
{
if (pe[i] != 0)
{
enan (y, NBITS);
return;
}
}
#else /* !IBMPC */
for (i = 1; i < 8; i++)
{
if (pe[i] != 0)
{
enan (y, NBITS);
return;
}
}
#endif /* !IBMPC */
#endif /* NANS */
eclear (y);
einfin (y, ldp);
if (*e & 0x8000)
eneg (y);
return;
}
#endif /* INFINITY */
yy[E] = r;
p = &yy[M + 1];
#ifdef IBMPC
for (i = 0; i < 7; i++)
*p++ = *(--e);
#else /* IBMPC */
++e;
for (i = 0; i < 7; i++)
*p++ = *e++;
#endif /* IBMPC */
/* If denormal, remove the implied bit; else shift down 1. */
if (r == 0)
{
yy[M] = 0;
}
else
{
yy[M] = 1;
eshift (yy, -1);
}
emovo (yy, y, ldp);
}
/* move out internal format to ieee long double */
static void
__attribute__ ((__unused__))
toe113 (short unsigned int *a, short unsigned int *b)
{
register unsigned short *p, *q;
unsigned short i;
#ifdef NANS
if (eiisnan (a))
{
enan (b, 113);
return;
}
#endif
p = a;
#ifdef MIEEE
q = b;
#else
q = b + 7; /* point to output exponent */
#endif
/* If not denormal, delete the implied bit. */
if (a[E] != 0)
{
eshup1 (a);
}
/* combine sign and exponent */
i = *p++;
#ifdef MIEEE
if (i)
*q++ = *p++ | 0x8000;
else
*q++ = *p++;
#else
if (i)
*q-- = *p++ | 0x8000;
else
*q-- = *p++;
#endif
/* skip over guard word */
++p;
/* move the significand */
#ifdef MIEEE
for (i = 0; i < 7; i++)
*q++ = *p++;
#else
for (i = 0; i < 7; i++)
*q-- = *p++;
#endif
}
#endif /* LDBL_MANT_DIG > 64 */
#if LDBL_MANT_DIG == 64
static void
e64toe (short unsigned int *pe, short unsigned int *y, LDPARMS * ldp)
{
unsigned short yy[NI];
unsigned short *p, *q, *e;
int i;
e = pe;
p = yy;
for (i = 0; i < NE - 5; i++)
*p++ = 0;
#ifdef IBMPC
for (i = 0; i < 5; i++)
*p++ = *e++;
#endif
#ifdef DEC
for (i = 0; i < 5; i++)
*p++ = *e++;
#endif
#ifdef MIEEE
p = &yy[0] + (NE - 1);
*p-- = *e++;
++e; /* MIEEE skips over 2nd short */
for (i = 0; i < 4; i++)
*p-- = *e++;
#endif
#ifdef IBMPC
/* For Intel long double, shift denormal significand up 1
-- but only if the top significand bit is zero. */
if ((yy[NE - 1] & 0x7fff) == 0 && (yy[NE - 2] & 0x8000) == 0)
{
unsigned short temp[NI + 1];
emovi (yy, temp);
eshup1 (temp);
emovo (temp, y, ldp);
return;
}
#endif
#ifdef USE_INFINITY
/* Point to the exponent field. */
p = &yy[NE - 1];
if ((*p & 0x7fff) == 0x7fff)
{
#ifdef NANS
#ifdef IBMPC
for (i = 0; i < 4; i++)
{
if ((i != 3 && pe[i] != 0)
/* Check for Intel long double infinity pattern. */
|| (i == 3 && pe[i] != 0x8000))
{
enan (y, NBITS);
return;
}
}
#endif
#ifdef MIEEE
for (i = 2; i <= 5; i++)
{
if (pe[i] != 0)
{
enan (y, NBITS);
return;
}
}
#endif
#endif /* NANS */
eclear (y);
einfin (y, ldp);
if (*p & 0x8000)
eneg (y);
return;
}
#endif /* USE_INFINITY */
p = yy;
q = y;
for (i = 0; i < NE; i++)
*q++ = *p++;
}
/* move out internal format to ieee long double */
static void
__attribute__ ((__unused__))
toe64 (short unsigned int *a, short unsigned int *b)
{
register unsigned short *p, *q;
unsigned short i;
#ifdef NANS
if (eiisnan (a))
{
enan (b, 64);
return;
}
#endif
#ifdef IBMPC
/* Shift Intel denormal significand down 1. */
if (a[E] == 0)
eshdn1 (a);
#endif
p = a;
#ifdef MIEEE
q = b;
#else
q = b + 4; /* point to output exponent */
/* NOTE: Intel data type is 96 bits wide, clear the last word here. */
*(q + 1) = 0;
#endif
/* combine sign and exponent */
i = *p++;
#ifdef MIEEE
if (i)
*q++ = *p++ | 0x8000;
else
*q++ = *p++;
*q++ = 0; /* leave 2nd short blank */
#else
if (i)
*q-- = *p++ | 0x8000;
else
*q-- = *p++;
#endif
/* skip over guard word */
++p;
/* move the significand */
#ifdef MIEEE
for (i = 0; i < 4; i++)
*q++ = *p++;
#else
#ifdef USE_INFINITY
#ifdef IBMPC
if (eiisinf (a))
{
/* Intel long double infinity. */
*q-- = 0x8000;
*q-- = 0;
*q-- = 0;
*q = 0;
return;
}
#endif /* IBMPC */
#endif /* USE_INFINITY */
for (i = 0; i < 4; i++)
*q-- = *p++;
#endif
}
#endif /* LDBL_MANT_DIG == 64 */
#if LDBL_MANT_DIG == 53
/*
; Convert IEEE double precision to e type
; double d;
; unsigned short x[N+2];
; e53toe( &d, x );
*/
static void
e53toe (short unsigned int *pe, short unsigned int *y, LDPARMS * ldp)
{
#ifdef DEC
dectoe (pe, y); /* see etodec.c */
#else
register unsigned short r;
register unsigned short *p, *e;
unsigned short yy[NI];
int denorm, k;
e = pe;
denorm = 0; /* flag if denormalized number */
ecleaz (yy);
#ifdef IBMPC
e += 3;
#endif
#ifdef DEC
e += 3;
#endif
r = *e;
yy[0] = 0;
if (r & 0x8000)
yy[0] = 0xffff;
yy[M] = (r & 0x0f) | 0x10;
r &= ~0x800f; /* strip sign and 4 significand bits */
#ifdef USE_INFINITY
if (r == 0x7ff0)
{
#ifdef NANS
#ifdef IBMPC
if (((pe[3] & 0xf) != 0) || (pe[2] != 0)
|| (pe[1] != 0) || (pe[0] != 0))
{
enan (y, NBITS);
return;
}
#else /* !IBMPC */
if (((pe[0] & 0xf) != 0) || (pe[1] != 0)
|| (pe[2] != 0) || (pe[3] != 0))
{
enan (y, NBITS);
return;
}
#endif /* !IBMPC */
#endif /* NANS */
eclear (y);
einfin (y, ldp);
if (yy[0])
eneg (y);
return;
}
#endif
r >>= 4;
/* If zero exponent, then the significand is denormalized.
* So, take back the understood high significand bit. */
if (r == 0)
{
denorm = 1;
yy[M] &= ~0x10;
}
r += EXONE - 01777;
yy[E] = r;
p = &yy[M + 1];
#ifdef IBMPC
*p++ = *(--e);
*p++ = *(--e);
*p++ = *(--e);
#else /* !IBMPC */
++e;
*p++ = *e++;
*p++ = *e++;
*p++ = *e++;
#endif /* !IBMPC */
(void) eshift (yy, -5);
if (denorm)
{ /* if zero exponent, then normalize the significand */
if ((k = enormlz (yy)) > NBITS)
ecleazs (yy);
else
yy[E] -= (unsigned short) (k - 1);
}
emovo (yy, y, ldp);
#endif /* !DEC */
}
/*
; e type to IEEE double precision
; double d;
; unsigned short x[NE];
; etoe53( x, &d );
*/
#ifdef DEC
static void
etoe53 (x, e)
unsigned short *x, *e;
{
etodec (x, e); /* see etodec.c */
}
static void
__attribute__ ((__unused__))
toe53 (x, y)
unsigned short *x, *y;
{
todec (x, y);
}
#else
static void
__attribute__ ((__unused__))
toe53 (short unsigned int *x, short unsigned int *y)
{
unsigned short i;
unsigned short *p;
#ifdef NANS
if (eiisnan (x))
{
enan (y, 53);
return;
}
#endif
p = &x[0];
#ifdef IBMPC
y += 3;
#endif
#ifdef DEC
y += 3;
#endif
*y = 0; /* output high order */
if (*p++)
*y = 0x8000; /* output sign bit */
i = *p++;
if (i >= (unsigned int) 2047)
{ /* Saturate at largest number less than infinity. */
#ifdef USE_INFINITY
*y |= 0x7ff0;
#ifdef IBMPC
*(--y) = 0;
*(--y) = 0;
*(--y) = 0;
#else /* !IBMPC */
++y;
*y++ = 0;
*y++ = 0;
*y++ = 0;
#endif /* IBMPC */
#else /* !USE_INFINITY */
*y |= (unsigned short) 0x7fef;
#ifdef IBMPC
*(--y) = 0xffff;
*(--y) = 0xffff;
*(--y) = 0xffff;
#else /* !IBMPC */
++y;
*y++ = 0xffff;
*y++ = 0xffff;
*y++ = 0xffff;
#endif
#endif /* !USE_INFINITY */
return;
}
if (i == 0)
{
(void) eshift (x, 4);
}
else
{
i <<= 4;
(void) eshift (x, 5);
}
i |= *p++ & (unsigned short) 0x0f; /* *p = xi[M] */
*y |= (unsigned short) i; /* high order output already has sign bit set */
#ifdef IBMPC
*(--y) = *p++;
*(--y) = *p++;
*(--y) = *p;
#else /* !IBMPC */
++y;
*y++ = *p++;
*y++ = *p++;
*y++ = *p++;
#endif /* !IBMPC */
}
#endif /* not DEC */
#endif /* LDBL_MANT_DIG == 53 */
#if LDBL_MANT_DIG == 24
/*
; Convert IEEE single precision to e type
; float d;
; unsigned short x[N+2];
; dtox( &d, x );
*/
void
e24toe (short unsigned int *pe, short unsigned int *y, LDPARMS * ldp)
{
register unsigned short r;
register unsigned short *p, *e;
unsigned short yy[NI];
int denorm, k;
e = pe;
denorm = 0; /* flag if denormalized number */
ecleaz (yy);
#ifdef IBMPC
e += 1;
#endif
#ifdef DEC
e += 1;
#endif
r = *e;
yy[0] = 0;
if (r & 0x8000)
yy[0] = 0xffff;
yy[M] = (r & 0x7f) | 0200;
r &= ~0x807f; /* strip sign and 7 significand bits */
#ifdef USE_INFINITY
if (r == 0x7f80)
{
#ifdef NANS
#ifdef MIEEE
if (((pe[0] & 0x7f) != 0) || (pe[1] != 0))
{
enan (y, NBITS);
return;
}
#else /* !MIEEE */
if (((pe[1] & 0x7f) != 0) || (pe[0] != 0))
{
enan (y, NBITS);
return;
}
#endif /* !MIEEE */
#endif /* NANS */
eclear (y);
einfin (y, ldp);
if (yy[0])
eneg (y);
return;
}
#endif
r >>= 7;
/* If zero exponent, then the significand is denormalized.
* So, take back the understood high significand bit. */
if (r == 0)
{
denorm = 1;
yy[M] &= ~0200;
}
r += EXONE - 0177;
yy[E] = r;
p = &yy[M + 1];
#ifdef IBMPC
*p++ = *(--e);
#endif
#ifdef DEC
*p++ = *(--e);
#endif
#ifdef MIEEE
++e;
*p++ = *e++;
#endif
(void) eshift (yy, -8);
if (denorm)
{ /* if zero exponent, then normalize the significand */
if ((k = enormlz (yy)) > NBITS)
ecleazs (yy);
else
yy[E] -= (unsigned short) (k - 1);
}
emovo (yy, y, ldp);
}
static void
__attribute__ ((__unused__))
toe24 (short unsigned int *x, short unsigned int *y)
{
unsigned short i;
unsigned short *p;
#ifdef NANS
if (eiisnan (x))
{
enan (y, 24);
return;
}
#endif
p = &x[0];
#ifdef IBMPC
y += 1;
#endif
#ifdef DEC
y += 1;
#endif
*y = 0; /* output high order */
if (*p++)
*y = 0x8000; /* output sign bit */
i = *p++;
if (i >= 255)
{ /* Saturate at largest number less than infinity. */
#ifdef USE_INFINITY
*y |= (unsigned short) 0x7f80;
#ifdef IBMPC
*(--y) = 0;
#endif
#ifdef DEC
*(--y) = 0;
#endif
#ifdef MIEEE
++y;
*y = 0;
#endif
#else /* !USE_INFINITY */
*y |= (unsigned short) 0x7f7f;
#ifdef IBMPC
*(--y) = 0xffff;
#endif
#ifdef DEC
*(--y) = 0xffff;
#endif
#ifdef MIEEE
++y;
*y = 0xffff;
#endif
#endif /* !USE_INFINITY */
return;
}
if (i == 0)
{
(void) eshift (x, 7);
}
else
{
i <<= 7;
(void) eshift (x, 8);
}
i |= *p++ & (unsigned short) 0x7f; /* *p = xi[M] */
*y |= i; /* high order output already has sign bit set */
#ifdef IBMPC
*(--y) = *p;
#endif
#ifdef DEC
*(--y) = *p;
#endif
#ifdef MIEEE
++y;
*y = *p;
#endif
}
#endif /* LDBL_MANT_DIG == 24 */
/* Compare two e type numbers.
*
* unsigned short a[NE], b[NE];
* ecmp( a, b );
*
* returns +1 if a > b
* 0 if a == b
* -1 if a < b
* -2 if either a or b is a NaN.
*/
static int
ecmp (const short unsigned int *a, const short unsigned int *b)
{
unsigned short ai[NI], bi[NI];
register unsigned short *p, *q;
register int i;
int msign;
#ifdef NANS
if (eisnan (a) || eisnan (b))
return (-2);
#endif
emovi (a, ai);
p = ai;
emovi (b, bi);
q = bi;
if (*p != *q)
{ /* the signs are different */
/* -0 equals + 0 */
for (i = 1; i < NI - 1; i++)
{
if (ai[i] != 0)
goto nzro;
if (bi[i] != 0)
goto nzro;
}
return (0);
nzro:
if (*p == 0)
return (1);
else
return (-1);
}
/* both are the same sign */
if (*p == 0)
msign = 1;
else
msign = -1;
i = NI - 1;
do
{
if (*p++ != *q++)
{
goto diff;
}
}
while (--i > 0);
return (0); /* equality */
diff:
if (*(--p) > *(--q))
return (msign); /* p is bigger */
else
return (-msign); /* p is littler */
}
/*
; Shift significand
;
; Shifts significand area up or down by the number of bits
; given by the variable sc.
*/
static int
eshift (short unsigned int *x, int sc)
{
unsigned short lost;
unsigned short *p;
if (sc == 0)
return (0);
lost = 0;
p = x + NI - 1;
if (sc < 0)
{
sc = -sc;
while (sc >= 16)
{
lost |= *p; /* remember lost bits */
eshdn6 (x);
sc -= 16;
}
while (sc >= 8)
{
lost |= *p & 0xff;
eshdn8 (x);
sc -= 8;
}
while (sc > 0)
{
lost |= *p & 1;
eshdn1 (x);
sc -= 1;
}
}
else
{
while (sc >= 16)
{
eshup6 (x);
sc -= 16;
}
while (sc >= 8)
{
eshup8 (x);
sc -= 8;
}
while (sc > 0)
{
eshup1 (x);
sc -= 1;
}
}
if (lost)
lost = 1;
return ((int) lost);
}
/*
; normalize
;
; Shift normalizes the significand area pointed to by argument
; shift count (up = positive) is returned.
*/
static int
enormlz (short unsigned int *x)
{
register unsigned short *p;
int sc;
sc = 0;
p = &x[M];
if (*p != 0)
goto normdn;
++p;
if (*p & 0x8000)
return (0); /* already normalized */
while (*p == 0)
{
eshup6 (x);
sc += 16;
/* With guard word, there are NBITS+16 bits available.
* return true if all are zero.
*/
if (sc > NBITS)
return (sc);
}
/* see if high byte is zero */
while ((*p & 0xff00) == 0)
{
eshup8 (x);
sc += 8;
}
/* now shift 1 bit at a time */
while ((*p & 0x8000) == 0)
{
eshup1 (x);
sc += 1;
if (sc > (NBITS + 16))
{
mtherr ("enormlz", UNDERFLOW);
return (sc);
}
}
return (sc);
/* Normalize by shifting down out of the high guard word
of the significand */
normdn:
if (*p & 0xff00)
{
eshdn8 (x);
sc -= 8;
}
while (*p != 0)
{
eshdn1 (x);
sc -= 1;
if (sc < -NBITS)
{
mtherr ("enormlz", OVERFLOW);
return (sc);
}
}
return (sc);
}
/* Convert e type number to decimal format ASCII string.
* The constants are for 64 bit precision.
*/
#define NTEN 12
#define MAXP 4096
#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 */
};
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 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 */
};
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
/* ASCII string outputs for unix */
#if 0
void
_IO_ldtostr (x, string, ndigs, flags, fmt)
long double *x;
char *string;
int ndigs;
int flags;
char fmt;
{
unsigned short w[NI];
char *t, *u;
LDPARMS rnd;
LDPARMS *ldp = &rnd;
rnd.rlast = -1;
rnd.rndprc = NBITS;
if (sizeof (long double) == 16)
e113toe ((unsigned short *) x, w, ldp);
else
e64toe ((unsigned short *) x, w, ldp);
etoasc (w, string, ndigs, -1, ldp);
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';
}
}
#endif
/* This routine will not return more than NDEC+1 digits. */
char *
_ldtoa_r (struct _reent *ptr, long double d, int mode, int ndigits,
int *decpt, int *sign, char **rve)
{
unsigned short e[NI];
char *s, *p;
int i, j, k;
int orig_ndigits;
LDPARMS rnd;
LDPARMS *ldp = &rnd;
char *outstr;
union uconv du;
du.d = d;
orig_ndigits = ndigits;
rnd.rlast = -1;
rnd.rndprc = NBITS;
_REENT_CHECK_MP (ptr);
/* reentrancy addition to use mprec storage pool */
if (_REENT_MP_RESULT (ptr))
{
_REENT_MP_RESULT (ptr)->_k = _REENT_MP_RESULT_K (ptr);
_REENT_MP_RESULT (ptr)->_maxwds = 1 << _REENT_MP_RESULT_K (ptr);
Bfree (ptr, _REENT_MP_RESULT (ptr));
_REENT_MP_RESULT (ptr) = 0;
}
#if LDBL_MANT_DIG == 24
e24toe (&du.pe, e, ldp);
#elif LDBL_MANT_DIG == 53
e53toe (&du.pe, e, ldp);
#elif LDBL_MANT_DIG == 64
e64toe (&du.pe, e, ldp);
#else
e113toe (&du.pe, e, ldp);
#endif
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;
char outbuf[ndigits + MAX_EXP_DIGITS + 10];
etoasc (e, outbuf, ndigits, mode, ldp);
s = outbuf;
if (eisinf (e) || eisnan (e))
{
*decpt = 9999;
goto stripspaces;
}
*decpt = ldp->outexpon + 1;
/* Transform the string returned by etoasc into what the caller wants. */
/* Look for decimal point and delete it from the string. */
s = outbuf;
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 > outbuf)
--s;
*s = '\0';
stripspaces:
/* Strip leading spaces and sign. */
p = outbuf;
while (*p == ' ' || *p == '-')
++p;
/* Find new end of string. */
s = outbuf;
while ((*s++ = *p++) != '\0')
;
--s;
/* Strip trailing zeros. */
if (mode == 2)
k = 1;
else if (ndigits > ldp->outexpon)
k = ndigits;
else
k = ldp->outexpon;
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 = 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) = eBalloc (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 = outstr + (s - outbuf);
return outstr;
}
/* Routine used to tell if long double is NaN or Infinity or regular number.
Returns: 0 = regular number
1 = Nan
2 = Infinity
*/
int
_ldcheck (long double *d)
{
unsigned short e[NI];
LDPARMS rnd;
LDPARMS *ldp = &rnd;
union uconv du;
rnd.rlast = -1;
rnd.rndprc = NBITS;
du.d = *d;
#if LDBL_MANT_DIG == 24
e24toe (&du.pe, e, ldp);
#elif LDBL_MANT_DIG == 53
e53toe (&du.pe, e, ldp);
#elif LDBL_MANT_DIG == 64
e64toe (&du.pe, e, ldp);
#else
e113toe (&du.pe, e, ldp);
#endif
if ((e[NE - 1] & 0x7fff) == 0x7fff)
{
#ifdef NANS
if (eisnan (e))
return (1);
#endif
return (2);
}
else
return (0);
} /* _ldcheck */
static void
etoasc (short unsigned int *x, char *string, int ndigits, int outformat,
LDPARMS * ldp)
{
long digit;
unsigned short y[NI], t[NI], u[NI], w[NI];
const unsigned short *p, *r, *ten;
unsigned short sign;
int i, j, k, expon, rndsav, ndigs;
char *s, *ss;
unsigned short m;
unsigned short *equot = ldp->equot;
ndigs = ndigits;
rndsav = ldp->rndprc;
#ifdef NANS
if (eisnan (x))
{
sprintf (string, " NaN ");
expon = 9999;
goto bxit;
}
#endif
ldp->rndprc = NBITS; /* set to full precision */
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, ldp);
efloor (t, w, ldp);
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, ldp);
emul (p, t, t, ldp);
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, ldp);
expon -= 1;
}
}
else
{
emovi (y, w);
/* Note that this loop does not access the incoming string array,
* which may be shorter than NDEC + 1 bytes! */
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, ldp);
}
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, ldp);
emul (r, t, t, ldp);
expon += k;
}
k /= 2;
if (k == 0)
break;
p += NE;
r += NE;
}
ediv (t, eone, t, ldp);
}
isone:
/* Find the first (leading) digit. */
emovi (t, w);
emovz (w, t);
emovi (y, w);
emovz (w, y);
eiremain (t, y, ldp);
digit = equot[NI - 1];
while ((digit == 0) && (ecmp (y, ezero) != 0))
{
eshup1 (y);
emovz (y, u);
eshup1 (u);
eshup1 (u);
eaddm (u, y);
eiremain (t, y, ldp);
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, ldp);
*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)
{
emovo (y, t, ldp);
if (ecmp (t, ezero) != 0)
goto roun; /* round to nearest */
if (ndigs < 0 || (*(s - 1 - (*(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:
#ifdef __GO32__
if (expon >= 0)
sprintf (ss, "e+%02d", expon);
else
sprintf (ss, "e-%02d", -expon);
#else
sprintf (ss, "E%d", expon);
#endif
bxit:
ldp->rndprc = rndsav;
ldp->outexpon = expon;
}
#if 0 /* Broken, unusable implementation of strtold */
/*
; ASCTOQ
; ASCTOQ.MAC LATEST REV: 11 JAN 84
; SLM, 3 JAN 78
;
; Convert ASCII string to quadruple precision floating point
;
; Numeric input is free field decimal number
; with max of 15 digits with or without
; decimal point entered as ASCII from teletype.
; Entering E after the number followed by a second
; number causes the second number to be interpreted
; as a power of 10 to be multiplied by the first number
; (i.e., "scientific" notation).
;
; Usage:
; asctoq( string, q );
*/
long double
_strtold (char *s, char **se)
{
union uconv x;
LDPARMS rnd;
LDPARMS *ldp = &rnd;
int lenldstr;
rnd.rlast = -1;
rnd.rndprc = NBITS;
lenldstr = asctoeg (s, &x.pe, LDBL_MANT_DIG, ldp);
if (se)
*se = s + lenldstr;
return x.d;
}
#define REASONABLE_LEN 200
static int
asctoeg (char *ss, short unsigned int *y, int oprec, LDPARMS * ldp)
{
unsigned short yy[NI], xt[NI], tt[NI];
int esign, decflg, sgnflg, nexp, exp, prec, lost;
int k, trail, c, rndsav;
long lexp;
unsigned short nsign;
const unsigned short *p;
char *sp, *s, *lstr;
int lenldstr;
int mflag = 0;
char tmpstr[REASONABLE_LEN];
/* Copy the input string. */
c = strlen (ss) + 2;
if (c <= REASONABLE_LEN)
lstr = tmpstr;
else
{
lstr = (char *) calloc (c, 1);
mflag = 1;
}
s = ss;
lenldstr = 0;
while (*s == ' ') /* skip leading spaces */
{
++s;
++lenldstr;
}
sp = lstr;
for (k = 0; k < c; k++)
{
if ((*sp++ = *s++) == '\0')
break;
}
*sp = '\0';
s = lstr;
rndsav = ldp->rndprc;
ldp->rndprc = NBITS; /* Set to full precision */
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;
/* Check for syntax error */
c = *sp & 0x7f;
if ((c != 'e') && (c != 'E') && (c != '\0')
&& (c != '\n') && (c != '\r') && (c != ' ') && (c != ','))
goto error;
--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;
}
switch (*s)
{
case 'z':
break;
case 'E':
case 'e':
goto expnt;
case '.': /* decimal point */
if (decflg)
goto error;
++decflg;
break;
case '-':
nsign = 0xffff;
if (sgnflg)
goto error;
++sgnflg;
break;
case '+':
if (sgnflg)
goto error;
++sgnflg;
break;
case ',':
case ' ':
case '\0':
case '\n':
case '\r':
goto daldone;
case 'i':
case 'I':
goto infinite;
default:
error:
#ifdef NANS
enan (yy, NI * 16);
#else
mtherr ("asctoe", DOMAIN);
ecleaz (yy);
#endif
goto aexit;
}
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 *= 10;
exp += *s++ - '0';
if (exp > 4977)
{
if (esign < 0)
goto zero;
else
goto infinite;
}
}
if (esign < 0)
exp = -exp;
if (exp > 4932)
{
infinite:
ecleaz (yy);
yy[E] = 0x7fff; /* infinity */
goto aexit;
}
if (exp < -4977)
{
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, ldp);
/* 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, ldp);
lexp += EXONE;
nexp -= 4096;
}
}
p = &etens[NTEN][0];
emov (eone, xt);
exp = 1;
do
{
if (exp & nexp)
emul (p, xt, xt, ldp);
p -= NE;
exp = exp + exp;
}
while (exp <= MAXP);
emovi (xt, tt);
if (esign < 0)
{
lexp -= tt[E];
k = edivm (tt, yy, ldp);
lexp += EXONE;
}
else
{
lexp += tt[E];
k = emulm (tt, yy, ldp);
lexp -= EXONE - 1;
}
expdon:
/* Round and convert directly to the destination type */
if (oprec == 53)
lexp -= EXONE - 0x3ff;
else if (oprec == 24)
lexp -= EXONE - 0177;
#ifdef DEC
else if (oprec == 56)
lexp -= EXONE - 0201;
#endif
ldp->rndprc = oprec;
emdnorm (yy, k, 0, lexp, 64, ldp);
aexit:
ldp->rndprc = rndsav;
yy[0] = nsign;
switch (oprec)
{
#ifdef DEC
case 56:
todec (yy, y); /* see etodec.c */
break;
#endif
#if LDBL_MANT_DIG == 53
case 53:
toe53 (yy, y);
break;
#elif LDBL_MANT_DIG == 24
case 24:
toe24 (yy, y);
break;
#elif LDBL_MANT_DIG == 64
case 64:
toe64 (yy, y);
break;
#elif LDBL_MANT_DIG == 113
case 113:
toe113 (yy, y);
break;
#else
case NBITS:
emovo (yy, y, ldp);
break;
#endif
}
lenldstr += s - lstr;
if (mflag)
free (lstr);
return lenldstr;
}
#endif
/* y = largest integer not greater than x
* (truncated toward minus infinity)
*
* unsigned short x[NE], y[NE]
* LDPARMS *ldp
*
* efloor( x, y, ldp );
*/
static const unsigned short bmask[] = {
0xffff,
0xfffe,
0xfffc,
0xfff8,
0xfff0,
0xffe0,
0xffc0,
0xff80,
0xff00,
0xfe00,
0xfc00,
0xf800,
0xf000,
0xe000,
0xc000,
0x8000,
0x0000,
};
static void
efloor (short unsigned int *x, short unsigned int *y, LDPARMS * ldp)
{
register unsigned short *p;
int e, expon, i;
unsigned short f[NE];
emov (x, f); /* leave in external format */
expon = (int) f[NE - 1];
e = (expon & 0x7fff) - (EXONE - 1);
if (e <= 0)
{
eclear (y);
goto isitneg;
}
/* number of bits to clear out */
e = NBITS - e;
emov (f, y);
if (e <= 0)
return;
p = &y[0];
while (e >= 16)
{
*p++ = 0;
e -= 16;
}
/* clear the remaining bits */
*p &= bmask[e];
/* truncate negatives toward minus infinity */
isitneg:
if ((unsigned short) expon & (unsigned short) 0x8000)
{
for (i = 0; i < NE - 1; i++)
{
if (f[i] != y[i])
{
esub (eone, y, y, ldp);
break;
}
}
}
}
static void
eiremain (short unsigned int *den, short unsigned int *num, LDPARMS * ldp)
{
long ld, ln;
unsigned short j;
unsigned short *equot = ldp->equot;
ld = den[E];
ld -= enormlz (den);
ln = num[E];
ln -= enormlz (num);
ecleaz (equot);
while (ln >= ld)
{
if (ecmpm (den, num) <= 0)
{
esubm (den, num);
j = 1;
}
else
{
j = 0;
}
eshup1 (equot);
equot[NI - 1] |= j;
eshup1 (num);
ln -= 1;
}
emdnorm (num, 0, 0, ln, 0, ldp);
}
/* NaN bit patterns
*/
#ifdef MIEEE
#if !defined(__mips)
static const unsigned short nan113[8] = {
0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff
};
static const unsigned short nan64[6] = {
0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff
};
static const unsigned short nan53[4] = { 0x7fff, 0xffff, 0xffff, 0xffff };
static const unsigned short nan24[2] = { 0x7fff, 0xffff };
#elif defined(__mips_nan2008) /* __mips */
static const unsigned short nan113[8] = { 0x7fff, 0x8000, 0, 0, 0, 0, 0, 0 };
static const unsigned short nan64[6] = { 0x7fff, 0xc000, 0, 0, 0, 0 };
static const unsigned short nan53[4] = { 0x7ff8, 0, 0, 0 };
static const unsigned short nan24[2] = { 0x7fc0, 0 };
#else /* __mips && !__mips_nan2008 */
static const unsigned short nan113[8] = {
0x7fff, 0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff
};
static const unsigned short nan64[6] = {
0x7fff, 0xbfff, 0xffff, 0xffff, 0xffff, 0xffff
};
static const unsigned short nan53[4] = { 0x7ff7, 0xffff, 0xffff, 0xffff };
static const unsigned short nan24[2] = { 0x7fbf, 0xffff };
#endif /* __mips && !__mips_nan2008 */
#else /* !MIEEE */
#if !defined(__mips) || defined(__mips_nan2008)
static const unsigned short nan113[8] = { 0, 0, 0, 0, 0, 0, 0x8000, 0x7fff };
static const unsigned short nan64[6] = { 0, 0, 0, 0, 0xc000, 0x7fff };
static const unsigned short nan53[4] = { 0, 0, 0, 0x7ff8 };
static const unsigned short nan24[2] = { 0, 0x7fc0 };
#else /* __mips && !__mips_nan2008 */
static const unsigned short nan113[8] = {
0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0x7fff, 0x7fff
};
static const unsigned short nan64[6] = {
0xffff, 0xffff, 0xffff, 0xffff, 0xbfff, 0x7fff
};
static const unsigned short nan53[4] = { 0xffff, 0xffff, 0xffff, 0x7ff7 };
static const unsigned short nan24[2] = { 0xffff, 0x7fbf };
#endif /* __mips && !__mips_nan2008 */
#endif /* !MIEEE */
static void
enan (short unsigned int *nan, int size)
{
int i, n;
const unsigned short *p;
switch (size)
{
#ifndef DEC
case 113:
n = 8;
p = nan113;
break;
case 64:
n = 6;
p = nan64;
break;
case 53:
n = 4;
p = nan53;
break;
case 24:
n = 2;
p = nan24;
break;
case NBITS:
#if !defined(__mips) || defined(__mips_nan2008)
for (i = 0; i < NE - 2; i++)
*nan++ = 0;
*nan++ = 0xc000;
#else /* __mips && !__mips_nan2008 */
for (i = 0; i < NE - 2; i++)
*nan++ = 0xffff;
*nan++ = 0xbfff;
#endif /* __mips && !__mips_nan2008 */
*nan++ = 0x7fff;
return;
case NI * 16:
*nan++ = 0;
*nan++ = 0x7fff;
*nan++ = 0;
#if !defined(__mips) || defined(__mips_nan2008)
*nan++ = 0xc000;
for (i = 4; i < NI - 1; i++)
*nan++ = 0;
#else /* __mips && !__mips_nan2008 */
*nan++ = 0xbfff;
for (i = 4; i < NI - 1; i++)
*nan++ = 0xffff;
#endif /* __mips && !__mips_nan2008 */
*nan++ = 0;
return;
#endif
default:
mtherr ("enan", DOMAIN);
return;
}
for (i = 0; i < n; i++)
*nan++ = *p++;
}
Reference in New Issue View Git Blame Copy Permalink