From 14faef01e4db4ab03a12ed08d6cee91f1f5a0af5 Mon Sep 17 00:00:00 2001 From: Danny Smith Date: Tue, 10 Dec 2002 10:48:18 +0000 Subject: Merge from mingw trunk (changes since 2002-11-19). --- winsup/mingw/ChangeLog | 50 +- winsup/mingw/include/math.h | 63 +- winsup/mingw/include/stdlib.h | 2 + winsup/mingw/include/wchar.h | 1 + winsup/mingw/mingwex/Makefile.in | 42 +- winsup/mingw/mingwex/ldtoa.c | 614 ++++++++++++++ winsup/mingw/mingwex/math/cephes_emath.c | 1318 ++++++++++++++++++++++++++++++ winsup/mingw/mingwex/math/cephes_emath.h | 710 ++++++++++++++++ winsup/mingw/mingwex/math/cephes_mconf.h | 98 +++ winsup/mingw/mingwex/math/lgamma.c | 359 ++++++++ winsup/mingw/mingwex/math/lgammaf.c | 253 ++++++ winsup/mingw/mingwex/math/lgammal.c | 416 ++++++++++ winsup/mingw/mingwex/math/powif2.c | 164 ++++ winsup/mingw/mingwex/math/s_erf.c | 342 ++++++++ winsup/mingw/mingwex/math/sf_erf.c | 259 ++++++ winsup/mingw/mingwex/math/tgamma.c | 385 +++++++++ winsup/mingw/mingwex/math/tgammaf.c | 265 ++++++ winsup/mingw/mingwex/math/tgammal.c | 501 ++++++++++++ winsup/mingw/mingwex/strtold.c | 384 +++++++++ winsup/mingw/mingwex/wcstold.c | 76 ++ 20 files changed, 6264 insertions(+), 38 deletions(-) create mode 100644 winsup/mingw/mingwex/ldtoa.c create mode 100644 winsup/mingw/mingwex/math/cephes_emath.c create mode 100644 winsup/mingw/mingwex/math/cephes_emath.h create mode 100644 winsup/mingw/mingwex/math/lgamma.c create mode 100644 winsup/mingw/mingwex/math/lgammaf.c create mode 100644 winsup/mingw/mingwex/math/lgammal.c create mode 100644 winsup/mingw/mingwex/math/powif2.c create mode 100644 winsup/mingw/mingwex/math/s_erf.c create mode 100644 winsup/mingw/mingwex/math/sf_erf.c create mode 100644 winsup/mingw/mingwex/math/tgamma.c create mode 100644 winsup/mingw/mingwex/math/tgammaf.c create mode 100644 winsup/mingw/mingwex/math/tgammal.c create mode 100644 winsup/mingw/mingwex/strtold.c create mode 100644 winsup/mingw/mingwex/wcstold.c diff --git a/winsup/mingw/ChangeLog b/winsup/mingw/ChangeLog index 0482b380e23..41aa680ac3f 100644 --- a/winsup/mingw/ChangeLog +++ b/winsup/mingw/ChangeLog @@ -1,6 +1,54 @@ +2002-12-10 Danny Smith + + Merge from mingw trunk (changes since 2002-11-19). + + 2002-12-08 Danny Smith + + * mingwex/math/s_erf.c: New file. + * mingwex/math/sf_erf.c: New file. + * mingwex/Makefile.in (MATH_DISTFILES): Add new files. + (MATH_OBJS): Add new objects. + * include/math.h (erf[f]): Add prototypes. + (erfc[f]): Add prototypes. + + 2002-12-07 Danny Smith + + * include/math.h: Add traditional/XOPEN math constants. + + 2002-11-27 Danny Smith + + * mingwex/math/lgamma.c: New file. + * mingwex/math/lgammaf.c: New file. + * mingwex/math/lgammal.c: New file. + * mingwex/math/tgamma.c: New file. + * mingwex/math/tgammaf.c: New file. + * mingwex/math/tgammal.c: New file. + * mingwex/math/cephes_mconf (polevlf): Add float version. + (p1evlf): Likewise. + Define _CEPHES_USE_ERRNO. + * mingwex/Makefile.in (MATH_DISTFILES): Add new files. + (MATH_OBJS): Add new objects. + * include/math.h (lgamma[fl]): Add prototypes. + (tgamma[fl]): Add prototypes. + + 2002-11-26 Danny Smith + + * mingwex/strtold.c: New file. + * mingwex/wcstold.c: New file. + * mingwex/ldtoa.c: New file. + * mingwex/math/cephes_emath.h: New file. + * mingwex/math/cephes_emath.c: New file. + * mingwex/Makefile.in (DISTFILES): Add new files. + (MATH_DISTFILES): Ditto. + (STDLIB_OBJS): New. Define as strtold.c wcstold.c. + (MATH_OBJS): Add cephes_emath.o. + (LIB_OBJS): Add $(STDLIB_OBJS). + * include/stdlib.h (strtold, wcstold): Add prototypes. + * include/wchar.h (wcstold): Add prototype. + 2002-11-19 Danny Smith - Merge from mingw trunk, . + Merge from mingw trunk (changes since 2002-10-04). 2002-11-09 Danny Smith diff --git a/winsup/mingw/include/math.h b/winsup/mingw/include/math.h index bf34ce89227..bba364da1b8 100644 --- a/winsup/mingw/include/math.h +++ b/winsup/mingw/include/math.h @@ -58,6 +58,22 @@ #endif /* Not _NO_OLDNAMES */ #endif /* Not __STRICT_ANSI__ */ +/* Traditional/XOPEN math constants (double precison) */ +#ifndef __STRICT_ANSI__ +#define M_E 2.7182818284590452354 +#define M_LOG2E 1.4426950408889634074 +#define M_LOG10E 0.43429448190325182765 +#define M_LN2 0.69314718055994530942 +#define M_LN10 2.30258509299404568402 +#define M_PI 3.14159265358979323846 +#define M_PI_2 1.57079632679489661923 +#define M_PI_4 0.78539816339744830962 +#define M_1_PI 0.31830988618379067154 +#define M_2_PI 0.63661977236758134308 +#define M_2_SQRTPI 1.12837916709551257390 +#define M_SQRT2 1.41421356237309504880 +#define M_SQRT1_2 0.70710678118654752440 +#endif /* These are also defined in Mingw float.h; needed here as well to work around GCC build issues. */ @@ -143,12 +159,6 @@ double log (double); double log10 (double); double pow (double, double); double sqrt (double); -extern __inline__ double sqrt (double __x) -{ - double res; - __asm__ ("fsqrt;" : "=t" (res) : "0" (__x)); - return res; -} double ceil (double); double floor (double); double fabs (double); @@ -504,21 +514,32 @@ extern __inline__ float powf (float __x, float __y) extern long double powl (long double, long double); /* 7.12.7.5 The sqrt functions. Double in C89. */ -extern __inline__ float sqrtf (float __x) -{ - float res; - __asm__ ("fsqrt" : "=t" (res) : "0" (__x)); - return res; -} - -extern __inline__ long double sqrtl (long double __x) -{ - long double res; - __asm__ ("fsqrt" : "=t" (res) : "0" (__x)); - return res; -} - -/* 7.12.8 Error and gamma functions: TODO */ +extern float sqrtf (float); +extern long double sqrtl (long double); + +/* 7.12.8.1 The erf functions */ +extern double erf (double); +extern float erff (float); +/* TODO +extern long double erfl (long double); +*/ + +/* 7.12.8.2 The erfc functions */ +extern double erfc (double); +extern float erfcf (float); +/* TODO +extern long double erfcl (long double); +*/ + +/* 7.12.8.3 The lgamma functions */ +extern double lgamma (double); +extern float lgammaf (float); +extern long double lgammal (long double); + +/* 77.12.8.4 The tgamma functions */ +extern double tgamma (double); +extern float tgammaf (float); +extern long double tgammal (long double); /* 7.12.9.1 Double in C89 */ extern float ceilf (float); diff --git a/winsup/mingw/include/stdlib.h b/winsup/mingw/include/stdlib.h index 6f0ee4e49de..577ee11a1bb 100644 --- a/winsup/mingw/include/stdlib.h +++ b/winsup/mingw/include/stdlib.h @@ -323,6 +323,7 @@ double strtod (const char*, char**); #if !defined __NO_ISOCEXT /* extern stubs in static libmingwex.a */ extern __inline__ float strtof (const char *__nptr, char **__endptr) { return (strtod (__nptr, __endptr));} +long double strtold (const char * __restrict__, char ** __restrict__); #endif /* __NO_ISOCEXT */ long strtol (const char*, char**, int); @@ -343,6 +344,7 @@ double wcstod (const wchar_t*, wchar_t**); #if !defined __NO_ISOCEXT /* extern stub in static libmingwex.a */ extern __inline__ float wcstof( const wchar_t *__nptr, wchar_t **__endptr) { return (wcstod(__nptr, __endptr)); } +long double wcstold (const wchar_t * __restrict__, wchar_t ** __restrict__); #endif /* __NO_ISOCEXT */ long wcstol (const wchar_t*, wchar_t**, int); diff --git a/winsup/mingw/include/wchar.h b/winsup/mingw/include/wchar.h index bf493bc1976..9c62139537f 100644 --- a/winsup/mingw/include/wchar.h +++ b/winsup/mingw/include/wchar.h @@ -407,6 +407,7 @@ double wcstod (const wchar_t*, wchar_t**); #if !defined __NO_ISOCEXT /* extern stub in static libmingwex.a */ extern __inline__ float wcstof( const wchar_t *__nptr, wchar_t **__endptr) { return (wcstod(__nptr, __endptr)); } +long double wcstold (const wchar_t * __restrict__, wchar_t ** __restrict__); #endif /* __NO_ISOCEXT */ __END_CSTD_NAMESPACE #define _WSTDLIB_DEFINED diff --git a/winsup/mingw/mingwex/Makefile.in b/winsup/mingw/mingwex/Makefile.in index 8ede9f9e97f..bb96cac7173 100644 --- a/winsup/mingw/mingwex/Makefile.in +++ b/winsup/mingw/mingwex/Makefile.in @@ -29,15 +29,16 @@ DISTFILES = Makefile.in configure configure.in \ _Exit.c atoll.c dirent.c feclearexcept.c fegetenv.c \ fegetexceptflag.c fegetround.c feholdexcept.c feraiseexcept.c \ fesetenv.c fesetexceptflag.c fesetround.c fetestexcept.c \ - feupdateenv.c fwide.c imaxabs.c imaxdiv.c lltoa.c lltow.c \ + feupdateenv.c fwide.c imaxabs.c imaxdiv.c ldtoa.c lltoa.c lltow.c \ mbsinit.c mingw-fseek.c sitest.c snprintf.c snwprintf.c \ - strtof.c strtoimax.c strtoumax.c testwmem.c ulltoa.c ulltow.c \ - vsnprintf.c vsnwprintf.c wcstof.c wcstoimax.c wcstoumax.c \ - wdirent.c wmemchr.c wmemcmp.c wmemcpy.c wmemmove.c wmemset.c \ - wtoll.c + strtof.c strtoimax.c strtold.c strtoumax.c testwmem.c \ + ulltoa.c ulltow.c vsnprintf.c vsnwprintf.c wcstof.c \ + wcstoimax.c wcstold.c wcstoumax.c wdirent.c wmemchr.c \ + wmemcmp.c wmemcpy.c wmemmove.c wmemset.c wtoll.c MATH_DISTFILES = \ acosf.c acosl.c asinf.c asinl.c atan2f.c atan2l.c \ - atanf.c atanl.c cbrt.c cbrtf.c cbrtl.c ceilf.S ceill.S cephes_mconf.h \ + atanf.c atanl.c cbrt.c cbrtf.c cbrtl.c ceilf.S ceill.S \ + cephes_emath.h cephes_emath.c cephes_mconf.h \ copysign.S copysignf.S copysignl.S cosf.S coshf.c coshl.c cosl.S \ exp2.S exp2f.S exp2l.S expf.c expl.c fabs.c fabsf.c fabsl.c \ fdim.c fdimf.c fdiml.c floorf.S floorl.S fma.S fmaf.S fmal.c \ @@ -45,7 +46,8 @@ MATH_DISTFILES = \ fmodl.c fp_consts.c fp_consts.h fp_constsf.c fp_constsl.c \ fpclassify.c fpclassifyf.c fpclassifyl.c \ frexpf.c frexpl.S fucom.c hypotf.c hypotl.c ilogb.S ilogbf.S \ - ilogbl.S isnan.c isnanf.c isnanl.c ldexpf.c ldexpl.c llrint.c \ + ilogbl.S isnan.c isnanf.c isnanl.c ldexpf.c ldexpl.c \ + lgamma.c lgammaf.c lgammal.c llrint.c \ llrintf.c llrintl.c llround.c llroundf.c llroundl.c \ log10f.S log10l.S log1p.S log1pf.S log1pl.S log2.S log2f.S \ log2l.S logb.c logbf.c logbl.c logf.S logl.S lrint.c lrintf.c \ @@ -54,9 +56,10 @@ MATH_DISTFILES = \ pow.c powf.c powi.c powif.c powil.c powl.c \ remainder.S remainderf.S remainderl.S remquo.S \ remquof.S remquol.S rint.c rintf.c rintl.c round.c roundf.c \ - roundl.c scalbn.S scalbnf.S scalbnl.S signbit.c signbitf.c \ - signbitl.c sinf.S sinhf.c sinhl.c sinl.S sqrtf.c sqrtl.c \ - tanf.S tanhf.c tanhl.c tanl.S trunc.c truncf.c truncl.c + roundl.c scalbn.S scalbnf.S scalbnl.S s_erf.c sf_erf.c \ + signbit.c signbitf.c signbitl.c sinf.S sinhf.c sinhl.c sinl.S \ + sqrtf.c sqrtl.c tanf.S tanhf.c tanhl.c tanl.S tgamma.c \ + tgammaf.c tgammal.c trunc.c truncf.c truncl.c CC = @CC@ # FIXME: Which is it, CC or CC_FOR_TARGET? @@ -93,6 +96,8 @@ Q8_OBJS = \ fwide.o imaxabs.o imaxdiv.o mbsinit.o \ strtoimax.o strtoumax.o wcstoimax.o wcstoumax.o \ wmemchr.o wmemcmp.o wmemcpy.o wmemmove.o wmemset.o +STDLIB_OBJS = \ + strtold.o wcstold.o STDLIB_STUB_OBJS = \ lltoa.o ulltoa.o \ lltow.o ulltow.o \ @@ -104,6 +109,7 @@ STDIO_STUB_OBJS = \ MATH_OBJS = \ acosf.o acosl.o asinf.o asinl.o atan2f.o atan2l.o \ atanf.o atanl.o cbrt.o cbrtf.o cbrtl.o ceilf.o ceill.o \ + cephes_emath.o \ copysign.o copysignf.o copysignl.o cosf.o coshf.o coshl.o cosl.o \ exp2.o exp2f.o exp2l.o expf.o expl.o fabs.o fabsf.o fabsl.o \ fdim.o fdimf.o fdiml.o floorf.o floorl.o fma.o fmaf.o fmal.o \ @@ -111,7 +117,8 @@ MATH_OBJS = \ fmodl.o fp_consts.o fp_constsf.o fp_constsl.o \ fpclassify.o fpclassifyf.o fpclassifyl.o \ frexpf.o frexpl.o fucom.o hypotf.o hypotl.o ilogb.o ilogbf.o \ - ilogbl.o isnan.o isnanf.o isnanl.o ldexpf.o ldexpl.o llrint.o \ + ilogbl.o isnan.o isnanf.o isnanl.o ldexpf.o ldexpl.o \ + lgamma.o lgammaf.o lgammal.o llrint.o \ llrintf.o llrintl.o llround.o llroundf.o llroundl.o \ log10f.o log10l.o log1p.o log1pf.o log1pl.o log2.o log2f.o \ log2l.o logb.o logbf.o logbl.o logf.o logl.o lrint.o lrintf.o \ @@ -120,9 +127,10 @@ MATH_OBJS = \ pow.o powf.o powi.o powif.o powil.o powl.o \ remainder.o remainderf.o remainderl.o remquo.o \ remquof.o remquol.o rint.o rintf.o rintl.o round.o roundf.o \ - roundl.o scalbn.o scalbnf.o scalbnl.o signbit.o signbitf.o \ - signbitl.o sinf.o sinhf.o sinhl.o sinl.o sqrtf.o sqrtl.o \ - tanf.o tanhf.o tanhl.o tanl.o trunc.o truncf.o truncl.o + roundl.o scalbn.o scalbnf.o scalbnl.o s_erf.o sf_erf.o \ + signbit.o signbitf.o signbitl.o sinf.o sinhf.o sinhl.o sinl.o \ + sqrtf.o sqrtl.o tanf.o tanhf.o tanhl.o tanl.o tgamma.o \ + tgammaf.o tgammal.o trunc.o truncf.o truncl.o FENV_OBJS = fesetround.o fegetround.o \ fegetenv.o fesetenv.o feupdateenv.o \ feclearexcept.o feholdexcept.o fegetexceptflag.o \ @@ -132,8 +140,8 @@ POSIX_OBJS = \ REPLACE_OBJS = \ mingw-fseek.o -LIB_OBJS = $(Q8_OBJS) $(STDLIB_STUB_OBJS) $(STDIO_STUB_OBJS) \ - $(MATH_OBJS) $(FENV_OBJS) $(POSIX_OBJS) \ +LIB_OBJS = $(Q8_OBJS) $(STDLIB_OBJS) $(STDLIB_STUB_OBJS) \ + $(STDIO_STUB_OBJS) $(MATH_OBJS) $(FENV_OBJS) $(POSIX_OBJS) \ $(REPLACE_OBJS) LIBS = $(LIBMINGWEX_A) @@ -184,6 +192,8 @@ distclean: # Dependancies # wdirent.o: $(srcdir)/dirent.c $(srcdir)/wdirent.c +strtold.o: $(srcdir)/strtold.c $(srcdir)/math/cephes_emath.h +wcstold.o: $(srcdir)/wcstold.c $(srcdir)/math/cephes_emath.h dist: diff --git a/winsup/mingw/mingwex/ldtoa.c b/winsup/mingw/mingwex/ldtoa.c new file mode 100644 index 00000000000..09dc7911b3c --- /dev/null +++ b/winsup/mingw/mingwex/ldtoa.c @@ -0,0 +1,614 @@ +/* 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 + * + */ + + +#ifdef USE_LDTOA + +#include "math/cephes_emath.h" + +#if NE == 10 + +/* 1.0E0 */ +static const unsigned short __eone[NE] = + {0x0000, 0x0000, 0x0000, 0x0000, + 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3fff,}; + +#else + +static const unsigned short __eone[NE] = { +0, 0000000,0000000,0000000,0100000,0x3fff,}; +#endif + + +#if NE == 10 +static const unsigned short __etens[NTEN + 1][NE] = +{ + {0x6576, 0x4a92, 0x804a, 0x153f, + 0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,}, /* 10**4096 */ + {0x6a32, 0xce52, 0x329a, 0x28ce, + 0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,}, /* 10**2048 */ + {0x526c, 0x50ce, 0xf18b, 0x3d28, + 0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,}, + {0x9c66, 0x58f8, 0xbc50, 0x5c54, + 0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,}, + {0x851e, 0xeab7, 0x98fe, 0x901b, + 0xddbb, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,}, + {0x0235, 0x0137, 0x36b1, 0x336c, + 0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,}, + {0x50f8, 0x25fb, 0xc76b, 0x6b71, + 0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,}, + {0x0000, 0x0000, 0x0000, 0x0000, + 0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,}, + {0x0000, 0x0000, 0x0000, 0x0000, + 0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,}, + {0x0000, 0x0000, 0x0000, 0x0000, + 0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,}, + {0x0000, 0x0000, 0x0000, 0x0000, + 0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,}, + {0x0000, 0x0000, 0x0000, 0x0000, + 0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,}, + {0x0000, 0x0000, 0x0000, 0x0000, + 0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,}, /* 10**1 */ +}; + +#else +static const unsigned short __etens[NTEN+1][NE] = { +{0xc94c,0x979a,0x8a20,0x5202,0xc460,0x7525,},/* 10**4096 */ +{0xa74d,0x5de4,0xc53d,0x3b5d,0x9e8b,0x5a92,},/* 10**2048 */ +{0x650d,0x0c17,0x8175,0x7586,0xc976,0x4d48,}, +{0xcc65,0x91c6,0xa60e,0xa0ae,0xe319,0x46a3,}, +{0xddbc,0xde8d,0x9df9,0xebfb,0xaa7e,0x4351,}, +{0xc66f,0x8cdf,0x80e9,0x47c9,0x93ba,0x41a8,}, +{0x3cbf,0xa6d5,0xffcf,0x1f49,0xc278,0x40d3,}, +{0xf020,0xb59d,0x2b70,0xada8,0x9dc5,0x4069,}, +{0x0000,0x0000,0x0400,0xc9bf,0x8e1b,0x4034,}, +{0x0000,0x0000,0x0000,0x2000,0xbebc,0x4019,}, +{0x0000,0x0000,0x0000,0x0000,0x9c40,0x400c,}, +{0x0000,0x0000,0x0000,0x0000,0xc800,0x4005,}, +{0x0000,0x0000,0x0000,0x0000,0xa000,0x4002,}, /* 10**1 */ +}; +#endif + +#if NE == 10 +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 __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 + +/* This routine will not return more than NDEC+1 digits. */ +void __etoasc(short unsigned int * __restrict__ x, + char * __restrict__ string, + const int ndigits, const int outformat, + int* outexp) +{ +long digit; +unsigned short y[NI], t[NI], u[NI], w[NI], equot[NI]; +const unsigned short *r, *p; +const unsigned short *ten; +unsigned short sign; +int i, j, k, expon, ndigs; +char *s, *ss; +unsigned short m; + +ndigs = ndigits; +#ifdef NANS +if( __eisnan(x) ) + { + sprintf( string, " NaN " ); + expon = 9999; + goto bxit; + } +#endif +__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>= 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 ); + __emul( p, t, t ); + 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 ); + expon -= 1; + } + } + else + { + __emovi( y, w ); + for( i=0; i 0 ) + { + if( __ecmp( p, w ) >= 0 ) + { + __emul( r, w, w ); + __emul( r, t, t ); + expon += k; + } + k /= 2; + if( k == 0 ) + break; + p += NE; + r += NE; + } + __ediv( t, __eone, t ); + } +isone: +/* Find the first (leading) digit. */ +__emovi( t, w ); +__emovz( w, t ); +__emovi( y, w ); +__emovz( w, y ); +__eiremain( t, y, equot); +digit = equot[NI-1]; +while( (digit == 0) && (__eiszero(y) == 0) ) + { + __eshup1( y ); + __emovz( y, u ); + __eshup1( u ); + __eshup1( u ); + __eaddm( u, y ); + __eiremain( t, y, equot); + 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, equot); + *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 ) + { + if( __eiiszero(y) == 0 ) + goto roun; /* round to nearest */ + if( (*(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: +#if defined (__GO32__) || defined (__MINGW32__) +if( expon >= 0 ) + sprintf( ss, "e+%02d", expon ); +else + sprintf( ss, "e-%02d", -expon ); +#else + sprintf( ss, "E%d", expon ); +#endif +bxit: + +if (outexp) + *outexp = expon; +} + + +/* FIXME: Not thread safe */ +static char outstr[128]; + + char * +_IO_ldtoa(long double d, int mode, int ndigits, int *decpt, + int *sign, char **rve) +{ +unsigned short e[NI]; +char *s, *p; +int k; +int outexpon = 0; + +union + { + unsigned short int us[6]; + long double ld; + } xx; +xx.ld = d; +__e64toe(xx.us, e ); +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; + +__etoasc( e, outstr, ndigits, mode, &outexpon ); +s = outstr; +if( __eisinf(e) || __eisnan(e) ) + { + *decpt = 9999; + goto stripspaces; + } +*decpt = outexpon + 1; + +/* Transform the string returned by etoasc into what the caller wants. */ + +/* Look for decimal point and delete it from the string. */ +s = outstr; +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 != 'e' && s > outstr) + --s; +*s = '\0'; + +stripspaces: + +/* Strip leading spaces and sign. */ +p = outstr; +while( *p == ' ' || *p == '-') + ++p; + +/* Find new end of string. */ +s = outstr; +while( (*s++ = *p++) != '\0' ) + ; +--s; + +/* Strip trailing zeros. */ +if( mode == 2 ) + k = 1; +else if( ndigits > outexpon ) + k = ndigits; +else + k = outexpon; + +while( *(s-1) == '0' && ((s - outstr) > 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 + outexpon) < 0)) + { + s = outstr; + *s = '\0'; + *decpt = 0; + } + +if( rve ) + *rve = s; +return outstr; +} + +void +_IO_ldtostr(long double *x, char *string, int ndigs, int flags, char fmt) +{ +unsigned short w[NI]; +char *t, *u; +int outexpon = 0; +int outformat = -1; +char dec_sym = *(localeconv()->decimal_point); + +__e64toe( (unsigned short *)x, w ); +__etoasc( w, string, ndigs, outformat, &outexpon ); + +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'; + } +if (dec_sym != '.') + { + t = string; + while (*t != '.') + ++t; + *t = dec_sym; + } +} + +#endif /* USE_LDTOA */ diff --git a/winsup/mingw/mingwex/math/cephes_emath.c b/winsup/mingw/mingwex/math/cephes_emath.c new file mode 100644 index 00000000000..ab798a2d255 --- /dev/null +++ b/winsup/mingw/mingwex/math/cephes_emath.c @@ -0,0 +1,1318 @@ +/* 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 + * + */ + + +#include "cephes_emath.h" + +/* + * The constants are for 64 bit precision. + */ + + +/* Move in external format number, + * converting it to internal format. + */ +void __emovi(const short unsigned int * __restrict__ a, + short unsigned int * __restrict__ 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 INFINITY +if( (*(q-1) & 0x7fff) == 0x7fff ) + { +#ifdef NANS + if( __eisnan(a) ) + { + *q++ = 0; + for( i=3; i> 16) + (m >> 16) + *pp; + *pp = (unsigned short )carry; + *(pp-1) = carry >> 16; + } + } +for( i=M; i 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 ); + } + } + __esubm( tprod, num ); + equot[i] = tquot; + __eshup6(num); + } +/* test for nonzero remainder after roundoff bit */ +p = &num[M]; +j = 0; +for( i=M; i NBITS ) + { + __ecleazs( s ); + return; + } +#endif +exp -= j; +#ifndef 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; +if (rndprc == 64) + { + rw = 7; + re = 6; + rbit[NI-2] = 0; + rbit[6] = 1; + } + +/* 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) && (rndprc != NBITS) ) +#else +if( (exp <= 0) && (rndprc != 64) && (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[rw] & rmsk; +if( rndprc < NBITS ) + { + i = rw + 1; + while( i < NI ) + { + if( s[i] ) + r |= 1; + s[i] = 0; + ++i; + } + } +s[rw] &= ~rmsk; +if( (r & rmbit) != 0 ) + { + if( r == rmbit ) + { + if( lost == 0 ) + { /* round to even */ + if( (s[re] & 1) == 0 ) + goto mddone; + } + else + { + if( subflg != 0 ) + goto mddone; + } + } + __eaddm( rbit, s ); + } +mddone: +#if IBMPC +if( (exp <= 0) && (rndprc != NBITS) ) +#else +if( (exp <= 0) && (rndprc != 64) && (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 INFINITY +overf: +#endif +#ifdef INFINITY + s[1] = 32767; + for( i=2; i b + * 0 if a == b + * -1 if a < b + * -2 if either a or b is a NaN. + */ +int __ecmp(const short unsigned int * __restrict__ a, + const short unsigned int * __restrict__ 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 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. +*/ +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. +*/ +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 ); +} + + +/* Move internal format number out, + * converting it to external format. + */ +void __emovo(const short unsigned int * __restrict__ a, + short unsigned int * __restrict__ b) +{ +register const unsigned short *p; +register unsigned short *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 INFINITY +if( *(p-1) == 0x7fff ) + { +#ifdef NANS + if( __eiisnan(a) ) + { + __enan_NBITS( b ); + return; + } +#endif + __einfin(b); + return; + } +#endif +/* skip over guard word */ +++p; +/* move the significand */ +for( i=0; i= 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, NBITS ); +} + + +void __eadd1(const short unsigned int * __restrict__ a, + const short unsigned int * __restrict__ b, + short unsigned int * __restrict__ c, + int subflg) +{ +unsigned short ai[NI], bi[NI], ci[NI]; +int i, lost, j, k; +long lt, lta, ltb; + +#ifdef 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( sub ) + 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 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, NBITS); + +done: +__emovo( bi, c ); +} + + +/* y = largest integer not greater than x + * (truncated toward minus infinity) + * + * unsigned short x[NE], y[NE] + * + * efloor( x, y ); + */ + + +void __efloor(short unsigned int *x, short unsigned int *y) +{ +register unsigned short *p; +int e, expon, i; +unsigned short f[NE]; +const unsigned short bmask[] = { +0xffff, +0xfffe, +0xfffc, +0xfff8, +0xfff0, +0xffe0, +0xffc0, +0xff80, +0xff00, +0xfe00, +0xfc00, +0xf800, +0xf000, +0xe000, +0xc000, +0x8000, +0x0000, +}; + +__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 + * + */ + + +/* 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 +#include +#include +#include +#include +#include +#include + +#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 +#define INFINITY + +/* 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 + +extern const unsigned short __etens[NTEN + 1][NE]; + +/* +; 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 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 *(--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>= 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>= 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> 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 + */ + + +#ifndef __MINGW32__ +#include "mconf.h" +#ifdef ANSIPROT +extern double pow ( double, double ); +extern double log ( double ); +extern double exp ( double ); +extern double sin ( double ); +extern double polevl ( double, void *, int ); +extern double p1evl ( double, void *, int ); +extern double floor ( double ); +extern double fabs ( double ); +extern int isnan ( double ); +extern int isfinite ( double ); +#else +double pow(), log(), exp(), sin(), polevl(), p1evl(), floor(), fabs(); +int isnan(), isfinite(); +#endif +#ifdef INFINITIES +extern double INFINITY; +#endif +#ifdef NANS +extern double NAN; +#endif +#else /* __MINGW32__ */ +#include "cephes_mconf.h" +#endif /* __MINGW32__ */ + + +/* A[]: Stirling's formula expansion of log gamma + * B[], C[]: log gamma function between 2 and 3 + */ +#ifdef UNK +static double A[] = { + 8.11614167470508450300E-4, +-5.95061904284301438324E-4, + 7.93650340457716943945E-4, +-2.77777777730099687205E-3, + 8.33333333333331927722E-2 +}; +static double B[] = { +-1.37825152569120859100E3, +-3.88016315134637840924E4, +-3.31612992738871184744E5, +-1.16237097492762307383E6, +-1.72173700820839662146E6, +-8.53555664245765465627E5 +}; +static double C[] = { +/* 1.00000000000000000000E0, */ +-3.51815701436523470549E2, +-1.70642106651881159223E4, +-2.20528590553854454839E5, +-1.13933444367982507207E6, +-2.53252307177582951285E6, +-2.01889141433532773231E6 +}; +/* log( sqrt( 2*pi ) ) */ +static double LS2PI = 0.91893853320467274178; +#define MAXLGM 2.556348e305 +static double LOGPI = 1.14472988584940017414; +#endif + +#ifdef DEC +static const unsigned short A[] = { +0035524,0141201,0034633,0031405, +0135433,0176755,0126007,0045030, +0035520,0006371,0003342,0172730, +0136066,0005540,0132605,0026407, +0037252,0125252,0125252,0125132 +}; +static const unsigned short B[] = { +0142654,0044014,0077633,0035410, +0144027,0110641,0125335,0144760, +0144641,0165637,0142204,0047447, +0145215,0162027,0146246,0155211, +0145322,0026110,0010317,0110130, +0145120,0061472,0120300,0025363 +}; +static const unsigned short C[] = { +/*0040200,0000000,0000000,0000000*/ +0142257,0164150,0163630,0112622, +0143605,0050153,0156116,0135272, +0144527,0056045,0145642,0062332, +0145213,0012063,0106250,0001025, +0145432,0111254,0044577,0115142, +0145366,0071133,0050217,0005122 +}; +/* log( sqrt( 2*pi ) ) */ +static const unsigned short LS2P[] = {040153,037616,041445,0172645,}; +#define LS2PI *(double *)LS2P +#define MAXLGM 2.035093e36 +static const unsigned short LPI[4] = { +0040222,0103202,0043475,0006750, +}; +#define LOGPI *(double *)LPI + +#endif + +#ifdef IBMPC +static const unsigned short A[] = { +0x6661,0x2733,0x9850,0x3f4a, +0xe943,0xb580,0x7fbd,0xbf43, +0x5ebb,0x20dc,0x019f,0x3f4a, +0xa5a1,0x16b0,0xc16c,0xbf66, +0x554b,0x5555,0x5555,0x3fb5 +}; +static const unsigned short B[] = { +0x6761,0x8ff3,0x8901,0xc095, +0xb93e,0x355b,0xf234,0xc0e2, +0x89e5,0xf890,0x3d73,0xc114, +0xdb51,0xf994,0xbc82,0xc131, +0xf20b,0x0219,0x4589,0xc13a, +0x055e,0x5418,0x0c67,0xc12a +}; +static const unsigned short C[] = { +/*0x0000,0x0000,0x0000,0x3ff0,*/ +0x12b2,0x1cf3,0xfd0d,0xc075, +0xd757,0x7b89,0xaa0d,0xc0d0, +0x4c9b,0xb974,0xeb84,0xc10a, +0x0043,0x7195,0x6286,0xc131, +0xf34c,0x892f,0x5255,0xc143, +0xe14a,0x6a11,0xce4b,0xc13e +}; +/* log( sqrt( 2*pi ) ) */ +static const unsigned short LS2P[] = { +0xbeb5,0xc864,0x67f1,0x3fed +}; +#define LS2PI *(double *)LS2P +#define MAXLGM 2.556348e305 +static const unsigned short LPI[4] = { +0xa1bd,0x48e7,0x50d0,0x3ff2, +}; +#define LOGPI *(double *)LPI +#endif + +#ifdef MIEEE +static const unsigned short A[] = { +0x3f4a,0x9850,0x2733,0x6661, +0xbf43,0x7fbd,0xb580,0xe943, +0x3f4a,0x019f,0x20dc,0x5ebb, +0xbf66,0xc16c,0x16b0,0xa5a1, +0x3fb5,0x5555,0x5555,0x554b +}; +static const unsigned short B[] = { +0xc095,0x8901,0x8ff3,0x6761, +0xc0e2,0xf234,0x355b,0xb93e, +0xc114,0x3d73,0xf890,0x89e5, +0xc131,0xbc82,0xf994,0xdb51, +0xc13a,0x4589,0x0219,0xf20b, +0xc12a,0x0c67,0x5418,0x055e +}; +static const unsigned short C[] = { +0xc075,0xfd0d,0x1cf3,0x12b2, +0xc0d0,0xaa0d,0x7b89,0xd757, +0xc10a,0xeb84,0xb974,0x4c9b, +0xc131,0x6286,0x7195,0x0043, +0xc143,0x5255,0x892f,0xf34c, +0xc13e,0xce4b,0x6a11,0xe14a +}; +/* log( sqrt( 2*pi ) ) */ +static const unsigned short LS2P[] = { +0x3fed,0x67f1,0xc864,0xbeb5 +}; +#define LS2PI *(double *)LS2P +#define MAXLGM 2.556348e305 +static unsigned short LPI[4] = { +0x3ff2,0x50d0,0x48e7,0xa1bd, +}; +#define LOGPI *(double *)LPI +#endif + + +/* Logarithm of gamma function */ +/* Reentrant version */ + +double __lgamma_r(double x, int* sgngam) +{ +double p, q, u, w, z; +int i; + +*sgngam = 1; +#ifdef NANS +if( isnan(x) ) + return(x); +#endif + +#ifdef INFINITIES +if( !isfinite(x) ) + return(INFINITY); +#endif + +if( x < -34.0 ) + { + q = -x; + w = __lgamma_r(q, sgngam); /* note this modifies sgngam! */ + p = floor(q); + if( p == q ) + { +lgsing: + _SET_ERRNO(EDOM); + mtherr( "lgam", SING ); +#ifdef INFINITIES + return (INFINITY); +#else + return (MAXNUM); +#endif + } + i = p; + if( (i & 1) == 0 ) + *sgngam = -1; + else + *sgngam = 1; + z = q - p; + if( z > 0.5 ) + { + p += 1.0; + z = p - q; + } + z = q * sin( PI * z ); + if( z == 0.0 ) + goto lgsing; +/* z = log(PI) - log( z ) - w;*/ + z = LOGPI - log( z ) - w; + return( z ); + } + +if( x < 13.0 ) + { + z = 1.0; + p = 0.0; + u = x; + while( u >= 3.0 ) + { + p -= 1.0; + u = x + p; + z *= u; + } + while( u < 2.0 ) + { + if( u == 0.0 ) + goto lgsing; + z /= u; + p += 1.0; + u = x + p; + } + if( z < 0.0 ) + { + *sgngam = -1; + z = -z; + } + else + *sgngam = 1; + if( u == 2.0 ) + return( log(z) ); + p -= 2.0; + x = x + p; + p = x * polevl( x, B, 5 ) / p1evl( x, C, 6); + return( log(z) + p ); + } + +if( x > MAXLGM ) + { + _SET_ERRNO(ERANGE); + mtherr( "lgamma", OVERFLOW ); +#ifdef INFINITIES + return( *sgngam * INFINITY ); +#else + return( *sgngam * MAXNUM ); +#endif + } + +q = ( x - 0.5 ) * log(x) - x + LS2PI; +if( x > 1.0e8 ) + return( q ); + +p = 1.0/(x*x); +if( x >= 1000.0 ) + q += (( 7.9365079365079365079365e-4 * p + - 2.7777777777777777777778e-3) *p + + 0.0833333333333333333333) / x; +else + q += polevl( p, A, 4 ) / x; +return( q ); +} + +/* This is the C99 version */ + +double lgamma(double x) +{ + int local_sgngam=0; + return (__lgamma_r(x, &local_sgngam)); +} diff --git a/winsup/mingw/mingwex/math/lgammaf.c b/winsup/mingw/mingwex/math/lgammaf.c new file mode 100644 index 00000000000..20982f999fc --- /dev/null +++ b/winsup/mingw/mingwex/math/lgammaf.c @@ -0,0 +1,253 @@ +/* lgamf() + * + * Natural logarithm of gamma function + * + * + * + * SYNOPSIS: + * + * float x, y, __lgammaf_r(); + * int* sgngamf; + * y = __lgammaf_r( x, sgngamf ); + * + * float x, y, lgammaf(); + * y = lgammaf( x); + * + * + * + * DESCRIPTION: + * + * Returns the base e (2.718...) logarithm of the absolute + * value of the gamma function of the argument. In the reentrant + * version the sign (+1 or -1) of the gamma function is returned in + * variable referenced by sgngamf. + * + * For arguments greater than 6.5, the logarithm of the gamma + * function is approximated by the logarithmic version of + * Stirling's formula. Arguments between 0 and +6.5 are reduced by + * by recurrence to the interval [.75,1.25] or [1.5,2.5] of a rational + * approximation. The cosecant reflection formula is employed for + * arguments less than zero. + * + * Arguments greater than MAXLGM = 2.035093e36 return MAXNUM and an + * error message. + * + * + * + * ACCURACY: + * + * + * + * arithmetic domain # trials peak rms + * IEEE -100,+100 500,000 7.4e-7 6.8e-8 + * The error criterion was relative when the function magnitude + * was greater than one but absolute when it was less than one. + * The routine has low relative error for positive arguments. + * + * The following test used the relative error criterion. + * IEEE -2, +3 100000 4.0e-7 5.6e-8 + * + */ + + +/* + Cephes Math Library Release 2.7: July, 1998 + Copyright 1984, 1987, 1989, 1992, 1998 by Stephen L. Moshier +*/ + +/* + 26-11-2002 Modified for mingw. + Danny Smith +*/ + + +/* log gamma(x+2), -.5 < x < .5 */ +static const float B[] = { + 6.055172732649237E-004, +-1.311620815545743E-003, + 2.863437556468661E-003, +-7.366775108654962E-003, + 2.058355474821512E-002, +-6.735323259371034E-002, + 3.224669577325661E-001, + 4.227843421859038E-001 +}; + +/* log gamma(x+1), -.25 < x < .25 */ +static const float C[] = { + 1.369488127325832E-001, +-1.590086327657347E-001, + 1.692415923504637E-001, +-2.067882815621965E-001, + 2.705806208275915E-001, +-4.006931650563372E-001, + 8.224670749082976E-001, +-5.772156501719101E-001 +}; + +/* log( sqrt( 2*pi ) ) */ +static const float LS2PI = 0.91893853320467274178; +#define MAXLGM 2.035093e36 +static const float PIINV = 0.318309886183790671538; + +#ifndef __MINGW32__ +#include "mconf.h" +float floorf(float); +float polevlf( float, float *, int ); +float p1evlf( float, float *, int ); +#else +#include "cephes_mconf.h" +#endif + +/* Reentrant version */ +/* Logarithm of gamma function */ + +float __lgammaf_r( float x, int* sgngamf ) +{ +float p, q, w, z; +float nx, tx; +int i, direction; + +*sgngamf = 1; +#ifdef NANS +if( isnan(x) ) + return(x); +#endif + +#ifdef INFINITIES +if( !isfinite(x) ) + return(x); +#endif + + +if( x < 0.0 ) + { + q = -x; + w = __lgammaf_r(q, sgngamf); /* note this modifies sgngam! */ + p = floorf(q); + if( p == q ) + { +lgsing: + _SET_ERRNO(EDOM); + mtherr( "lgamf", SING ); +#ifdef INFINITIES + return (INFINITYF); +#else + return( *sgngamf * MAXNUMF ); +#endif + } + i = p; + if( (i & 1) == 0 ) + *sgngamf = -1; + else + *sgngamf = 1; + z = q - p; + if( z > 0.5 ) + { + p += 1.0; + z = p - q; + } + z = q * sinf( PIF * z ); + if( z == 0.0 ) + goto lgsing; + z = -logf( PIINV*z ) - w; + return( z ); + } + +if( x < 6.5 ) + { + direction = 0; + z = 1.0; + tx = x; + nx = 0.0; + if( x >= 1.5 ) + { + while( tx > 2.5 ) + { + nx -= 1.0; + tx = x + nx; + z *=tx; + } + x += nx - 2.0; +iv1r5: + p = x * polevlf( x, B, 7 ); + goto cont; + } + if( x >= 1.25 ) + { + z *= x; + x -= 1.0; /* x + 1 - 2 */ + direction = 1; + goto iv1r5; + } + if( x >= 0.75 ) + { + x -= 1.0; + p = x * polevlf( x, C, 7 ); + q = 0.0; + goto contz; + } + while( tx < 1.5 ) + { + if( tx == 0.0 ) + goto lgsing; + z *=tx; + nx += 1.0; + tx = x + nx; + } + direction = 1; + x += nx - 2.0; + p = x * polevlf( x, B, 7 ); + +cont: + if( z < 0.0 ) + { + *sgngamf = -1; + z = -z; + } + else + { + *sgngamf = 1; + } + q = logf(z); + if( direction ) + q = -q; +contz: + return( p + q ); + } + +if( x > MAXLGM ) + { + _SET_ERRNO(ERANGE); + mtherr( "lgamf", OVERFLOW ); +#ifdef INFINITIES + return( *sgngamf * INFINITYF ); +#else + return( *sgngamf * MAXNUMF ); +#endif + + } + +/* Note, though an asymptotic formula could be used for x >= 3, + * there is cancellation error in the following if x < 6.5. */ +q = LS2PI - x; +q += ( x - 0.5 ) * logf(x); + +if( x <= 1.0e4 ) + { + z = 1.0/x; + p = z * z; + q += (( 6.789774945028216E-004 * p + - 2.769887652139868E-003 ) * p + + 8.333316229807355E-002 ) * z; + } +return( q ); +} + +/* This is the C99 version */ + +float lgammaf(float x) +{ + int local_sgngamf=0; + return (__lgammaf_r(x, &local_sgngamf)); +} diff --git a/winsup/mingw/mingwex/math/lgammal.c b/winsup/mingw/mingwex/math/lgammal.c new file mode 100644 index 00000000000..d2b306afd77 --- /dev/null +++ b/winsup/mingw/mingwex/math/lgammal.c @@ -0,0 +1,416 @@ +/* lgaml() + * + * Natural logarithm of gamma function + * + * + * + * SYNOPSIS: + * + * long double x, y, __lgammal_r(); + * int* sgngaml; + * y = __lgammal_r( x, sgngaml ); + * + * long double x, y, lgammal(); + * y = lgammal( x); + * + * + * + * DESCRIPTION: + * + * Returns the base e (2.718...) logarithm of the absolute + * value of the gamma function of the argument. In the reentrant + * version, the sign (+1 or -1) of the gamma function is returned + * in the variable referenced by sgngaml. + * + * For arguments greater than 33, the logarithm of the gamma + * function is approximated by the logarithmic version of + * Stirling's formula using a polynomial approximation of + * degree 4. Arguments between -33 and +33 are reduced by + * recurrence to the interval [2,3] of a rational approximation. + * The cosecant reflection formula is employed for arguments + * less than -33. + * + * Arguments greater than MAXLGML (10^4928) return MAXNUML. + * + * + * + * ACCURACY: + * + * + * arithmetic domain # trials peak rms + * IEEE -40, 40 100000 2.2e-19 4.6e-20 + * IEEE 10^-2000,10^+2000 20000 1.6e-19 3.3e-20 + * The error criterion was relative when the function magnitude + * was greater than one but absolute when it was less than one. + * + */ + +/* + * Copyright 1994 by Stephen L. Moshier + */ + +/* + * 26-11-2002 Modified for mingw. + * Danny Smith + */ + +#ifndef __MINGW32__ +#include "mconf.h" +#ifdef ANSIPROT +extern long double fabsl ( long double ); +extern long double lgaml ( long double ); +extern long double logl ( long double ); +extern long double expl ( long double ); +extern long double gammal ( long double ); +extern long double sinl ( long double ); +extern long double floorl ( long double ); +extern long double powl ( long double, long double ); +extern long double polevll ( long double, void *, int ); +extern long double p1evll ( long double, void *, int ); +extern int isnanl ( long double ); +extern int isfinitel ( long double ); +#else +long double fabsl(), lgaml(), logl(), expl(), gammal(), sinl(); +long double floorl(), powl(), polevll(), p1evll(), isnanl(), isfinitel(); +#endif +#ifdef INFINITIES +extern long double INFINITYL; +#endif +#ifdef NANS +extern long double NANL; +#endif +#else /* __MINGW32__ */ +#include "cephes_mconf.h" +#endif /* __MINGW32__ */ + +#if UNK +static long double S[9] = { +-1.193945051381510095614E-3L, + 7.220599478036909672331E-3L, +-9.622023360406271645744E-3L, +-4.219773360705915470089E-2L, + 1.665386113720805206758E-1L, +-4.200263503403344054473E-2L, +-6.558780715202540684668E-1L, + 5.772156649015328608253E-1L, + 1.000000000000000000000E0L, +}; +#endif +#if IBMPC +static const unsigned short S[] = { +0xbaeb,0xd6d3,0x25e5,0x9c7e,0xbff5, XPD +0xfe9a,0xceb4,0xc74e,0xec9a,0x3ff7, XPD +0x9225,0xdfef,0xb0e9,0x9da5,0xbff8, XPD +0x10b0,0xec17,0x87dc,0xacd7,0xbffa, XPD +0x6b8d,0x7515,0x1905,0xaa89,0x3ffc, XPD +0xf183,0x126b,0xf47d,0xac0a,0xbffa, XPD +0x7bf6,0x57d1,0xa013,0xa7e7,0xbffe, XPD +0xc7a9,0x7db0,0x67e3,0x93c4,0x3ffe, XPD +0x0000,0x0000,0x0000,0x8000,0x3fff, XPD +}; +#endif +#if MIEEE +static long S[27] = { +0xbff50000,0x9c7e25e5,0xd6d3baeb, +0x3ff70000,0xec9ac74e,0xceb4fe9a, +0xbff80000,0x9da5b0e9,0xdfef9225, +0xbffa0000,0xacd787dc,0xec1710b0, +0x3ffc0000,0xaa891905,0x75156b8d, +0xbffa0000,0xac0af47d,0x126bf183, +0xbffe0000,0xa7e7a013,0x57d17bf6, +0x3ffe0000,0x93c467e3,0x7db0c7a9, +0x3fff0000,0x80000000,0x00000000, +}; +#endif + +#if UNK +static long double SN[9] = { + 1.133374167243894382010E-3L, + 7.220837261893170325704E-3L, + 9.621911155035976733706E-3L, +-4.219773343731191721664E-2L, +-1.665386113944413519335E-1L, +-4.200263503402112910504E-2L, + 6.558780715202536547116E-1L, + 5.772156649015328608727E-1L, +-1.000000000000000000000E0L, +}; +#endif +#if IBMPC +static const unsigned SN[] = { +0x5dd1,0x02de,0xb9f7,0x948d,0x3ff5, XPD +0x989b,0xdd68,0xc5f1,0xec9c,0x3ff7, XPD +0x2ca1,0x18f0,0x386f,0x9da5,0x3ff8, XPD +0x783f,0x41dd,0x87d1,0xacd7,0xbffa, XPD +0x7a5b,0xd76d,0x1905,0xaa89,0xbffc, XPD +0x7f64,0x1234,0xf47d,0xac0a,0xbffa, XPD +0x5e26,0x57d1,0xa013,0xa7e7,0x3ffe, XPD +0xc7aa,0x7db0,0x67e3,0x93c4,0x3ffe, XPD +0x0000,0x0000,0x0000,0x8000,0xbfff, XPD +}; +#endif +#if MIEEE +static long SN[27] = { +0x3ff50000,0x948db9f7,0x02de5dd1, +0x3ff70000,0xec9cc5f1,0xdd68989b, +0x3ff80000,0x9da5386f,0x18f02ca1, +0xbffa0000,0xacd787d1,0x41dd783f, +0xbffc0000,0xaa891905,0xd76d7a5b, +0xbffa0000,0xac0af47d,0x12347f64, +0x3ffe0000,0xa7e7a013,0x57d15e26, +0x3ffe0000,0x93c467e3,0x7db0c7aa, +0xbfff0000,0x80000000,0x00000000, +}; +#endif + + +/* A[]: Stirling's formula expansion of log gamma + * B[], C[]: log gamma function between 2 and 3 + */ + + +/* log gamma(x) = ( x - 0.5 ) * log(x) - x + LS2PI + 1/x A(1/x^2) + * x >= 8 + * Peak relative error 1.51e-21 + * Relative spread of error peaks 5.67e-21 + */ +#if UNK +static long double A[7] = { + 4.885026142432270781165E-3L, +-1.880801938119376907179E-3L, + 8.412723297322498080632E-4L, +-5.952345851765688514613E-4L, + 7.936507795855070755671E-4L, +-2.777777777750349603440E-3L, + 8.333333333333331447505E-2L, +}; +#endif +#if IBMPC +static const unsigned short A[] = { +0xd984,0xcc08,0x91c2,0xa012,0x3ff7, XPD +0x3d91,0x0304,0x3da1,0xf685,0xbff5, XPD +0x3bdc,0xaad1,0xd492,0xdc88,0x3ff4, XPD +0x8b20,0x9fce,0x844e,0x9c09,0xbff4, XPD +0xf8f2,0x30e5,0x0092,0xd00d,0x3ff4, XPD +0x4d88,0x03a8,0x60b6,0xb60b,0xbff6, XPD +0x9fcc,0xaaaa,0xaaaa,0xaaaa,0x3ffb, XPD +}; +#endif +#if MIEEE +static long A[21] = { +0x3ff70000,0xa01291c2,0xcc08d984, +0xbff50000,0xf6853da1,0x03043d91, +0x3ff40000,0xdc88d492,0xaad13bdc, +0xbff40000,0x9c09844e,0x9fce8b20, +0x3ff40000,0xd00d0092,0x30e5f8f2, +0xbff60000,0xb60b60b6,0x03a84d88, +0x3ffb0000,0xaaaaaaaa,0xaaaa9fcc, +}; +#endif + +/* log gamma(x+2) = x B(x)/C(x) + * 0 <= x <= 1 + * Peak relative error 7.16e-22 + * Relative spread of error peaks 4.78e-20 + */ +#if UNK +static long double B[7] = { +-2.163690827643812857640E3L, +-8.723871522843511459790E4L, +-1.104326814691464261197E6L, +-6.111225012005214299996E6L, +-1.625568062543700591014E7L, +-2.003937418103815175475E7L, +-8.875666783650703802159E6L, +}; +static long double C[7] = { +/* 1.000000000000000000000E0L,*/ +-5.139481484435370143617E2L, +-3.403570840534304670537E4L, +-6.227441164066219501697E5L, +-4.814940379411882186630E6L, +-1.785433287045078156959E7L, +-3.138646407656182662088E7L, +-2.099336717757895876142E7L, +}; +#endif +#if IBMPC +static const unsigned short B[] = { +0x9557,0x4995,0x0da1,0x873b,0xc00a, XPD +0xfe44,0x9af8,0x5b8c,0xaa63,0xc00f, XPD +0x5aa8,0x7cf5,0x3684,0x86ce,0xc013, XPD +0x259a,0x258c,0xf206,0xba7f,0xc015, XPD +0xbe18,0x1ca3,0xc0a0,0xf80a,0xc016, XPD +0x168f,0x2c42,0x6717,0x98e3,0xc017, XPD +0x2051,0x9d55,0x92c8,0x876e,0xc016, XPD +}; +static const unsigned short C[] = { +/*0x0000,0x0000,0x0000,0x8000,0x3fff, XPD*/ +0xaa77,0xcf2f,0xae76,0x807c,0xc008, XPD +0xb280,0x0d74,0xb55a,0x84f3,0xc00e, XPD +0xa505,0xcd30,0x81dc,0x9809,0xc012, XPD +0x3369,0x4246,0xb8c2,0x92f0,0xc015, XPD +0x63cf,0x6aee,0xbe6f,0x8837,0xc017, XPD +0x26bb,0xccc7,0xb009,0xef75,0xc017, XPD +0x462b,0xbae8,0xab96,0xa02a,0xc017, XPD +}; +#endif +#if MIEEE +static long B[21] = { +0xc00a0000,0x873b0da1,0x49959557, +0xc00f0000,0xaa635b8c,0x9af8fe44, +0xc0130000,0x86ce3684,0x7cf55aa8, +0xc0150000,0xba7ff206,0x258c259a, +0xc0160000,0xf80ac0a0,0x1ca3be18, +0xc0170000,0x98e36717,0x2c42168f, +0xc0160000,0x876e92c8,0x9d552051, +}; +static long C[21] = { +/*0x3fff0000,0x80000000,0x00000000,*/ +0xc0080000,0x807cae76,0xcf2faa77, +0xc00e0000,0x84f3b55a,0x0d74b280, +0xc0120000,0x980981dc,0xcd30a505, +0xc0150000,0x92f0b8c2,0x42463369, +0xc0170000,0x8837be6f,0x6aee63cf, +0xc0170000,0xef75b009,0xccc726bb, +0xc0170000,0xa02aab96,0xbae8462b, +}; +#endif + +/* log( sqrt( 2*pi ) ) */ +static const long double LS2PI = 0.91893853320467274178L; +#define MAXLGM 1.04848146839019521116e+4928L + + +/* Logarithm of gamma function */ +/* Reentrant version */ + +long double __lgammal_r(long double x, int* sgngaml) +{ +long double p, q, w, z, f, nx; +int i; + +*sgngaml = 1; +#ifdef NANS +if( isnanl(x) ) + return(NANL); +#endif +#ifdef INFINITIES +if( !isfinitel(x) ) + return(INFINITYL); +#endif +if( x < -34.0L ) + { + q = -x; + w = __lgammal_r(q, sgngaml); /* note this modifies sgngam! */ + p = floorl(q); + if( p == q ) + { +lgsing: + _SET_ERRNO(EDOM); + mtherr( "lgammal", SING ); +#ifdef INFINITIES + return (INFINITYL); +#else + return (MAXNUML); +#endif + } + i = p; + if( (i & 1) == 0 ) + *sgngaml = -1; + else + *sgngaml = 1; + z = q - p; + if( z > 0.5L ) + { + p += 1.0L; + z = p - q; + } + z = q * sinl( PIL * z ); + if( z == 0.0L ) + goto lgsing; +/* z = LOGPI - logl( z ) - w; */ + z = logl( PIL/z ) - w; + return( z ); + } + +if( x < 13.0L ) + { + z = 1.0L; + nx = floorl( x + 0.5L ); + f = x - nx; + while( x >= 3.0L ) + { + nx -= 1.0L; + x = nx + f; + z *= x; + } + while( x < 2.0L ) + { + if( fabsl(x) <= 0.03125 ) + goto lsmall; + z /= nx + f; + nx += 1.0L; + x = nx + f; + } + if( z < 0.0L ) + { + *sgngaml = -1; + z = -z; + } + else + *sgngaml = 1; + if( x == 2.0L ) + return( logl(z) ); + x = (nx - 2.0L) + f; + p = x * polevll( x, B, 6 ) / p1evll( x, C, 7); + return( logl(z) + p ); + } + +if( x > MAXLGM ) + { + _SET_ERRNO(ERANGE); + mtherr( "lgammal", OVERFLOW ); +#ifdef INFINITIES + return( *sgngaml * INFINITYL ); +#else + return( *sgngaml * MAXNUML ); +#endif + } + +q = ( x - 0.5L ) * logl(x) - x + LS2PI; +if( x > 1.0e10L ) + return(q); +p = 1.0L/(x*x); +q += polevll( p, A, 6 ) / x; +return( q ); + + +lsmall: +if( x == 0.0L ) + goto lgsing; +if( x < 0.0L ) + { + x = -x; + q = z / (x * polevll( x, SN, 8 )); + } +else + q = z / (x * polevll( x, S, 8 )); +if( q < 0.0L ) + { + *sgngaml = -1; + q = -q; + } +else + *sgngaml = 1; +q = logl( q ); +return(q); +} + +/* This is the C99 version */ + +long double lgammal(long double x) +{ + int local_sgngaml=0; + return (__lgammal_r(x, &local_sgngaml)); +} diff --git a/winsup/mingw/mingwex/math/powif2.c b/winsup/mingw/mingwex/math/powif2.c new file mode 100644 index 00000000000..a181d57e292 --- /dev/null +++ b/winsup/mingw/mingwex/math/powif2.c @@ -0,0 +1,164 @@ +/* powif.c + * + * Real raised to integer power + * + * + * + * SYNOPSIS: + * + * float x, y, powif(); + * int n; + * + * y = powif( x, n ); + * + * + * + * DESCRIPTION: + * + * Returns argument x raised to the nth power. + * The routine efficiently decomposes n as a sum of powers of + * two. The desired power is a product of two-to-the-kth + * powers of x. Thus to compute the 32767 power of x requires + * 28 multiplications instead of 32767 multiplications. + * + * + * + * ACCURACY: + * + * + * Relative error: + * arithmetic x domain n domain # trials peak rms + * IEEE .04,26 -26,26 100000 1.1e-6 2.0e-7 + * IEEE 1,2 -128,128 100000 1.1e-5 1.0e-6 + * + * Returns MAXNUMF on overflow, zero on underflow. + * + */ + +/* powi.c */ + +/* +Cephes Math Library Release 2.2: June, 1992 +Copyright 1984, 1987, 1989 by Stephen L. Moshier +Direct inquiries to 30 Frost Street, Cambridge, MA 02140 +*/ + +#include "mconf.h" +extern float MAXNUMF, MAXLOGF, MINLOGF, LOGE2F; + +#ifdef ANSIC +float frexpf( float, int * ); + +float powif( float x, int nn ) +#else +float frexpf(); + +float powif( x, nn ) +double x; +int nn; +#endif +{ +int n, e, sign, asign, lx; +float w, y, s; + +if( x == 0.0 ) + { + if( nn == 0 ) + return( 1.0 ); + else if( nn < 0 ) + return( MAXNUMF ); + else + return( 0.0 ); + } + +if( nn == 0 ) + return( 1.0 ); + + +if( x < 0.0 ) + { + asign = -1; + x = -x; + } +else + asign = 0; + + +if( nn < 0 ) + { + sign = -1; + n = -nn; +/* + x = 1.0/x; +*/ + } +else + { + sign = 0; + n = nn; + } + +/* Overflow detection */ + +/* Calculate approximate logarithm of answer */ +s = frexpf( x, &lx ); +e = (lx - 1)*n; +if( (e == 0) || (e > 64) || (e < -64) ) + { + s = (s - 7.0710678118654752e-1) / (s + 7.0710678118654752e-1); + s = (2.9142135623730950 * s - 0.5 + lx) * nn * LOGE2F; + } +else + { + s = LOGE2F * e; + } + +if( s > MAXLOGF ) + { + mtherr( "powi", OVERFLOW ); + y = MAXNUMF; + goto done; + } + +if( s < MINLOGF ) + return(0.0); + +/* Handle tiny denormal answer, but with less accuracy + * since roundoff error in 1.0/x will be amplified. + * The precise demarcation should be the gradual underflow threshold. + */ +if( s < (-MAXLOGF+2.0) ) + { + x = 1.0/x; + sign = 0; + } + +/* First bit of the power */ +if( n & 1 ) + y = x; + +else + { + y = 1.0; + asign = 0; + } + +w = x; +n >>= 1; +while( n ) + { + w = w * w; /* arg to the 2-to-the-kth power */ + if( n & 1 ) /* if that bit is set, then include in product */ + y *= w; + n >>= 1; + } + + +done: + +if( asign ) + y = -y; /* odd power of negative number */ +if( sign ) + y = 1.0/y; +return(y); +} diff --git a/winsup/mingw/mingwex/math/s_erf.c b/winsup/mingw/mingwex/math/s_erf.c new file mode 100644 index 00000000000..4673f48b372 --- /dev/null +++ b/winsup/mingw/mingwex/math/s_erf.c @@ -0,0 +1,342 @@ + +/* @(#)s_erf.c 1.3 95/01/18 */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* double erf(double x) + * double erfc(double x) + * x + * 2 |\ + * erf(x) = --------- | exp(-t*t)dt + * sqrt(pi) \| + * 0 + * + * erfc(x) = 1-erf(x) + * Note that + * erf(-x) = -erf(x) + * erfc(-x) = 2 - erfc(x) + * + * Method: + * 1. For |x| in [0, 0.84375] + * erf(x) = x + x*R(x^2) + * erfc(x) = 1 - erf(x) if x in [-.84375,0.25] + * = 0.5 + ((0.5-x)-x*R) if x in [0.25,0.84375] + * where R = P/Q where P is an odd poly of degree 8 and + * Q is an odd poly of degree 10. + * -57.90 + * | R - (erf(x)-x)/x | <= 2 + * + * + * Remark. The formula is derived by noting + * erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....) + * and that + * 2/sqrt(pi) = 1.128379167095512573896158903121545171688 + * is close to one. The interval is chosen because the fix + * point of erf(x) is near 0.6174 (i.e., erf(x)=x when x is + * near 0.6174), and by some experiment, 0.84375 is chosen to + * guarantee the error is less than one ulp for erf. + * + * 2. For |x| in [0.84375,1.25], let s = |x| - 1, and + * c = 0.84506291151 rounded to single (24 bits) + * erf(x) = sign(x) * (c + P1(s)/Q1(s)) + * erfc(x) = (1-c) - P1(s)/Q1(s) if x > 0 + * 1+(c+P1(s)/Q1(s)) if x < 0 + * |P1/Q1 - (erf(|x|)-c)| <= 2**-59.06 + * Remark: here we use the taylor series expansion at x=1. + * erf(1+s) = erf(1) + s*Poly(s) + * = 0.845.. + P1(s)/Q1(s) + * That is, we use rational approximation to approximate + * erf(1+s) - (c = (single)0.84506291151) + * Note that |P1/Q1|< 0.078 for x in [0.84375,1.25] + * where + * P1(s) = degree 6 poly in s + * Q1(s) = degree 6 poly in s + * + * 3. For x in [1.25,1/0.35(~2.857143)], + * erfc(x) = (1/x)*exp(-x*x-0.5625+R1/S1) + * erf(x) = 1 - erfc(x) + * where + * R1(z) = degree 7 poly in z, (z=1/x^2) + * S1(z) = degree 8 poly in z + * + * 4. For x in [1/0.35,28] + * erfc(x) = (1/x)*exp(-x*x-0.5625+R2/S2) if x > 0 + * = 2.0 - (1/x)*exp(-x*x-0.5625+R2/S2) if -6 x >= 28 + * erf(x) = sign(x) *(1 - tiny) (raise inexact) + * erfc(x) = tiny*tiny (raise underflow) if x > 0 + * = 2 - tiny if x<0 + * + * 7. Special case: + * erf(0) = 0, erf(inf) = 1, erf(-inf) = -1, + * erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2, + * erfc/erf(NaN) is NaN + */ + + +/* #include "fdlibm.h" */ + +#include +#include + +#define __ieee754_exp exp + +typedef union +{ + double value; + struct + { + uint32_t lsw; + uint32_t msw; + } parts; +} ieee_double_shape_type; + + +static inline int __get_hi_word(const double x) +{ + ieee_double_shape_type u; + u.value = x; + return u.parts.msw; +} + +static inline void __trunc_lo_word(double *x) +{ + ieee_double_shape_type u; + u.value = *x; + u.parts.lsw = 0; + *x = u.value; +} + + +#ifdef __STDC__ +static const double +#else +static double +#endif +tiny = 1e-300, +half= 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */ +one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ +two = 2.00000000000000000000e+00, /* 0x40000000, 0x00000000 */ + /* c = (float)0.84506291151 */ +erx = 8.45062911510467529297e-01, /* 0x3FEB0AC1, 0x60000000 */ +/* + * Coefficients for approximation to erf on [0,0.84375] + */ +efx = 1.28379167095512586316e-01, /* 0x3FC06EBA, 0x8214DB69 */ +efx8= 1.02703333676410069053e+00, /* 0x3FF06EBA, 0x8214DB69 */ +pp0 = 1.28379167095512558561e-01, /* 0x3FC06EBA, 0x8214DB68 */ +pp1 = -3.25042107247001499370e-01, /* 0xBFD4CD7D, 0x691CB913 */ +pp2 = -2.84817495755985104766e-02, /* 0xBF9D2A51, 0xDBD7194F */ +pp3 = -5.77027029648944159157e-03, /* 0xBF77A291, 0x236668E4 */ +pp4 = -2.37630166566501626084e-05, /* 0xBEF8EAD6, 0x120016AC */ +qq1 = 3.97917223959155352819e-01, /* 0x3FD97779, 0xCDDADC09 */ +qq2 = 6.50222499887672944485e-02, /* 0x3FB0A54C, 0x5536CEBA */ +qq3 = 5.08130628187576562776e-03, /* 0x3F74D022, 0xC4D36B0F */ +qq4 = 1.32494738004321644526e-04, /* 0x3F215DC9, 0x221C1A10 */ +qq5 = -3.96022827877536812320e-06, /* 0xBED09C43, 0x42A26120 */ +/* + * Coefficients for approximation to erf in [0.84375,1.25] + */ +pa0 = -2.36211856075265944077e-03, /* 0xBF6359B8, 0xBEF77538 */ +pa1 = 4.14856118683748331666e-01, /* 0x3FDA8D00, 0xAD92B34D */ +pa2 = -3.72207876035701323847e-01, /* 0xBFD7D240, 0xFBB8C3F1 */ +pa3 = 3.18346619901161753674e-01, /* 0x3FD45FCA, 0x805120E4 */ +pa4 = -1.10894694282396677476e-01, /* 0xBFBC6398, 0x3D3E28EC */ +pa5 = 3.54783043256182359371e-02, /* 0x3FA22A36, 0x599795EB */ +pa6 = -2.16637559486879084300e-03, /* 0xBF61BF38, 0x0A96073F */ +qa1 = 1.06420880400844228286e-01, /* 0x3FBB3E66, 0x18EEE323 */ +qa2 = 5.40397917702171048937e-01, /* 0x3FE14AF0, 0x92EB6F33 */ +qa3 = 7.18286544141962662868e-02, /* 0x3FB2635C, 0xD99FE9A7 */ +qa4 = 1.26171219808761642112e-01, /* 0x3FC02660, 0xE763351F */ +qa5 = 1.36370839120290507362e-02, /* 0x3F8BEDC2, 0x6B51DD1C */ +qa6 = 1.19844998467991074170e-02, /* 0x3F888B54, 0x5735151D */ +/* + * Coefficients for approximation to erfc in [1.25,1/0.35] + */ +ra0 = -9.86494403484714822705e-03, /* 0xBF843412, 0x600D6435 */ +ra1 = -6.93858572707181764372e-01, /* 0xBFE63416, 0xE4BA7360 */ +ra2 = -1.05586262253232909814e+01, /* 0xC0251E04, 0x41B0E726 */ +ra3 = -6.23753324503260060396e+01, /* 0xC04F300A, 0xE4CBA38D */ +ra4 = -1.62396669462573470355e+02, /* 0xC0644CB1, 0x84282266 */ +ra5 = -1.84605092906711035994e+02, /* 0xC067135C, 0xEBCCABB2 */ +ra6 = -8.12874355063065934246e+01, /* 0xC0545265, 0x57E4D2F2 */ +ra7 = -9.81432934416914548592e+00, /* 0xC023A0EF, 0xC69AC25C */ +sa1 = 1.96512716674392571292e+01, /* 0x4033A6B9, 0xBD707687 */ +sa2 = 1.37657754143519042600e+02, /* 0x4061350C, 0x526AE721 */ +sa3 = 4.34565877475229228821e+02, /* 0x407B290D, 0xD58A1A71 */ +sa4 = 6.45387271733267880336e+02, /* 0x40842B19, 0x21EC2868 */ +sa5 = 4.29008140027567833386e+02, /* 0x407AD021, 0x57700314 */ +sa6 = 1.08635005541779435134e+02, /* 0x405B28A3, 0xEE48AE2C */ +sa7 = 6.57024977031928170135e+00, /* 0x401A47EF, 0x8E484A93 */ +sa8 = -6.04244152148580987438e-02, /* 0xBFAEEFF2, 0xEE749A62 */ +/* + * Coefficients for approximation to erfc in [1/.35,28] + */ +rb0 = -9.86494292470009928597e-03, /* 0xBF843412, 0x39E86F4A */ +rb1 = -7.99283237680523006574e-01, /* 0xBFE993BA, 0x70C285DE */ +rb2 = -1.77579549177547519889e+01, /* 0xC031C209, 0x555F995A */ +rb3 = -1.60636384855821916062e+02, /* 0xC064145D, 0x43C5ED98 */ +rb4 = -6.37566443368389627722e+02, /* 0xC083EC88, 0x1375F228 */ +rb5 = -1.02509513161107724954e+03, /* 0xC0900461, 0x6A2E5992 */ +rb6 = -4.83519191608651397019e+02, /* 0xC07E384E, 0x9BDC383F */ +sb1 = 3.03380607434824582924e+01, /* 0x403E568B, 0x261D5190 */ +sb2 = 3.25792512996573918826e+02, /* 0x40745CAE, 0x221B9F0A */ +sb3 = 1.53672958608443695994e+03, /* 0x409802EB, 0x189D5118 */ +sb4 = 3.19985821950859553908e+03, /* 0x40A8FFB7, 0x688C246A */ +sb5 = 2.55305040643316442583e+03, /* 0x40A3F219, 0xCEDF3BE6 */ +sb6 = 4.74528541206955367215e+02, /* 0x407DA874, 0xE79FE763 */ +sb7 = -2.24409524465858183362e+01; /* 0xC03670E2, 0x42712D62 */ + +#ifdef __STDC__ + double erf(double x) +#else + double erf(x) + double x; +#endif +{ + int hx,ix,i; + double R,S,P,Q,s,y,z,r; + hx = __get_hi_word(x); + ix = hx&0x7fffffff; + if(ix>=0x7ff00000) { /* erf(nan)=nan */ + i = ((unsigned)hx>>31)<<1; + return (double)(1-i)+one/x; /* erf(+-inf)=+-1 */ + } + + if(ix < 0x3feb0000) { /* |x|<0.84375 */ + if(ix < 0x3e300000) { /* |x|<2**-28 */ + if (ix < 0x00800000) + return 0.125*(8.0*x+efx8*x); /*avoid underflow */ + return x + efx*x; + } + z = x*x; + r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4))); + s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5)))); + y = r/s; + return x + x*y; + } + if(ix < 0x3ff40000) { /* 0.84375 <= |x| < 1.25 */ + s = fabs(x)-one; + P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6))))); + Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6))))); + if(hx>=0) return erx + P/Q; else return -erx - P/Q; + } + if (ix >= 0x40180000) { /* inf>|x|>=6 */ + if(hx>=0) return one-tiny; else return tiny-one; + } + x = fabs(x); + s = one/(x*x); + if(ix< 0x4006DB6E) { /* |x| < 1/0.35 */ + R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*( + ra5+s*(ra6+s*ra7)))))); + S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*( + sa5+s*(sa6+s*(sa7+s*sa8))))))); + } else { /* |x| >= 1/0.35 */ + R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*( + rb5+s*rb6))))); + S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*( + sb5+s*(sb6+s*sb7)))))); + } + z = x; + __trunc_lo_word(&z); + r = __ieee754_exp(-z*z-0.5625)*__ieee754_exp((z-x)*(z+x)+R/S); + if(hx>=0) return one-r/x; else return r/x-one; +} + +#ifdef __STDC__ + double erfc(double x) +#else + double erfc(x) + double x; +#endif +{ + int hx,ix; + double R,S,P,Q,s,y,z,r; + hx = __get_hi_word(x); + ix = hx&0x7fffffff; + if(ix>=0x7ff00000) { /* erfc(nan)=nan */ + /* erfc(+-inf)=0,2 */ + return (double)(((unsigned)hx>>31)<<1)+one/x; + } + + if(ix < 0x3feb0000) { /* |x|<0.84375 */ + if(ix < 0x3c700000) /* |x|<2**-56 */ + return one-x; + z = x*x; + r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4))); + s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5)))); + y = r/s; + if(hx < 0x3fd00000) { /* x<1/4 */ + return one-(x+x*y); + } else { + r = x*y; + r += (x-half); + return half - r ; + } + } + if(ix < 0x3ff40000) { /* 0.84375 <= |x| < 1.25 */ + s = fabs(x)-one; + P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6))))); + Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6))))); + if(hx>=0) { + z = one-erx; return z - P/Q; + } else { + z = erx+P/Q; return one+z; + } + } + if (ix < 0x403c0000) { /* |x|<28 */ + x = fabs(x); + s = one/(x*x); + if(ix< 0x4006DB6D) { /* |x| < 1/.35 ~ 2.857143*/ + R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*( + ra5+s*(ra6+s*ra7)))))); + S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*( + sa5+s*(sa6+s*(sa7+s*sa8))))))); + } else { /* |x| >= 1/.35 ~ 2.857143 */ + if(hx<0&&ix>=0x40180000) return two-tiny;/* x < -6 */ + R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*( + rb5+s*rb6))))); + S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*( + sb5+s*(sb6+s*sb7)))))); + } + z = x; + __trunc_lo_word(&z); + r = __ieee754_exp(-z*z-0.5625)* + __ieee754_exp((z-x)*(z+x)+R/S); + if(hx>0) return r/x; else return two-r/x; + } else { + if(hx>0) return tiny*tiny; else return two-tiny; + } +} diff --git a/winsup/mingw/mingwex/math/sf_erf.c b/winsup/mingw/mingwex/math/sf_erf.c new file mode 100644 index 00000000000..20a20fc25c8 --- /dev/null +++ b/winsup/mingw/mingwex/math/sf_erf.c @@ -0,0 +1,259 @@ +/* sf_erf.c -- float version of s_erf.c. + * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. + */ + +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* +#include "fdlibm.h" +*/ +#include +#define __ieee754_expf expf + +#include + +typedef union +{ + float value; + uint32_t word; +} ieee_float_shape_type; + +/* Get a 32 bit int from a float. */ + +static inline int +__get_float_word(float d) +{ + ieee_float_shape_type u; + u.value = d; + return u.word; +} + +/* Set a float from a 32 bit int. */ + +#define SET_FLOAT_WORD(d,i) \ +do { \ + ieee_float_shape_type sf_u; \ + sf_u.word = (i); \ + (d) = sf_u.value; \ +} while (0) + +static inline void __trunc_float_word(float * x) +{ + ieee_float_shape_type u; + u.value = * x; + u.word &= 0xfffff000; +} + +#ifdef __v810__ +#define const +#endif + +#ifdef __STDC__ +static const float +#else +static float +#endif +tiny = 1e-30, +half= 5.0000000000e-01, /* 0x3F000000 */ +one = 1.0000000000e+00, /* 0x3F800000 */ +two = 2.0000000000e+00, /* 0x40000000 */ + /* c = (subfloat)0.84506291151 */ +erx = 8.4506291151e-01, /* 0x3f58560b */ +/* + * Coefficients for approximation to erf on [0,0.84375] + */ +efx = 1.2837916613e-01, /* 0x3e0375d4 */ +efx8= 1.0270333290e+00, /* 0x3f8375d4 */ +pp0 = 1.2837916613e-01, /* 0x3e0375d4 */ +pp1 = -3.2504209876e-01, /* 0xbea66beb */ +pp2 = -2.8481749818e-02, /* 0xbce9528f */ +pp3 = -5.7702702470e-03, /* 0xbbbd1489 */ +pp4 = -2.3763017452e-05, /* 0xb7c756b1 */ +qq1 = 3.9791721106e-01, /* 0x3ecbbbce */ +qq2 = 6.5022252500e-02, /* 0x3d852a63 */ +qq3 = 5.0813062117e-03, /* 0x3ba68116 */ +qq4 = 1.3249473704e-04, /* 0x390aee49 */ +qq5 = -3.9602282413e-06, /* 0xb684e21a */ +/* + * Coefficients for approximation to erf in [0.84375,1.25] + */ +pa0 = -2.3621185683e-03, /* 0xbb1acdc6 */ +pa1 = 4.1485610604e-01, /* 0x3ed46805 */ +pa2 = -3.7220788002e-01, /* 0xbebe9208 */ +pa3 = 3.1834661961e-01, /* 0x3ea2fe54 */ +pa4 = -1.1089469492e-01, /* 0xbde31cc2 */ +pa5 = 3.5478305072e-02, /* 0x3d1151b3 */ +pa6 = -2.1663755178e-03, /* 0xbb0df9c0 */ +qa1 = 1.0642088205e-01, /* 0x3dd9f331 */ +qa2 = 5.4039794207e-01, /* 0x3f0a5785 */ +qa3 = 7.1828655899e-02, /* 0x3d931ae7 */ +qa4 = 1.2617121637e-01, /* 0x3e013307 */ +qa5 = 1.3637083583e-02, /* 0x3c5f6e13 */ +qa6 = 1.1984500103e-02, /* 0x3c445aa3 */ +/* + * Coefficients for approximation to erfc in [1.25,1/0.35] + */ +ra0 = -9.8649440333e-03, /* 0xbc21a093 */ +ra1 = -6.9385856390e-01, /* 0xbf31a0b7 */ +ra2 = -1.0558626175e+01, /* 0xc128f022 */ +ra3 = -6.2375331879e+01, /* 0xc2798057 */ +ra4 = -1.6239666748e+02, /* 0xc322658c */ +ra5 = -1.8460508728e+02, /* 0xc3389ae7 */ +ra6 = -8.1287437439e+01, /* 0xc2a2932b */ +ra7 = -9.8143291473e+00, /* 0xc11d077e */ +sa1 = 1.9651271820e+01, /* 0x419d35ce */ +sa2 = 1.3765776062e+02, /* 0x4309a863 */ +sa3 = 4.3456588745e+02, /* 0x43d9486f */ +sa4 = 6.4538726807e+02, /* 0x442158c9 */ +sa5 = 4.2900814819e+02, /* 0x43d6810b */ +sa6 = 1.0863500214e+02, /* 0x42d9451f */ +sa7 = 6.5702495575e+00, /* 0x40d23f7c */ +sa8 = -6.0424413532e-02, /* 0xbd777f97 */ +/* + * Coefficients for approximation to erfc in [1/.35,28] + */ +rb0 = -9.8649431020e-03, /* 0xbc21a092 */ +rb1 = -7.9928326607e-01, /* 0xbf4c9dd4 */ +rb2 = -1.7757955551e+01, /* 0xc18e104b */ +rb3 = -1.6063638306e+02, /* 0xc320a2ea */ +rb4 = -6.3756646729e+02, /* 0xc41f6441 */ +rb5 = -1.0250950928e+03, /* 0xc480230b */ +rb6 = -4.8351919556e+02, /* 0xc3f1c275 */ +sb1 = 3.0338060379e+01, /* 0x41f2b459 */ +sb2 = 3.2579251099e+02, /* 0x43a2e571 */ +sb3 = 1.5367296143e+03, /* 0x44c01759 */ +sb4 = 3.1998581543e+03, /* 0x4547fdbb */ +sb5 = 2.5530502930e+03, /* 0x451f90ce */ +sb6 = 4.7452853394e+02, /* 0x43ed43a7 */ +sb7 = -2.2440952301e+01; /* 0xc1b38712 */ + +#ifdef __STDC__ + float erff(float x) +#else + float erff(x) + float x; +#endif +{ + int32_t hx,ix,i; + float R,S,P,Q,s,y,z,r; + hx = __get_float_word(x); + ix = hx&0x7fffffff; + if(!(ix<0x7f800000L)) { /* erf(nan)=nan */ + i = ((uint32_t)hx>>31)<<1; + return (float)(1-i)+one/x; /* erf(+-inf)=+-1 */ + } + + if(ix < 0x3f580000) { /* |x|<0.84375 */ + if(ix < 0x31800000) { /* |x|<2**-28 */ + if (ix < 0x04000000) + /*avoid underflow */ + return (float)0.125*((float)8.0*x+efx8*x); + return x + efx*x; + } + z = x*x; + r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4))); + s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5)))); + y = r/s; + return x + x*y; + } + if(ix < 0x3fa00000) { /* 0.84375 <= |x| < 1.25 */ + s = fabsf(x)-one; + P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6))))); + Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6))))); + if(hx>=0) return erx + P/Q; else return -erx - P/Q; + } + if (ix >= 0x40c00000) { /* inf>|x|>=6 */ + if(hx>=0) return one-tiny; else return tiny-one; + } + x = fabsf(x); + s = one/(x*x); + if(ix< 0x4036DB6E) { /* |x| < 1/0.35 */ + R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*( + ra5+s*(ra6+s*ra7)))))); + S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*( + sa5+s*(sa6+s*(sa7+s*sa8))))))); + } else { /* |x| >= 1/0.35 */ + R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*( + rb5+s*rb6))))); + S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*( + sb5+s*(sb6+s*sb7)))))); + } + __trunc_float_word (&z); + r = __ieee754_expf(-z*z-(float)0.5625)*__ieee754_expf((z-x)*(z+x)+R/S); + if(hx>=0) return one-r/x; else return r/x-one; +} + +#ifdef __STDC__ + float erfcf(float x) +#else + float erfcf(x) + float x; +#endif +{ + int32_t hx,ix; + float R,S,P,Q,s,y,z,r; + hx = __get_float_word(x); + ix = hx&0x7fffffff; + if(!(ix<0x7f800000L)) { /* erfc(nan)=nan */ + /* erfc(+-inf)=0,2 */ + return (float)(((uint32_t)hx>>31)<<1)+one/x; + } + + if(ix < 0x3f580000) { /* |x|<0.84375 */ + if(ix < 0x23800000) /* |x|<2**-56 */ + return one-x; + z = x*x; + r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4))); + s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5)))); + y = r/s; + if(hx < 0x3e800000) { /* x<1/4 */ + return one-(x+x*y); + } else { + r = x*y; + r += (x-half); + return half - r ; + } + } + if(ix < 0x3fa00000) { /* 0.84375 <= |x| < 1.25 */ + s = fabsf(x)-one; + P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6))))); + Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6))))); + if(hx>=0) { + z = one-erx; return z - P/Q; + } else { + z = erx+P/Q; return one+z; + } + } + + if (ix < 0x41e00000) { /* |x|<28 */ + x = fabsf(x); + s = one/(x*x); + if(ix< 0x4036DB6D) { /* |x| < 1/.35 ~ 2.857143*/ + R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*( + ra5+s*(ra6+s*ra7)))))); + S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*( + sa5+s*(sa6+s*(sa7+s*sa8))))))); + } else { /* |x| >= 1/.35 ~ 2.857143 */ + if(hx<0&&ix>=0x40c00000) return two-tiny;/* x < -6 */ + R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*( + rb5+s*rb6))))); + S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*( + sb5+s*(sb6+s*sb7)))))); + } + __trunc_float_word (&z); + r = __ieee754_expf(-z*z-(float)0.5625)* + __ieee754_expf((z-x)*(z+x)+R/S); + if(hx>0) return r/x; else return two-r/x; + } else { + if(hx>0) return tiny*tiny; else return two-tiny; + } +} diff --git a/winsup/mingw/mingwex/math/tgamma.c b/winsup/mingw/mingwex/math/tgamma.c new file mode 100644 index 00000000000..c3912a890af --- /dev/null +++ b/winsup/mingw/mingwex/math/tgamma.c @@ -0,0 +1,385 @@ +/* gamma.c + * + * Gamma function + * + * + * + * SYNOPSIS: + * + * double x, y, __tgamma_r(); + * int* sgngam; + * y = __tgamma_r( x, sgngam ); + * + * double x, y, tgamma(); + * y = tgamma( x) + * + * + * + * DESCRIPTION: + * + * Returns gamma function of the argument. The result is + * correctly signed. In the reentrant version the sign (+1 or -1) + * is returned in the variable referenced by sgngam. + * + * Arguments |x| <= 34 are reduced by recurrence and the function + * approximated by a rational function of degree 6/7 in the + * interval (2,3). Large arguments are handled by Stirling's + * formula. Large negative arguments are made positive using + * a reflection formula. + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC -34, 34 10000 1.3e-16 2.5e-17 + * IEEE -170,-33 20000 2.3e-15 3.3e-16 + * IEEE -33, 33 20000 9.4e-16 2.2e-16 + * IEEE 33, 171.6 20000 2.3e-15 3.2e-16 + * + * Error for arguments outside the test range will be larger + * owing to error amplification by the exponential function. + * + */ + +/* +Cephes Math Library Release 2.8: June, 2000 +Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier +*/ + + +/* + * 26-11-2002 Modified for mingw. + * Danny Smith + */ + + +#ifndef __MINGW32__ +#include "mconf.h" +#else +#include "cephes_mconf.h" +#endif + +#ifdef UNK +static const double P[] = { + 1.60119522476751861407E-4, + 1.19135147006586384913E-3, + 1.04213797561761569935E-2, + 4.76367800457137231464E-2, + 2.07448227648435975150E-1, + 4.94214826801497100753E-1, + 9.99999999999999996796E-1 +}; +static const double Q[] = { +-2.31581873324120129819E-5, + 5.39605580493303397842E-4, +-4.45641913851797240494E-3, + 1.18139785222060435552E-2, + 3.58236398605498653373E-2, +-2.34591795718243348568E-1, + 7.14304917030273074085E-2, + 1.00000000000000000320E0 +}; +#define MAXGAM 171.624376956302725 +static const double LOGPI = 1.14472988584940017414; +#endif + +#ifdef DEC +static const unsigned short P[] = { +0035047,0162701,0146301,0005234, +0035634,0023437,0032065,0176530, +0036452,0137157,0047330,0122574, +0037103,0017310,0143041,0017232, +0037524,0066516,0162563,0164605, +0037775,0004671,0146237,0014222, +0040200,0000000,0000000,0000000 +}; +static const unsigned short Q[] = { +0134302,0041724,0020006,0116565, +0035415,0072121,0044251,0025634, +0136222,0003447,0035205,0121114, +0036501,0107552,0154335,0104271, +0037022,0135717,0014776,0171471, +0137560,0034324,0165024,0037021, +0037222,0045046,0047151,0161213, +0040200,0000000,0000000,0000000 +}; +#define MAXGAM 34.84425627277176174 +#endif + +#ifdef IBMPC +static const unsigned short P[] = { +0x2153,0x3998,0xfcb8,0x3f24, +0xbfab,0xe686,0x84e3,0x3f53, +0x14b0,0xe9db,0x57cd,0x3f85, +0x23d3,0x18c4,0x63d9,0x3fa8, +0x7d31,0xdcae,0x8da9,0x3fca, +0xe312,0x3993,0xa137,0x3fdf, +0x0000,0x0000,0x0000,0x3ff0 +}; +static const unsigned short Q[] = { +0xd3af,0x8400,0x487a,0xbef8, +0x2573,0x2915,0xae8a,0x3f41, +0xb44a,0xe750,0x40e4,0xbf72, +0xb117,0x5b1b,0x31ed,0x3f88, +0xde67,0xe33f,0x5779,0x3fa2, +0x87c2,0x9d42,0x071a,0xbfce, +0x3c51,0xc9cd,0x4944,0x3fb2, +0x0000,0x0000,0x0000,0x3ff0 +}; +#define MAXGAM 171.624376956302725 +#endif + +#ifdef MIEEE +static const unsigned short P[] = { +0x3f24,0xfcb8,0x3998,0x2153, +0x3f53,0x84e3,0xe686,0xbfab, +0x3f85,0x57cd,0xe9db,0x14b0, +0x3fa8,0x63d9,0x18c4,0x23d3, +0x3fca,0x8da9,0xdcae,0x7d31, +0x3fdf,0xa137,0x3993,0xe312, +0x3ff0,0x0000,0x0000,0x0000 +}; +static const unsigned short Q[] = { +0xbef8,0x487a,0x8400,0xd3af, +0x3f41,0xae8a,0x2915,0x2573, +0xbf72,0x40e4,0xe750,0xb44a, +0x3f88,0x31ed,0x5b1b,0xb117, +0x3fa2,0x5779,0xe33f,0xde67, +0xbfce,0x071a,0x9d42,0x87c2, +0x3fb2,0x4944,0xc9cd,0x3c51, +0x3ff0,0x0000,0x0000,0x0000 +}; +#define MAXGAM 171.624376956302725 +#endif + +/* Stirling's formula for the gamma function */ +#if UNK +static const double STIR[5] = { + 7.87311395793093628397E-4, +-2.29549961613378126380E-4, +-2.68132617805781232825E-3, + 3.47222221605458667310E-3, + 8.33333333333482257126E-2, +}; +#define MAXSTIR 143.01608 +static const double SQTPI = 2.50662827463100050242E0; +#endif +#if DEC +static const unsigned short STIR[20] = { +0035516,0061622,0144553,0112224, +0135160,0131531,0037460,0165740, +0136057,0134460,0037242,0077270, +0036143,0107070,0156306,0027751, +0037252,0125252,0125252,0146064, +}; +#define MAXSTIR 26.77 +static const unsigned short SQT[4] = { +0040440,0066230,0177661,0034055, +}; +#define SQTPI *(double *)SQT +#endif +#if IBMPC +static const unsigned short STIR[20] = { +0x7293,0x592d,0xcc72,0x3f49, +0x1d7c,0x27e6,0x166b,0xbf2e, +0x4fd7,0x07d4,0xf726,0xbf65, +0xc5fd,0x1b98,0x71c7,0x3f6c, +0x5986,0x5555,0x5555,0x3fb5, +}; +#define MAXSTIR 143.01608 +static const unsigned short SQT[4] = { +0x2706,0x1ff6,0x0d93,0x4004, +}; +#define SQTPI *(double *)SQT +#endif +#if MIEEE +static const unsigned short STIR[20] = { +0x3f49,0xcc72,0x592d,0x7293, +0xbf2e,0x166b,0x27e6,0x1d7c, +0xbf65,0xf726,0x07d4,0x4fd7, +0x3f6c,0x71c7,0x1b98,0xc5fd, +0x3fb5,0x5555,0x5555,0x5986, +}; +#define MAXSTIR 143.01608 +static const unsigned short SQT[4] = { +0x4004,0x0d93,0x1ff6,0x2706, +}; +#define SQTPI *(double *)SQT +#endif + +#ifndef __MINGW32__ +int sgngam = 0; +extern int sgngam; +extern double MAXLOG, MAXNUM, PI; +#ifdef ANSIPROT +extern double pow ( double, double ); +extern double log ( double ); +extern double exp ( double ); +extern double sin ( double ); +extern double polevl ( double, void *, int ); +extern double p1evl ( double, void *, int ); +extern double floor ( double ); +extern double fabs ( double ); +extern int isnan ( double ); +extern int isfinite ( double ); +static double stirf ( double ); +double lgam ( double ); +#else +double pow(), log(), exp(), sin(), polevl(), p1evl(), floor(), fabs(); +int isnan(), isfinite(); +static double stirf(); +double lgam(); +#endif +#ifdef INFINITIES +extern double INFINITY; +#endif +#ifdef NANS +extern double NAN; +#endif +#else /* __MINGW32__ */ +static double stirf ( double ); +#endif + +/* Gamma function computed by Stirling's formula. + * The polynomial STIR is valid for 33 <= x <= 172. + */ +static double stirf(x) +double x; +{ +double y, w, v; + +w = 1.0/x; +w = 1.0 + w * polevl( w, STIR, 4 ); +y = exp(x); +if( x > MAXSTIR ) + { /* Avoid overflow in pow() */ + v = pow( x, 0.5 * x - 0.25 ); + y = v * (v / y); + } +else + { + y = pow( x, x - 0.5 ) / y; + } +y = SQTPI * y * w; +return( y ); +} + + + +double __tgamma_r(double x, int* sgngam) +{ +double p, q, z; +int i; + +*sgngam = 1; +#ifdef NANS +if( isnan(x) ) + return(x); +#endif +#ifdef INFINITIES +#ifdef NANS +if( x == INFINITY ) + return(x); +if( x == -INFINITY ) + return(NAN); +#else +if( !isfinite(x) ) + return(x); +#endif +#endif +q = fabs(x); + +if( q > 33.0 ) + { + if( x < 0.0 ) + { + p = floor(q); + if( p == q ) + { +gsing: + _SET_ERRNO(EDOM); + mtherr( "tgamma", SING ); +#ifdef INFINITIES + return (INFINITY); +#else + return (MAXNUM); +#endif + } + i = p; + if( (i & 1) == 0 ) + *sgngam = -1; + z = q - p; + if( z > 0.5 ) + { + p += 1.0; + z = q - p; + } + z = q * sin( PI * z ); + if( z == 0.0 ) + { + _SET_ERRNO(ERANGE); + mtherr( "tgamma", OVERFLOW ); +#ifdef INFINITIES + return( *sgngam * INFINITY); +#else + return( *sgngam * MAXNUM); +#endif + } + z = fabs(z); + z = PI/(z * stirf(q) ); + } + else + { + z = stirf(x); + } + return( *sgngam * z ); + } + +z = 1.0; +while( x >= 3.0 ) + { + x -= 1.0; + z *= x; + } + +while( x < 0.0 ) + { + if( x > -1.E-9 ) + goto small; + z /= x; + x += 1.0; + } + +while( x < 2.0 ) + { + if( x < 1.e-9 ) + goto small; + z /= x; + x += 1.0; + } + +if( x == 2.0 ) + return(z); + +x -= 2.0; +p = polevl( x, P, 6 ); +q = polevl( x, Q, 7 ); +return( z * p / q ); + +small: +if( x == 0.0 ) + { + goto gsing; + } +else + return( z/((1.0 + 0.5772156649015329 * x) * x) ); +} + +/* This is the C99 version */ + +double tgamma(double x) +{ + int local_sgngam=0; + return (__tgamma_r(x, &local_sgngam)); +} diff --git a/winsup/mingw/mingwex/math/tgammaf.c b/winsup/mingw/mingwex/math/tgammaf.c new file mode 100644 index 00000000000..38e7d75a10a --- /dev/null +++ b/winsup/mingw/mingwex/math/tgammaf.c @@ -0,0 +1,265 @@ +/* gammaf.c + * + * Gamma function + * + * + * + * SYNOPSIS: + * + * float x, y, __tgammaf_r(); + * int* sgngamf; + * y = __tgammaf_r( x, sgngamf ); + * + * float x, y, tgammaf(); + * y = tgammaf( x); + * + * + * DESCRIPTION: + * + * Returns gamma function of the argument. The result is + * correctly signed. In the reentrant version the sign (+1 or -1) + * is returned in the variable referenced by sgngamf. + * + * Arguments between 0 and 10 are reduced by recurrence and the + * function is approximated by a polynomial function covering + * the interval (2,3). Large arguments are handled by Stirling's + * formula. Negative arguments are made positive using + * a reflection formula. + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE 0,-33 100,000 5.7e-7 1.0e-7 + * IEEE -33,0 100,000 6.1e-7 1.2e-7 + * + * + */ + +/* +Cephes Math Library Release 2.7: July, 1998 +Copyright 1984, 1987, 1989, 1992, 1998 by Stephen L. Moshier +*/ + + +/* + * 26-11-2002 Modified for mingw. + * Danny Smith + */ + + +#ifndef __MINGW32__ +#include "mconf.h" +#else +#include "cephes_mconf.h" +#endif + +/* define MAXGAM 34.84425627277176174 */ + +/* Stirling's formula for the gamma function + * gamma(x) = sqrt(2 pi) x^(x-.5) exp(-x) ( 1 + 1/x P(1/x) ) + * .028 < 1/x < .1 + * relative error < 1.9e-11 + */ +static const float STIR[] = { +-2.705194986674176E-003, + 3.473255786154910E-003, + 8.333331788340907E-002, +}; +static const float MAXSTIR = 26.77; +static const float SQTPIF = 2.50662827463100050242; /* sqrt( 2 pi ) */ + +#ifndef __MINGW32__ + +extern float MAXLOGF, MAXNUMF, PIF; + +#ifdef ANSIC +float expf(float); +float logf(float); +float powf( float, float ); +float sinf(float); +float gammaf(float); +float floorf(float); +static float stirf(float); +float polevlf( float, float *, int ); +float p1evlf( float, float *, int ); +#else +float expf(), logf(), powf(), sinf(), floorf(); +float polevlf(), p1evlf(); +static float stirf(); +#endif + +#else /* __MINGW32__ */ +static float stirf(float); +#endif + +/* Gamma function computed by Stirling's formula, + * sqrt(2 pi) x^(x-.5) exp(-x) (1 + 1/x P(1/x)) + * The polynomial STIR is valid for 33 <= x <= 172. + */ +static float stirf( float x ) +{ +float y, w, v; + +w = 1.0/x; +w = 1.0 + w * polevlf( w, STIR, 2 ); +y = expf( -x ); +if( x > MAXSTIR ) + { /* Avoid overflow in pow() */ + v = powf( x, 0.5 * x - 0.25 ); + y *= v; + y *= v; + } +else + { + y = powf( x, x - 0.5 ) * y; + } +y = SQTPIF * y * w; +return( y ); +} + + +/* gamma(x+2), 0 < x < 1 */ +static const float P[] = { + 1.536830450601906E-003, + 5.397581592950993E-003, + 4.130370201859976E-003, + 7.232307985516519E-002, + 8.203960091619193E-002, + 4.117857447645796E-001, + 4.227867745131584E-001, + 9.999999822945073E-001, +}; + +float __tgammaf_r( float x, int* sgngamf) +{ +float p, q, z, nz; +int i, direction, negative; + +#ifdef NANS +if( isnan(x) ) + return(x); +#endif +#ifdef INFINITIES +#ifdef NANS +if( x == INFINITYF ) + return(x); +if( x == -INFINITYF ) + return(NANF); +#else +if( !isfinite(x) ) + return(x); +#endif +#endif + +*sgngamf = 1; +negative = 0; +nz = 0.0; +if( x < 0.0 ) + { + negative = 1; + q = -x; + p = floorf(q); + if( p == q ) + { +gsing: + _SET_ERRNO(EDOM); + mtherr( "tgammaf", SING ); +#ifdef INFINITIES + return (INFINITYF); +#else + return (MAXNUMF); +#endif + } + i = p; + if( (i & 1) == 0 ) + *sgngamf = -1; + nz = q - p; + if( nz > 0.5 ) + { + p += 1.0; + nz = q - p; + } + nz = q * sinf( PIF * nz ); + if( nz == 0.0 ) + { + _SET_ERRNO(ERANGE); + mtherr( "tgamma", OVERFLOW ); +#ifdef INFINITIES + return( *sgngamf * INFINITYF); +#else + return( *sgngamf * MAXNUMF); +#endif + } + if( nz < 0 ) + nz = -nz; + x = q; + } +if( x >= 10.0 ) + { + z = stirf(x); + } +if( x < 2.0 ) + direction = 1; +else + direction = 0; +z = 1.0; +while( x >= 3.0 ) + { + x -= 1.0; + z *= x; + } +/* +while( x < 0.0 ) + { + if( x > -1.E-4 ) + goto small; + z *=x; + x += 1.0; + } +*/ +while( x < 2.0 ) + { + if( x < 1.e-4 ) + goto small; + z *=x; + x += 1.0; + } + +if( direction ) + z = 1.0/z; + +if( x == 2.0 ) + return(z); + +x -= 2.0; +p = z * polevlf( x, P, 7 ); + +gdone: + +if( negative ) + { + p = *sgngamf * PIF/(nz * p ); + } +return(p); + +small: +if( x == 0.0 ) + { + goto gsing; + } +else + { + p = z / ((1.0 + 0.5772156649015329 * x) * x); + goto gdone; + } +} + +/* This is the C99 version */ + +float tgammaf(float x) +{ + int local_sgngamf=0; + return (__tgammaf_r(x, &local_sgngamf)); +} diff --git a/winsup/mingw/mingwex/math/tgammal.c b/winsup/mingw/mingwex/math/tgammal.c new file mode 100644 index 00000000000..682a12e8e4b --- /dev/null +++ b/winsup/mingw/mingwex/math/tgammal.c @@ -0,0 +1,501 @@ +/* gammal.c + * + * Gamma function + * + * + * + * SYNOPSIS: + * + * long double x, y, __tgammal_r(); + * int* sgngaml; + * y = __tgammal_r( x, sgngaml ); + * + * long double x, y, tgammal(); + * y = tgammal( x); * + * + * + * DESCRIPTION: + * + * Returns gamma function of the argument. The result is + * correctly signed. In the reentrant version the sign (+1 or -1) + * is returned in the variable referenced by sgngamf. + * + * Arguments |x| <= 13 are reduced by recurrence and the function + * approximated by a rational function of degree 7/8 in the + * interval (2,3). Large arguments are handled by Stirling's + * formula. Large negative arguments are made positive using + * a reflection formula. + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -40,+40 10000 3.6e-19 7.9e-20 + * IEEE -1755,+1755 10000 4.8e-18 6.5e-19 + * + * Accuracy for large arguments is dominated by error in powl(). + * + */ + +/* +Copyright 1994 by Stephen L. Moshier +*/ + + +/* + * 26-11-2002 Modified for mingw. + * Danny Smith + */ + + +#ifndef __MINGW32__ +#include "mconf.h" +#else +#include "cephes_mconf.h" +#endif + +/* +gamma(x+2) = gamma(x+2) P(x)/Q(x) +0 <= x <= 1 +Relative error +n=7, d=8 +Peak error = 1.83e-20 +Relative error spread = 8.4e-23 +*/ + +#if UNK +static const long double P[8] = { + 4.212760487471622013093E-5L, + 4.542931960608009155600E-4L, + 4.092666828394035500949E-3L, + 2.385363243461108252554E-2L, + 1.113062816019361559013E-1L, + 3.629515436640239168939E-1L, + 8.378004301573126728826E-1L, + 1.000000000000000000009E0L, +}; +static const long double Q[9] = { +-1.397148517476170440917E-5L, + 2.346584059160635244282E-4L, +-1.237799246653152231188E-3L, +-7.955933682494738320586E-4L, + 2.773706565840072979165E-2L, +-4.633887671244534213831E-2L, +-2.243510905670329164562E-1L, + 4.150160950588455434583E-1L, + 9.999999999999999999908E-1L, +}; +#endif +#if IBMPC +static const short P[] = { +0x434a,0x3f22,0x2bda,0xb0b2,0x3ff0, XPD +0xf5aa,0xe82f,0x335b,0xee2e,0x3ff3, XPD +0xbe6c,0x3757,0xc717,0x861b,0x3ff7, XPD +0x7f43,0x5196,0xb166,0xc368,0x3ff9, XPD +0x9549,0x8eb5,0x8c3a,0xe3f4,0x3ffb, XPD +0x8d75,0x23af,0xc8e4,0xb9d4,0x3ffd, XPD +0x29cf,0x19b3,0x16c8,0xd67a,0x3ffe, XPD +0x0000,0x0000,0x0000,0x8000,0x3fff, XPD +}; +static const short Q[] = { +0x5473,0x2de8,0x1268,0xea67,0xbfee, XPD +0x334b,0xc2f0,0xa2dd,0xf60e,0x3ff2, XPD +0xbeed,0x1853,0xa691,0xa23d,0xbff5, XPD +0x296e,0x7cb1,0x5dfd,0xd08f,0xbff4, XPD +0x0417,0x7989,0xd7bc,0xe338,0x3ff9, XPD +0x3295,0x3698,0xd580,0xbdcd,0xbffa, XPD +0x75ef,0x3ab7,0x4ad3,0xe5bc,0xbffc, XPD +0xe458,0x2ec7,0xfd57,0xd47c,0x3ffd, XPD +0x0000,0x0000,0x0000,0x8000,0x3fff, XPD +}; +#endif +#if MIEEE +static const long P[24] = { +0x3ff00000,0xb0b22bda,0x3f22434a, +0x3ff30000,0xee2e335b,0xe82ff5aa, +0x3ff70000,0x861bc717,0x3757be6c, +0x3ff90000,0xc368b166,0x51967f43, +0x3ffb0000,0xe3f48c3a,0x8eb59549, +0x3ffd0000,0xb9d4c8e4,0x23af8d75, +0x3ffe0000,0xd67a16c8,0x19b329cf, +0x3fff0000,0x80000000,0x00000000, +}; +static const long Q[27] = { +0xbfee0000,0xea671268,0x2de85473, +0x3ff20000,0xf60ea2dd,0xc2f0334b, +0xbff50000,0xa23da691,0x1853beed, +0xbff40000,0xd08f5dfd,0x7cb1296e, +0x3ff90000,0xe338d7bc,0x79890417, +0xbffa0000,0xbdcdd580,0x36983295, +0xbffc0000,0xe5bc4ad3,0x3ab775ef, +0x3ffd0000,0xd47cfd57,0x2ec7e458, +0x3fff0000,0x80000000,0x00000000, +}; +#endif +/* +static const long double P[] = { +-3.01525602666895735709e0L, +-3.25157411956062339893e1L, +-2.92929976820724030353e2L, +-1.70730828800510297666e3L, +-7.96667499622741999770e3L, +-2.59780216007146401957e4L, +-5.99650230220855581642e4L, +-7.15743521530849602425e4L +}; +static const long double Q[] = { + 1.00000000000000000000e0L, +-1.67955233807178858919e1L, + 8.85946791747759881659e1L, + 5.69440799097468430177e1L, +-1.98526250512761318471e3L, + 3.31667508019495079814e3L, + 1.60577839621734713377e4L, +-2.97045081369399940529e4L, +-7.15743521530849602412e4L +}; +*/ +#define MAXGAML 1755.455L +/*static const long double LOGPI = 1.14472988584940017414L;*/ + +/* Stirling's formula for the gamma function +gamma(x) = sqrt(2 pi) x^(x-.5) exp(-x) (1 + 1/x P(1/x)) +z(x) = x +13 <= x <= 1024 +Relative error +n=8, d=0 +Peak error = 9.44e-21 +Relative error spread = 8.8e-4 +*/ +#if UNK +static const long double STIR[9] = { + 7.147391378143610789273E-4L, +-2.363848809501759061727E-5L, +-5.950237554056330156018E-4L, + 6.989332260623193171870E-5L, + 7.840334842744753003862E-4L, +-2.294719747873185405699E-4L, +-2.681327161876304418288E-3L, + 3.472222222230075327854E-3L, + 8.333333333333331800504E-2L, +}; +#endif +#if IBMPC +static const short STIR[] = { +0x6ede,0x69f7,0x54e3,0xbb5d,0x3ff4, XPD +0xc395,0x0295,0x4443,0xc64b,0xbfef, XPD +0xba6f,0x7c59,0x5e47,0x9bfb,0xbff4, XPD +0x5704,0x1a39,0xb11d,0x9293,0x3ff1, XPD +0x30b7,0x1a21,0x98b2,0xcd87,0x3ff4, XPD +0xbef3,0x7023,0x6a08,0xf09e,0xbff2, XPD +0x3a1c,0x5ac8,0x3478,0xafb9,0xbff6, XPD +0xc3c9,0x906e,0x38e3,0xe38e,0x3ff6, XPD +0xa1d5,0xaaaa,0xaaaa,0xaaaa,0x3ffb, XPD +}; +#endif +#if MIEEE +static const long STIR[27] = { +0x3ff40000,0xbb5d54e3,0x69f76ede, +0xbfef0000,0xc64b4443,0x0295c395, +0xbff40000,0x9bfb5e47,0x7c59ba6f, +0x3ff10000,0x9293b11d,0x1a395704, +0x3ff40000,0xcd8798b2,0x1a2130b7, +0xbff20000,0xf09e6a08,0x7023bef3, +0xbff60000,0xafb93478,0x5ac83a1c, +0x3ff60000,0xe38e38e3,0x906ec3c9, +0x3ffb0000,0xaaaaaaaa,0xaaaaa1d5, +}; +#endif +#define MAXSTIR 1024.0L +static const long double SQTPI = 2.50662827463100050242E0L; + +/* 1/gamma(x) = z P(z) + * z(x) = 1/x + * 0 < x < 0.03125 + * Peak relative error 4.2e-23 + */ +#if UNK +static const long double S[9] = { +-1.193945051381510095614E-3L, + 7.220599478036909672331E-3L, +-9.622023360406271645744E-3L, +-4.219773360705915470089E-2L, + 1.665386113720805206758E-1L, +-4.200263503403344054473E-2L, +-6.558780715202540684668E-1L, + 5.772156649015328608253E-1L, + 1.000000000000000000000E0L, +}; +#endif +#if IBMPC +static const unsigned short S[] = { +0xbaeb,0xd6d3,0x25e5,0x9c7e,0xbff5, XPD +0xfe9a,0xceb4,0xc74e,0xec9a,0x3ff7, XPD +0x9225,0xdfef,0xb0e9,0x9da5,0xbff8, XPD +0x10b0,0xec17,0x87dc,0xacd7,0xbffa, XPD +0x6b8d,0x7515,0x1905,0xaa89,0x3ffc, XPD +0xf183,0x126b,0xf47d,0xac0a,0xbffa, XPD +0x7bf6,0x57d1,0xa013,0xa7e7,0xbffe, XPD +0xc7a9,0x7db0,0x67e3,0x93c4,0x3ffe, XPD +0x0000,0x0000,0x0000,0x8000,0x3fff, XPD +}; +#endif +#if MIEEE +static const long S[27] = { +0xbff50000,0x9c7e25e5,0xd6d3baeb, +0x3ff70000,0xec9ac74e,0xceb4fe9a, +0xbff80000,0x9da5b0e9,0xdfef9225, +0xbffa0000,0xacd787dc,0xec1710b0, +0x3ffc0000,0xaa891905,0x75156b8d, +0xbffa0000,0xac0af47d,0x126bf183, +0xbffe0000,0xa7e7a013,0x57d17bf6, +0x3ffe0000,0x93c467e3,0x7db0c7a9, +0x3fff0000,0x80000000,0x00000000, +}; +#endif +/* 1/gamma(-x) = z P(z) + * z(x) = 1/x + * 0 < x < 0.03125 + * Peak relative error 5.16e-23 + * Relative error spread = 2.5e-24 + */ +#if UNK +static const long double SN[9] = { + 1.133374167243894382010E-3L, + 7.220837261893170325704E-3L, + 9.621911155035976733706E-3L, +-4.219773343731191721664E-2L, +-1.665386113944413519335E-1L, +-4.200263503402112910504E-2L, + 6.558780715202536547116E-1L, + 5.772156649015328608727E-1L, +-1.000000000000000000000E0L, +}; +#endif +#if IBMPC +static const unsigned short SN[] = { +0x5dd1,0x02de,0xb9f7,0x948d,0x3ff5, XPD +0x989b,0xdd68,0xc5f1,0xec9c,0x3ff7, XPD +0x2ca1,0x18f0,0x386f,0x9da5,0x3ff8, XPD +0x783f,0x41dd,0x87d1,0xacd7,0xbffa, XPD +0x7a5b,0xd76d,0x1905,0xaa89,0xbffc, XPD +0x7f64,0x1234,0xf47d,0xac0a,0xbffa, XPD +0x5e26,0x57d1,0xa013,0xa7e7,0x3ffe, XPD +0xc7aa,0x7db0,0x67e3,0x93c4,0x3ffe, XPD +0x0000,0x0000,0x0000,0x8000,0xbfff, XPD +}; +#endif +#if MIEEE +static const long SN[27] = { +0x3ff50000,0x948db9f7,0x02de5dd1, +0x3ff70000,0xec9cc5f1,0xdd68989b, +0x3ff80000,0x9da5386f,0x18f02ca1, +0xbffa0000,0xacd787d1,0x41dd783f, +0xbffc0000,0xaa891905,0xd76d7a5b, +0xbffa0000,0xac0af47d,0x12347f64, +0x3ffe0000,0xa7e7a013,0x57d15e26, +0x3ffe0000,0x93c467e3,0x7db0c7aa, +0xbfff0000,0x80000000,0x00000000, +}; +#endif + +#ifndef __MINGW32__ +extern long double MAXLOGL, MAXNUML, PIL; +/* #define PIL 3.14159265358979323846L */ +/* #define MAXNUML 1.189731495357231765021263853E4932L */ + +#ifdef ANSIPROT +extern long double fabsl ( long double ); +extern long double lgaml ( long double ); +extern long double logl ( long double ); +extern long double expl ( long double ); +extern long double gammal ( long double ); +extern long double sinl ( long double ); +extern long double floorl ( long double ); +extern long double powl ( long double, long double ); +extern long double polevll ( long double, void *, int ); +extern long double p1evll ( long double, void *, int ); +extern int isnanl ( long double ); +extern int isfinitel ( long double ); +static long double stirf ( long double ); +#else +long double fabsl(), lgaml(), logl(), expl(), gammal(), sinl(); +long double floorl(), powl(), polevll(), p1evll(), isnanl(), isfinitel(); +static long double stirf(); +#endif +#ifdef INFINITIES +extern long double INFINITYL; +#endif +#ifdef NANS +extern long double NANL; +#endif + +#else /* __MINGW32__ */ +static long double stirf ( long double ); +#endif + + +/* Gamma function computed by Stirling's formula. */ + +static long double stirf(x) +long double x; +{ +long double y, w, v; + +w = 1.0L/x; +/* For large x, use rational coefficients from the analytical expansion. */ +if( x > 1024.0L ) + w = (((((6.97281375836585777429E-5L * w + + 7.84039221720066627474E-4L) * w + - 2.29472093621399176955E-4L) * w + - 2.68132716049382716049E-3L) * w + + 3.47222222222222222222E-3L) * w + + 8.33333333333333333333E-2L) * w + + 1.0L; +else + w = 1.0L + w * polevll( w, STIR, 8 ); +y = expl(x); +if( x > MAXSTIR ) + { /* Avoid overflow in pow() */ + v = powl( x, 0.5L * x - 0.25L ); + y = v * (v / y); + } +else + { + y = powl( x, x - 0.5L ) / y; + } +y = SQTPI * y * w; +return( y ); +} + + +long double __tgammal_r(long double x, int* sgngaml) +{ +long double p, q, z; +int i; + +*sgngaml = 1; +#ifdef NANS +if( isnanl(x) ) + return(NANL); +#endif +#ifdef INFINITIES +#ifdef NANS +if( x == INFINITYL ) + return(x); +if( x == -INFINITYL ) + return(NANL); +#else +if( !isfinite(x) ) + return(x); +#endif +#endif +q = fabsl(x); + +if( q > 13.0L ) + { + if( q > MAXGAML ) + goto goverf; + if( x < 0.0L ) + { + p = floorl(q); + if( p == q ) + { +gsing: + _SET_ERRNO(EDOM); + mtherr( "tgammal", SING ); +#ifdef INFINITIES + return (INFINITYL); +#else + return( *sgngaml * MAXNUML); +#endif + } + i = p; + if( (i & 1) == 0 ) + *sgngaml = -1; + z = q - p; + if( z > 0.5L ) + { + p += 1.0L; + z = q - p; + } + z = q * sinl( PIL * z ); + z = fabsl(z) * stirf(q); + if( z <= PIL/MAXNUML ) + { +goverf: + _SET_ERRNO(ERANGE); + mtherr( "tgammal", OVERFLOW ); +#ifdef INFINITIES + return( *sgngaml * INFINITYL); +#else + return( *sgngaml * MAXNUML); +#endif + } + z = PIL/z; + } + else + { + z = stirf(x); + } + return( *sgngaml * z ); + } + +z = 1.0L; +while( x >= 3.0L ) + { + x -= 1.0L; + z *= x; + } + +while( x < -0.03125L ) + { + z /= x; + x += 1.0L; + } + +if( x <= 0.03125L ) + goto small; + +while( x < 2.0L ) + { + z /= x; + x += 1.0L; + } + +if( x == 2.0L ) + return(z); + +x -= 2.0L; +p = polevll( x, P, 7 ); +q = polevll( x, Q, 8 ); +return( z * p / q ); + +small: +if( x == 0.0L ) + { + goto gsing; + } +else + { + if( x < 0.0L ) + { + x = -x; + q = z / (x * polevll( x, SN, 8 )); + } + else + q = z / (x * polevll( x, S, 8 )); + } +return q; +} + + +/* This is the C99 version. */ + +long double tgammal(long double x) +{ + int local_sgngaml=0; + return (__tgammal_r(x, &local_sgngaml)); +} + diff --git a/winsup/mingw/mingwex/strtold.c b/winsup/mingw/mingwex/strtold.c new file mode 100644 index 00000000000..81db8a6f238 --- /dev/null +++ b/winsup/mingw/mingwex/strtold.c @@ -0,0 +1,384 @@ +/* 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 + */ + + +#include "math/cephes_emath.h" + +#if NE == 10 + +/* 1.0E0 */ +static const unsigned short __eone[NE] = + {0x0000, 0x0000, 0x0000, 0x0000, + 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3fff,}; +#else +static const unsigned short __eone[NE] = { +0, 0000000,0000000,0000000,0100000,0x3fff,}; +#endif + +#if NE == 10 +static const unsigned short __etens[NTEN + 1][NE] = +{ + {0x6576, 0x4a92, 0x804a, 0x153f, + 0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,}, /* 10**4096 */ + {0x6a32, 0xce52, 0x329a, 0x28ce, + 0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,}, /* 10**2048 */ + {0x526c, 0x50ce, 0xf18b, 0x3d28, + 0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,}, + {0x9c66, 0x58f8, 0xbc50, 0x5c54, + 0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,}, + {0x851e, 0xeab7, 0x98fe, 0x901b, + 0xddbb, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,}, + {0x0235, 0x0137, 0x36b1, 0x336c, + 0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,}, + {0x50f8, 0x25fb, 0xc76b, 0x6b71, + 0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,}, + {0x0000, 0x0000, 0x0000, 0x0000, + 0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,}, + {0x0000, 0x0000, 0x0000, 0x0000, + 0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,}, + {0x0000, 0x0000, 0x0000, 0x0000, + 0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,}, + {0x0000, 0x0000, 0x0000, 0x0000, + 0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,}, + {0x0000, 0x0000, 0x0000, 0x0000, + 0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,}, + {0x0000, 0x0000, 0x0000, 0x0000, + 0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,}, /* 10**1 */ +}; +#else +static const unsigned short __etens[NTEN+1][NE] = { +{0xc94c,0x979a,0x8a20,0x5202,0xc460,0x7525,},/* 10**4096 */ +{0xa74d,0x5de4,0xc53d,0x3b5d,0x9e8b,0x5a92,},/* 10**2048 */ +{0x650d,0x0c17,0x8175,0x7586,0xc976,0x4d48,}, +{0xcc65,0x91c6,0xa60e,0xa0ae,0xe319,0x46a3,}, +{0xddbc,0xde8d,0x9df9,0xebfb,0xaa7e,0x4351,}, +{0xc66f,0x8cdf,0x80e9,0x47c9,0x93ba,0x41a8,}, +{0x3cbf,0xa6d5,0xffcf,0x1f49,0xc278,0x40d3,}, +{0xf020,0xb59d,0x2b70,0xada8,0x9dc5,0x4069,}, +{0x0000,0x0000,0x0400,0xc9bf,0x8e1b,0x4034,}, +{0x0000,0x0000,0x0000,0x2000,0xbebc,0x4019,}, +{0x0000,0x0000,0x0000,0x0000,0x9c40,0x400c,}, +{0x0000,0x0000,0x0000,0x0000,0xc800,0x4005,}, +{0x0000,0x0000,0x0000,0x0000,0xa000,0x4002,}, /* 10**1 */ +}; +#endif + +int __asctoe64(const char * __restrict__ ss, short unsigned int * __restrict__ y) +{ +unsigned short yy[NI], xt[NI], tt[NI]; +int esign, decflg, sgnflg, nexp, exp, prec, lost; +int k, trail, c; +long lexp; +unsigned short nsign; +const unsigned short *p; +char *sp, *lstr; +char *s; + +char dec_sym = *(localeconv ()->decimal_point); + +int lenldstr = 0; + +/* Copy the input string. */ +c = strlen (ss) + 2; +lstr = (char *) alloca (c); +s = (char *) ss; +while( isspace ((int)(unsigned char)*s)) /* skip leading spaces */ + { + ++s; + ++lenldstr; + } +sp = lstr; +for( k=0; k= 0) && (k <= 9) ) + { +/* Ignore leading zeros */ + if( (prec == 0) && (decflg == 0) && (k == 0) ) + goto donchr; +/* Identify and strip trailing zeros after the decimal point. */ + if( (trail == 0) && (decflg != 0) ) + { + sp = s; + while( (*sp >= '0') && (*sp <= '9') ) + ++sp; + --sp; + while( *sp == '0' ) + *sp-- = 'z'; + trail = 1; + if( *s == 'z' ) + goto donchr; + } +/* If enough digits were given to more than fill up the yy register, + * continuing until overflow into the high guard word yy[2] + * guarantees that there will be a roundoff bit at the top + * of the low guard word after normalization. + */ + if( yy[2] == 0 ) + { + if( decflg ) + nexp += 1; /* count digits after decimal point */ + __eshup1( yy ); /* multiply current number by 10 */ + __emovz( yy, xt ); + __eshup1( xt ); + __eshup1( xt ); + __eaddm( xt, yy ); + __ecleaz( xt ); + xt[NI-2] = (unsigned short )k; + __eaddm( xt, yy ); + } + else + { + /* Mark any lost non-zero digit. */ + lost |= k; + /* Count lost digits before the decimal point. */ + if (decflg == 0) + nexp -= 1; + } + prec += 1; + goto donchr; + } +if (*s == dec_sym) + { + if( decflg ) + goto daldone; + ++decflg; + } +else + switch( *s ) + { + case 'z': + break; + case 'E': + case 'e': + goto expnt; + case '-': + nsign = 0xffff; + if( sgnflg ) + goto daldone; + ++sgnflg; + break; + case '+': + if( sgnflg ) + goto daldone; + ++sgnflg; + break; + case 'i': + case 'I': + { + s++; + if (*s != 'n' && *s != 'N') + goto zero; + s++; + if (*s != 'f' && *s != 'F') + goto zero; + s++; + if ((*s == 'i' || *s == 'I') && (s[1] == 'n' || s[1] == 'N') + && (s[2] == 'i' || s[2] == 'I') + && (s[3] == 't' || s[3] == 'T') + && (s[4] == 'y' || s[4] == 'Y')) + s += 5; + goto infinite; + } + case 'n': + case 'N': + { + s++; + if (*s != 'a' && *s != 'A') + goto zero; + s++; + if (*s != 'n' && *s != 'N') + goto zero; + s++; + __enan_NI16( yy ); + goto aexit; + } + default: + goto daldone; + } +donchr: +++s; +goto nxtcom; + +/* Exponent interpretation */ +expnt: + +esign = 1; +exp = 0; +++s; +/* check for + or - */ +if( *s == '-' ) + { + esign = -1; + ++s; + } +if( *s == '+' ) + ++s; +while( (*s >= '0') && (*s <= '9') && exp < 4978) + { + exp *= 10; + exp += *s++ - '0'; + } +if( esign < 0 ) + exp = -exp; +if( exp > 4932 ) + { + errno = ERANGE; +infinite: + __ecleaz(yy); + yy[E] = 0x7fff; /* infinity */ + goto aexit; + } +if( exp < -4977 ) + { + errno = ERANGE; +zero: + __ecleaz(yy); + goto aexit; + } + +daldone: +nexp = exp - nexp; +/* Pad trailing zeros to minimize power of 10, per IEEE spec. */ +while( (nexp > 0) && (yy[2] == 0) ) + { + __emovz( yy, xt ); + __eshup1( xt ); + __eshup1( xt ); + __eaddm( yy, xt ); + __eshup1( xt ); + if( xt[2] != 0 ) + break; + nexp -= 1; + __emovz( xt, yy ); + } +if( (k = __enormlz(yy)) > NBITS ) + { + __ecleaz(yy); + goto aexit; + } +lexp = (EXONE - 1 + NBITS) - k; +__emdnorm( yy, lost, 0, lexp, 64, NBITS ); +/* convert to external format */ + + +/* Multiply by 10**nexp. If precision is 64 bits, + * the maximum relative error incurred in forming 10**n + * for 0 <= n <= 324 is 8.2e-20, at 10**180. + * For 0 <= n <= 999, the peak relative error is 1.4e-19 at 10**947. + * For 0 >= n >= -999, it is -1.55e-19 at 10**-435. + */ +lexp = yy[E]; +if( nexp == 0 ) + { + k = 0; + goto expdon; + } +esign = 1; +if( nexp < 0 ) + { + nexp = -nexp; + esign = -1; + if( nexp > 4096 ) + { /* Punt. Can't handle this without 2 divides. */ + __emovi( __etens[0], tt ); + lexp -= tt[E]; + k = __edivm( tt, yy ); + lexp += EXONE; + nexp -= 4096; + } + } +p = &__etens[NTEN][0]; +__emov( __eone, xt ); +exp = 1; +do + { + if( exp & nexp ) + __emul( p, xt, xt ); + p -= NE; + exp = exp + exp; + } +while( exp <= MAXP ); + +__emovi( xt, tt ); +if( esign < 0 ) + { + lexp -= tt[E]; + k = __edivm( tt, yy ); + lexp += EXONE; + } +else + { + lexp += tt[E]; + k = __emulm( tt, yy ); + lexp -= EXONE - 1; + } + +expdon: + +/* Round and convert directly to the destination type */ + +__emdnorm( yy, k, 0, lexp, 64, 64 ); + +aexit: + +yy[0] = nsign; +__toe64( yy, y ); +return (lenldstr + s - lstr); +} + + +long double strtold (const char * __restrict__ s, char ** __restrict__ se) +{ + int lenldstr; + union + { + unsigned short int us[6]; + long double ld; + } xx = {{0}}; + + lenldstr = __asctoe64( s, xx.us); + if (se) + *se = (char*)s + lenldstr; + return xx.ld; +} diff --git a/winsup/mingw/mingwex/wcstold.c b/winsup/mingw/mingwex/wcstold.c new file mode 100644 index 00000000000..85298807c9a --- /dev/null +++ b/winsup/mingw/mingwex/wcstold.c @@ -0,0 +1,76 @@ +/* Wide char wrapper for strtold + * Revision history: + * 6 Nov 2002 Initial version. + * + * Contributor: Danny Smith + */ + + /* This routine has been placed in the public domain.*/ + +#define WIN32_LEAN_AND_MEAN +#include +#include +#include +#include +#include + +extern int __asctoe64(const char * __restrict__ ss, + short unsigned int * __restrict__ y); + + +static __inline__ unsigned int get_codepage (void) +{ + char* cp; + + /* + locale :: "lang[_country[.code_page]]" + | ".code_page" + */ + if ((cp = strchr(setlocale(LC_CTYPE, NULL), '.'))) + return atoi( cp + 1); + else + return 0; +} + +long double wcstold (const wchar_t * __restrict__ wcs, wchar_t ** __restrict__ wcse) +{ + char * cs; + int i; + int lenldstr; + union + { + unsigned short int us[6]; + long double ld; + } xx; + + unsigned int cp = get_codepage (); + + /* Allocate enough room for (possibly) mb chars */ + cs = (char *) malloc ((wcslen(wcs)+1) * MB_CUR_MAX); + + if (cp == 0) /* C locale */ + { + for (i = 0; (wcs[i] != 0) && wcs[i] <= 255; i++) + cs[i] = (char) wcs[i]; + cs[i] = '\0'; + } + else + { + int nbytes = -1; + int mb_len = 0; + /* loop through till we hit null or invalid character */ + for (i = 0; (wcs[i] != 0) && (nbytes != 0); i++) + { + nbytes = WideCharToMultiByte(cp, WC_COMPOSITECHECK | WC_SEPCHARS, + wcs + i, 1, cs + mb_len, MB_CUR_MAX, + NULL, NULL); + mb_len += nbytes; + } + cs[mb_len] = '\0'; + } + lenldstr = __asctoe64( cs, xx.us); + free (cs); + if (wcse) + *wcse = (wchar_t*) wcs + lenldstr; + return xx.ld; +} -- cgit v1.2.1