diff options
Diffstat (limited to 'winsup/mingw/mingwex/math/cephes_emath.h')
| -rw-r--r-- | winsup/mingw/mingwex/math/cephes_emath.h | 713 |
1 files changed, 0 insertions, 713 deletions
diff --git a/winsup/mingw/mingwex/math/cephes_emath.h b/winsup/mingw/mingwex/math/cephes_emath.h deleted file mode 100644 index 121937a..0000000 --- a/winsup/mingw/mingwex/math/cephes_emath.h +++ /dev/null @@ -1,713 +0,0 @@ -#ifndef _CEPHES_EMATH_H -#define _CEPHES_EMATH_H - -/* 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> - -#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, - long 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 *nan); -static __inline__ void __enan_NBITS (short unsigned int *nan); -static __inline__ void __enan_NI16 (short unsigned int *nan); -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) -{ -if (*((long double*) a) == 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* nan) -{ - - int i; - for( i=0; i<3; i++ ) - *nan++ = 0; - *nan++ = 0xc000; - *nan++ = 0x7fff; - *nan = 0; - return; -} - -static __inline__ void __enan_NBITS(unsigned short* nan) -{ - int i; - for( i=0; i<NE-2; i++ ) - *nan++ = 0; - *nan++ = 0xc000; - *nan = 0x7fff; - return; -} - -static __inline__ void __enan_NI16(unsigned short* nan) -{ - int i; - *nan++ = 0; - *nan++ = 0x7fff; - *nan++ = 0; - *nan++ = 0xc000; - for( i=4; i<NI; i++ ) - *nan++ = 0; - return; -} - -#endif /* _CEPHES_EMATH_H */ - |