diff options
Diffstat (limited to 'winsup/cygwin/math/cephes_emath.h')
| -rw-r--r-- | winsup/cygwin/math/cephes_emath.h | 719 |
1 files changed, 719 insertions, 0 deletions
diff --git a/winsup/cygwin/math/cephes_emath.h b/winsup/cygwin/math/cephes_emath.h new file mode 100644 index 0000000..b92d710 --- /dev/null +++ b/winsup/cygwin/math/cephes_emath.h @@ -0,0 +1,719 @@ +/** + * This file has no copyright assigned and is placed in the Public Domain. + * This file is part of the mingw-w64 runtime package. + * No warranty is given; refer to the file DISCLAIMER.PD within this package. + */ +#ifndef _CEPHES_EMATH_H +#define _CEPHES_EMATH_H + +/** + * This is a workaround for a gcc bug + */ +#define __restrict__ + +/* 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> + * + */ + + +/* 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 + * + * __asctoe64( string, &d ) ASCII string to long double + * __asctoeg( string, e, prec ) ASCII string to specified precision + * __e64toe( &d, e ) IEEE long double precision to e type + * __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 ) c = b / a + * __efloor( a, b ) 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 ) 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 ) c = b * a + * __eneg(e) e = -e + * __eround( a, b ) b = nearest integer value to a + * __esub( a, b, c ) 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 ) 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 + * __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 ) divide significands, bi = bi / ai + * __emdnorm(ai,l,s,exp) normalize and round off + * __emovi( a, ai ) convert external a to internal ai + * __emovo( ai, a ) convert internal ai to external a + * __emovz( ai, bi ) bi = ai, low guard word of bi = 0 + * __emulm( ai, bi ) 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 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). + */ + +#include <stdio.h> +#include <stdlib.h> +#include <string.h> +#include <errno.h> +#include <math.h> +#include <locale.h> +#include <ctype.h> + +#undef alloca +#define alloca __builtin_alloca + +/* Don't build non-ANSI _IO_ldtoa. It is not thread safe. */ +#ifndef USE_LDTOA +#define USE_LDTOA 0 +#endif + + + /* Number of 16 bit words in external x type format */ +#define NE 6 + + /* 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 + * = NBITS*log10(2) + */ +#define NDEC (NBITS*8/27) + + /* The exponent of 1.0 */ +#define EXONE (0x3fff) + + +#define mtherr(x,y) + + +extern long double strtold (const char * __restrict__ s, char ** __restrict__ se); +extern int __asctoe64(const char * __restrict__ ss, + short unsigned int * __restrict__ y); +extern void __emul(const short unsigned int * a, + const short unsigned int * b, + short unsigned int * c); +extern int __ecmp(const short unsigned int * __restrict__ a, + const short unsigned int * __restrict__ b); +extern int __enormlz(short unsigned int *x); +extern int __eshift(short unsigned int *x, int sc); +extern void __eaddm(const short unsigned int * __restrict__ x, + short unsigned int * __restrict__ y); +extern void __esubm(const short unsigned int * __restrict__ x, + short unsigned int * __restrict__ y); +extern void __emdnorm(short unsigned int *s, int lost, int subflg, + int exp, int rcntrl, const int rndprc); +extern void __toe64(short unsigned int * __restrict__ a, + short unsigned int * __restrict__ b); +extern int __edivm(short unsigned int * __restrict__ den, + short unsigned int * __restrict__ num); +extern int __emulm(const short unsigned int * __restrict__ a, + short unsigned int * __restrict__ b); +extern void __emovi(const short unsigned int * __restrict__ a, + short unsigned int * __restrict__ b); +extern void __emovo(const short unsigned int * __restrict__ a, + short unsigned int * __restrict__ b); + +#if USE_LDTOA + +extern char * _IO_ldtoa(long double, int, int, int *, int *, char **); +extern void _IO_ldtostr(long double *x, char *string, int ndigs, + int flags, char fmt); + +extern void __eiremain(short unsigned int * __restrict__ den, + short unsigned int *__restrict__ num, + short unsigned int *__restrict__ equot); +extern void __efloor(short unsigned int *x, short unsigned int *y); +extern void __eadd1(const short unsigned int * __restrict__ a, + const short unsigned int * __restrict__ b, + short unsigned int * __restrict__ c, + int subflg); +extern void __esub(const short unsigned int *a, const short unsigned int *b, + short unsigned int *c); +extern void __ediv(const short unsigned int *a, const short unsigned int *b, + short unsigned int *c); +extern void __e64toe(short unsigned int *pe, short unsigned int *y); + + +#endif + +static __inline__ int __eisneg(const short unsigned int *x); +static __inline__ int __eisinf(const short unsigned int *x); +static __inline__ int __eisnan(const short unsigned int *x); +static __inline__ int __eiszero(const short unsigned int *a); +static __inline__ void __emovz(register const short unsigned int * __restrict__ a, + register short unsigned int * __restrict__ b); +static __inline__ void __eclear(register short unsigned int *x); +static __inline__ void __ecleaz(register short unsigned int *xi); +static __inline__ void __ecleazs(register short unsigned int *xi); +static __inline__ int __eiisinf(const short unsigned int *x); +static __inline__ int __eiisnan(const short unsigned int *x); +static __inline__ int __eiiszero(const short unsigned int *x); +static __inline__ void __enan_64(short unsigned int *nanptr); +static __inline__ void __enan_NBITS (short unsigned int *nanptr); +static __inline__ void __enan_NI16 (short unsigned int *nanptr); +static __inline__ void __einfin(register short unsigned int *x); +static __inline__ void __eneg(short unsigned int *x); +static __inline__ void __eshup1(register short unsigned int *x); +static __inline__ void __eshup8(register short unsigned int *x); +static __inline__ void __eshup6(register short unsigned int *x); +static __inline__ void __eshdn1(register short unsigned int *x); +static __inline__ void __eshdn8(register short unsigned int *x); +static __inline__ void __eshdn6(register short unsigned int *x); + + + +/* Intel IEEE, low order words come first: + */ +#define IBMPC 1 + +/* Define 1 for ANSI C atan2() function + * See atan.c and clog.c. + */ +#define ANSIC 1 + +/*define VOLATILE volatile*/ +#define VOLATILE + +/* For 12-byte long doubles on an i386, pad a 16-bit short 0 + * to the end of real constants initialized by integer arrays. + * + * #define XPD 0, + * + * Otherwise, the type is 10 bytes long and XPD should be + * defined blank. + * + * #define XPD + */ +#define XPD 0, +/* #define XPD */ +#define NANS + +/* NaN's require infinity support. */ +#ifdef NANS +#ifndef INFINITY +#define INFINITY +#endif +#endif + +/* This handles 64-bit long ints. */ +#define LONGBITS (8 * sizeof(long)) + + +#define NTEN 12 +#define MAXP 4096 + +/* +; Clear out entire external format number. +; +; unsigned short x[]; +; eclear( x ); +*/ + +static __inline__ void __eclear(register short unsigned int *x) +{ + memset(x, 0, NE * sizeof(unsigned short)); +} + + +/* Move external format number from a to b. + * + * emov( a, b ); + */ + +static __inline__ void __emov(register const short unsigned int * __restrict__ a, + register short unsigned int * __restrict__ b) +{ + memcpy(b, a, NE * sizeof(unsigned short)); +} + + +/* +; Negate external format number +; +; unsigned short x[NE]; +; eneg( x ); +*/ + +static __inline__ 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 __inline__ 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 __inline__ 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 __inline__ int __eisnan(const short unsigned int *x) +{ +#ifdef NANS + int i; + /* NaN has maximum __exponent */ + if ((x[NE - 1] & 0x7fff) == 0x7fff) + /* ... 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 __inline__ void __einfin(register short unsigned int *x) +{ + register int i; +#ifdef INFINITY + for (i = 0; i < NE - 1; i++) + *x++ = 0; + *x |= 32767; +#else + for (i = 0; i < NE - 1; i++) + *x++ = 0xffff; + *x |= 32766; + *(x - 5) = 0; +#endif +} + +/* Clear out internal format number. + */ + +static __inline__ void __ecleaz(register short unsigned int *xi) +{ + memset(xi, 0, NI * sizeof(unsigned short)); +} + +/* same, but don't touch the sign. */ + +static __inline__ void __ecleazs(register short unsigned int *xi) +{ + ++xi; + memset(xi, 0, (NI-1) * sizeof(unsigned short)); +} + +/* Move internal format number from a to b. + */ +static __inline__ void __emovz(register const short unsigned int * __restrict__ a, + register short unsigned int * __restrict__ b) +{ + memcpy(b, a, (NI-1) * sizeof(unsigned short)); + b[NI - 1] = 0; +} + +/* Return nonzero if internal format number is a NaN. + */ + +static __inline__ int __eiisnan (const 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); +} + +/* Return nonzero if external format number is zero. */ + +static __inline__ int +__eiszero(const short unsigned int * a) +{ + union { + long double ld; + unsigned short sh[8]; + } av; + av.ld = 0.0; + memcpy (av.sh, a, 12); + if (av.ld == 0.0) + return (1); + return (0); +} + +/* Return nonzero if internal format number is zero. */ + +static __inline__ int +__eiiszero(const short unsigned int * ai) +{ + int i; + /* skip the sign word */ + for (i = 1; i < NI - 1; i++ ) + { + if (ai[i] != 0) + return (0); + } + return (1); +} + + +/* Return nonzero if internal format number is infinite. */ + +static __inline__ int +__eiisinf (const unsigned short *x) +{ +#ifdef NANS + if (__eiisnan (x)) + return (0); +#endif + if ((x[E] & 0x7fff) == 0x7fff) + return (1); + return (0); +} + +/* +; 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 __inline__ int __ecmpm(register const short unsigned int * __restrict__ a, + register const short unsigned int * __restrict__ 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 __inline__ 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 __inline__ 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 __inline__ 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 __inline__ 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 __inline__ 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 __inline__ 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 __inline__ void __enan_64(unsigned short* nanptr) +{ + int i; + for (i = 0; i < 3; i++) + *nanptr++ = 0; + *nanptr++ = 0xc000; + *nanptr++ = 0x7fff; + *nanptr = 0; + return; +} + +static __inline__ void __enan_NBITS(unsigned short* nanptr) +{ + int i; + for (i = 0; i < NE - 2; i++) + *nanptr++ = 0; + *nanptr++ = 0xc000; + *nanptr = 0x7fff; + return; +} + +static __inline__ void __enan_NI16(unsigned short* nanptr) +{ + int i; + *nanptr++ = 0; + *nanptr++ = 0x7fff; + *nanptr++ = 0; + *nanptr++ = 0xc000; + for (i = 4; i < NI; i++) + *nanptr++ = 0; + return; +} + +#endif /* _CEPHES_EMATH_H */ + |