diff options
author | Danny Smith <dannysmith@users.sourceforge.net> | 2002-07-29 07:00:10 +0400 |
---|---|---|
committer | Danny Smith <dannysmith@users.sourceforge.net> | 2002-07-29 07:00:10 +0400 |
commit | b8cdc234c67d883b7e9d6e4caf4ba15972ac7788 (patch) | |
tree | e45d96d2d976020ad018e6e919bbeae2500e884a /winsup/mingw/mingwex | |
parent | 840e611264ac6ad7d41bc0b275945f50067cda50 (diff) |
Add incomplet long double math support to libmingwex.a
Diffstat (limited to 'winsup/mingw/mingwex')
138 files changed, 4354 insertions, 219 deletions
diff --git a/winsup/mingw/mingwex/Makefile.in b/winsup/mingw/mingwex/Makefile.in index 59f8e40e4..4c346d346 100644 --- a/winsup/mingw/mingwex/Makefile.in +++ b/winsup/mingw/mingwex/Makefile.in @@ -3,8 +3,8 @@ # # This makefile requires GNU make. -VPATH = @srcdir@ srcdir = @srcdir@ +VPATH = $(srcdir):$(srcdir)/math objdir = . target_alias = @target_alias@ @@ -26,92 +26,43 @@ INSTALL_DATA = @INSTALL_DATA@ INSTALL_PROGRAM = @INSTALL_PROGRAM@ mkinstalldirs = $(SHELL) $(srcdir)/../mkinstalldirs DISTFILES = Makefile.in configure configure.in \ - mingw-fseek.c \ - _Exit.c \ - atoll.c \ - copysignl.S \ - dirent.c \ - fdim.c \ - fdimf.c \ - fdiml.c \ - feclearexcept.c \ - fegetenv.c \ - fegetexceptflag.c \ - fegetround.c \ - feholdexcept.c \ - feraiseexcept.c \ - fesetenv.c \ - fesetround.c \ - fetestexcept.c \ - fesetexceptflag.c \ - feupdateenv.c \ - fma.S \ - fmaf.S \ - fmal.c \ - fmax.c \ - fmaxf.c \ - fmaxl.c \ - fmin.c \ - fminf.c \ - fminl.c \ - fp_consts.c \ - fpclassify.c \ - fpclassifyf.c \ - fpclassifyl.c \ - fucom.c \ - fwide.c \ - imaxabs.c \ - imaxdiv.c \ - isnan.c \ - isnanf.c \ - isnanl.c \ - lltoa.c \ - lltow.c \ - log2.c \ - log2f.c \ - log2l.c \ - math_stubs.c \ - mbsinit.c \ - rint.c \ - rintf.c \ - rintl.c \ - round.c \ - roundf.c \ - roundl.c \ - signbit.c \ - signbitf.c \ - signbitl.c \ - sitest.c \ - snprintf.c \ - snwprintf.c \ - strtof.c \ - strtoimax.c \ - strtoumax.c \ - testwmem.c \ - trunc.c \ - truncf.c \ - truncl.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 \ + _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 \ + 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 +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 \ + 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 \ + fmax.c fmaxf.c fmaxl.c fmin.c fminf.c fminl.c fmodf.c \ + 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 \ + 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 \ + lrintl.c lround.c lroundf.c lroundl.c modff.c modfl.c \ + nearbyint.S nearbyintf.S nearbyintl.S nextafterf.c powf.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 CC = @CC@ # FIXME: Which is it, CC or CC_FOR_TARGET? CC_FOR_TARGET = $(CC) AS_FOR_TARGET = $(AS) CFLAGS = @CFLAGS@ -Wall -CXXFLAGS = @CXXFLAGS@ +CXXFLAGS = @CXXFLAGS@ OPTFLAGS= -fomit-frame-pointer # compiling with Cygwin? @@ -150,22 +101,26 @@ STDLIB_STUB_OBJS = \ STDIO_STUB_OBJS = \ snprintf.o vsnprintf.o snwprintf.o vsnwprintf.o 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 \ + 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 \ + fmax.o fmaxf.o fmaxl.o fmin.o fminf.o fminl.o fmodf.o \ + fmodl.o fp_consts.o fp_constsf.o fp_constsl.o \ fpclassify.o fpclassifyf.o fpclassifyl.o \ - fucom.o \ - round.o roundf.o roundl.o \ - rint.o rintf.o rintl.o \ - signbit.o signbitf.o signbitl.o \ - trunc.o truncf.o truncl.o \ - isnan.o isnanf.o isnanl.o \ - fp_consts.o \ - fdim.o fdimf.o fdiml.o \ - fmax.o fmaxf.o fmaxl.o \ - fmin.o fminf.o fminl.o \ - fma.o fmaf.o fmal.o \ - log2.o log2f.o log2l.o \ - copysignl.o -MATH_STUB_OBJS = \ - math_stubs.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 \ + 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 \ + lrintl.o lround.o lroundf.o lroundl.o modff.o modfl.o \ + nearbyint.o nearbyintf.o nearbyintl.o nextafterf.o powf.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 FENV_OBJS = fesetround.o fegetround.o \ fegetenv.o fesetenv.o feupdateenv.o \ feclearexcept.o feholdexcept.o fegetexceptflag.o \ @@ -176,7 +131,7 @@ REPLACE_OBJS = \ mingw-fseek.o LIB_OBJS = $(Q8_OBJS) $(STDLIB_STUB_OBJS) $(STDIO_STUB_OBJS) \ - $(MATH_OBJS) $(MATH_STUB_OBJS) $(FENV_OBJS) $(POSIX_OBJS) \ + $(MATH_OBJS) $(FENV_OBJS) $(POSIX_OBJS) \ $(REPLACE_OBJS) LIBS = $(LIBMINGWEX_A) @@ -185,6 +140,7 @@ DLLS = all: $(LIBMINGWEX_A) $(LIBMINGWEX_A): $(LIB_OBJS) + rm -f $(LIBMINGWEX_A) $(AR) $(ARFLAGS) $@ $(LIB_OBJS) $(RANLIB) $@ @@ -234,3 +190,9 @@ dist: @for i in $(DISTFILES); do\ cp -p $(srcdir)/$$i $(distdir)/mingwex/$$i ; \ done + mkdir $(distdir)/mingwex/math + chmod 755 $(distdir)//mingwex/math + @for i in $(MATH_DISTFILES); do\ + cp -p $(srcdir)/math/$$i $(distdir)/mingwex/math/$$i ; \ + done + diff --git a/winsup/mingw/mingwex/fmal.c b/winsup/mingw/mingwex/fmal.c deleted file mode 100644 index b827875bf..000000000 --- a/winsup/mingw/mingwex/fmal.c +++ /dev/null @@ -1,4 +0,0 @@ -long double -fmal ( long double _x, long double _y, long double _z){ -return ((_x * _y) + _z); -} diff --git a/winsup/mingw/mingwex/fp_consts.c b/winsup/mingw/mingwex/fp_consts.c deleted file mode 100644 index 9293eeda4..000000000 --- a/winsup/mingw/mingwex/fp_consts.c +++ /dev/null @@ -1,81 +0,0 @@ -/* Floating point consts needed by STL class mumeric_limits<float> - and numeric_limits<double>. Also used as return values by nan, inf */ - -#include <math.h> -/* -According to IEEE 754 a QNaN has exponent bits of all 1 values and -initial significand bit of 1. A SNaN has has an exponent of all 1 -values and initial significand bit of 0 (with one or more other -significand bits of 1). An Inf has significand of 0 and -exponent of all 1 values. A denormal value has all exponent bits of 0. - -The following does _not_ follow those rules, but uses values -equal to those exported from MS C++ runtime lib, msvcprt.dll -for float and double. MSVC however, does not have long doubles. -*/ - - -#define __DOUBLE_INF_REP { 0, 0, 0, 0x7ff0 } -#define __DOUBLE_QNAN_REP { 0, 0, 0, 0xfff8 } /* { 0, 0, 0, 0x7ff8 } */ -#define __DOUBLE_SNAN_REP { 0, 0, 0, 0xfff0 } /* { 1, 0, 0, 0x7ff0 } */ -#define __DOUBLE_DENORM_REP {1, 0, 0, 0} -#define D_NAN_MASK 0x7ff0000000000000LL /* this will mask NaN's and Inf's */ - -#define __FLOAT_INF_REP { 0, 0x7f80 } -#define __FLOAT_QNAN_REP { 0, 0xffc0 } /* { 0, 0x7fc0 } */ -#define __FLOAT_SNAN_REP { 0, 0xff80 } /* { 1, 0x7f80 } */ -#define __FLOAT_DENORM_REP {1,0} -#define F_NAN_MASK 0x7f800000 - -/* This assumes no implicit (hidden) bit in extended mode */ -#define __LONG_DOUBLE_INF_REP { 0, 0, 0, 0x8000, 0x7fff } -#define __LONG_DOUBLE_QNAN_REP { 0, 0, 0, 0xc000, 0xffff } -#define __LONG_DOUBLE_SNAN_REP { 0, 0, 0, 0x8000, 0xffff } -#define __LONG_DOUBLE_DENORM_REP {1, 0, 0, 0, 0} - -union _ieee_rep -{ - unsigned short rep[5]; - float float_val; - double double_val; - long double ldouble_val; -} ; - -const union _ieee_rep __QNAN = { __DOUBLE_QNAN_REP }; -/* -const union _ieee_rep __SNAN = { __DOUBLE_SNAN_REP }; -const union _ieee_rep __INF = { __DOUBLE_INF_REP }; -const union _ieee_rep __DENORM = { __DOUBLE_DENORM_REP }; -*/ -/* ISO C99 */ - -#undef nan -/* FIXME */ -double nan (const char * tagp __attribute__((unused)) ) - { return __QNAN.double_val; } - - -const union _ieee_rep __QNANF = { __FLOAT_QNAN_REP }; -/* -const union _ieee_rep __SNANF = { __FLOAT_SNAN_REP }; -const union _ieee_rep __INFF = { __FLOAT_INF_REP }; -const union _ieee_rep __DENORMF = { __FLOAT_DENORM_REP }; -*/ - -#undef nanf -/* FIXME */ -float nanf(const char * tagp __attribute__((unused)) ) - { return __QNANF.float_val;} - - -const union _ieee_rep __QNANL = { __LONG_DOUBLE_QNAN_REP }; -/* -const union _ieee_rep __SNANL = { __LONG_DOUBLE_SNAN_REP }; -const union _ieee_rep __INFL = { __LONG_DOUBLE_INF_REP }; -const union _ieee_rep __DENORML = { __LONG_DOUBLE_DENORM_REP }; -*/ - -#undef nanl -/* FIXME */ -long double nanl (const char * tagp __attribute__((unused)) ) - { return __QNANL.ldouble_val;} diff --git a/winsup/mingw/mingwex/log2.c b/winsup/mingw/mingwex/log2.c deleted file mode 100644 index 4f14f2614..000000000 --- a/winsup/mingw/mingwex/log2.c +++ /dev/null @@ -1,8 +0,0 @@ -#include <math.h> -double -log2 (double _x) -{ - double retval; - __asm__ ("fyl2x;" : "=t" (retval) : "0" (_x), "u" (1.0L) : "st(1)"); - return retval; -} diff --git a/winsup/mingw/mingwex/log2f.c b/winsup/mingw/mingwex/log2f.c deleted file mode 100644 index 576396c59..000000000 --- a/winsup/mingw/mingwex/log2f.c +++ /dev/null @@ -1,8 +0,0 @@ -#include <math.h> -float -log2f (float _x) -{ - float retval; - __asm__ ("fyl2x;" : "=t" (retval) : "0" (_x), "u" (1.0L) : "st(1)"); - return retval; -} diff --git a/winsup/mingw/mingwex/log2l.c b/winsup/mingw/mingwex/log2l.c deleted file mode 100644 index c67d7701d..000000000 --- a/winsup/mingw/mingwex/log2l.c +++ /dev/null @@ -1,8 +0,0 @@ -#include <math.h> -long double -log2l (long double _x) -{ - long double retval; - __asm__ ("fyl2x;" : "=t" (retval) : "0" (_x), "u" (1.0L) : "st(1)"); - return retval; -} diff --git a/winsup/mingw/mingwex/math/acosf.c b/winsup/mingw/mingwex/math/acosf.c new file mode 100644 index 000000000..364f6a90c --- /dev/null +++ b/winsup/mingw/mingwex/math/acosf.c @@ -0,0 +1,23 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + +#include <math.h> + +float +acosf (float x) +{ + float res; + + /* acosl = atanl (sqrtl(1 - x^2) / x) */ + asm ( "fld %%st\n\t" + "fmul %%st(0)\n\t" /* x^2 */ + "fld1\n\t" + "fsubp\n\t" /* 1 - x^2 */ + "fsqrt\n\t" /* sqrtl (1 - x^2) */ + "fxch %%st(1)\n\t" + "fpatan" + : "=t" (res) : "0" (x) : "st(1)"); + return res; +} diff --git a/winsup/mingw/mingwex/math/acosl.c b/winsup/mingw/mingwex/math/acosl.c new file mode 100644 index 000000000..f98d2cdc1 --- /dev/null +++ b/winsup/mingw/mingwex/math/acosl.c @@ -0,0 +1,25 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * + * Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>. + */ + +#include <math.h> + +long double +acosl (long double x) +{ + long double res; + + /* acosl = atanl (sqrtl(1 - x^2) / x) */ + asm ( "fld %%st\n\t" + "fmul %%st(0)\n\t" /* x^2 */ + "fld1\n\t" + "fsubp\n\t" /* 1 - x^2 */ + "fsqrt\n\t" /* sqrtl (1 - x^2) */ + "fxch %%st(1)\n\t" + "fpatan" + : "=t" (res) : "0" (x) : "st(1)"); + return res; +} diff --git a/winsup/mingw/mingwex/math/asinf.c b/winsup/mingw/mingwex/math/asinf.c new file mode 100644 index 000000000..e79429ec8 --- /dev/null +++ b/winsup/mingw/mingwex/math/asinf.c @@ -0,0 +1,20 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + +/* asin = atan (x / sqrt(1 - x^2)) */ + +float asinf (float x) +{ + float res; + + asm ( "fld %%st\n\t" + "fmul %%st(0)\n\t" /* x^2 */ + "fld1\n\t" + "fsubp\n\t" /* 1 - x^2 */ + "fsqrt\n\t" /* sqrt (1 - x^2) */ + "fpatan" + : "=t" (res) : "0" (x) : "st(1)"); + return res; +} diff --git a/winsup/mingw/mingwex/math/asinl.c b/winsup/mingw/mingwex/math/asinl.c new file mode 100644 index 000000000..a2ac32b39 --- /dev/null +++ b/winsup/mingw/mingwex/math/asinl.c @@ -0,0 +1,21 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * Adapted for long double type by Danny Smith <dannysmith@users.sourceforge.net>. + */ + +/* asin = atan (x / sqrt(1 - x^2)) */ + +long double asinl (long double x) +{ + long double res; + + asm ( "fld %%st\n\t" + "fmul %%st(0)\n\t" /* x^2 */ + "fld1\n\t" + "fsubp\n\t" /* 1 - x^2 */ + "fsqrt\n\t" /* sqrt (1 - x^2) */ + "fpatan" + : "=t" (res) : "0" (x) : "st(1)"); + return res; +} diff --git a/winsup/mingw/mingwex/math/atan2f.c b/winsup/mingw/mingwex/math/atan2f.c new file mode 100644 index 000000000..52ec6f672 --- /dev/null +++ b/winsup/mingw/mingwex/math/atan2f.c @@ -0,0 +1,15 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * + */ + +#include <math.h> + +float +atan2f (float y, float x) +{ + float res; + asm ("fpatan" : "=t" (res) : "u" (y), "0" (x) : "st(1)"); + return res; +} diff --git a/winsup/mingw/mingwex/math/atan2l.c b/winsup/mingw/mingwex/math/atan2l.c new file mode 100644 index 000000000..efd62c1ec --- /dev/null +++ b/winsup/mingw/mingwex/math/atan2l.c @@ -0,0 +1,16 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * + * Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>. + */ + +#include <math.h> + +long double +atan2l (long double y, long double x) +{ + long double res; + asm ("fpatan" : "=t" (res) : "u" (y), "0" (x) : "st(1)"); + return res; +} diff --git a/winsup/mingw/mingwex/math/atanf.c b/winsup/mingw/mingwex/math/atanf.c new file mode 100644 index 000000000..ae70d5daa --- /dev/null +++ b/winsup/mingw/mingwex/math/atanf.c @@ -0,0 +1,17 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * + */ + +#include <math.h> + +float +atanf (float x) +{ + float res; + + asm ("fld1\n\t" + "fpatan" : "=t" (res) : "0" (x)); + return res; +} diff --git a/winsup/mingw/mingwex/math/atanl.c b/winsup/mingw/mingwex/math/atanl.c new file mode 100644 index 000000000..5de06d35b --- /dev/null +++ b/winsup/mingw/mingwex/math/atanl.c @@ -0,0 +1,19 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * + * Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>. + */ + +#include <math.h> + +long double +atanl (long double x) +{ + long double res; + + asm ("fld1\n\t" + "fpatan" + : "=t" (res) : "0" (x)); + return res; +} diff --git a/winsup/mingw/mingwex/math/cbrt.c b/winsup/mingw/mingwex/math/cbrt.c new file mode 100644 index 000000000..93f5c819c --- /dev/null +++ b/winsup/mingw/mingwex/math/cbrt.c @@ -0,0 +1,162 @@ +/* cbrt.c + * + * Cube root + * + * + * + * SYNOPSIS: + * + * double x, y, cbrt(); + * + * y = cbrt( x ); + * + * + * + * DESCRIPTION: + * + * Returns the cube root of the argument, which may be negative. + * + * Range reduction involves determining the power of 2 of + * the argument. A polynomial of degree 2 applied to the + * mantissa, and multiplication by the cube root of 1, 2, or 4 + * approximates the root to within about 0.1%. Then Newton's + * iteration is used three times to converge to an accurate + * result. + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * DEC -10,10 200000 1.8e-17 6.2e-18 + * IEEE 0,1e308 30000 1.5e-16 5.0e-17 + * + */ +/* cbrt.c */ + +/* +Cephes Math Library Release 2.2: January, 1991 +Copyright 1984, 1991 by Stephen L. Moshier +Direct inquiries to 30 Frost Street, Cambridge, MA 02140 +*/ + +/* + Modified for mingwex.a + 2002-07-01 Danny Smith <dannysmith@users.sourceforge.net> + */ +#ifdef __MINGW32__ +#include <math.h> +#include "cephes_mconf.h" +#else +#include "mconf.h" +#endif + + +static const double CBRT2 = 1.2599210498948731647672; +static const double CBRT4 = 1.5874010519681994747517; +static const double CBRT2I = 0.79370052598409973737585; +static const double CBRT4I = 0.62996052494743658238361; + +#ifndef __MINGW32__ +#ifdef ANSIPROT +extern double frexp ( double, int * ); +extern double ldexp ( double, int ); +extern int isnan ( double ); +extern int isfinite ( double ); +#else +double frexp(), ldexp(); +int isnan(), isfinite(); +#endif +#endif + +double cbrt(x) +double x; +{ +int e, rem, sign; +double z; + +#ifdef __MINGW32__ +if (!isfinite (x) || x == 0 ) + return x; +#else + +#ifdef NANS +if( isnan(x) ) + return x; +#endif +#ifdef INFINITIES +if( !isfinite(x) ) + return x; +#endif +if( x == 0 ) + return( x ); + +#endif /* __MINGW32__ */ + +if( x > 0 ) + sign = 1; +else + { + sign = -1; + x = -x; + } + +z = x; +/* extract power of 2, leaving + * mantissa between 0.5 and 1 + */ +x = frexp( x, &e ); + +/* Approximate cube root of number between .5 and 1, + * peak relative error = 9.2e-6 + */ +x = (((-1.3466110473359520655053e-1 * x + + 5.4664601366395524503440e-1) * x + - 9.5438224771509446525043e-1) * x + + 1.1399983354717293273738e0 ) * x + + 4.0238979564544752126924e-1; + +/* exponent divided by 3 */ +if( e >= 0 ) + { + rem = e; + e /= 3; + rem -= 3*e; + if( rem == 1 ) + x *= CBRT2; + else if( rem == 2 ) + x *= CBRT4; + } + + +/* argument less than 1 */ + +else + { + e = -e; + rem = e; + e /= 3; + rem -= 3*e; + if( rem == 1 ) + x *= CBRT2I; + else if( rem == 2 ) + x *= CBRT4I; + e = -e; + } + +/* multiply by power of 2 */ +x = ldexp( x, e ); + +/* Newton iteration */ +x -= ( x - (z/(x*x)) )*0.33333333333333333333; +#ifdef DEC +x -= ( x - (z/(x*x)) )/3.0; +#else +x -= ( x - (z/(x*x)) )*0.33333333333333333333; +#endif + +if( sign < 0 ) + x = -x; +return(x); +} diff --git a/winsup/mingw/mingwex/math/cbrtf.c b/winsup/mingw/mingwex/math/cbrtf.c new file mode 100644 index 000000000..537cf8d98 --- /dev/null +++ b/winsup/mingw/mingwex/math/cbrtf.c @@ -0,0 +1,147 @@ +/* cbrtf.c + * + * Cube root + * + * + * + * SYNOPSIS: + * + * float x, y, cbrtf(); + * + * y = cbrtf( x ); + * + * + * + * DESCRIPTION: + * + * Returns the cube root of the argument, which may be negative. + * + * Range reduction involves determining the power of 2 of + * the argument. A polynomial of degree 2 applied to the + * mantissa, and multiplication by the cube root of 1, 2, or 4 + * approximates the root to within about 0.1%. Then Newton's + * iteration is used to converge to an accurate result. + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE 0,1e38 100000 7.6e-8 2.7e-8 + * + */ +/* cbrt.c */ + +/* +Cephes Math Library Release 2.2: June, 1992 +Copyright 1984, 1987, 1988, 1992 by Stephen L. Moshier +Direct inquiries to 30 Frost Street, Cambridge, MA 02140 +*/ + +/* + Modified for mingwex.a + 2002-07-01 Danny Smith <dannysmith@users.sourceforge.net> + */ +#ifdef __MINGW32__ +#include <math.h> +#include "cephes_mconf.h" +#else +#include "mconf.h" +#endif + +static const float CBRT2 = 1.25992104989487316477; +static const float CBRT4 = 1.58740105196819947475; + +#ifndef __MINGW32__ +#ifdef ANSIC +float frexpf(float, int *), ldexpf(float, int); + +float cbrtf( float xx ) +#else +float frexpf(), ldexpf(); + +float cbrtf(xx) +double xx; +#endif +{ +int e, rem, sign; +float x, z; + +x = xx; + +#else /* __MINGW32__ */ +float cbrtf (float x) +{ +int e, rem, sign; +float z; +#endif /* __MINGW32__ */ + +#ifdef __MINGW32__ +if (!isfinite (x) || x == 0.0F ) + return x; +#else +if( x == 0 ) + return( 0.0 ); +#endif +if( x > 0 ) + sign = 1; +else + { + sign = -1; + x = -x; + } + +z = x; +/* extract power of 2, leaving + * mantissa between 0.5 and 1 + */ +x = frexpf( x, &e ); + +/* Approximate cube root of number between .5 and 1, + * peak relative error = 9.2e-6 + */ +x = (((-0.13466110473359520655053 * x + + 0.54664601366395524503440 ) * x + - 0.95438224771509446525043 ) * x + + 1.1399983354717293273738 ) * x + + 0.40238979564544752126924; + +/* exponent divided by 3 */ +if( e >= 0 ) + { + rem = e; + e /= 3; + rem -= 3*e; + if( rem == 1 ) + x *= CBRT2; + else if( rem == 2 ) + x *= CBRT4; + } + + +/* argument less than 1 */ + +else + { + e = -e; + rem = e; + e /= 3; + rem -= 3*e; + if( rem == 1 ) + x /= CBRT2; + else if( rem == 2 ) + x /= CBRT4; + e = -e; + } + +/* multiply by power of 2 */ +x = ldexpf( x, e ); + +/* Newton iteration */ +x -= ( x - (z/(x*x)) ) * 0.333333333333; + +if( sign < 0 ) + x = -x; +return(x); +} diff --git a/winsup/mingw/mingwex/math/cbrtl.c b/winsup/mingw/mingwex/math/cbrtl.c new file mode 100644 index 000000000..36bd48f70 --- /dev/null +++ b/winsup/mingw/mingwex/math/cbrtl.c @@ -0,0 +1,161 @@ +/* cbrtl.c + * + * Cube root, long double precision + * + * + * + * SYNOPSIS: + * + * long double x, y, cbrtl(); + * + * y = cbrtl( x ); + * + * + * + * DESCRIPTION: + * + * Returns the cube root of the argument, which may be negative. + * + * Range reduction involves determining the power of 2 of + * the argument. A polynomial of degree 2 applied to the + * mantissa, and multiplication by the cube root of 1, 2, or 4 + * approximates the root to within about 0.1%. Then Newton's + * iteration is used three times to converge to an accurate + * result. + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE .125,8 80000 7.0e-20 2.2e-20 + * IEEE exp(+-707) 100000 7.0e-20 2.4e-20 + * + */ + + +/* +Cephes Math Library Release 2.2: January, 1991 +Copyright 1984, 1991 by Stephen L. Moshier +Direct inquiries to 30 Frost Street, Cambridge, MA 02140 +*/ + +/* + Modified for mingwex.a + 2002-07-01 Danny Smith <dannysmith@users.sourceforge.net> + */ +#ifdef __MINGW32__ +#include "cephes_mconf.h" +#else +#include "mconf.h" +#endif + +static const long double CBRT2 = 1.2599210498948731647672L; +static const long double CBRT4 = 1.5874010519681994747517L; +static const long double CBRT2I = 0.79370052598409973737585L; +static const long double CBRT4I = 0.62996052494743658238361L; + +#ifndef __MINGW32__ + +#ifdef ANSIPROT +extern long double frexpl ( long double, int * ); +extern long double ldexpl ( long double, int ); +extern int isnanl ( long double ); +#else +long double frexpl(), ldexpl(); +extern int isnanl(); +#endif + +#ifdef INFINITIES +extern long double INFINITYL; +#endif + +#endif /* __MINGW32__ */ + +long double cbrtl(x) +long double x; +{ +int e, rem, sign; +long double z; + +#ifdef __MINGW32__ +if (!isfinite (x) || x == 0.0L) + return(x); +#else + +#ifdef NANS +if(isnanl(x)) + return(x); +#endif +#ifdef INFINITIES +if( x == INFINITYL) + return(x); +if( x == -INFINITYL) + return(x); +#endif +if( x == 0 ) + return( x ); + +#endif /* __MINGW32__ */ + +if( x > 0 ) + sign = 1; +else + { + sign = -1; + x = -x; + } + +z = x; +/* extract power of 2, leaving + * mantissa between 0.5 and 1 + */ +x = frexpl( x, &e ); + +/* Approximate cube root of number between .5 and 1, + * peak relative error = 1.2e-6 + */ +x = (((( 1.3584464340920900529734e-1L * x + - 6.3986917220457538402318e-1L) * x + + 1.2875551670318751538055e0L) * x + - 1.4897083391357284957891e0L) * x + + 1.3304961236013647092521e0L) * x + + 3.7568280825958912391243e-1L; + +/* exponent divided by 3 */ +if( e >= 0 ) + { + rem = e; + e /= 3; + rem -= 3*e; + if( rem == 1 ) + x *= CBRT2; + else if( rem == 2 ) + x *= CBRT4; + } +else + { /* argument less than 1 */ + e = -e; + rem = e; + e /= 3; + rem -= 3*e; + if( rem == 1 ) + x *= CBRT2I; + else if( rem == 2 ) + x *= CBRT4I; + e = -e; + } + +/* multiply by power of 2 */ +x = ldexpl( x, e ); + +/* Newton iteration */ + +x -= ( x - (z/(x*x)) )*0.3333333333333333333333L; +x -= ( x - (z/(x*x)) )*0.3333333333333333333333L; + +if( sign < 0 ) + x = -x; +return(x); +} diff --git a/winsup/mingw/mingwex/math/ceilf.S b/winsup/mingw/mingwex/math/ceilf.S new file mode 100644 index 000000000..ffcdfc687 --- /dev/null +++ b/winsup/mingw/mingwex/math/ceilf.S @@ -0,0 +1,31 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + + .file "ceilf.S" + .text + .align 4 +.globl _ceilf + .def _ceilf; .scl 2; .type 32; .endef +_ceilf: + flds 4(%esp) + subl $8,%esp + + fstcw 4(%esp) /* store fpu control word */ + + /* We use here %edx although only the low 1 bits are defined. + But none of the operations should care and they are faster + than the 16 bit operations. */ + movl $0x0800,%edx /* round towards +oo */ + orl 4(%esp),%edx + andl $0xfbff,%edx + movl %edx,(%esp) + fldcw (%esp) /* load modified control word */ + + frndint /* round */ + + fldcw 4(%esp) /* restore original control word */ + + addl $8,%esp + ret diff --git a/winsup/mingw/mingwex/math/ceill.S b/winsup/mingw/mingwex/math/ceill.S new file mode 100644 index 000000000..29cb27a62 --- /dev/null +++ b/winsup/mingw/mingwex/math/ceill.S @@ -0,0 +1,33 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * Changes for long double by Ulrich Drepper <drepper@cygnus.com> + */ + + + .file "ceill.S" + .text + .align 4 +.globl _ceill + .def _ceill; .scl 2; .type 32; .endef +_ceill: + fldt 4(%esp) + subl $8,%esp + + fstcw 4(%esp) /* store fpu control word */ + + /* We use here %edx although only the low 1 bits are defined. + But none of the operations should care and they are faster + than the 16 bit operations. */ + movl $0x0800,%edx /* round towards +oo */ + orl 4(%esp),%edx + andl $0xfbff,%edx + movl %edx,(%esp) + fldcw (%esp) /* load modified control word */ + + frndint /* round */ + + fldcw 4(%esp) /* restore original control word */ + + addl $8,%esp + ret diff --git a/winsup/mingw/mingwex/math/cephes_mconf.h b/winsup/mingw/mingwex/math/cephes_mconf.h new file mode 100644 index 000000000..ba8400f34 --- /dev/null +++ b/winsup/mingw/mingwex/math/cephes_mconf.h @@ -0,0 +1,134 @@ +#include <math.h> +#include <errno.h> + +/* constants used by cephes functions */ + +#define MAXNUML 1.189731495357231765021263853E4932L +#define MAXLOGL 1.1356523406294143949492E4L +#define MINLOGL -1.13994985314888605586758E4L +#define LOGE2L 6.9314718055994530941723E-1L +#define LOG2EL 1.4426950408889634073599E0L +#define PIL 3.1415926535897932384626L +#define PIO2L 1.5707963267948966192313L +#define PIO4L 7.8539816339744830961566E-1L + +#define isfinitel isfinite +#define isinfl isinf +#define isnanl isnan +#define signbitl signbit + + +#define IBMPC 1 +#define ANSIPROT 1 +#define MINUSZERO 1 +#define INFINITIES 1 +#define NANS 1 +#define DENORMAL 1 +#define NEGZEROL (-0.0L) +extern long double __INFL; +#define INFINITYL (__INFL) +extern long double __QNANL; +#define NANL (__QNANL) +#define VOLATILE +#define mtherr(fname, code) +#define XPD 0, + +#ifdef _CEPHES_USE_ERRNO +#define _SET_ERRNO(x) errno = (x) +#else +#define _SET_ERRNO(x) +#endif +/* +Cephes Math Library Release 2.2: July, 1992 +Copyright 1984, 1987, 1988, 1992 by Stephen L. Moshier +Direct inquiries to 30 Frost Street, Cambridge, MA 02140 +*/ + + +/* polevll.c + * p1evll.c + * + * Evaluate polynomial + * + * + * + * SYNOPSIS: + * + * int N; + * long double x, y, coef[N+1], polevl[]; + * + * y = polevll( x, coef, N ); + * + * + * + * DESCRIPTION: + * + * Evaluates polynomial of degree N: + * + * 2 N + * y = C + C x + C x +...+ C x + * 0 1 2 N + * + * Coefficients are stored in reverse order: + * + * coef[0] = C , ..., coef[N] = C . + * N 0 + * + * The function p1evll() assumes that coef[N] = 1.0 and is + * omitted from the array. Its calling arguments are + * otherwise the same as polevll(). + * + * + * SPEED: + * + * In the interest of speed, there are no checks for out + * of bounds arithmetic. This routine is used by most of + * the functions in the library. Depending on available + * equipment features, the user may wish to rewrite the + * program in microcode or assembly language. + * + */ + +/* Polynomial evaluator: + * P[0] x^n + P[1] x^(n-1) + ... + P[n] + */ +static __inline__ long double polevll( x, p, n ) +long double x; +const void *p; +int n; +{ +register long double y; +register long double *P = (long double *)p; + +y = *P++; +do + { + y = y * x + *P++; + } +while( --n ); +return(y); +} + + + +/* Polynomial evaluator: + * x^n + P[0] x^(n-1) + P[1] x^(n-2) + ... + P[n] + */ +static __inline__ long double p1evll( x, p, n ) +long double x; +const void *p; +int n; +{ +register long double y; +register long double *P = (long double *)p; + +n -= 1; +y = x + *P++; +do + { + y = y * x + *P++; + } +while( --n ); +return( y ); +} + diff --git a/winsup/mingw/mingwex/math/copysign.S b/winsup/mingw/mingwex/math/copysign.S new file mode 100644 index 000000000..60d6c72db --- /dev/null +++ b/winsup/mingw/mingwex/math/copysign.S @@ -0,0 +1,19 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + + .file "copysign.S" + .text + .align 4 +.globl _copysign + .def _copysign; .scl 2; .type 32; .endef +_copysign: + movl 16(%esp),%edx + movl 8(%esp),%eax + andl $0x80000000,%edx + andl $0x7fffffff,%eax + orl %edx,%eax + movl %eax,8(%esp) + fldl 4(%esp) + ret diff --git a/winsup/mingw/mingwex/math/copysignf.S b/winsup/mingw/mingwex/math/copysignf.S new file mode 100644 index 000000000..8a60c463c --- /dev/null +++ b/winsup/mingw/mingwex/math/copysignf.S @@ -0,0 +1,19 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + + .file "copysignf.S" + .text + .align 4 +.globl _copysignf + .def _copysignf; .scl 2; .type 32; .endef +_copysignf: + movl 8(%esp),%edx + movl 4(%esp),%eax + andl $0x80000000,%edx + andl $0x7fffffff,%eax + orl %edx,%eax + movl %eax,4(%esp) + flds 4(%esp) + ret diff --git a/winsup/mingw/mingwex/copysignl.S b/winsup/mingw/mingwex/math/copysignl.S index 500607e24..4143b37f8 100644 --- a/winsup/mingw/mingwex/copysignl.S +++ b/winsup/mingw/mingwex/math/copysignl.S @@ -6,8 +6,7 @@ .file "copysignl.S" .text - .align 2 - .p2align 4,,15 + .align 4 .globl _copysignl .def _copysignl; .scl 2; .type 32; .endef _copysignl: diff --git a/winsup/mingw/mingwex/math/cosf.S b/winsup/mingw/mingwex/math/cosf.S new file mode 100644 index 000000000..862f6ce6c --- /dev/null +++ b/winsup/mingw/mingwex/math/cosf.S @@ -0,0 +1,29 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * + * Removed glibc header dependancy by Danny Smith + * <dannysmith@users.sourceforge.net> + */ + .file "cosf.S" + .text + .align 4 +.globl _cosl + .def _cosf; .scl 2; .type 32; .endef +_cosf: + flds 4(%esp) + fcos + fnstsw %ax + testl $0x400,%eax + jnz 1f + ret +1: fldpi + fadd %st(0) + fxch %st(1) +2: fprem1 + fnstsw %ax + testl $0x400,%eax + jnz 2b + fstp %st(1) + fcos + ret diff --git a/winsup/mingw/mingwex/math/coshf.c b/winsup/mingw/mingwex/math/coshf.c new file mode 100644 index 000000000..4e44f0811 --- /dev/null +++ b/winsup/mingw/mingwex/math/coshf.c @@ -0,0 +1,3 @@ +#include <math.h> +float coshf (float x) + {return (float) cosh (x);} diff --git a/winsup/mingw/mingwex/math/coshl.c b/winsup/mingw/mingwex/math/coshl.c new file mode 100644 index 000000000..c698e50c0 --- /dev/null +++ b/winsup/mingw/mingwex/math/coshl.c @@ -0,0 +1,110 @@ +/* coshl.c + * + * Hyperbolic cosine, long double precision + * + * + * + * SYNOPSIS: + * + * long double x, y, coshl(); + * + * y = coshl( x ); + * + * + * + * DESCRIPTION: + * + * Returns hyperbolic cosine of argument in the range MINLOGL to + * MAXLOGL. + * + * cosh(x) = ( exp(x) + exp(-x) )/2. + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE +-10000 30000 1.1e-19 2.8e-20 + * + * + * ERROR MESSAGES: + * + * message condition value returned + * cosh overflow |x| > MAXLOGL+LOGE2L INFINITYL + * + * + */ + + +/* +Cephes Math Library Release 2.7: May, 1998 +Copyright 1985, 1991, 1998 by Stephen L. Moshier +*/ + +/* +Modified for mingw +2002-07-22 Danny Smith <dannysmith@users.sourceforge.net> +*/ + +#ifdef __MINGW32__ +#include "cephes_mconf.h" +#else +#include "mconf.h" +#endif + +#ifndef _SET_ERRNO +#define _SET_ERRNO(x) +#endif + + +#ifndef __MINGW32__ +extern long double MAXLOGL, MAXNUML, LOGE2L; +#ifdef ANSIPROT +extern long double expl ( long double ); +extern int isnanl ( long double ); +#else +long double expl(), isnanl(); +#endif +#ifdef INFINITIES +extern long double INFINITYL; +#endif +#ifdef NANS +extern long double NANL; +#endif +#endif /* __MINGW32__ */ + +long double coshl(x) +long double x; +{ +long double y; + +#ifdef NANS +if( isnanl(x) ) + { + _SET_ERRNO(EDOM); + return(x); + } +#endif +if( x < 0 ) + x = -x; +if( x > (MAXLOGL + LOGE2L) ) + { + mtherr( "coshl", OVERFLOW ); + _SET_ERRNO(ERANGE); +#ifdef INFINITIES + return( INFINITYL ); +#else + return( MAXNUML ); +#endif + } +if( x >= (MAXLOGL - LOGE2L) ) + { + y = expl(0.5L * x); + y = (0.5L * y) * y; + return(y); + } +y = expl(x); +y = 0.5L * (y + 1.0L / y); +return( y ); +} diff --git a/winsup/mingw/mingwex/math/cosl.S b/winsup/mingw/mingwex/math/cosl.S new file mode 100644 index 000000000..59d9858b3 --- /dev/null +++ b/winsup/mingw/mingwex/math/cosl.S @@ -0,0 +1,30 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * + * Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>. + * Removed glibc header dependancy by Danny Smith + * <dannysmith@users.sourceforge.net> + */ + .file "cosl.S" + .text + .align 4 +.globl _cosl + .def _cosl; .scl 2; .type 32; .endef +_cosl: + fldt 4(%esp) + fcos + fnstsw %ax + testl $0x400,%eax + jnz 1f + ret +1: fldpi + fadd %st(0) + fxch %st(1) +2: fprem1 + fnstsw %ax + testl $0x400,%eax + jnz 2b + fstp %st(1) + fcos + ret diff --git a/winsup/mingw/mingwex/math/exp2.S b/winsup/mingw/mingwex/math/exp2.S new file mode 100644 index 000000000..320065726 --- /dev/null +++ b/winsup/mingw/mingwex/math/exp2.S @@ -0,0 +1,39 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Adapted for exp2 by Ulrich Drepper <drepper@cygnus.com>. + * Public domain. + */ + + .file "exp2.S" + .text + .align 4 +.globl _exp2 + .def _exp2; .scl 2; .type 32; .endef +_exp2: + fldl 4(%esp) +/* I added the following ugly construct because exp(+-Inf) resulted + in NaN. The ugliness results from the bright minds at Intel. + For the i686 the code can be written better. + -- drepper@cygnus.com. */ + fxam /* Is NaN or +-Inf? */ + fstsw %ax + movb $0x45, %dh + andb %ah, %dh + cmpb $0x05, %dh + je 1f /* Is +-Inf, jump. */ + fld %st + frndint /* int(x) */ + fsubr %st,%st(1) /* fract(x) */ + fxch + f2xm1 /* 2^(fract(x)) - 1 */ + fld1 + faddp /* 2^(fract(x)) */ + fscale /* e^x */ + fstp %st(1) + ret + +1: testl $0x200, %eax /* Test sign. */ + jz 2f /* If positive, jump. */ + fstp %st + fldz /* Set result to 0. */ +2: ret diff --git a/winsup/mingw/mingwex/math/exp2f.S b/winsup/mingw/mingwex/math/exp2f.S new file mode 100644 index 000000000..0707a0cc6 --- /dev/null +++ b/winsup/mingw/mingwex/math/exp2f.S @@ -0,0 +1,39 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Adapted for exp2 by Ulrich Drepper <drepper@cygnus.com>. + * Public domain. + */ + + .file "exp2f.S" + .text + .align 4 +.globl _exp2f + .def _exp2f; .scl 2; .type 32; .endef +_exp2f: + flds 4(%esp) +/* I added the following ugly construct because exp(+-Inf) resulted + in NaN. The ugliness results from the bright minds at Intel. + For the i686 the code can be written better. + -- drepper@cygnus.com. */ + fxam /* Is NaN or +-Inf? */ + fstsw %ax + movb $0x45, %dh + andb %ah, %dh + cmpb $0x05, %dh + je 1f /* Is +-Inf, jump. */ + fld %st + frndint /* int(x) */ + fsubr %st,%st(1) /* fract(x) */ + fxch + f2xm1 /* 2^(fract(x)) - 1 */ + fld1 + faddp /* 2^(fract(x)) */ + fscale /* e^x */ + fstp %st(1) + ret + +1: testl $0x200, %eax /* Test sign. */ + jz 2f /* If positive, jump. */ + fstp %st + fldz /* Set result to 0. */ +2: ret diff --git a/winsup/mingw/mingwex/math/exp2l.S b/winsup/mingw/mingwex/math/exp2l.S new file mode 100644 index 000000000..2457c26f4 --- /dev/null +++ b/winsup/mingw/mingwex/math/exp2l.S @@ -0,0 +1,39 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Adapted for exp2 by Ulrich Drepper <drepper@cygnus.com>. + * Public domain. + */ + + .file "exp2l.S" + .text + .align 4 +.globl _exp2l + .def _exp2l; .scl 2; .type 32; .endef +_exp2l: + fldt 4(%esp) +/* I added the following ugly construct because exp(+-Inf) resulted + in NaN. The ugliness results from the bright minds at Intel. + For the i686 the code can be written better. + -- drepper@cygnus.com. */ + fxam /* Is NaN or +-Inf? */ + fstsw %ax + movb $0x45, %dh + andb %ah, %dh + cmpb $0x05, %dh + je 1f /* Is +-Inf, jump. */ + fld %st + frndint /* int(x) */ + fsubr %st,%st(1) /* fract(x) */ + fxch + f2xm1 /* 2^(fract(x)) - 1 */ + fld1 + faddp /* 2^(fract(x)) */ + fscale /* e^x */ + fstp %st(1) + ret + +1: testl $0x200, %eax /* Test sign. */ + jz 2f /* If positive, jump. */ + fstp %st + fldz /* Set result to 0. */ +2: ret diff --git a/winsup/mingw/mingwex/math/expf.c b/winsup/mingw/mingwex/math/expf.c new file mode 100644 index 000000000..e56e0bc6e --- /dev/null +++ b/winsup/mingw/mingwex/math/expf.c @@ -0,0 +1,3 @@ +#include <math.h> +float expf (float x) + {return (float) exp (x);} diff --git a/winsup/mingw/mingwex/math/expl.c b/winsup/mingw/mingwex/math/expl.c new file mode 100644 index 000000000..5f8face2c --- /dev/null +++ b/winsup/mingw/mingwex/math/expl.c @@ -0,0 +1,77 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * + * Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>. + */ + +/* + * The 8087 method for the exponential function is to calculate + * exp(x) = 2^(x log2(e)) + * after separating integer and fractional parts + * x log2(e) = i + f, |f| <= .5 + * 2^i is immediate but f needs to be precise for long double accuracy. + * Suppress range reduction error in computing f by the following. + * Separate x into integer and fractional parts + * x = xi + xf, |xf| <= .5 + * Separate log2(e) into the sum of an exact number c0 and small part c1. + * c0 + c1 = log2(e) to extra precision + * Then + * f = (c0 xi - i) + c0 xf + c1 x + * where c0 xi is exact and so also is (c0 xi - i). + * -- moshier@na-net.ornl.gov + */ + +#include <math.h> + +static long double c0 = 1.44268798828125L; +static long double c1 = 7.05260771340735992468e-6L; + +long double +expl (long double x) +{ + long double res; + +/* I added the following ugly construct because expl(+-Inf) resulted + in NaN. The ugliness results from the bright minds at Intel. + For the i686 the code can be written better. + -- drepper@cygnus.com. */ + asm ("fxam\n\t" /* Is NaN or +-Inf? */ + "fstsw %%ax\n\t" + "movb $0x45, %%dh\n\t" + "andb %%ah, %%dh\n\t" + "cmpb $0x05, %%dh\n\t" + "je 1f\n\t" /* Is +-Inf, jump. */ + "fldl2e\n\t" /* 1 log2(e) */ + "fmul %%st(1),%%st\n\t" /* 1 x log2(e) */ + "frndint\n\t" /* 1 i */ + "fld %%st(1)\n\t" /* 2 x */ + "frndint\n\t" /* 2 xi */ + "fld %%st(1)\n\t" /* 3 i */ + "fldt %2\n\t" /* 4 c0 */ + "fld %%st(2)\n\t" /* 5 xi */ + "fmul %%st(1),%%st\n\t" /* 5 c0 xi */ + "fsubp %%st,%%st(2)\n\t" /* 4 f = c0 xi - i */ + "fld %%st(4)\n\t" /* 5 x */ + "fsub %%st(3),%%st\n\t" /* 5 xf = x - xi */ + "fmulp %%st,%%st(1)\n\t" /* 4 c0 xf */ + "faddp %%st,%%st(1)\n\t" /* 3 f = f + c0 xf */ + "fldt %3\n\t" /* 4 */ + "fmul %%st(4),%%st\n\t" /* 4 c1 * x */ + "faddp %%st,%%st(1)\n\t" /* 3 f = f + c1 * x */ + "f2xm1\n\t" /* 3 2^(fract(x * log2(e))) - 1 */ + "fld1\n\t" /* 4 1.0 */ + "faddp\n\t" /* 3 2^(fract(x * log2(e))) */ + "fstp %%st(1)\n\t" /* 2 */ + "fscale\n\t" /* 2 scale factor is st(1); e^x */ + "fstp %%st(1)\n\t" /* 1 */ + "fstp %%st(1)\n\t" /* 0 */ + "jmp 2f\n\t" + "1:\ttestl $0x200, %%eax\n\t" /* Test sign. */ + "jz 2f\n\t" /* If positive, jump. */ + "fstp %%st\n\t" + "fldz\n\t" /* Set result to 0. */ + "2:\t\n" + : "=t" (res) : "0" (x), "m" (c0), "m" (c1) : "ax", "dx"); + return res; +} diff --git a/winsup/mingw/mingwex/math/fabs.c b/winsup/mingw/mingwex/math/fabs.c new file mode 100644 index 000000000..c2074e8cb --- /dev/null +++ b/winsup/mingw/mingwex/math/fabs.c @@ -0,0 +1,10 @@ +#include <math.h> + +double +fabs (double x) +{ + double res; + + asm ("fabs;" : "=t" (res) : "0" (x)); + return res; +} diff --git a/winsup/mingw/mingwex/math/fabsf.c b/winsup/mingw/mingwex/math/fabsf.c new file mode 100644 index 000000000..6580f955c --- /dev/null +++ b/winsup/mingw/mingwex/math/fabsf.c @@ -0,0 +1,9 @@ +#include <math.h> + +float +fabsf (float x) +{ + float res; + asm ("fabs;" : "=t" (res) : "0" (x)); + return res; +} diff --git a/winsup/mingw/mingwex/math/fabsl.c b/winsup/mingw/mingwex/math/fabsl.c new file mode 100644 index 000000000..eead724d4 --- /dev/null +++ b/winsup/mingw/mingwex/math/fabsl.c @@ -0,0 +1,9 @@ +#include <math.h> + +long double +fabsl (long double x) +{ + long double res; + asm ("fabs;" : "=t" (res) : "0" (x)); + return res; +} diff --git a/winsup/mingw/mingwex/fdim.c b/winsup/mingw/mingwex/math/fdim.c index 330b09241..330b09241 100644 --- a/winsup/mingw/mingwex/fdim.c +++ b/winsup/mingw/mingwex/math/fdim.c diff --git a/winsup/mingw/mingwex/fdimf.c b/winsup/mingw/mingwex/math/fdimf.c index 02bfc6e5e..02bfc6e5e 100644 --- a/winsup/mingw/mingwex/fdimf.c +++ b/winsup/mingw/mingwex/math/fdimf.c diff --git a/winsup/mingw/mingwex/fdiml.c b/winsup/mingw/mingwex/math/fdiml.c index 1c3d0aaaa..1c3d0aaaa 100644 --- a/winsup/mingw/mingwex/fdiml.c +++ b/winsup/mingw/mingwex/math/fdiml.c diff --git a/winsup/mingw/mingwex/math/files.txt b/winsup/mingw/mingwex/math/files.txt new file mode 100644 index 000000000..b72dcb1d4 --- /dev/null +++ b/winsup/mingw/mingwex/math/files.txt @@ -0,0 +1 @@ +cvs -z9 add 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 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 files.txt floorf.S floorl.S fma.S fmaf.S fmal.c fmax.c fmaxf.c fmaxl.c fmin.c fminf.c fminl.c fmodf.c 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 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 lrintl.c lround.c lroundf.c lroundl.c modff.c modfl.c nearbyint.S nearbyintf.S nearbyintl.S nextafterf.c powf.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 diff --git a/winsup/mingw/mingwex/math/floorf.S b/winsup/mingw/mingwex/math/floorf.S new file mode 100644 index 000000000..8ae8100a7 --- /dev/null +++ b/winsup/mingw/mingwex/math/floorf.S @@ -0,0 +1,35 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * + * Changes for long double by Ulrich Drepper <drepper@cygnus.com> + * + * Removed header file dependency for use in libmingwex.a by + * Danny Smith <dannysmith@users.sourceforge.net> + */ + .file "floorf.S" + .text + .align 4 +.globl _floorf + .def _floorf; .scl 2; .type 32; .endef +_floorf: + flds 4(%esp) + subl $8,%esp + + fstcw 4(%esp) /* store fpu control word */ + + /* We use here %edx although only the low 1 bits are defined. + But none of the operations should care and they are faster + than the 16 bit operations. */ + movl $0x400,%edx /* round towards -oo */ + orl 4(%esp),%edx + andl $0xf7ff,%edx + movl %edx,(%esp) + fldcw (%esp) /* load modified control word */ + + frndint /* round */ + + fldcw 4(%esp) /* restore original control word */ + + addl $8,%esp + ret diff --git a/winsup/mingw/mingwex/math/floorl.S b/winsup/mingw/mingwex/math/floorl.S new file mode 100644 index 000000000..5ab9214b5 --- /dev/null +++ b/winsup/mingw/mingwex/math/floorl.S @@ -0,0 +1,33 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * + * Changes for long double by Ulrich Drepper <drepper@cygnus.com> + * + */ + .file "floorl.S" + .text + .align 4 +.globl _floorl + .def _floorl; .scl 2; .type 32; .endef +_floorl: + fldt 4(%esp) + subl $8,%esp + + fstcw 4(%esp) /* store fpu control word */ + + /* We use here %edx although only the low 1 bits are defined. + But none of the operations should care and they are faster + than the 16 bit operations. */ + movl $0x400,%edx /* round towards -oo */ + orl 4(%esp),%edx + andl $0xf7ff,%edx + movl %edx,(%esp) + fldcw (%esp) /* load modified control word */ + + frndint /* round */ + + fldcw 4(%esp) /* restore original control word */ + + addl $8,%esp + ret diff --git a/winsup/mingw/mingwex/fma.S b/winsup/mingw/mingwex/math/fma.S index d6226653c..d6226653c 100644 --- a/winsup/mingw/mingwex/fma.S +++ b/winsup/mingw/mingwex/math/fma.S diff --git a/winsup/mingw/mingwex/fmaf.S b/winsup/mingw/mingwex/math/fmaf.S index 0d64ac2f1..0d64ac2f1 100644 --- a/winsup/mingw/mingwex/fmaf.S +++ b/winsup/mingw/mingwex/math/fmaf.S diff --git a/winsup/mingw/mingwex/math/fmal.c b/winsup/mingw/mingwex/math/fmal.c new file mode 100644 index 000000000..1fbd41d28 --- /dev/null +++ b/winsup/mingw/mingwex/math/fmal.c @@ -0,0 +1,5 @@ +long double +fmal ( long double _x, long double _y, long double _z) +{ + return ((_x * _y) + _z); +} diff --git a/winsup/mingw/mingwex/fmax.c b/winsup/mingw/mingwex/math/fmax.c index 35c1f45e5..35c1f45e5 100644 --- a/winsup/mingw/mingwex/fmax.c +++ b/winsup/mingw/mingwex/math/fmax.c diff --git a/winsup/mingw/mingwex/fmaxf.c b/winsup/mingw/mingwex/math/fmaxf.c index 079a7e746..079a7e746 100644 --- a/winsup/mingw/mingwex/fmaxf.c +++ b/winsup/mingw/mingwex/math/fmaxf.c diff --git a/winsup/mingw/mingwex/fmaxl.c b/winsup/mingw/mingwex/math/fmaxl.c index 4e38da476..4e38da476 100644 --- a/winsup/mingw/mingwex/fmaxl.c +++ b/winsup/mingw/mingwex/math/fmaxl.c diff --git a/winsup/mingw/mingwex/fmin.c b/winsup/mingw/mingwex/math/fmin.c index 96a6ed111..96a6ed111 100644 --- a/winsup/mingw/mingwex/fmin.c +++ b/winsup/mingw/mingwex/math/fmin.c diff --git a/winsup/mingw/mingwex/fminf.c b/winsup/mingw/mingwex/math/fminf.c index f3d71480d..f3d71480d 100644 --- a/winsup/mingw/mingwex/fminf.c +++ b/winsup/mingw/mingwex/math/fminf.c diff --git a/winsup/mingw/mingwex/fminl.c b/winsup/mingw/mingwex/math/fminl.c index d8a3fea2c..d8a3fea2c 100644 --- a/winsup/mingw/mingwex/fminl.c +++ b/winsup/mingw/mingwex/math/fminl.c diff --git a/winsup/mingw/mingwex/math/fmodf.c b/winsup/mingw/mingwex/math/fmodf.c new file mode 100644 index 000000000..6405d725f --- /dev/null +++ b/winsup/mingw/mingwex/math/fmodf.c @@ -0,0 +1,23 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * + * Adapted for float type by Danny Smith + * <dannysmith@users.sourceforge.net>. + */ + +#include <math.h> + +float +fmodf (float x, float y) +{ + float res; + + asm ("1:\tfprem\n\t" + "fstsw %%ax\n\t" + "sahf\n\t" + "jp 1b\n\t" + "fstp %%st(1)" + : "=t" (res) : "0" (x), "u" (y) : "ax", "st(1)"); + return res; +} diff --git a/winsup/mingw/mingwex/math/fmodl.c b/winsup/mingw/mingwex/math/fmodl.c new file mode 100644 index 000000000..f1c97f10b --- /dev/null +++ b/winsup/mingw/mingwex/math/fmodl.c @@ -0,0 +1,22 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * + * Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>. + */ + +#include <math.h> + +long double +fmodl (long double x, long double y) +{ + long double res; + + asm ("1:\tfprem\n\t" + "fstsw %%ax\n\t" + "sahf\n\t" + "jp 1b\n\t" + "fstp %%st(1)" + : "=t" (res) : "0" (x), "u" (y) : "ax", "st(1)"); + return res; +} diff --git a/winsup/mingw/mingwex/math/fp_consts.c b/winsup/mingw/mingwex/math/fp_consts.c new file mode 100644 index 000000000..285c9d7dc --- /dev/null +++ b/winsup/mingw/mingwex/math/fp_consts.c @@ -0,0 +1,14 @@ + +#include "fp_consts.h" +const union _ieee_rep __QNAN = { __DOUBLE_QNAN_REP }; +const union _ieee_rep __SNAN = { __DOUBLE_SNAN_REP }; +const union _ieee_rep __INF = { __DOUBLE_INF_REP }; +const union _ieee_rep __DENORM = { __DOUBLE_DENORM_REP }; + +/* ISO C99 */ +#undef nan +/* FIXME */ +double nan (const char * tagp __attribute__((unused)) ) + { return __QNAN.double_val; } + + diff --git a/winsup/mingw/mingwex/math/fp_consts.h b/winsup/mingw/mingwex/math/fp_consts.h new file mode 100644 index 000000000..249339501 --- /dev/null +++ b/winsup/mingw/mingwex/math/fp_consts.h @@ -0,0 +1,48 @@ +#ifndef _FP_CONSTS_H +#define _FP_CONSTS_H + +/* +According to IEEE 754 a QNaN has exponent bits of all 1 values and +initial significand bit of 1. A SNaN has has an exponent of all 1 +values and initial significand bit of 0 (with one or more other +significand bits of 1). An Inf has significand of 0 and +exponent of all 1 values. A denormal value has all exponent bits of 0. + +The following does _not_ follow those rules, but uses values +equal to those exported from MS C++ runtime lib, msvcprt.dll +for float and double. MSVC however, does not have long doubles. +*/ + + +#define __DOUBLE_INF_REP { 0, 0, 0, 0x7ff0 } +#define __DOUBLE_QNAN_REP { 0, 0, 0, 0xfff8 } /* { 0, 0, 0, 0x7ff8 } */ +#define __DOUBLE_SNAN_REP { 0, 0, 0, 0xfff0 } /* { 1, 0, 0, 0x7ff0 } */ +#define __DOUBLE_DENORM_REP {1, 0, 0, 0} + +#define D_NAN_MASK 0x7ff0000000000000LL /* this will mask NaN's and Inf's */ + +#define __FLOAT_INF_REP { 0, 0x7f80 } +#define __FLOAT_QNAN_REP { 0, 0xffc0 } /* { 0, 0x7fc0 } */ +#define __FLOAT_SNAN_REP { 0, 0xff80 } /* { 1, 0x7f80 } */ +#define __FLOAT_DENORM_REP {1,0} + +#define F_NAN_MASK 0x7f800000 + +/* + This assumes no implicit (hidden) bit in extended mode. + Padded to 96 bits + */ +#define __LONG_DOUBLE_INF_REP { 0, 0, 0, 0x8000, 0x7fff, 0 } +#define __LONG_DOUBLE_QNAN_REP { 0, 0, 0, 0xc000, 0xffff, 0 } +#define __LONG_DOUBLE_SNAN_REP { 0, 0, 0, 0x8000, 0xffff, 0 } +#define __LONG_DOUBLE_DENORM_REP {1, 0, 0, 0, 0, 0} + +union _ieee_rep +{ + unsigned short rep[6]; + float float_val; + double double_val; + long double ldouble_val; +} ; + +#endif diff --git a/winsup/mingw/mingwex/math/fp_constsf.c b/winsup/mingw/mingwex/math/fp_constsf.c new file mode 100644 index 000000000..5a4afef2b --- /dev/null +++ b/winsup/mingw/mingwex/math/fp_constsf.c @@ -0,0 +1,12 @@ +#include "fp_consts.h" + +const union _ieee_rep __QNANF = { __FLOAT_QNAN_REP }; +const union _ieee_rep __SNANF = { __FLOAT_SNAN_REP }; +const union _ieee_rep __INFF = { __FLOAT_INF_REP }; +const union _ieee_rep __DENORMF = { __FLOAT_DENORM_REP }; + +/* ISO C99 */ +#undef nanf +/* FIXME */ +float nanf(const char * tagp __attribute__((unused)) ) + { return __QNANF.float_val;} diff --git a/winsup/mingw/mingwex/math/fp_constsl.c b/winsup/mingw/mingwex/math/fp_constsl.c new file mode 100644 index 000000000..44fdb7fd3 --- /dev/null +++ b/winsup/mingw/mingwex/math/fp_constsl.c @@ -0,0 +1,12 @@ +#include "fp_consts.h" + +const union _ieee_rep __QNANL = { __LONG_DOUBLE_QNAN_REP }; +const union _ieee_rep __SNANL = { __LONG_DOUBLE_SNAN_REP }; +const union _ieee_rep __INFL = { __LONG_DOUBLE_INF_REP }; +const union _ieee_rep __DENORML = { __LONG_DOUBLE_DENORM_REP }; + + +#undef nanl +/* FIXME */ +long double nanl (const char * tagp __attribute__((unused)) ) + { return __QNANL.ldouble_val; } diff --git a/winsup/mingw/mingwex/fpclassify.c b/winsup/mingw/mingwex/math/fpclassify.c index f8cd8cb44..f8cd8cb44 100644 --- a/winsup/mingw/mingwex/fpclassify.c +++ b/winsup/mingw/mingwex/math/fpclassify.c diff --git a/winsup/mingw/mingwex/fpclassifyf.c b/winsup/mingw/mingwex/math/fpclassifyf.c index aca4e59f1..aca4e59f1 100644 --- a/winsup/mingw/mingwex/fpclassifyf.c +++ b/winsup/mingw/mingwex/math/fpclassifyf.c diff --git a/winsup/mingw/mingwex/fpclassifyl.c b/winsup/mingw/mingwex/math/fpclassifyl.c index 9979d6278..9979d6278 100644 --- a/winsup/mingw/mingwex/fpclassifyl.c +++ b/winsup/mingw/mingwex/math/fpclassifyl.c diff --git a/winsup/mingw/mingwex/math/frexpf.c b/winsup/mingw/mingwex/math/frexpf.c new file mode 100644 index 000000000..df262abc5 --- /dev/null +++ b/winsup/mingw/mingwex/math/frexpf.c @@ -0,0 +1,3 @@ +#include <math.h> +float frexpf (float x, int* expn) + {return (float)frexp(x, expn);} diff --git a/winsup/mingw/mingwex/math/frexpl.S b/winsup/mingw/mingwex/math/frexpl.S new file mode 100644 index 000000000..2b691c87f --- /dev/null +++ b/winsup/mingw/mingwex/math/frexpl.S @@ -0,0 +1,71 @@ +/* + Cephes Math Library Release 2.7: May, 1998 + Copyright 1984, 1987, 1988, 1992, 1998 by Stephen L. Moshier + + Extracted from floorl.387 for use in libmingwex.a by + Danny Smith <dannysmith@users.sourceforge.net> + 2002-06-20 +*/ + +/* + * frexpl(long double x, int* expnt) extracts the exponent from x. + * It returns an integer power of two to expnt and the significand + * between 0.5 and 1 to y. Thus x = y * 2**expn. + */ + .align 2 +.globl _frexpl +_frexpl: + pushl %ebp + movl %esp,%ebp + subl $24,%esp + pushl %esi + pushl %ebx + fldt 8(%ebp) + movl 20(%ebp),%ebx + fld %st(0) + fstpt -12(%ebp) + leal -4(%ebp),%ecx + movw -4(%ebp),%dx + andl $32767,%edx + jne L25 + fldz + fucompp + fnstsw %ax + andb $68,%ah + xorb $64,%ah + jne L21 + movl $0,(%ebx) + fldz + jmp L24 + .align 2,0x90 + .align 2,0x90 +L21: + fldt -12(%ebp) + fadd %st(0),%st + fstpt -12(%ebp) + decl %edx + movw (%ecx),%si + andl $32767,%esi + jne L22 + cmpl $-66,%edx + jg L21 +L22: + addl %esi,%edx + jmp L19 + .align 2,0x90 +L25: + fstp %st(0) +L19: + addl $-16382,%edx + movl %edx,(%ebx) + movw (%ecx),%ax + andl $-32768,%eax + orl $16382,%eax + movw %ax,(%ecx) + fldt -12(%ebp) +L24: + leal -32(%ebp),%esp + popl %ebx + popl %esi + leave + ret diff --git a/winsup/mingw/mingwex/fucom.c b/winsup/mingw/mingwex/math/fucom.c index 80c937262..80c937262 100644 --- a/winsup/mingw/mingwex/fucom.c +++ b/winsup/mingw/mingwex/math/fucom.c diff --git a/winsup/mingw/mingwex/math/hypotf.c b/winsup/mingw/mingwex/math/hypotf.c new file mode 100644 index 000000000..ee67a45dc --- /dev/null +++ b/winsup/mingw/mingwex/math/hypotf.c @@ -0,0 +1,4 @@ +#include <math.h> + +float hypotf (float x, float y) + { return (float) _hypot (x, y);} diff --git a/winsup/mingw/mingwex/math/hypotl.c b/winsup/mingw/mingwex/math/hypotl.c new file mode 100644 index 000000000..2bcfc8638 --- /dev/null +++ b/winsup/mingw/mingwex/math/hypotl.c @@ -0,0 +1,58 @@ +#include <math.h> + +typedef union +{ + long double value; + struct + { + unsigned mantissa[2]; + unsigned sign_exponent : 16; + unsigned empty : 16; + } parts; +} ieee_long_double_shape_type; + + + +/* Get int from the exponent of a long double. */ +static __inline__ void +get_ld_exp (unsigned exp,long double d) +{ + ieee_long_double_shape_type u; + u.value = d; + exp = u.parts.sign_exponent; +} + +/* Set exponent of a long double from an int. */ +static __inline__ void +set_ld_exp (long double d,unsigned exp) +{ + ieee_long_double_shape_type u; + u.value = d; + u.parts.sign_exponent = exp; + d = u.value; +} + +long double +hypotl (long double x, long double y) +{ + unsigned exx = 0U; + unsigned eyy = 0U; + unsigned scale; + long double xx =fabsl(x); + long double yy =fabsl(y); + if (!isfinite(xx) || !isfinite(yy)) + return xx + yy; /* Return INF or NAN. */ + + /* Scale to avoid overflow.*/ + get_ld_exp (exx, xx); + get_ld_exp (eyy, yy); + scale = (exx > eyy ? exx : eyy); + if (scale == 0) + return 0.0L; + set_ld_exp (xx, exx - scale); + set_ld_exp (yy, eyy - scale); + xx = sqrtl(xx * xx + yy * yy); + get_ld_exp (exx,xx); + set_ld_exp (xx, exx + scale); + return xx; +} diff --git a/winsup/mingw/mingwex/math/ilogb.S b/winsup/mingw/mingwex/math/ilogb.S new file mode 100644 index 000000000..2335b5146 --- /dev/null +++ b/winsup/mingw/mingwex/math/ilogb.S @@ -0,0 +1,37 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + + + .file "ilogb.S" + .text + .align 4 +.globl _ilogb + .def _ilogb; .scl 2; .type 32; .endef +_ilogb: + + fldl 4(%esp) +/* I added the following ugly construct because ilogb(+-Inf) is + required to return INT_MAX in ISO C99. + -- jakub@redhat.com. */ + fxam /* Is NaN or +-Inf? */ + fstsw %ax + movb $0x45, %dh + andb %ah, %dh + cmpb $0x05, %dh + je 1f /* Is +-Inf, jump. */ + + fxtract + pushl %eax + fstp %st + + fistpl (%esp) + fwait + popl %eax + + ret + +1: fstp %st + movl $0x7fffffff, %eax + ret diff --git a/winsup/mingw/mingwex/math/ilogbf.S b/winsup/mingw/mingwex/math/ilogbf.S new file mode 100644 index 000000000..fa3e78e84 --- /dev/null +++ b/winsup/mingw/mingwex/math/ilogbf.S @@ -0,0 +1,35 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + + .file "ilogbf.S" + .text + .align 4 +.globl _ilogbf + .def _ilogbf; .scl 2; .type 32; .endef +_ilogbf: + flds 4(%esp) +/* I added the following ugly construct because ilogb(+-Inf) is + required to return INT_MAX in ISO C99. + -- jakub@redhat.com. */ + fxam /* Is NaN or +-Inf? */ + fstsw %ax + movb $0x45, %dh + andb %ah, %dh + cmpb $0x05, %dh + je 1f /* Is +-Inf, jump. */ + + fxtract + pushl %eax + fstp %st + + fistpl (%esp) + fwait + popl %eax + + ret + +1: fstp %st + movl $0x7fffffff, %eax + ret diff --git a/winsup/mingw/mingwex/math/ilogbl.S b/winsup/mingw/mingwex/math/ilogbl.S new file mode 100644 index 000000000..b9dc6ea72 --- /dev/null +++ b/winsup/mingw/mingwex/math/ilogbl.S @@ -0,0 +1,36 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Changes for long double by Ulrich Drepper <drepper@cygnus.com> + * Public domain. + */ + + .file "ilogbl.S" + .text + .align 4 +.globl _ilogbl + .def _ilogbl; .scl 2; .type 32; .endef +_ilogbl: + fldt 4(%esp) +/* I added the following ugly construct because ilogb(+-Inf) is + required to return INT_MAX in ISO C99. + -- jakub@redhat.com. */ + fxam /* Is NaN or +-Inf? */ + fstsw %ax + movb $0x45, %dh + andb %ah, %dh + cmpb $0x05, %dh + je 1f /* Is +-Inf, jump. */ + + fxtract + pushl %eax + fstp %st + + fistpl (%esp) + fwait + popl %eax + + ret + +1: fstp %st + movl $0x7fffffff, %eax + ret diff --git a/winsup/mingw/mingwex/isnan.c b/winsup/mingw/mingwex/math/isnan.c index b38bc290e..b38bc290e 100644 --- a/winsup/mingw/mingwex/isnan.c +++ b/winsup/mingw/mingwex/math/isnan.c diff --git a/winsup/mingw/mingwex/isnanf.c b/winsup/mingw/mingwex/math/isnanf.c index 73fe0eb02..73fe0eb02 100644 --- a/winsup/mingw/mingwex/isnanf.c +++ b/winsup/mingw/mingwex/math/isnanf.c diff --git a/winsup/mingw/mingwex/isnanl.c b/winsup/mingw/mingwex/math/isnanl.c index 86d0088b4..86d0088b4 100644 --- a/winsup/mingw/mingwex/isnanl.c +++ b/winsup/mingw/mingwex/math/isnanl.c diff --git a/winsup/mingw/mingwex/math/ldexpf.c b/winsup/mingw/mingwex/math/ldexpf.c new file mode 100644 index 000000000..5d01a0184 --- /dev/null +++ b/winsup/mingw/mingwex/math/ldexpf.c @@ -0,0 +1,3 @@ +#include <math.h> +float ldexpf (float x, int expn) + {return (float) ldexp (x, expn);} diff --git a/winsup/mingw/mingwex/math/ldexpl.c b/winsup/mingw/mingwex/math/ldexpl.c new file mode 100644 index 000000000..e4477bf6c --- /dev/null +++ b/winsup/mingw/mingwex/math/ldexpl.c @@ -0,0 +1,14 @@ +#include <math.h> +#include <errno.h> + + +long double ldexpl(long double x, int expn) +{ + if (isfinite (x) && x != 0.0L) + { + x = scalbnl (x , expn); + if (!isfinite (x) || x == 0.0) errno = ERANGE; + } + return x; +} + diff --git a/winsup/mingw/mingwex/math/llrint.c b/winsup/mingw/mingwex/math/llrint.c new file mode 100644 index 000000000..b6d9f3273 --- /dev/null +++ b/winsup/mingw/mingwex/math/llrint.c @@ -0,0 +1,10 @@ +#include <math.h> + +long long llrint (double x) +{ + long long retval; + __asm__ __volatile__ \ + ("fistpll %0" : "=m" (retval) : "t" (x) : "st"); \ + return retval; +} + diff --git a/winsup/mingw/mingwex/math/llrintf.c b/winsup/mingw/mingwex/math/llrintf.c new file mode 100644 index 000000000..7fa67dbdf --- /dev/null +++ b/winsup/mingw/mingwex/math/llrintf.c @@ -0,0 +1,9 @@ +#include <math.h> + +long long llrintf (float x) +{ + long long retval; + __asm__ __volatile__ \ + ("fistpll %0" : "=m" (retval) : "t" (x) : "st"); \ + return retval; +} diff --git a/winsup/mingw/mingwex/math/llrintl.c b/winsup/mingw/mingwex/math/llrintl.c new file mode 100644 index 000000000..948d96265 --- /dev/null +++ b/winsup/mingw/mingwex/math/llrintl.c @@ -0,0 +1,10 @@ +#include <math.h> + +long long llrintl (long double x) +{ + long long retval; + __asm__ __volatile__ \ + ("fistpll %0" : "=m" (retval) : "t" (x) : "st"); \ + return retval; +} + diff --git a/winsup/mingw/mingwex/math/llround.c b/winsup/mingw/mingwex/math/llround.c new file mode 100644 index 000000000..8b9fe9781 --- /dev/null +++ b/winsup/mingw/mingwex/math/llround.c @@ -0,0 +1,24 @@ +#include <fenv.h> +#include <math.h> + +long +lround (double x) { + long retval; + unsigned short saved_cw, _cw; + __asm__ ( + "fnstcw %0;" : "=m" (saved_cw) + ); /* save control word */ + _cw = ~(FE_TONEAREST | FE_DOWNWARD | FE_UPWARD | FE_TOWARDZERO) + | (x > 0.0 ? FE_UPWARD : FE_DOWNWARD); /* round away from zero */ + __asm__ ( + "fldcw %0;" : : "m" (_cw) + ); /* load the rounding control */ + __asm__ __volatile__ ( + "fistpll %0" : "=m" (retval) : "t" (x) : "st" + ); + __asm__ ( + "fldcw %0;" : : "m" (saved_cw) + ); /* restore control word */ + return retval; +} + diff --git a/winsup/mingw/mingwex/math/llroundf.c b/winsup/mingw/mingwex/math/llroundf.c new file mode 100644 index 000000000..b6c431f9d --- /dev/null +++ b/winsup/mingw/mingwex/math/llroundf.c @@ -0,0 +1,23 @@ +#include <fenv.h> +#include <math.h> + +long long +llroundf (float x) { + long long retval; + unsigned short saved_cw, _cw; + __asm__ ( + "fnstcw %0;" : "=m" (saved_cw) + ); /* save control word */ + _cw = ~(FE_TONEAREST | FE_DOWNWARD | FE_UPWARD | FE_TOWARDZERO) + | (x > 0.0 ? FE_UPWARD : FE_DOWNWARD); /* round away from zero */ + __asm__ ( + "fldcw %0;" : : "m" (_cw) + ); /* load the rounding control */ + __asm__ __volatile__ ( + "fistpll %0" : "=m" (retval) : "t" (x) : "st" + ); + __asm__ ( + "fldcw %0;" : : "m" (saved_cw) + ); /* restore control word */ + return retval; +} diff --git a/winsup/mingw/mingwex/math/llroundl.c b/winsup/mingw/mingwex/math/llroundl.c new file mode 100644 index 000000000..ad334fb2c --- /dev/null +++ b/winsup/mingw/mingwex/math/llroundl.c @@ -0,0 +1,22 @@ +#include <fenv.h> +#include <math.h> + +long long +llroundl (long double x) { + long long retval; + unsigned short saved_cw, _cw; + __asm__ ( + "fnstcw %0;" : "=m" (saved_cw) + ); /* save control word */ + _cw = ~(FE_TONEAREST | FE_DOWNWARD | FE_UPWARD | FE_TOWARDZERO) + | (x > 0.0 ? FE_UPWARD : FE_DOWNWARD); /* round away from zero */ + __asm__ ( + "fldcw %0;" : : "m" (_cw) + ); /* load the rounding control */ + __asm__ __volatile__ ( + "fistpll %0" : "=m" (retval) : "t" (x) : "st"); + __asm__ ( + "fldcw %0;" : : "m" (saved_cw) + ); /* restore control word */ + return retval; +} diff --git a/winsup/mingw/mingwex/math/log10f.S b/winsup/mingw/mingwex/math/log10f.S new file mode 100644 index 000000000..90fc9af92 --- /dev/null +++ b/winsup/mingw/mingwex/math/log10f.S @@ -0,0 +1,48 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * Adapted for float type by Ulrich Drepper <drepper@cygnus.com>. + * + * Changed to use fyl2xp1 for values near 1, <drepper@cygnus.com>. + */ + + .file "log10f.S" + .text + .align 4 +one: .double 1.0 + /* It is not important that this constant is precise. It is only + a value which is known to be on the safe side for using the + fyl2xp1 instruction. */ +limit: .double 0.29 + + .text + .align 4 +.globl _log10f + .def _log10f; .scl 2; .type 32; .endef +_log10f: + fldlg2 // log10(2) + flds 4(%esp) // x : log10(2) + fxam + fnstsw + fld %st // x : x : log10(2) + sahf + jc 3f // in case x is NaN or ±Inf +4: fsubl one // x-1 : x : log10(2) + fld %st // x-1 : x-1 : x : log10(2) + fabs // |x-1| : x-1 : x : log10(2) + fcompl limit // x-1 : x : log10(2) + fnstsw // x-1 : x : log10(2) + andb $0x45, %ah + jz 2f + fstp %st(1) // x-1 : log10(2) + fyl2xp1 // log10(x) + ret + +2: fstp %st(0) // x : log10(2) + fyl2x // log10(x) + ret + +3: jp 4b // in case x is ±Inf + fstp %st(1) + fstp %st(1) + ret diff --git a/winsup/mingw/mingwex/math/log10l.S b/winsup/mingw/mingwex/math/log10l.S new file mode 100644 index 000000000..8c046a09d --- /dev/null +++ b/winsup/mingw/mingwex/math/log10l.S @@ -0,0 +1,52 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * + * Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>. + * + * Changed to use fyl2xp1 for values near 1, <drepper@cygnus.com>. + * + * Removed header file dependency for use in libmingwex.a by + * Danny Smith <dannysmith@users.sourceforge.net> + */ + + .file "log10l.S" + .text + .align 4 +one: .double 1.0 + /* It is not important that this constant is precise. It is only + a value which is known to be on the safe side for using the + fyl2xp1 instruction. */ +limit: .double 0.29 + + .text + .align 4 +.globl _log10l + .def _log10l; .scl 2; .type 32; .endef +_log10l: + fldlg2 // log10(2) + fldt 4(%esp) // x : log10(2) + fxam + fnstsw + fld %st // x : x : log10(2) + sahf + jc 3f // in case x is NaN or ±Inf +4: fsubl one // x-1 : x : log10(2) + fld %st // x-1 : x-1 : x : log10(2) + fabs // |x-1| : x-1 : x : log10(2) + fcompl limit // x-1 : x : log10(2) + fnstsw // x-1 : x : log10(2) + andb $0x45, %ah + jz 2f + fstp %st(1) // x-1 : log10(2) + fyl2xp1 // log10(x) + ret + +2: fstp %st(0) // x : log10(2) + fyl2x // log10(x) + ret + +3: jp 4b // in case x is ±Inf + fstp %st(1) + fstp %st(1) + ret diff --git a/winsup/mingw/mingwex/math/log1p.S b/winsup/mingw/mingwex/math/log1p.S new file mode 100644 index 000000000..a38816cb3 --- /dev/null +++ b/winsup/mingw/mingwex/math/log1p.S @@ -0,0 +1,47 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * Removed header file dependency for use in libmingwex.a by + * Danny Smith <dannysmith@users.sourceforge.net> + */ + + .file "log1p.S" + .text + .align 4 + /* The fyl2xp1 can only be used for values in + -1 + sqrt(2) / 2 <= x <= 1 - sqrt(2) / 2 + 0.29 is a safe value. + */ +limit: .double 0.29 +one: .double 1.0 +/* + * Use the fyl2xp1 function when the argument is in the range -0.29 to 0.29, + * otherwise fyl2x with the needed extra computation. + */ +.globl _log1p; + .def _log1p; .scl 2; .type 32; .endef +_log1p: + fldln2 + fldl 4(%esp) + fxam + fnstsw + fld %st + sahf + jc 3f // in case x is NaN or ±Inf + +4: fabs + fcompl limit + fnstsw + sahf + jc 2f + faddl one + fyl2x + ret + +2: fyl2xp1 + ret + +3: jp 4b // in case x is ±Inf + fstp %st(1) + fstp %st(1) + ret diff --git a/winsup/mingw/mingwex/math/log1pf.S b/winsup/mingw/mingwex/math/log1pf.S new file mode 100644 index 000000000..1d9949f2a --- /dev/null +++ b/winsup/mingw/mingwex/math/log1pf.S @@ -0,0 +1,47 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * Removed header file dependency for use in libmingwex.a by + * Danny Smith <dannysmith@users.sourceforge.net> + */ + + .file "log1pf.S" + .text + .align 4 + /* The fyl2xp1 can only be used for values in + -1 + sqrt(2) / 2 <= x <= 1 - sqrt(2) / 2 + 0.29 is a safe value. + */ +limit: .float 0.29 +one: .float 1.0 +/* + * Use the fyl2xp1 function when the argument is in the range -0.29 to 0.29, + * otherwise fyl2x with the needed extra computation. + */ +.globl _log1pf; + .def _log1pf; .scl 2; .type 32; .endef +_log1pf: + fldln2 + flds 4(%esp) + fxam + fnstsw + fld %st + sahf + jc 3f // in case x is NaN or ±Inf + +4: fabs + fcomps limit + fnstsw + sahf + jc 2f + fadds one + fyl2x + ret + +2: fyl2xp1 + ret + +3: jp 4b // in case x is ±Inf + fstp %st(1) + fstp %st(1) + ret diff --git a/winsup/mingw/mingwex/math/log1pl.S b/winsup/mingw/mingwex/math/log1pl.S new file mode 100644 index 000000000..5ce4fbaaa --- /dev/null +++ b/winsup/mingw/mingwex/math/log1pl.S @@ -0,0 +1,54 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * + * Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>. +* Removed header file dependency for use in libmingwex.a by + * Danny Smith <dannysmith@users.sourceforge.net> + */ + + .file "log1pl.S" + .text + .align 4 + /* The fyl2xp1 can only be used for values in + -1 + sqrt(2) / 2 <= x <= 1 - sqrt(2) / 2 + 0.29 is a safe value. + */ +limit: .tfloat 0.29 + /* Please note: we use a double value here. Since 1.0 has + an exact representation this does not effect the accuracy + but it helps to optimize the code. */ +one: .double 1.0 + +/* + * Use the fyl2xp1 function when the argument is in the range -0.29 to 0.29, + * otherwise fyl2x with the needed extra computation. + */ +.globl _log1pl; + .def _log1pl; .scl 2; .type 32; .endef +_log1pl: + fldln2 + fldt 4(%esp) + fxam + fnstsw + fld %st + sahf + jc 3f // in case x is NaN or ±Inf +4: + fabs + fldt limit + fcompp + fnstsw + sahf + jnc 2f + faddl one + fyl2x + ret + +2: fyl2xp1 + ret + +3: jp 4b // in case x is ±Inf + fstp %st(1) + fstp %st(1) + ret diff --git a/winsup/mingw/mingwex/math/log2.S b/winsup/mingw/mingwex/math/log2.S new file mode 100644 index 000000000..08f008310 --- /dev/null +++ b/winsup/mingw/mingwex/math/log2.S @@ -0,0 +1,51 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Adapted for use as log2 by Ulrich Drepper <drepper@cygnus.com>. + * Public domain. + * + * Changed to use fyl2xp1 for values near 1, <drepper@cygnus.com>. + * + * Removed header file dependency for use in libmingwex.a by + * Danny Smith <dannysmith@users.sourceforge.net> + */ + + .file "log2.S" + .text + .align 4 +one: .double 1.0 + /* It is not important that this constant is precise. It is only + a value which is known to be on the safe side for using the + fyl2xp1 instruction. */ +limit: .double 0.29 + + .text + .align 4 +.globl _log2 + .def _log2; .scl 2; .type 32; .endef +_log2: + fldl one + fldl 4(%esp) // x : 1 + fxam + fnstsw + fld %st // x : x : 1 + sahf + jc 3f // in case x is NaN or ±Inf +4: fsub %st(2), %st // x-1 : x : 1 + fld %st // x-1 : x-1 : x : 1 + fabs // |x-1| : x-1 : x : 1 + fcompl limit // x-1 : x : 1 + fnstsw // x-1 : x : 1 + andb $0x45, %ah + jz 2f + fstp %st(1) // x-1 : 1 + fyl2xp1 // log(x) + ret + +2: fstp %st(0) // x : 1 + fyl2x // log(x) + ret + +3: jp 4b // in case x is ±Inf + fstp %st(1) + fstp %st(1) + ret diff --git a/winsup/mingw/mingwex/math/log2f.S b/winsup/mingw/mingwex/math/log2f.S new file mode 100644 index 000000000..211abba3d --- /dev/null +++ b/winsup/mingw/mingwex/math/log2f.S @@ -0,0 +1,51 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Adapted for use as log2 by Ulrich Drepper <drepper@cygnus.com>. + * Public domain. + * + * Changed to use fyl2xp1 for values near 1, <drepper@cygnus.com>. + * + * Removed header file dependency for use in libmingwex.a by + * Danny Smith <dannysmith@users.sourceforge.net> + */ + + .file "log2f.S" + .text + .align 4 +one: .double 1.0 + /* It is not important that this constant is precise. It is only + a value which is known to be on the safe side for using the + fyl2xp1 instruction. */ +limit: .double 0.29 + + .text + .align 4 +.globl _log2f + .def _log2f; .scl 2; .type 32; .endef +_log2f: + fldl one + flds 4(%esp) // x : 1 + fxam + fnstsw + fld %st // x : x : 1 + sahf + jc 3f // in case x is NaN or ±Inf +4: fsub %st(2), %st // x-1 : x : 1 + fld %st // x-1 : x-1 : x : 1 + fabs // |x-1| : x-1 : x : 1 + fcompl limit // x-1 : x : 1 + fnstsw // x-1 : x : 1 + andb $0x45, %ah + jz 2f + fstp %st(1) // x-1 : 1 + fyl2xp1 // log(x) + ret + +2: fstp %st(0) // x : 1 + fyl2x // log(x) + ret + +3: jp 4b // in case x is ±Inf + fstp %st(1) + fstp %st(1) + ret diff --git a/winsup/mingw/mingwex/math/log2l.S b/winsup/mingw/mingwex/math/log2l.S new file mode 100644 index 000000000..52503fc52 --- /dev/null +++ b/winsup/mingw/mingwex/math/log2l.S @@ -0,0 +1,48 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Adapted for use as log2 by Ulrich Drepper <drepper@cygnus.com>. + * Public domain. + * + * Changed to use fyl2xp1 for values near 1, <drepper@cygnus.com>. + */ + + .file "log2l.S" + .text + .align 4 +one: .double 1.0 + /* It is not important that this constant is precise. It is only + a value which is known to be on the safe side for using the + fyl2xp1 instruction. */ +limit: .double 0.29 + + .text + .align 4 +.globl _log2l + .def _log2l; .scl 2; .type 32; .endef +_log2l: + fldl one + fldt 4(%esp) // x : 1 + fxam + fnstsw + fld %st // x : x : 1 + sahf + jc 3f // in case x is NaN or ±Inf +4: fsub %st(2), %st // x-1 : x : 1 + fld %st // x-1 : x-1 : x : 1 + fabs // |x-1| : x-1 : x : 1 + fcompl limit // x-1 : x : 1 + fnstsw // x-1 : x : 1 + andb $0x45, %ah + jz 2f + fstp %st(1) // x-1 : 1 + fyl2xp1 // log(x) + ret + +2: fstp %st(0) // x : 1 + fyl2x // log(x) + ret + +3: jp 4b // in case x is ±Inf + fstp %st(1) + fstp %st(1) + ret diff --git a/winsup/mingw/mingwex/math/logb.c b/winsup/mingw/mingwex/math/logb.c new file mode 100644 index 000000000..cdff13647 --- /dev/null +++ b/winsup/mingw/mingwex/math/logb.c @@ -0,0 +1,16 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Changes for long double by Ulrich Drepper <drepper@cygnus.com> + * Public domain. + */ + +#include <math.h> + +double +logb (double x) +{ + double res; + asm ("fxtract\n\t" + "fstp %%st" : "=t" (res) : "0" (x)); + return res; +} diff --git a/winsup/mingw/mingwex/math/logbf.c b/winsup/mingw/mingwex/math/logbf.c new file mode 100644 index 000000000..b5f57d2e1 --- /dev/null +++ b/winsup/mingw/mingwex/math/logbf.c @@ -0,0 +1,16 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Changes for long double by Ulrich Drepper <drepper@cygnus.com> + * Public domain. + */ + +#include <math.h> + +float +logbf (float x) +{ + float res; + asm ("fxtract\n\t" + "fstp %%st" : "=t" (res) : "0" (x)); + return res; +} diff --git a/winsup/mingw/mingwex/math/logbl.c b/winsup/mingw/mingwex/math/logbl.c new file mode 100644 index 000000000..f1448eb99 --- /dev/null +++ b/winsup/mingw/mingwex/math/logbl.c @@ -0,0 +1,17 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Changes for long double by Ulrich Drepper <drepper@cygnus.com> + * Public domain. + */ + +#include <math.h> + +long double +logbl (long double x) +{ + long double res; + + asm ("fxtract\n\t" + "fstp %%st" : "=t" (res) : "0" (x)); + return res; +} diff --git a/winsup/mingw/mingwex/math/logf.S b/winsup/mingw/mingwex/math/logf.S new file mode 100644 index 000000000..32119ecde --- /dev/null +++ b/winsup/mingw/mingwex/math/logf.S @@ -0,0 +1,39 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * Adapted for float by Ulrich Drepper <drepper@cygnus.com>. + * + * Changed to use fyl2xp1 for values near 1, <drepper@cygnus.com>. + */ + + .file "logf.S" + .text + .align 4 +one: .double 1.0 + /* It is not important that this constant is precise. It is only + a value which is known to be on the safe side for using the + fyl2xp1 instruction. */ +limit: .double 0.29 + + .text + .align 4 +.globl _logf + .def _logf; .scl 2; .type 32; .endef +_logf: + fldln2 // log(2) + flds 4(%esp) // x : log(2) + fld %st // x : x : log(2) + fsubl one // x-1 : x : log(2) + fld %st // x-1 : x-1 : x : log(2) + fabs // |x-1| : x-1 : x : log(2) + fcompl limit // x-1 : x : log(2) + fnstsw // x-1 : x : log(2) + andb $0x45, %ah + jz 2f + fstp %st(1) // x-1 : log(2) + fyl2xp1 // log(x) + ret + +2: fstp %st(0) // x : log(2) + fyl2x // log(x) + ret diff --git a/winsup/mingw/mingwex/math/logl.S b/winsup/mingw/mingwex/math/logl.S new file mode 100644 index 000000000..8dc144915 --- /dev/null +++ b/winsup/mingw/mingwex/math/logl.S @@ -0,0 +1,40 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * + * Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>. + * + * Removed header file dependency for use in libmingwex.a by + * Danny Smith <dannysmith@users.sourceforge.net> + */ + .file "logl.S" + .text + .align 4 +one: .double 1.0 + /* It is not important that this constant is precise. It is only + a value which is known to be on the safe side for using the + fyl2xp1 instruction. */ +limit: .double 0.29 + + .text + .align 4 +.globl _logl + .def _logl; .scl 2; .type 32; .endef +_logl: + fldln2 // log(2) + fldt 4(%esp) // x : log(2) + fld %st // x : x : log(2) + fsubl one // x-1 : x : log(2) + fld %st // x-1 : x-1 : x : log(2) + fabs // |x-1| : x-1 : x : log(2) + fcompl limit // x-1 : x : log(2) + fnstsw // x-1 : x : log(2) + andb $0x45, %ah + jz 2f + fstp %st(1) // x-1 : log(2) + fyl2xp1 // log(x) + ret + +2: fstp %st(0) // x : log(2) + fyl2x // log(x) + ret diff --git a/winsup/mingw/mingwex/math/lrint.c b/winsup/mingw/mingwex/math/lrint.c new file mode 100644 index 000000000..7dfa233a8 --- /dev/null +++ b/winsup/mingw/mingwex/math/lrint.c @@ -0,0 +1,9 @@ +#include <math.h> + +long lrint (double x) +{ + long retval; + __asm__ __volatile__ \ + ("fistpl %0" : "=m" (retval) : "t" (x) : "st"); \ + return retval; +} diff --git a/winsup/mingw/mingwex/math/lrintf.c b/winsup/mingw/mingwex/math/lrintf.c new file mode 100644 index 000000000..24b7a7d07 --- /dev/null +++ b/winsup/mingw/mingwex/math/lrintf.c @@ -0,0 +1,9 @@ +#include <math.h> + +long lrintf (float x) +{ + long retval; + __asm__ __volatile__ \ + ("fistpl %0" : "=m" (retval) : "t" (x) : "st"); \ + return retval; +} diff --git a/winsup/mingw/mingwex/math/lrintl.c b/winsup/mingw/mingwex/math/lrintl.c new file mode 100644 index 000000000..f55599332 --- /dev/null +++ b/winsup/mingw/mingwex/math/lrintl.c @@ -0,0 +1,10 @@ +#include <math.h> + +long lrintl (long double x) +{ + long retval; + __asm__ __volatile__ \ + ("fistpl %0" : "=m" (retval) : "t" (x) : "st"); \ + return retval; +} + diff --git a/winsup/mingw/mingwex/math/lround.c b/winsup/mingw/mingwex/math/lround.c new file mode 100644 index 000000000..4b0b0d047 --- /dev/null +++ b/winsup/mingw/mingwex/math/lround.c @@ -0,0 +1,24 @@ +#include <fenv.h> +#include <math.h> + +long +lround (double x) { + long retval; + unsigned short saved_cw, _cw; + __asm__ ( + "fnstcw %0;" : "=m" (saved_cw) + ); /* save control word */ + _cw = ~(FE_TONEAREST | FE_DOWNWARD | FE_UPWARD | FE_TOWARDZERO) + | (x > 0.0 ? FE_UPWARD : FE_DOWNWARD); /* round away from zero */ + __asm__ ( + "fldcw %0;" : : "m" (_cw) + ); /* load the rounding control */ + __asm__ __volatile__ ( + "fistpl %0" : "=m" (retval) : "t" (x) : "st" + ); + __asm__ ( + "fldcw %0;" : : "m" (saved_cw) + ); /* restore control word */ + return retval; +} + diff --git a/winsup/mingw/mingwex/math/lroundf.c b/winsup/mingw/mingwex/math/lroundf.c new file mode 100644 index 000000000..40051bfcb --- /dev/null +++ b/winsup/mingw/mingwex/math/lroundf.c @@ -0,0 +1,23 @@ +#include <fenv.h> +#include <math.h> + +long +lroundf (float x) { + long retval; + unsigned short saved_cw, _cw; + __asm__ ( + "fnstcw %0;" : "=m" (saved_cw) + ); /* save control word */ + _cw = ~(FE_TONEAREST | FE_DOWNWARD | FE_UPWARD | FE_TOWARDZERO) + | (x > 0.0 ? FE_UPWARD : FE_DOWNWARD); /* round away from zero */ + __asm__ ( + "fldcw %0;" : : "m" (_cw) + ); /* load the rounding control */ + __asm__ __volatile__ ( + "fistpl %0" : "=m" (retval) : "t" (x) : "st" + ); + __asm__ ( + "fldcw %0;" : : "m" (saved_cw) + ); /* restore control word */ + return retval; +} diff --git a/winsup/mingw/mingwex/math/lroundl.c b/winsup/mingw/mingwex/math/lroundl.c new file mode 100644 index 000000000..6e69c2fdd --- /dev/null +++ b/winsup/mingw/mingwex/math/lroundl.c @@ -0,0 +1,22 @@ +#include <fenv.h> +#include <math.h> + +long +lroundl (long double x) { + long retval; + unsigned short saved_cw, _cw; + __asm__ ( + "fnstcw %0;" : "=m" (saved_cw) + ); /* save control word */ + _cw = ~(FE_TONEAREST | FE_DOWNWARD | FE_UPWARD | FE_TOWARDZERO) + | (x > 0.0 ? FE_UPWARD : FE_DOWNWARD); /* round away from zero */ + __asm__ ( + "fldcw %0;" : : "m" (_cw) + ); /* load the rounding control */ + __asm__ __volatile__ ( + "fistpl %0" : "=m" (retval) : "t" (x) : "st"); + __asm__ ( + "fldcw %0;" : : "m" (saved_cw) + ); /* restore control word */ + return retval; +} diff --git a/winsup/mingw/mingwex/math/modff.c b/winsup/mingw/mingwex/math/modff.c new file mode 100644 index 000000000..9e34c7c95 --- /dev/null +++ b/winsup/mingw/mingwex/math/modff.c @@ -0,0 +1,21 @@ +#include <fenv.h> +#include <math.h> +#include <errno.h> +#define FE_ROUNDING_MASK \ + (FE_TONEAREST | FE_DOWNWARD | FE_UPWARD | FE_TOWARDZERO) + +float +modff (float value, float* iptr) +{ + float int_part; + unsigned short saved_cw; + /* truncate */ + asm ("fnstcw %0;" : "=m" (saved_cw)); /* save control word */ + asm ("fldcw %0;" : : "m" ((saved_cw & ~FE_ROUNDING_MASK) + | FE_TOWARDZERO)); + asm ("frndint;" : "=t" (int_part) : "0" (value)); /* round */ + asm ("fldcw %0;" : : "m" (saved_cw)); /* restore saved cw */ + if (iptr) + *iptr = int_part; + return (isinf (value) ? 0.0F : value - int_part); +} diff --git a/winsup/mingw/mingwex/math/modfl.c b/winsup/mingw/mingwex/math/modfl.c new file mode 100644 index 000000000..5663956f0 --- /dev/null +++ b/winsup/mingw/mingwex/math/modfl.c @@ -0,0 +1,21 @@ +#include <fenv.h> +#include <math.h> +#include <errno.h> +#define FE_ROUNDING_MASK \ + (FE_TONEAREST | FE_DOWNWARD | FE_UPWARD | FE_TOWARDZERO) + +long double +modfl (long double value, long double* iptr) +{ + long double int_part; + unsigned short saved_cw; + /* truncate */ + asm ("fnstcw %0;" : "=m" (saved_cw)); /* save control word */ + asm ("fldcw %0;" : : "m" ((saved_cw & ~FE_ROUNDING_MASK) + | FE_TOWARDZERO)); + asm ("frndint;" : "=t" (int_part) : "0" (value)); /* round */ + asm ("fldcw %0;" : : "m" (saved_cw)); /* restore saved cw */ + if (iptr) + *iptr = int_part; + return (isinf (value) ? 0.0L : value - int_part); +} diff --git a/winsup/mingw/mingwex/math/nearbyint.S b/winsup/mingw/mingwex/math/nearbyint.S new file mode 100644 index 000000000..9730aeebf --- /dev/null +++ b/winsup/mingw/mingwex/math/nearbyint.S @@ -0,0 +1,30 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * + * Adapted for use as nearbyint by Ulrich Drepper <drepper@cygnus.com>. + * + * Removed header file dependency for use in libmingwex.a by + * Danny Smith <dannysmith@users.sourceforge.net> + */ + + .file "nearbyint.S" + .text + .align 4 +.globl _nearbyint + .def _nearbyint; .scl 2; .type 32; .endef +_nearbyint: + fldl 4(%esp) + pushl %eax + pushl %ecx + fnstcw (%esp) + movl (%esp), %eax + orl $0x20, %eax + movl %eax, 4(%esp) + fldcw 4(%esp) + frndint + fclex + fldcw (%esp) + popl %ecx + popl %eax + ret diff --git a/winsup/mingw/mingwex/math/nearbyintf.S b/winsup/mingw/mingwex/math/nearbyintf.S new file mode 100644 index 000000000..1c5734084 --- /dev/null +++ b/winsup/mingw/mingwex/math/nearbyintf.S @@ -0,0 +1,29 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * Adapted for use as nearbyint by Ulrich Drepper <drepper@cygnus.com>. + * + * Removed header file dependency for use in libmingwex.a by + * Danny Smith <dannysmith@users.sourceforge.net> + */ + + .file "nearbyintf.S" + .text + .align 4 +.globl _nearbyintf + .def _nearbyintf; .scl 2; .type 32; .endef +_nearbyintf: + flds 4(%esp) + pushl %eax + pushl %ecx + fnstcw (%esp) + movl (%esp), %eax + orl $0x20, %eax + movl %eax, 4(%esp) + fldcw 4(%esp) + frndint + fclex + fldcw (%esp) + popl %ecx + popl %eax + ret diff --git a/winsup/mingw/mingwex/math/nearbyintl.S b/winsup/mingw/mingwex/math/nearbyintl.S new file mode 100644 index 000000000..7dbc2a8b7 --- /dev/null +++ b/winsup/mingw/mingwex/math/nearbyintl.S @@ -0,0 +1,30 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * + * Adaptedfor use as nearbyint by Ulrich Drepper <drepper@cygnus.com>. + * + * Removed header file dependency for use in libmingwex.a by + * Danny Smith <dannysmith@users.sourceforge.net> + */ + + .file "nearbyintl.S" + .text + .align 4 +.globl _nearbyintl + .def _nearbyintl; .scl 2; .type 32; .endef +_nearbyintl: + fldt 4(%esp) + pushl %eax + pushl %ecx + fnstcw (%esp) + movl (%esp), %eax + orl $0x20, %eax + movl %eax, 4(%esp) + fldcw 4(%esp) + frndint + fclex + fldcw (%esp) + popl %ecx + popl %eax + ret diff --git a/winsup/mingw/mingwex/math/nextafterf.c b/winsup/mingw/mingwex/math/nextafterf.c new file mode 100644 index 000000000..ac1fbee26 --- /dev/null +++ b/winsup/mingw/mingwex/math/nextafterf.c @@ -0,0 +1,29 @@ +#include <math.h> + +float +nextafterf (float x, float y) +{ + union + { + float f; + unsigned int i; + } u; + if (isnan (y) || isnan (x)) + return x + y; + if (x == y ) + return x; + u.f = x; + if (u.i == 0u) + { + if (y > 0.0F) + u.i = 1; + else + u.i = 0x80000001; + return u.f; + } + if (((x > 0.0F) ^ (y > x)) == 0) + u.i++; + else + u.i--; + return u.f; +} diff --git a/winsup/mingw/mingwex/math/powf.c b/winsup/mingw/mingwex/math/powf.c new file mode 100644 index 000000000..1af4d2d8f --- /dev/null +++ b/winsup/mingw/mingwex/math/powf.c @@ -0,0 +1,3 @@ +#include <math.h> +float powf (float x, float y) + {return (float) pow (x, y);} diff --git a/winsup/mingw/mingwex/math/powil.c b/winsup/mingw/mingwex/math/powil.c new file mode 100644 index 000000000..e900dd881 --- /dev/null +++ b/winsup/mingw/mingwex/math/powil.c @@ -0,0 +1,179 @@ +/* powil.c + * + * Real raised to integer power, long double precision + * + * + * + * SYNOPSIS: + * + * long double x, y, powil(); + * int n; + * + * y = powil( 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 .001,1000 -1022,1023 50000 4.3e-17 7.8e-18 + * IEEE 1,2 -1022,1023 20000 3.9e-17 7.6e-18 + * IEEE .99,1.01 0,8700 10000 3.6e-16 7.2e-17 + * + * Returns INFINITY on overflow, zero on underflow. + * + */ + +/* powil.c */ + +/* +Cephes Math Library Release 2.2: December, 1990 +Copyright 1984, 1990 by Stephen L. Moshier +Direct inquiries to 30 Frost Street, Cambridge, MA 02140 +*/ + +/* +Modified for mingw +2002-07-22 Danny Smith <dannysmith@users.sourceforge.net> +*/ + +#ifdef __MINGW32__ +#include "cephes_mconf.h" +#else +#include "mconf.h" +extern long double MAXNUML, MAXLOGL, MINLOGL; +extern long double LOGE2L; +#ifdef ANSIPROT +extern long double frexpl ( long double, int * ); +#else +long double frexpl(); +#endif +#endif /* __MINGW32__ */ + +#ifndef _SET_ERRNO +#define _SET_ERRNO(x) +#endif + +long double powil( x, nn ) +long double x; +int nn; +{ +long double w, y; +long double s; +int n, e, sign, asign, lx; + +if( x == 0.0L ) + { + if( nn == 0 ) + return( 1.0L ); + else if( nn < 0 ) + return( INFINITYL ); + else + return( 0.0L ); + } + +if( nn == 0 ) + return( 1.0L ); + + +if( x < 0.0L ) + { + asign = -1; + x = -x; + } +else + asign = 0; + + +if( nn < 0 ) + { + sign = -1; + n = -nn; + } +else + { + sign = 1; + n = nn; + } + +/* Overflow detection */ + +/* Calculate approximate logarithm of answer */ +s = x; +s = frexpl( s, &lx ); +e = (lx - 1)*n; +if( (e == 0) || (e > 64) || (e < -64) ) + { + s = (s - 7.0710678118654752e-1L) / (s + 7.0710678118654752e-1L); + s = (2.9142135623730950L * s - 0.5L + lx) * nn * LOGE2L; + } +else + { + s = LOGE2L * e; + } + +if( s > MAXLOGL ) + { + mtherr( "powil", OVERFLOW ); + _SET_ERRNO(ERANGE); + y = INFINITYL; + goto done; + } + +if( s < MINLOGL ) + { + mtherr( "powil", UNDERFLOW ); + _SET_ERRNO(ERANGE); + return(0.0L); + } +/* 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 < (-MAXLOGL+2.0L) ) + { + x = 1.0L/x; + sign = -sign; + } + +/* First bit of the power */ +if( n & 1 ) + y = x; + +else + { + y = 1.0L; + 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 < 0 ) + y = 1.0L/y; +return(y); +} diff --git a/winsup/mingw/mingwex/math/powl.c b/winsup/mingw/mingwex/math/powl.c new file mode 100644 index 000000000..a94ede965 --- /dev/null +++ b/winsup/mingw/mingwex/math/powl.c @@ -0,0 +1,803 @@ +/* powl.c + * + * Power function, long double precision + * + * + * + * SYNOPSIS: + * + * long double x, y, z, powl(); + * + * z = powl( x, y ); + * + * + * + * DESCRIPTION: + * + * Computes x raised to the yth power. Analytically, + * + * x**y = exp( y log(x) ). + * + * Following Cody and Waite, this program uses a lookup table + * of 2**-i/32 and pseudo extended precision arithmetic to + * obtain several extra bits of accuracy in both the logarithm + * and the exponential. + * + * + * + * ACCURACY: + * + * The relative error of pow(x,y) can be estimated + * by y dl ln(2), where dl is the absolute error of + * the internally computed base 2 logarithm. At the ends + * of the approximation interval the logarithm equal 1/32 + * and its relative error is about 1 lsb = 1.1e-19. Hence + * the predicted relative error in the result is 2.3e-21 y . + * + * Relative error: + * arithmetic domain # trials peak rms + * + * IEEE +-1000 40000 2.8e-18 3.7e-19 + * .001 < x < 1000, with log(x) uniformly distributed. + * -1000 < y < 1000, y uniformly distributed. + * + * IEEE 0,8700 60000 6.5e-18 1.0e-18 + * 0.99 < x < 1.01, 0 < y < 8700, uniformly distributed. + * + * + * ERROR MESSAGES: + * + * message condition value returned + * pow overflow x**y > MAXNUM INFINITY + * pow underflow x**y < 1/MAXNUM 0.0 + * pow domain x<0 and y noninteger 0.0 + * + */ + +/* +Cephes Math Library Release 2.7: May, 1998 +Copyright 1984, 1991, 1998 by Stephen L. Moshier +*/ + +/* +Modified for mingw +2002-07-22 Danny Smith <dannysmith@users.sourceforge.net> +*/ + +#ifdef __MINGW32__ +#include "cephes_mconf.h" +#else +#include "mconf.h" + +static char fname[] = {"powl"}; +#endif + +#ifndef _SET_ERRNO +#define _SET_ERRNO(x) +#endif + + +/* Table size */ +#define NXT 32 +/* log2(Table size) */ +#define LNXT 5 + +#ifdef UNK +/* log(1+x) = x - .5x^2 + x^3 * P(z)/Q(z) + * on the domain 2^(-1/32) - 1 <= x <= 2^(1/32) - 1 + */ +static long double P[] = { + 8.3319510773868690346226E-4L, + 4.9000050881978028599627E-1L, + 1.7500123722550302671919E0L, + 1.4000100839971580279335E0L, +}; +static long double Q[] = { +/* 1.0000000000000000000000E0L,*/ + 5.2500282295834889175431E0L, + 8.4000598057587009834666E0L, + 4.2000302519914740834728E0L, +}; +/* A[i] = 2^(-i/32), rounded to IEEE long double precision. + * If i is even, A[i] + B[i/2] gives additional accuracy. + */ +static long double A[33] = { + 1.0000000000000000000000E0L, + 9.7857206208770013448287E-1L, + 9.5760328069857364691013E-1L, + 9.3708381705514995065011E-1L, + 9.1700404320467123175367E-1L, + 8.9735453750155359320742E-1L, + 8.7812608018664974155474E-1L, + 8.5930964906123895780165E-1L, + 8.4089641525371454301892E-1L, + 8.2287773907698242225554E-1L, + 8.0524516597462715409607E-1L, + 7.8799042255394324325455E-1L, + 7.7110541270397041179298E-1L, + 7.5458221379671136985669E-1L, + 7.3841307296974965571198E-1L, + 7.2259040348852331001267E-1L, + 7.0710678118654752438189E-1L, + 6.9195494098191597746178E-1L, + 6.7712777346844636413344E-1L, + 6.6261832157987064729696E-1L, + 6.4841977732550483296079E-1L, + 6.3452547859586661129850E-1L, + 6.2092890603674202431705E-1L, + 6.0762367999023443907803E-1L, + 5.9460355750136053334378E-1L, + 5.8186242938878875689693E-1L, + 5.6939431737834582684856E-1L, + 5.5719337129794626814472E-1L, + 5.4525386633262882960438E-1L, + 5.3357020033841180906486E-1L, + 5.2213689121370692017331E-1L, + 5.1094857432705833910408E-1L, + 5.0000000000000000000000E-1L, +}; +static long double B[17] = { + 0.0000000000000000000000E0L, + 2.6176170809902549338711E-20L, +-1.0126791927256478897086E-20L, + 1.3438228172316276937655E-21L, + 1.2207982955417546912101E-20L, +-6.3084814358060867200133E-21L, + 1.3164426894366316434230E-20L, +-1.8527916071632873716786E-20L, + 1.8950325588932570796551E-20L, + 1.5564775779538780478155E-20L, + 6.0859793637556860974380E-21L, +-2.0208749253662532228949E-20L, + 1.4966292219224761844552E-20L, + 3.3540909728056476875639E-21L, +-8.6987564101742849540743E-22L, +-1.2327176863327626135542E-20L, + 0.0000000000000000000000E0L, +}; + +/* 2^x = 1 + x P(x), + * on the interval -1/32 <= x <= 0 + */ +static long double R[] = { + 1.5089970579127659901157E-5L, + 1.5402715328927013076125E-4L, + 1.3333556028915671091390E-3L, + 9.6181291046036762031786E-3L, + 5.5504108664798463044015E-2L, + 2.4022650695910062854352E-1L, + 6.9314718055994530931447E-1L, +}; + +#define douba(k) A[k] +#define doubb(k) B[k] +#define MEXP (NXT*16384.0L) +/* The following if denormal numbers are supported, else -MEXP: */ +#ifdef DENORMAL +#define MNEXP (-NXT*(16384.0L+64.0L)) +#else +#define MNEXP (-NXT*16384.0L) +#endif +/* log2(e) - 1 */ +#define LOG2EA 0.44269504088896340735992L +#endif + + +#ifdef IBMPC +static const unsigned short P[] = { +0xb804,0xa8b7,0xc6f4,0xda6a,0x3ff4, XPD +0x7de9,0xcf02,0x58c0,0xfae1,0x3ffd, XPD +0x405a,0x3722,0x67c9,0xe000,0x3fff, XPD +0xcd99,0x6b43,0x87ca,0xb333,0x3fff, XPD +}; +static const unsigned short Q[] = { +/* 0x0000,0x0000,0x0000,0x8000,0x3fff, */ +0x6307,0xa469,0x3b33,0xa800,0x4001, XPD +0xfec2,0x62d7,0xa51c,0x8666,0x4002, XPD +0xda32,0xd072,0xa5d7,0x8666,0x4001, XPD +}; +static const unsigned short A[] = { +0x0000,0x0000,0x0000,0x8000,0x3fff, XPD +0x033a,0x722a,0xb2db,0xfa83,0x3ffe, XPD +0xcc2c,0x2486,0x7d15,0xf525,0x3ffe, XPD +0xf5cb,0xdcda,0xb99b,0xefe4,0x3ffe, XPD +0x392f,0xdd24,0xc6e7,0xeac0,0x3ffe, XPD +0x48a8,0x7c83,0x06e7,0xe5b9,0x3ffe, XPD +0xe111,0x2a94,0xdeec,0xe0cc,0x3ffe, XPD +0x3755,0xdaf2,0xb797,0xdbfb,0x3ffe, XPD +0x6af4,0xd69d,0xfcca,0xd744,0x3ffe, XPD +0xe45a,0xf12a,0x1d91,0xd2a8,0x3ffe, XPD +0x80e4,0x1f84,0x8c15,0xce24,0x3ffe, XPD +0x27a3,0x6e2f,0xbd86,0xc9b9,0x3ffe, XPD +0xdadd,0x5506,0x2a11,0xc567,0x3ffe, XPD +0x9456,0x6670,0x4cca,0xc12c,0x3ffe, XPD +0x36bf,0x580c,0xa39f,0xbd08,0x3ffe, XPD +0x9ee9,0x62fb,0xaf47,0xb8fb,0x3ffe, XPD +0x6484,0xf9de,0xf333,0xb504,0x3ffe, XPD +0x2590,0xd2ac,0xf581,0xb123,0x3ffe, XPD +0x4ac6,0x42a1,0x3eea,0xad58,0x3ffe, XPD +0x0ef8,0xea7c,0x5ab4,0xa9a1,0x3ffe, XPD +0x38ea,0xb151,0xd6a9,0xa5fe,0x3ffe, XPD +0x6819,0x0c49,0x4303,0xa270,0x3ffe, XPD +0x11ae,0x91a1,0x3260,0x9ef5,0x3ffe, XPD +0x5539,0xd54e,0x39b9,0x9b8d,0x3ffe, XPD +0xa96f,0x8db8,0xf051,0x9837,0x3ffe, XPD +0x0961,0xfef7,0xefa8,0x94f4,0x3ffe, XPD +0xc336,0xab11,0xd373,0x91c3,0x3ffe, XPD +0x53c0,0x45cd,0x398b,0x8ea4,0x3ffe, XPD +0xd6e7,0xea8b,0xc1e3,0x8b95,0x3ffe, XPD +0x8527,0x92da,0x0e80,0x8898,0x3ffe, XPD +0x7b15,0xcc48,0xc367,0x85aa,0x3ffe, XPD +0xa1d7,0xac2b,0x8698,0x82cd,0x3ffe, XPD +0x0000,0x0000,0x0000,0x8000,0x3ffe, XPD +}; +static const unsigned short B[] = { +0x0000,0x0000,0x0000,0x0000,0x0000, XPD +0x1f87,0xdb30,0x18f5,0xf73a,0x3fbd, XPD +0xac15,0x3e46,0x2932,0xbf4a,0xbfbc, XPD +0x7944,0xba66,0xa091,0xcb12,0x3fb9, XPD +0xff78,0x40b4,0x2ee6,0xe69a,0x3fbc, XPD +0xc895,0x5069,0xe383,0xee53,0xbfbb, XPD +0x7cde,0x9376,0x4325,0xf8ab,0x3fbc, XPD +0xa10c,0x25e0,0xc093,0xaefd,0xbfbd, XPD +0x7d3e,0xea95,0x1366,0xb2fb,0x3fbd, XPD +0x5d89,0xeb34,0x5191,0x9301,0x3fbd, XPD +0x80d9,0xb883,0xfb10,0xe5eb,0x3fbb, XPD +0x045d,0x288c,0xc1ec,0xbedd,0xbfbd, XPD +0xeded,0x5c85,0x4630,0x8d5a,0x3fbd, XPD +0x9d82,0xe5ac,0x8e0a,0xfd6d,0x3fba, XPD +0x6dfd,0xeb58,0xaf14,0x8373,0xbfb9, XPD +0xf938,0x7aac,0x91cf,0xe8da,0xbfbc, XPD +0x0000,0x0000,0x0000,0x0000,0x0000, XPD +}; +static const unsigned short R[] = { +0xa69b,0x530e,0xee1d,0xfd2a,0x3fee, XPD +0xc746,0x8e7e,0x5960,0xa182,0x3ff2, XPD +0x63b6,0xadda,0xfd6a,0xaec3,0x3ff5, XPD +0xc104,0xfd99,0x5b7c,0x9d95,0x3ff8, XPD +0xe05e,0x249d,0x46b8,0xe358,0x3ffa, XPD +0x5d1d,0x162c,0xeffc,0xf5fd,0x3ffc, XPD +0x79aa,0xd1cf,0x17f7,0xb172,0x3ffe, XPD +}; + +/* 10 byte sizes versus 12 byte */ +#define douba(k) (*(long double *)(&A[(sizeof( long double )/2)*(k)])) +#define doubb(k) (*(long double *)(&B[(sizeof( long double )/2)*(k)])) +#define MEXP (NXT*16384.0L) +#ifdef DENORMAL +#define MNEXP (-NXT*(16384.0L+64.0L)) +#else +#define MNEXP (-NXT*16384.0L) +#endif +static const unsigned short L[] = {0xc2ef,0x705f,0xeca5,0xe2a8,0x3ffd, XPD}; +#define LOG2EA (*(long double *)(&L[0])) +#endif + +#ifdef MIEEE +static long P[] = { +0x3ff40000,0xda6ac6f4,0xa8b7b804, +0x3ffd0000,0xfae158c0,0xcf027de9, +0x3fff0000,0xe00067c9,0x3722405a, +0x3fff0000,0xb33387ca,0x6b43cd99, +}; +static long Q[] = { +/* 0x3fff0000,0x80000000,0x00000000, */ +0x40010000,0xa8003b33,0xa4696307, +0x40020000,0x8666a51c,0x62d7fec2, +0x40010000,0x8666a5d7,0xd072da32, +}; +static long A[] = { +0x3fff0000,0x80000000,0x00000000, +0x3ffe0000,0xfa83b2db,0x722a033a, +0x3ffe0000,0xf5257d15,0x2486cc2c, +0x3ffe0000,0xefe4b99b,0xdcdaf5cb, +0x3ffe0000,0xeac0c6e7,0xdd24392f, +0x3ffe0000,0xe5b906e7,0x7c8348a8, +0x3ffe0000,0xe0ccdeec,0x2a94e111, +0x3ffe0000,0xdbfbb797,0xdaf23755, +0x3ffe0000,0xd744fcca,0xd69d6af4, +0x3ffe0000,0xd2a81d91,0xf12ae45a, +0x3ffe0000,0xce248c15,0x1f8480e4, +0x3ffe0000,0xc9b9bd86,0x6e2f27a3, +0x3ffe0000,0xc5672a11,0x5506dadd, +0x3ffe0000,0xc12c4cca,0x66709456, +0x3ffe0000,0xbd08a39f,0x580c36bf, +0x3ffe0000,0xb8fbaf47,0x62fb9ee9, +0x3ffe0000,0xb504f333,0xf9de6484, +0x3ffe0000,0xb123f581,0xd2ac2590, +0x3ffe0000,0xad583eea,0x42a14ac6, +0x3ffe0000,0xa9a15ab4,0xea7c0ef8, +0x3ffe0000,0xa5fed6a9,0xb15138ea, +0x3ffe0000,0xa2704303,0x0c496819, +0x3ffe0000,0x9ef53260,0x91a111ae, +0x3ffe0000,0x9b8d39b9,0xd54e5539, +0x3ffe0000,0x9837f051,0x8db8a96f, +0x3ffe0000,0x94f4efa8,0xfef70961, +0x3ffe0000,0x91c3d373,0xab11c336, +0x3ffe0000,0x8ea4398b,0x45cd53c0, +0x3ffe0000,0x8b95c1e3,0xea8bd6e7, +0x3ffe0000,0x88980e80,0x92da8527, +0x3ffe0000,0x85aac367,0xcc487b15, +0x3ffe0000,0x82cd8698,0xac2ba1d7, +0x3ffe0000,0x80000000,0x00000000, +}; +static long B[51] = { +0x00000000,0x00000000,0x00000000, +0x3fbd0000,0xf73a18f5,0xdb301f87, +0xbfbc0000,0xbf4a2932,0x3e46ac15, +0x3fb90000,0xcb12a091,0xba667944, +0x3fbc0000,0xe69a2ee6,0x40b4ff78, +0xbfbb0000,0xee53e383,0x5069c895, +0x3fbc0000,0xf8ab4325,0x93767cde, +0xbfbd0000,0xaefdc093,0x25e0a10c, +0x3fbd0000,0xb2fb1366,0xea957d3e, +0x3fbd0000,0x93015191,0xeb345d89, +0x3fbb0000,0xe5ebfb10,0xb88380d9, +0xbfbd0000,0xbeddc1ec,0x288c045d, +0x3fbd0000,0x8d5a4630,0x5c85eded, +0x3fba0000,0xfd6d8e0a,0xe5ac9d82, +0xbfb90000,0x8373af14,0xeb586dfd, +0xbfbc0000,0xe8da91cf,0x7aacf938, +0x00000000,0x00000000,0x00000000, +}; +static long R[] = { +0x3fee0000,0xfd2aee1d,0x530ea69b, +0x3ff20000,0xa1825960,0x8e7ec746, +0x3ff50000,0xaec3fd6a,0xadda63b6, +0x3ff80000,0x9d955b7c,0xfd99c104, +0x3ffa0000,0xe35846b8,0x249de05e, +0x3ffc0000,0xf5fdeffc,0x162c5d1d, +0x3ffe0000,0xb17217f7,0xd1cf79aa, +}; + +#define douba(k) (*(long double *)&A[3*(k)]) +#define doubb(k) (*(long double *)&B[3*(k)]) +#define MEXP (NXT*16384.0L) +#ifdef DENORMAL +#define MNEXP (-NXT*(16384.0L+64.0L)) +#else +#define MNEXP (-NXT*16382.0L) +#endif +static long L[3] = {0x3ffd0000,0xe2a8eca5,0x705fc2ef}; +#define LOG2EA (*(long double *)(&L[0])) +#endif + + +#define F W +#define Fa Wa +#define Fb Wb +#define G W +#define Ga Wa +#define Gb u +#define H W +#define Ha Wb +#define Hb Wb + +#ifndef __MINGW32__ +extern long double MAXNUML; +#endif + +static VOLATILE long double z; +static long double w, W, Wa, Wb, ya, yb, u; + +#ifdef __MINGW32__ +static __inline__ long double reducl( long double ); +extern long double powil ( long double, int ); +extern long double powl ( long double x, long double y); +#else +#ifdef ANSIPROT +extern long double floorl ( long double ); +extern long double fabsl ( long double ); +extern long double frexpl ( long double, int * ); +extern long double ldexpl ( long double, int ); +extern long double polevll ( long double, void *, int ); +extern long double p1evll ( long double, void *, int ); +extern long double powil ( long double, int ); +extern int isnanl ( long double ); +extern int isfinitel ( long double ); +static long double reducl( long double ); +extern int signbitl ( long double ); +#else +long double floorl(), fabsl(), frexpl(), ldexpl(); +long double polevll(), p1evll(), powil(); +static long double reducl(); +int isnanl(), isfinitel(), signbitl(); +#endif + +#ifdef INFINITIES +extern long double INFINITYL; +#else +#define INFINITYL MAXNUML +#endif + +#ifdef NANS +extern long double NANL; +#endif +#ifdef MINUSZERO +extern long double NEGZEROL; +#endif + +#endif /* __MINGW32__ */ + +long double powl( x, y ) +long double x, y; +{ +/* double F, Fa, Fb, G, Ga, Gb, H, Ha, Hb */ +int i, nflg, iyflg, yoddint; +long e; + +if( y == 0.0L ) + return( 1.0L ); + +#ifdef NANS +if( isnanl(x) ) + { + _SET_ERRNO (EDOM); + return( x ); + } +if( isnanl(y) ) + { + _SET_ERRNO (EDOM); + return( y ); + } +#endif + +if( y == 1.0L ) + return( x ); + +if( isinfl(y) && (x == -1.0L || x == 1.0L) ) + return( y ); + +if( x == 1.0L ) + return( 1.0L ); + +if( y >= MAXNUML ) + { + _SET_ERRNO (ERANGE); +#ifdef INFINITIES + if( x > 1.0L ) + return( INFINITYL ); +#else + if( x > 1.0L ) + return( MAXNUML ); +#endif + if( x > 0.0L && x < 1.0L ) + return( 0.0L ); +#ifdef INFINITIES + if( x < -1.0L ) + return( INFINITYL ); +#else + if( x < -1.0L ) + return( MAXNUML ); +#endif + if( x > -1.0L && x < 0.0L ) + return( 0.0L ); + } +if( y <= -MAXNUML ) + { + _SET_ERRNO (ERANGE); + if( x > 1.0L ) + return( 0.0L ); +#ifdef INFINITIES + if( x > 0.0L && x < 1.0L ) + return( INFINITYL ); +#else + if( x > 0.0L && x < 1.0L ) + return( MAXNUML ); +#endif + if( x < -1.0L ) + return( 0.0L ); +#ifdef INFINITIES + if( x > -1.0L && x < 0.0L ) + return( INFINITYL ); +#else + if( x > -1.0L && x < 0.0L ) + return( MAXNUML ); +#endif + } +if( x >= MAXNUML ) + { +#if INFINITIES + if( y > 0.0L ) + return( INFINITYL ); +#else + if( y > 0.0L ) + return( MAXNUML ); +#endif + return( 0.0L ); + } + +w = floorl(y); +/* Set iyflg to 1 if y is an integer. */ +iyflg = 0; +if( w == y ) + iyflg = 1; + +/* Test for odd integer y. */ +yoddint = 0; +if( iyflg ) + { + ya = fabsl(y); + ya = floorl(0.5L * ya); + yb = 0.5L * fabsl(w); + if( ya != yb ) + yoddint = 1; + } + +if( x <= -MAXNUML ) + { + if( y > 0.0L ) + { +#ifdef INFINITIES + if( yoddint ) + return( -INFINITYL ); + return( INFINITYL ); +#else + if( yoddint ) + return( -MAXNUML ); + return( MAXNUML ); +#endif + } + if( y < 0.0L ) + { +#ifdef MINUSZERO + if( yoddint ) + return( NEGZEROL ); +#endif + return( 0.0 ); + } + } + + +nflg = 0; /* flag = 1 if x<0 raised to integer power */ +if( x <= 0.0L ) + { + if( x == 0.0L ) + { + if( y < 0.0 ) + { +#ifdef MINUSZERO + if( signbitl(x) && yoddint ) + return( -INFINITYL ); +#endif +#ifdef INFINITIES + return( INFINITYL ); +#else + return( MAXNUML ); +#endif + } + if( y > 0.0 ) + { +#ifdef MINUSZERO + if( signbitl(x) && yoddint ) + return( NEGZEROL ); +#endif + return( 0.0 ); + } + if( y == 0.0L ) + return( 1.0L ); /* 0**0 */ + else + return( 0.0L ); /* 0**y */ + } + else + { + if( iyflg == 0 ) + { /* noninteger power of negative number */ + mtherr( fname, DOMAIN ); + _SET_ERRNO (EDOM); +#ifdef NANS + return(NANL); +#else + return(0.0L); +#endif + } + nflg = 1; + } + } + +/* Integer power of an integer. */ + +if( iyflg ) + { + i = w; + w = floorl(x); + if( (w == x) && (fabsl(y) < 32768.0) ) + { + w = powil( x, (int) y ); + return( w ); + } + } + + +if( nflg ) + x = fabsl(x); + +/* separate significand from exponent */ +x = frexpl( x, &i ); +e = i; + +/* find significand in antilog table A[] */ +i = 1; +if( x <= douba(17) ) + i = 17; +if( x <= douba(i+8) ) + i += 8; +if( x <= douba(i+4) ) + i += 4; +if( x <= douba(i+2) ) + i += 2; +if( x >= douba(1) ) + i = -1; +i += 1; + + +/* Find (x - A[i])/A[i] + * in order to compute log(x/A[i]): + * + * log(x) = log( a x/a ) = log(a) + log(x/a) + * + * log(x/a) = log(1+v), v = x/a - 1 = (x-a)/a + */ +x -= douba(i); +x -= doubb(i/2); +x /= douba(i); + + +/* rational approximation for log(1+v): + * + * log(1+v) = v - v**2/2 + v**3 P(v) / Q(v) + */ +z = x*x; +w = x * ( z * polevll( x, P, 3 ) / p1evll( x, Q, 3 ) ); +w = w - ldexpl( z, -1 ); /* w - 0.5 * z */ + +/* Convert to base 2 logarithm: + * multiply by log2(e) = 1 + LOG2EA + */ +z = LOG2EA * w; +z += w; +z += LOG2EA * x; +z += x; + +/* Compute exponent term of the base 2 logarithm. */ +w = -i; +w = ldexpl( w, -LNXT ); /* divide by NXT */ +w += e; +/* Now base 2 log of x is w + z. */ + +/* Multiply base 2 log by y, in extended precision. */ + +/* separate y into large part ya + * and small part yb less than 1/NXT + */ +ya = reducl(y); +yb = y - ya; + +/* (w+z)(ya+yb) + * = w*ya + w*yb + z*y + */ +F = z * y + w * yb; +Fa = reducl(F); +Fb = F - Fa; + +G = Fa + w * ya; +Ga = reducl(G); +Gb = G - Ga; + +H = Fb + Gb; +Ha = reducl(H); +w = ldexpl( Ga+Ha, LNXT ); + +/* Test the power of 2 for overflow */ +if( w > MEXP ) + { +/* printf( "w = %.4Le ", w ); */ + _SET_ERRNO (ERANGE); + mtherr( fname, OVERFLOW ); + return( MAXNUML ); + } + +if( w < MNEXP ) + { +/* printf( "w = %.4Le ", w ); */ + _SET_ERRNO (ERANGE); + mtherr( fname, UNDERFLOW ); + return( 0.0L ); + } + +e = w; +Hb = H - Ha; + +if( Hb > 0.0L ) + { + e += 1; + Hb -= (1.0L/NXT); /*0.0625L;*/ + } + +/* Now the product y * log2(x) = Hb + e/NXT. + * + * Compute base 2 exponential of Hb, + * where -0.0625 <= Hb <= 0. + */ +z = Hb * polevll( Hb, R, 6 ); /* z = 2**Hb - 1 */ + +/* Express e/NXT as an integer plus a negative number of (1/NXT)ths. + * Find lookup table entry for the fractional power of 2. + */ +if( e < 0 ) + i = 0; +else + i = 1; +i = e/NXT + i; +e = NXT*i - e; +w = douba( e ); +z = w * z; /* 2**-e * ( 1 + (2**Hb-1) ) */ +z = z + w; +z = ldexpl( z, i ); /* multiply by integer power of 2 */ + +if( nflg ) + { +/* For negative x, + * find out if the integer exponent + * is odd or even. + */ + w = ldexpl( y, -1 ); + w = floorl(w); + w = ldexpl( w, 1 ); + if( w != y ) + z = -z; /* odd exponent */ + } + +return( z ); +} + + +/* Find a multiple of 1/NXT that is within 1/NXT of x. */ +static __inline__ long double reducl(x) +long double x; +{ +long double t; + +t = ldexpl( x, LNXT ); +t = floorl( t ); +t = ldexpl( t, -LNXT ); +return(t); +} + +/* powil.c + * + * Real raised to integer power, long double precision + * + * + * + * SYNOPSIS: + * + * long double x, y, powil(); + * int n; + * + * y = powil( 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 .001,1000 -1022,1023 50000 4.3e-17 7.8e-18 + * IEEE 1,2 -1022,1023 20000 3.9e-17 7.6e-18 + * IEEE .99,1.01 0,8700 10000 3.6e-16 7.2e-17 + * + * Returns INFINITY on overflow, zero on underflow. + * + */ diff --git a/winsup/mingw/mingwex/math/remainder.S b/winsup/mingw/mingwex/math/remainder.S new file mode 100644 index 000000000..01930d3ba --- /dev/null +++ b/winsup/mingw/mingwex/math/remainder.S @@ -0,0 +1,19 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + + .file "remainder.S" + .text + .align 4 +.globl _remainder + .def _remainder; .scl 2; .type 32; .endef +_remainder: + fldl 12(%esp) + fldl 4(%esp) +1: fprem1 + fstsw %ax + sahf + jp 1b + fstp %st(1) + ret diff --git a/winsup/mingw/mingwex/math/remainderf.S b/winsup/mingw/mingwex/math/remainderf.S new file mode 100644 index 000000000..81e78415a --- /dev/null +++ b/winsup/mingw/mingwex/math/remainderf.S @@ -0,0 +1,19 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + + .file "remainderf.S" + .text + .align 4 +.globl _remainder + .def _remainderf; .scl 2; .type 32; .endef +_remainderf: + flds 8(%esp) + flds 4(%esp) +1: fprem1 + fstsw %ax + sahf + jp 1b + fstp %st(1) + ret diff --git a/winsup/mingw/mingwex/math/remainderl.S b/winsup/mingw/mingwex/math/remainderl.S new file mode 100644 index 000000000..b5ce3736d --- /dev/null +++ b/winsup/mingw/mingwex/math/remainderl.S @@ -0,0 +1,22 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>. + * Removed header file dependency for use in libmingwex.a by + * Danny Smith <dannysmith@users.sourceforge.net> + */ + + .file "remainderl.S" + .text + .align 4 +.globl _remainderl + .def _remainderl; .scl 2; .type 32; .endef +_remainderl: + fldt 16(%esp) + fldt 4(%esp) +1: fprem1 + fstsw %ax + sahf + jp 1b + fstp %st(1) + ret diff --git a/winsup/mingw/mingwex/math/remquo.S b/winsup/mingw/mingwex/math/remquo.S new file mode 100644 index 000000000..987c37ca5 --- /dev/null +++ b/winsup/mingw/mingwex/math/remquo.S @@ -0,0 +1,38 @@ +/* + * Written by Ulrich Drepper <drepper@cygnus.com>. + * Based on e_remainder by J.T. Conklin <jtc@netbsd.org>. + * Removed header file dependency for use in libmingwex.a by + * Danny Smith <dannysmith@users.sourceforge.ne + * Public domain. + */ + + .file "remquo.S" + .text + .align 4; +.globl _remquo; +_remquo: + fldl 4 +8(%esp) + fldl 4(%esp) +1: fprem1 + fstsw %ax + sahf + jp 1b + fstp %st(1) + movl %eax, %ecx + shrl $8, %eax + shrl $12, %ecx + andl $4, %ecx + andl $3, %eax + orl %eax, %ecx + movl $0xef2960, %eax + shrl %cl, %eax + andl $3, %eax + movl 4 +8 +8(%esp), %ecx + movl 4 +4(%esp), %edx + xorl 4 +8 +4(%esp), %edx + testl $0x80000000, %edx + jz 1f + negl %eax +1: movl %eax, (%ecx) + + ret diff --git a/winsup/mingw/mingwex/math/remquof.S b/winsup/mingw/mingwex/math/remquof.S new file mode 100644 index 000000000..af540ef5b --- /dev/null +++ b/winsup/mingw/mingwex/math/remquof.S @@ -0,0 +1,38 @@ +/* + * Written by Ulrich Drepper <drepper@cygnus.com>. + * Based on e_remainder by J.T. Conklin <jtc@netbsd.org>. + * Removed header file dependency for use in libmingwex.a by + * Danny Smith <dannysmith@users.sourceforge.ne + * Public domain. + */ + + .file "remquo.S" + .text + .align 4; +.globl _remquof; +_remquof: + flds 4 +4(%esp) + flds 4(%esp) +1: fprem1 + fstsw %ax + sahf + jp 1b + fstp %st(1) + movl %eax, %ecx + shrl $8, %eax + shrl $12, %ecx + andl $4, %ecx + andl $3, %eax + orl %eax, %ecx + movl $0xef2960, %eax + shrl %cl, %eax + andl $3, %eax + movl 4 +4 +4(%esp), %ecx + movl 4(%esp), %edx + xorl 4 +4(%esp), %edx + testl $0x80000000, %edx + jz 1f + negl %eax +1: movl %eax, (%ecx) + + ret diff --git a/winsup/mingw/mingwex/math/remquol.S b/winsup/mingw/mingwex/math/remquol.S new file mode 100644 index 000000000..e6f1b5420 --- /dev/null +++ b/winsup/mingw/mingwex/math/remquol.S @@ -0,0 +1,36 @@ +/* + * Written by Ulrich Drepper <drepper@cygnus.com>. + * Based on e_remainder by J.T. Conklin <jtc@netbsd.org>. + * Removed header file dependency for use in libmingwex.a by + * Danny Smith <dannysmith@users.sourceforge.net> + * Public domain. + */ + .text + .align 4; +.globl _remquol; + _remquol: + fldt 4 +12(%esp) + fldt 4(%esp) +1: fprem1 + fstsw %ax + sahf + jp 1b + fstp %st(1) + movl %eax, %ecx + shrl $8, %eax + shrl $12, %ecx + andl $4, %ecx + andl $3, %eax + orl %eax, %ecx + movl $0xef2960, %eax + shrl %cl, %eax + andl $3, %eax + movl 4 +12 +12(%esp), %ecx + movl 4 +8(%esp), %edx + xorl 4 +12 +8(%esp), %edx + testl $0x8000, %edx + jz 1f + negl %eax +1: movl %eax, (%ecx) + + ret diff --git a/winsup/mingw/mingwex/rint.c b/winsup/mingw/mingwex/math/rint.c index 3198f4b26..3198f4b26 100644 --- a/winsup/mingw/mingwex/rint.c +++ b/winsup/mingw/mingwex/math/rint.c diff --git a/winsup/mingw/mingwex/rintf.c b/winsup/mingw/mingwex/math/rintf.c index 0b05e8f89..0b05e8f89 100644 --- a/winsup/mingw/mingwex/rintf.c +++ b/winsup/mingw/mingwex/math/rintf.c diff --git a/winsup/mingw/mingwex/rintl.c b/winsup/mingw/mingwex/math/rintl.c index ffc9d1107..ffc9d1107 100644 --- a/winsup/mingw/mingwex/rintl.c +++ b/winsup/mingw/mingwex/math/rintl.c diff --git a/winsup/mingw/mingwex/round.c b/winsup/mingw/mingwex/math/round.c index 9d8e949e4..9d8e949e4 100644 --- a/winsup/mingw/mingwex/round.c +++ b/winsup/mingw/mingwex/math/round.c diff --git a/winsup/mingw/mingwex/roundf.c b/winsup/mingw/mingwex/math/roundf.c index 6ae81bdd8..6ae81bdd8 100644 --- a/winsup/mingw/mingwex/roundf.c +++ b/winsup/mingw/mingwex/math/roundf.c diff --git a/winsup/mingw/mingwex/roundl.c b/winsup/mingw/mingwex/math/roundl.c index de3334a62..de3334a62 100644 --- a/winsup/mingw/mingwex/roundl.c +++ b/winsup/mingw/mingwex/math/roundl.c diff --git a/winsup/mingw/mingwex/math/scalbn.S b/winsup/mingw/mingwex/math/scalbn.S new file mode 100644 index 000000000..76e2d396e --- /dev/null +++ b/winsup/mingw/mingwex/math/scalbn.S @@ -0,0 +1,19 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + + .file "scalbn.S" + .text + .align 4 +.globl _scalbn + .def _scalbn; .scl 2; .type 32; .endef +_scalbn: + fildl 12(%esp) + fldl 4(%esp) + fscale + fstp %st(1) + ret + +.globl _scalbln + .set _scalbln,_scalbn diff --git a/winsup/mingw/mingwex/math/scalbnf.S b/winsup/mingw/mingwex/math/scalbnf.S new file mode 100644 index 000000000..1fe42a3de --- /dev/null +++ b/winsup/mingw/mingwex/math/scalbnf.S @@ -0,0 +1,19 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + */ + + .file "scalbnf.S" + .text + .align 4 +.globl _scalbnf + .def _scalbnf; .scl 2; .type 32; .endef +_scalbnf: + fildl 8(%esp) + flds 4(%esp) + fscale + fstp %st(1) + ret + +.globl _scalblnf + .set _scalblnf,_scalbnf diff --git a/winsup/mingw/mingwex/math/scalbnl.S b/winsup/mingw/mingwex/math/scalbnl.S new file mode 100644 index 000000000..77eaff7be --- /dev/null +++ b/winsup/mingw/mingwex/math/scalbnl.S @@ -0,0 +1,20 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Changes for long double by Ulrich Drepper <drepper@cygnus.com> + * Public domain. + */ + + .file "scalbnl.S" + .text + .align 4 +.globl _scalbnl + .def _scalbnl; .scl 2; .type 32; .endef +_scalbnl: + fildl 16(%esp) + fldt 4(%esp) + fscale + fstp %st(1) + ret + +.globl _scalblnl + .set _scalblnl,_scalbnl diff --git a/winsup/mingw/mingwex/signbit.c b/winsup/mingw/mingwex/math/signbit.c index 7f86c86a3..7f86c86a3 100644 --- a/winsup/mingw/mingwex/signbit.c +++ b/winsup/mingw/mingwex/math/signbit.c diff --git a/winsup/mingw/mingwex/signbitf.c b/winsup/mingw/mingwex/math/signbitf.c index 5bbf675ad..5bbf675ad 100644 --- a/winsup/mingw/mingwex/signbitf.c +++ b/winsup/mingw/mingwex/math/signbitf.c diff --git a/winsup/mingw/mingwex/signbitl.c b/winsup/mingw/mingwex/math/signbitl.c index 78f990350..78f990350 100644 --- a/winsup/mingw/mingwex/signbitl.c +++ b/winsup/mingw/mingwex/math/signbitl.c diff --git a/winsup/mingw/mingwex/math/sinf.S b/winsup/mingw/mingwex/math/sinf.S new file mode 100644 index 000000000..23e986d11 --- /dev/null +++ b/winsup/mingw/mingwex/math/sinf.S @@ -0,0 +1,32 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * + * Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>. + * + * Removed header file dependency for use in libmingwex.a by + * Danny Smith <dannysmith@users.sourceforge.net> + */ + + .file "sinf.S" + .text + .align 4 +.globl _sinf + .def _sinf; .scl 2; .type 32; .endef +_sinf: + flds 4(%esp) + fsin + fnstsw %ax + testl $0x400,%eax + jnz 1f + ret +1: fldpi + fadd %st(0) + fxch %st(1) +2: fprem1 + fnstsw %ax + testl $0x400,%eax + jnz 2b + fstp %st(1) + fsin + ret diff --git a/winsup/mingw/mingwex/math/sinhf.c b/winsup/mingw/mingwex/math/sinhf.c new file mode 100644 index 000000000..3d6bcff41 --- /dev/null +++ b/winsup/mingw/mingwex/math/sinhf.c @@ -0,0 +1,3 @@ +#include <math.h> +float sinhf (float x) + {return (float) sinh (x);} diff --git a/winsup/mingw/mingwex/math/sinhl.c b/winsup/mingw/mingwex/math/sinhl.c new file mode 100644 index 000000000..ca6a370b9 --- /dev/null +++ b/winsup/mingw/mingwex/math/sinhl.c @@ -0,0 +1,172 @@ +/* sinhl.c + * + * Hyperbolic sine, long double precision + * + * + * + * SYNOPSIS: + * + * long double x, y, sinhl(); + * + * y = sinhl( x ); + * + * + * + * DESCRIPTION: + * + * Returns hyperbolic sine of argument in the range MINLOGL to + * MAXLOGL. + * + * The range is partitioned into two segments. If |x| <= 1, a + * rational function of the form x + x**3 P(x)/Q(x) is employed. + * Otherwise the calculation is sinh(x) = ( exp(x) - exp(-x) )/2. + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -2,2 10000 1.5e-19 3.9e-20 + * IEEE +-10000 30000 1.1e-19 2.8e-20 + * + */ + +/* +Cephes Math Library Release 2.7: January, 1998 +Copyright 1984, 1991, 1998 by Stephen L. Moshier +*/ + +/* +Modified for mingw +2002-07-22 Danny Smith <dannysmith@users.sourceforge.net> +*/ + +#ifdef __MINGW32__ +#include "cephes_mconf.h" +#else +#include "mconf.h" +#endif + +#ifndef _SET_ERRNO +#define _SET_ERRNO(x) +#endif + +#ifdef UNK +static long double P[] = { + 1.7550769032975377032681E-6L, + 4.1680702175874268714539E-4L, + 3.0993532520425419002409E-2L, + 9.9999999999999999998002E-1L, +}; +static long double Q[] = { + 1.7453965448620151484660E-8L, +-5.9116673682651952419571E-6L, + 1.0599252315677389339530E-3L, +-1.1403880487744749056675E-1L, + 6.0000000000000000000200E0L, +}; +#endif + +#ifdef IBMPC +static const unsigned short P[] = { +0xec6a,0xd942,0xfbb3,0xeb8f,0x3feb, XPD +0x365e,0xb30a,0xe437,0xda86,0x3ff3, XPD +0x8890,0x01f6,0x2612,0xfde6,0x3ff9, XPD +0x0000,0x0000,0x0000,0x8000,0x3fff, XPD +}; +static const unsigned short Q[] = { +0x4edd,0x4c21,0xad09,0x95ed,0x3fe5, XPD +0x4376,0x9b70,0xd605,0xc65c,0xbfed, XPD +0xc8ad,0x5d21,0x3069,0x8aed,0x3ff5, XPD +0x9c32,0x6374,0x2d4b,0xe98d,0xbffb, XPD +0x0000,0x0000,0x0000,0xc000,0x4001, XPD +}; +#endif + +#ifdef MIEEE +static long P[] = { +0x3feb0000,0xeb8ffbb3,0xd942ec6a, +0x3ff30000,0xda86e437,0xb30a365e, +0x3ff90000,0xfde62612,0x01f68890, +0x3fff0000,0x80000000,0x00000000, +}; +static long Q[] = { +0x3fe50000,0x95edad09,0x4c214edd, +0xbfed0000,0xc65cd605,0x9b704376, +0x3ff50000,0x8aed3069,0x5d21c8ad, +0xbffb0000,0xe98d2d4b,0x63749c32, +0x40010000,0xc0000000,0x00000000, +}; +#endif + +#ifndef __MINGW32__ +extern long double MAXNUML, MAXLOGL, MINLOGL, LOGE2L; +#ifdef ANSIPROT +extern long double fabsl ( long double ); +extern long double expl ( long double ); +extern long double polevll ( long double, void *, int ); +extern long double p1evll ( long double, void *, int ); +#else +long double fabsl(), expl(), polevll(), p1evll(); +#endif +#ifdef INFINITIES +extern long double INFINITYL; +#endif +#ifdef NANS +extern long double NANL; +#endif +#endif /* __MINGW32__ */ + +long double sinhl(x) +long double x; +{ +long double a; + +#ifdef MINUSZERO +if( x == 0.0 ) + return(x); +#endif +#ifdef NANS +if (isnanl(x)) + { + _SET_ERRNO(EDOM); + } +#endif +a = fabsl(x); +if( (x > (MAXLOGL + LOGE2L)) || (x > -(MINLOGL-LOGE2L) ) ) + { + mtherr( "sinhl", DOMAIN ); + _SET_ERRNO(ERANGE); +#ifdef INFINITIES + if( x > 0.0L ) + return( INFINITYL ); + else + return( -INFINITYL ); +#else + if( x > 0.0L ) + return( MAXNUML ); + else + return( -MAXNUML ); +#endif + } +if( a > 1.0L ) + { + if( a >= (MAXLOGL - LOGE2L) ) + { + a = expl(0.5L*a); + a = (0.5L * a) * a; + if( x < 0.0L ) + a = -a; + return(a); + } + a = expl(a); + a = 0.5L*a - (0.5L/a); + if( x < 0.0L ) + a = -a; + return(a); + } + +a *= a; +return( x + x * a * (polevll(a,P,3)/polevll(a,Q,4)) ); +} diff --git a/winsup/mingw/mingwex/math/sinl.S b/winsup/mingw/mingwex/math/sinl.S new file mode 100644 index 000000000..16b2d9e50 --- /dev/null +++ b/winsup/mingw/mingwex/math/sinl.S @@ -0,0 +1,32 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * + * Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>. + * + * Removed header file dependency for use in libmingwex.a by + * Danny Smith <dannysmith@users.sourceforge.net> + */ + + .file "sinl.S" + .text + .align 4 +.globl _sinl + .def _sinl; .scl 2; .type 32; .endef +_sinl: + fldt 4(%esp) + fsin + fnstsw %ax + testl $0x400,%eax + jnz 1f + ret +1: fldpi + fadd %st(0) + fxch %st(1) +2: fprem1 + fnstsw %ax + testl $0x400,%eax + jnz 2b + fstp %st(1) + fsin + ret diff --git a/winsup/mingw/mingwex/math/sqrtf.c b/winsup/mingw/mingwex/math/sqrtf.c new file mode 100644 index 000000000..55ca39dbe --- /dev/null +++ b/winsup/mingw/mingwex/math/sqrtf.c @@ -0,0 +1,9 @@ +#include <math.h> + +float +sqrtf (float x) +{ + float res; + asm ("fsqrt" : "=t" (res) : "0" (x)); + return res; +} diff --git a/winsup/mingw/mingwex/math/sqrtl.c b/winsup/mingw/mingwex/math/sqrtl.c new file mode 100644 index 000000000..0bd301390 --- /dev/null +++ b/winsup/mingw/mingwex/math/sqrtl.c @@ -0,0 +1,8 @@ +#include <math.h> +long double +sqrtl (long double x) +{ + long double res; + asm ("fsqrt" : "=t" (res) : "0" (x)); + return res; +} diff --git a/winsup/mingw/mingwex/math/tanf.S b/winsup/mingw/mingwex/math/tanf.S new file mode 100644 index 000000000..540fc6836 --- /dev/null +++ b/winsup/mingw/mingwex/math/tanf.S @@ -0,0 +1,31 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * + * Removed header file dependency for use in libmingwex.a by + * Danny Smith <dannysmith@users.sourceforge.net> + */ + .file "tanf.S" + .text + .align 4 +.globl _tanf + .def _tanf; .scl 2; .type 32; .endef +_tanf: + flds 4(%esp) + fptan + fnstsw %ax + testl $0x400,%eax + jnz 1f + fstp %st(0) + ret +1: fldpi + fadd %st(0) + fxch %st(1) +2: fprem1 + fstsw %ax + testl $0x400,%eax + jnz 2b + fstp %st(1) + fptan + fstp %st(0) + ret diff --git a/winsup/mingw/mingwex/math/tanhf.c b/winsup/mingw/mingwex/math/tanhf.c new file mode 100644 index 000000000..b7c56f05c --- /dev/null +++ b/winsup/mingw/mingwex/math/tanhf.c @@ -0,0 +1,3 @@ +#include <math.h> +float tanhf (float x) + {return (float) tanh (x);} diff --git a/winsup/mingw/mingwex/math/tanhl.c b/winsup/mingw/mingwex/math/tanhl.c new file mode 100644 index 000000000..3727b3990 --- /dev/null +++ b/winsup/mingw/mingwex/math/tanhl.c @@ -0,0 +1,151 @@ +/* tanhl.c + * + * Hyperbolic tangent, long double precision + * + * + * + * SYNOPSIS: + * + * long double x, y, tanhl(); + * + * y = tanhl( x ); + * + * + * + * DESCRIPTION: + * + * Returns hyperbolic tangent of argument in the range MINLOGL to + * MAXLOGL. + * + * A rational function is used for |x| < 0.625. The form + * x + x**3 P(x)/Q(x) of Cody _& Waite is employed. + * Otherwise, + * tanh(x) = sinh(x)/cosh(x) = 1 - 2/(exp(2x) + 1). + * + * + * + * ACCURACY: + * + * Relative error: + * arithmetic domain # trials peak rms + * IEEE -2,2 30000 1.3e-19 2.4e-20 + * + */ + +/* +Cephes Math Library Release 2.7: May, 1998 +Copyright 1984, 1987, 1989, 1998 by Stephen L. Moshier +*/ + +/* +Modified for mingw +2002-07-22 Danny Smith <dannysmith@users.sourceforge.net> +*/ + +#ifdef __MINGW32__ +#include "cephes_mconf.h" +#else +#include "mconf.h" +#endif + +#ifndef _SET_ERRNO +#define _SET_ERRNO(x) +#endif + +#ifdef UNK +static long double P[] = { +-6.8473739392677100872869E-5L, +-9.5658283111794641589011E-1L, +-8.4053568599672284488465E1L, +-1.3080425704712825945553E3L, +}; +static long double Q[] = { +/* 1.0000000000000000000000E0L,*/ + 9.6259501838840336946872E1L, + 1.8218117903645559060232E3L, + 3.9241277114138477845780E3L, +}; +#endif + +#ifdef IBMPC +static short P[] = { +0xd2a4,0x1b0c,0x8f15,0x8f99,0xbff1, XPD +0x5959,0x9111,0x9cc7,0xf4e2,0xbffe, XPD +0xb576,0xef5e,0x6d57,0xa81b,0xc005, XPD +0xe3be,0xbfbd,0x5cbc,0xa381,0xc009, XPD +}; +static short Q[] = { +/*0x0000,0x0000,0x0000,0x8000,0x3fff,*/ +0x687f,0xce24,0xdd6c,0xc084,0x4005, XPD +0x3793,0xc95f,0xfa2f,0xe3b9,0x4009, XPD +0xd5a2,0x1f9c,0x0b1b,0xf542,0x400a, XPD +}; +#endif + +#ifdef MIEEE +static long P[] = { +0xbff10000,0x8f998f15,0x1b0cd2a4, +0xbffe0000,0xf4e29cc7,0x91115959, +0xc0050000,0xa81b6d57,0xef5eb576, +0xc0090000,0xa3815cbc,0xbfbde3be, +}; +static long Q[] = { +/*0x3fff0000,0x80000000,0x00000000,*/ +0x40050000,0xc084dd6c,0xce24687f, +0x40090000,0xe3b9fa2f,0xc95f3793, +0x400a0000,0xf5420b1b,0x1f9cd5a2, +}; +#endif + +#ifndef __MINGW32__ +extern long double MAXLOGL; +#ifdef ANSIPROT +extern long double fabsl ( long double ); +extern long double expl ( long double ); +extern long double polevll ( long double, void *, int ); +extern long double p1evll ( long double, void *, int ); +#else +long double fabsl(), expl(), polevll(), p1evll(); +#endif +#endif /* __MINGW32__ */ + +long double tanhl(x) +long double x; +{ +long double s, z; + +#ifdef MINUSZERO +if( x == 0.0L ) + return(x); +#endif +if (isnanl(x)) + { + _SET_ERRNO (EDOM); + return x; + } + +z = fabsl(x); +if( z > 0.5L * MAXLOGL ) + { + _SET_ERRNO (ERANGE); + if( x > 0 ) + return( 1.0L ); + else + return( -1.0L ); + } +if( z >= 0.625L ) + { + s = expl(2.0*z); + z = 1.0L - 2.0/(s + 1.0L); + if( x < 0 ) + z = -z; + } +else + { + s = x * x; + z = polevll( s, P, 3 )/p1evll(s, Q, 3); + z = x * s * z; + z = x + z; + } +return( z ); +} diff --git a/winsup/mingw/mingwex/math/tanl.S b/winsup/mingw/mingwex/math/tanl.S new file mode 100644 index 000000000..fd30019a8 --- /dev/null +++ b/winsup/mingw/mingwex/math/tanl.S @@ -0,0 +1,33 @@ +/* + * Written by J.T. Conklin <jtc@netbsd.org>. + * Public domain. + * + * Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>. + * + * Removed header file dependency for use in libmingwex.a by + * Danny Smith <dannysmith@users.sourceforge.net> + */ + .file "tanl.S" + .text + .align 4 +.globl _tanl + .def _tanl; .scl 2; .type 32; .endef +_tanl: + fldt 4(%esp) + fptan + fnstsw %ax + testl $0x400,%eax + jnz 1f + fstp %st(0) + ret +1: fldpi + fadd %st(0) + fxch %st(1) +2: fprem1 + fstsw %ax + testl $0x400,%eax + jnz 2b + fstp %st(1) + fptan + fstp %st(0) + ret diff --git a/winsup/mingw/mingwex/trunc.c b/winsup/mingw/mingwex/math/trunc.c index 2b9931255..2b9931255 100644 --- a/winsup/mingw/mingwex/trunc.c +++ b/winsup/mingw/mingwex/math/trunc.c diff --git a/winsup/mingw/mingwex/truncf.c b/winsup/mingw/mingwex/math/truncf.c index 53fccb153..53fccb153 100644 --- a/winsup/mingw/mingwex/truncf.c +++ b/winsup/mingw/mingwex/math/truncf.c diff --git a/winsup/mingw/mingwex/truncl.c b/winsup/mingw/mingwex/math/truncl.c index 908197acc..908197acc 100644 --- a/winsup/mingw/mingwex/truncl.c +++ b/winsup/mingw/mingwex/math/truncl.c diff --git a/winsup/mingw/mingwex/math_stubs.c b/winsup/mingw/mingwex/math_stubs.c deleted file mode 100644 index 225439846..000000000 --- a/winsup/mingw/mingwex/math_stubs.c +++ /dev/null @@ -1,12 +0,0 @@ -
-#include <math.h>
-
-double copysign (double x, double y) {return _copysign(x, y);}
-float copysignf (float x, float y) {return _copysign(x, y);}
-double logb (double x) {return _logb(x);}
-float logbf (float x) {return _logb( x );}
-double nextafter(double x, double y) {return _nextafter(x, y);}
-float nextafterf(float x, float y) {return _nextafter(x, y);}
-double scalb (double x, long i) {return _scalb (x, i);}
-float scalbf (float x, long i) {return _scalb(x, i);}
-
|