Am 21.06.2014 17:10, schrieb André Hentschel:
> Hi,
> this time i just put some commits on github for review [1],
> so please review, then i'll commit it except of the temporary tests.
> (btw, this is the work of more than two weeks...)
>
> [1] https://github.com/AndreRH/mingw-w64/commits/master
>
Here it is as patch (without Makefile.in)
diff --git a/mingw-w64-crt/Makefile.am b/mingw-w64-crt/Makefile.am
index b530a7b..5379a1d 100644
--- a/mingw-w64-crt/Makefile.am
+++ b/mingw-w64-crt/Makefile.am
@@ -195,6 +195,7 @@ src_msvcrt32=\
$(src_msvcrt) \
misc/lc_locale_func.c
+# These mingwex sources are target independent:
src_libmingwex=\
crt/dllentry.c crt/dllmain.c \
\
@@ -219,34 +220,34 @@ src_libmingwex=\
gdtoa/strtopx.c gdtoa/sum.c gdtoa/ulp.c \
\
math/_chgsignl.S math/ceil.S math/ceilf.S math/ceill.S math/copysignl.S math/cosl_internal.S \
- math/exp2.S math/exp2f.S math/exp2l.S math/floor.S math/floorf.S math/floorl.S \
- math/fma.S math/fmaf.S math/frexpl.S math/ilogb.S math/ilogbf.S math/ilogbl.S \
- math/internal_logl.S math/log10l.S math/log1p.S math/log1pf.S math/log1pl.S math/log2.S \
- math/log2f.S math/log2l.S math/nearbyint.S math/nearbyintf.S math/nearbyintl.S math/remainder.S \
+ math/exp2.S math/exp2l.S math/floor.S math/floorf.S math/floorl.S \
+ math/fma.S math/fmaf.S math/frexpl.S math/ilogb.S math/ilogbl.S \
+ math/internal_logl.S math/log10l.S math/log1p.S math/log1pl.S math/log2.S \
+ math/log2l.S math/nearbyint.S math/nearbyintf.S math/nearbyintl.S math/remainder.S \
math/remainderf.S math/remainderl.S math/remquo.S math/remquof.S math/remquol.S math/scalbn.S \
- math/scalbnf.S math/scalbnl.S math/sinl_internal.S math/tanl.S math/trunc.S math/truncf.S \
- math/acosf.c math/acosh.c math/acoshf.c math/acoshl.c math/acosl.c math/asinf.c \
- math/asinh.c math/asinhf.c math/asinhl.c math/asinl.c math/atan2f.c math/atan2l.c \
- math/atanf.c math/atanh.c math/atanhf.c math/atanhl.c math/atanl.c math/cbrt.c \
+ math/scalbnl.S math/sinl_internal.S math/tanl.S math/trunc.S math/truncf.S \
+ math/acosh.c math/acoshl.c math/acosl.c \
+ math/asinh.c math/asinhl.c math/asinl.c math/atan2l.c \
+ math/atanh.c math/atanhl.c math/atanl.c math/cbrt.c \
math/cbrtf.c math/cbrtl.c math/cephes_emath.c math/copysign.c math/copysignf.c math/cos.c \
- math/cosf.c math/coshf.c math/coshl.c math/cosl.c math/cossin.c math/erfl.c \
- math/exp.c math/expf.c math/expl.c math/expm1.c math/expm1f.c math/expm1l.c \
+ math/coshf.c math/coshl.c math/cosl.c math/cossin.c math/erfl.c \
+ math/exp.c math/expf.c math/expl.c math/expm1.c math/expm1l.c \
math/fabs.c math/fabsf.c math/fabsl.c math/fdim.c math/fdimf.c math/fdiml.c \
math/fmal.c math/fmax.c math/fmaxf.c math/fmaxl.c math/fmin.c math/fminf.c \
- math/fminl.c math/fmod.c math/fmodf.c math/fmodl.c math/fp_consts.c math/fp_constsf.c \
+ math/fminl.c math/fmod.c math/fmodl.c math/fp_consts.c math/fp_constsf.c \
math/fp_constsl.c math/fpclassify.c math/fpclassifyf.c math/fpclassifyl.c math/frexpf.c math/fucom.c \
math/hypotf.c math/hypotl.c math/isnan.c math/isnanf.c math/isnanl.c math/ldexp.c \
math/ldexpf.c math/ldexpl.c math/lgamma.c math/lgammaf.c math/lgammal.c math/llrint.c \
math/signgam.c \
math/llrintf.c math/llrintl.c math/llround.c math/llroundf.c math/llroundl.c math/log.c \
- math/log10f.c math/logb.c math/logbf.c math/logbl.c math/logf.c math/logl.c \
+ math/log10f.c math/logb.c math/logbl.c math/logf.c math/logl.c \
math/lrint.c math/lrintf.c math/lrintl.c math/lround.c math/lroundf.c math/lroundl.c \
math/modf.c math/modff.c math/modfl.c math/nextafterf.c math/nextafterl.c math/nexttoward.c \
math/nexttowardf.c math/pow.c math/powf.c math/powi.c math/powif.c math/powil.c \
math/powl.c math/rint.c math/rintf.c math/rintl.c math/round.c math/roundf.c \
math/roundl.c math/s_erf.c math/sf_erf.c math/signbit.c math/signbitf.c math/signbitl.c \
- math/sin.c math/sinf.c math/sinhf.c math/sinhl.c math/sinl.c math/sqrt.c \
- math/sqrtf.c math/sqrtl.c math/tanf.c math/tanhf.c math/tanhl.c math/tgamma.c \
+ math/sin.c math/sinhf.c math/sinhl.c math/sinl.c math/sqrt.c \
+ math/sqrtf.c math/sqrtl.c math/tanhf.c math/tanhl.c math/tgamma.c \
math/tgammaf.c math/tgammal.c math/truncl.c \
math/acosh.def.h math/cos.def.h math/exp.def.h math/expm1.def.h math/log.def.h math/pow.def.h \
math/powi.def.h math/sin.def.h math/sqrt.def.h \
@@ -292,6 +293,16 @@ src_libmingwex=\
stdio/vscanf.c stdio/vsnprintf.c stdio/vsnprintf_s.c stdio/vsnwprintf.c stdio/vsscanf.c \
stdio/vswscanf.c stdio/vwscanf.c stdio/wtoll.c stdio/mingw_asprintf.c stdio/mingw_vasprintf.c
+# these only go into the 64 bit version:
+src_libmingwex64=math/acoshf.c math/asinhf.c math/atanhf.c math/cosf.c math/exp2f.S math/expm1f.c math/fmodf.c math/ilogbf.S math/log1pf.S math/log2f.S math/logbf.c math/scalbnf.S math/sinf.c math/tanf.c math/asinf.c math/acosf.c math/atan2f.c math/atanf.c
+
+# these only go into the 32 bit version:
+src_libmingwex32=math/acoshf.c math/asinhf.c math/atanhf.c math/cosf.c math/exp2f.S math/expm1f.c math/fmodf.c math/ilogbf.S math/log1pf.S math/log2f.S math/logbf.c math/scalbnf.S math/sinf.c math/tanf.c math/asinf.c math/acosf.c math/atan2f.c math/atanf.c
+
+# these only go into the ARM32 version:
+src_libmingwexarm32=math/softmath/e_fmodf.c math/softmath/e_powf.c math/softmath/softmath.c
+
+
# These intrinsics are target independent:
src_intrincs= \
intrincs/__movsb.c intrincs/__movsd.c intrincs/__movsw.c intrincs/__stosb.c intrincs/__stosd.c \
@@ -513,7 +524,7 @@ lib32_libmingw32_a_SOURCES = $(src_libmingw32)
lib32_LIBRARIES += lib32/libmingwex.a
lib32_libmingwex_a_CPPFLAGS=$(CPPFLAGS32) $(extra_include) $(AM_CPPFLAGS)
-lib32_libmingwex_a_SOURCES = $(src_libmingwex) $(src_dfp_math)
+lib32_libmingwex_a_SOURCES = $(src_libmingwex) $(src_libmingwex32) $(src_dfp_math)
lib32_libmingwex_a_LIBADD = $(obj_dfpsrc32)
lib32_LIBRARIES += lib32/libmoldname.a
@@ -758,7 +769,7 @@ lib64_libmingw32_a_SOURCES = $(src_libmingw32)
lib64_LIBRARIES += lib64/libmingwex.a
lib64_libmingwex_a_CPPFLAGS=$(CPPFLAGS64) $(extra_include) $(AM_CPPFLAGS)
-lib64_libmingwex_a_SOURCES = $(src_libmingwex) $(src_dfp_math)
+lib64_libmingwex_a_SOURCES = $(src_libmingwex) $(src_libmingwex64) $(src_dfp_math)
lib64_libmingwex_a_LIBADD = $(obj_dfpsrc64)
lib64_LIBRARIES += lib64/libmoldname.a
@@ -1269,7 +1280,7 @@ libarm32_libmingw32_a_SOURCES = $(src_libmingw32)
libarm32_LIBRARIES += libarm32/libmingwex.a
libarm32_libmingwex_a_CPPFLAGS=$(CPPFLAGSARM32) $(extra_include) $(AM_CPPFLAGS)
-libarm32_libmingwex_a_SOURCES = $(src_libmingwex) $(src_dfp_math)
+libarm32_libmingwex_a_SOURCES = $(src_libmingwex) $(src_libmingwexarm32) $(src_dfp_math)
libarm32_libmingwex_a_LIBADD = $(obj_dfpsrc32)
libarm32_LIBRARIES += libarm32/libmoldname.a
diff --git a/mingw-w64-crt/math/softmath/e_fmodf.c b/mingw-w64-crt/math/softmath/e_fmodf.c
new file mode 100644
index 0000000..c57daf4
--- /dev/null
+++ b/mingw-w64-crt/math/softmath/e_fmodf.c
@@ -0,0 +1,94 @@
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+/*
+ * __ieee754_fmodf(x,y)
+ * Return x mod y in exact arithmetic
+ * Method: shift and subtract
+ */
+
+#include "math.h"
+#include "softmath_private.h"
+
+float bsd__ieee754_fmodf(float x, float y)
+{
+ int32_t n,hx,hy,hz,ix,iy,sx,i;
+
+ GET_FLOAT_WORD(hx,x);
+ GET_FLOAT_WORD(hy,y);
+ sx = hx&0x80000000; /* sign of x */
+ hx ^=sx; /* |x| */
+ hy &= 0x7fffffff; /* |y| */
+
+ /* purge off exception values */
+ if(hy==0||(hx>=0x7f800000)|| /* y=0,or x not finite */
+ (hy>0x7f800000)) /* or y is NaN */
+ return (x*y)/(x*y);
+ if(hx<hy) return x; /* |x|<|y| return x */
+ if(hx==hy)
+ return Zero[(u_int32_t)sx>>31]; /* |x|=|y| return x*0*/
+
+ /* determine ix = ilogb(x) */
+ if(hx<0x00800000) { /* subnormal x */
+ for (ix = -126,i=(hx<<8); i>0; i<<=1) ix -=1;
+ } else ix = (hx>>23)-127;
+
+ /* determine iy = ilogb(y) */
+ if(hy<0x00800000) { /* subnormal y */
+ for (iy = -126,i=(hy<<8); i>=0; i<<=1) iy -=1;
+ } else iy = (hy>>23)-127;
+
+ /* set up {hx,lx}, {hy,ly} and align y to x */
+ if(ix >= -126)
+ hx = 0x00800000|(0x007fffff&hx);
+ else { /* subnormal x, shift x to normal */
+ n = -126-ix;
+ hx = hx<<n;
+ }
+ if(iy >= -126)
+ hy = 0x00800000|(0x007fffff&hy);
+ else { /* subnormal y, shift y to normal */
+ n = -126-iy;
+ hy = hy<<n;
+ }
+
+ /* fix point fmod */
+ n = ix - iy;
+ while(n--) {
+ hz=hx-hy;
+ if(hz<0){hx = hx+hx;}
+ else {
+ if(hz==0) /* return sign(x)*0 */
+ return Zero[(u_int32_t)sx>>31];
+ hx = hz+hz;
+ }
+ }
+ hz=hx-hy;
+ if(hz>=0) {hx=hz;}
+
+ /* convert back to floating value and restore the sign */
+ if(hx==0) /* return sign(x)*0 */
+ return Zero[(u_int32_t)sx>>31];
+ while(hx<0x00800000) { /* normalize x */
+ hx = hx+hx;
+ iy -= 1;
+ }
+ if(iy>= -126) { /* normalize output */
+ hx = ((hx-0x00800000)|((iy+127)<<23));
+ SET_FLOAT_WORD(x,hx|sx);
+ } else { /* subnormal output */
+ n = -126 - iy;
+ hx >>= n;
+ SET_FLOAT_WORD(x,hx|sx);
+ x *= one; /* create necessary signal */
+ }
+ return x; /* exact output */
+}
diff --git a/mingw-w64-crt/math/softmath/e_powf.c b/mingw-w64-crt/math/softmath/e_powf.c
new file mode 100644
index 0000000..faaf3dd
--- /dev/null
+++ b/mingw-w64-crt/math/softmath/e_powf.c
@@ -0,0 +1,208 @@
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#include "math.h"
+#include "softmath_private.h"
+
+float bsd__ieee754_powf(float x, float y)
+{
+ float z,ax,z_h,z_l,p_h,p_l;
+ float y1,t1,t2,r,s,sn,t,u,v,w;
+ int32_t i,j,k,yisint,n;
+ int32_t hx,hy,ix,iy,is;
+
+ GET_FLOAT_WORD(hx,x);
+ GET_FLOAT_WORD(hy,y);
+ ix = hx&0x7fffffff; iy = hy&0x7fffffff;
+
+ /* y==zero: x**0 = 1 */
+ if(iy==0) return one;
+
+ /* x==1: 1**y = 1, even if y is NaN */
+ if (hx==0x3f800000) return one;
+
+ /* y!=zero: result is NaN if either arg is NaN */
+ if(ix > 0x7f800000 ||
+ iy > 0x7f800000)
+ return (x+0.0F)+(y+0.0F);
+
+ /* determine if y is an odd int when x < 0
+ * yisint = 0 ... y is not an integer
+ * yisint = 1 ... y is an odd int
+ * yisint = 2 ... y is an even int
+ */
+ yisint = 0;
+ if(hx<0) {
+ if(iy>=0x4b800000) yisint = 2; /* even integer y */
+ else if(iy>=0x3f800000) {
+ k = (iy>>23)-0x7f; /* exponent */
+ j = iy>>(23-k);
+ if((j<<(23-k))==iy) yisint = 2-(j&1);
+ }
+ }
+
+ /* special value of y */
+ if (iy==0x7f800000) { /* y is +-inf */
+ if (ix==0x3f800000)
+ return one; /* (-1)**+-inf is NaN */
+ else if (ix > 0x3f800000)/* (|x|>1)**+-inf = inf,0 */
+ return (hy>=0)? y: zero;
+ else /* (|x|<1)**-,+inf = inf,0 */
+ return (hy<0)?-y: zero;
+ }
+ if(iy==0x3f800000) { /* y is +-1 */
+ if(hy<0) return one/x; else return x;
+ }
+ if(hy==0x40000000) return x*x; /* y is 2 */
+ if(hy==0x3f000000) { /* y is 0.5 */
+ if(hx>=0) /* x >= +0 */
+ return sqrtf(x);
+ }
+
+ ax = fabsf(x);
+ /* special value of x */
+ if(ix==0x7f800000||ix==0||ix==0x3f800000){
+ z = ax; /*x is +-0,+-inf,+-1*/
+ if(hy<0) z = one/z; /* z = (1/|x|) */
+ if(hx<0) {
+ if(((ix-0x3f800000)|yisint)==0) {
+ z = (z-z)/(z-z); /* (-1)**non-int is NaN */
+ } else if(yisint==1)
+ z = -z; /* (x<0)**odd = -(|x|**odd) */
+ }
+ return z;
+ }
+
+ n = ((u_int32_t)hx>>31)-1;
+
+ /* (x<0)**(non-int) is NaN */
+ if((n|yisint)==0) return (x-x)/(x-x);
+
+ sn = one; /* s (sign of result -ve**odd) = -1 else = 1 */
+ if((n|(yisint-1))==0) sn = -one; /* (-ve)**(odd int) */
+
+ /* |y| is huge */
+ if(iy>0x4d000000) { /* if |y| > 2**27 */
+ /* over/underflow if x is not close to one */
+ if(ix<0x3f7ffff8) return (hy<0)? sn*huge*huge:sn*tiny*tiny;
+ if(ix>0x3f800007) return (hy>0)? sn*huge*huge:sn*tiny*tiny;
+ /* now |1-x| is tiny <= 2**-20, suffice to compute
+ log(x) by x-x^2/2+x^3/3-x^4/4 */
+ t = ax-1; /* t has 20 trailing zeros */
+ w = (t*t)*((float)0.5-t*((float)0.333333333333-t*(float)0.25));
+ u = ivln2_h*t; /* ivln2_h has 16 sig. bits */
+ v = t*ivln2_l-w*ivln2;
+ t1 = u+v;
+ GET_FLOAT_WORD(is,t1);
+ SET_FLOAT_WORD(t1,is&0xfffff000);
+ t2 = v-(t1-u);
+ } else {
+ float s2,s_h,s_l,t_h,t_l;
+ n = 0;
+ /* take care subnormal number */
+ if(ix<0x00800000)
+ {ax *= two24; n -= 24; GET_FLOAT_WORD(ix,ax); }
+ n += ((ix)>>23)-0x7f;
+ j = ix&0x007fffff;
+ /* determine interval */
+ ix = j|0x3f800000; /* normalize ix */
+ if(j<=0x1cc471) k=0; /* |x|<sqrt(3/2) */
+ else if(j<0x5db3d7) k=1; /* |x|<sqrt(3) */
+ else {k=0;n+=1;ix -= 0x00800000;}
+ SET_FLOAT_WORD(ax,ix);
+
+ /* compute s = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
+ u = ax-bp[k]; /* bp[0]=1.0, bp[1]=1.5 */
+ v = one/(ax+bp[k]);
+ s = u*v;
+ s_h = s;
+ GET_FLOAT_WORD(is,s_h);
+ SET_FLOAT_WORD(s_h,is&0xfffff000);
+ /* t_h=ax+bp[k] High */
+ is = ((ix>>1)&0xfffff000)|0x20000000;
+ SET_FLOAT_WORD(t_h,is+0x00400000+(k<<21));
+ t_l = ax - (t_h-bp[k]);
+ s_l = v*((u-s_h*t_h)-s_h*t_l);
+ /* compute log(ax) */
+ s2 = s*s;
+ r = s2*s2*(L1+s2*(L2+s2*(L3+s2*(L4+s2*(L5+s2*L6)))));
+ r += s_l*(s_h+s);
+ s2 = s_h*s_h;
+ t_h = (float)3.0+s2+r;
+ GET_FLOAT_WORD(is,t_h);
+ SET_FLOAT_WORD(t_h,is&0xfffff000);
+ t_l = r-((t_h-(float)3.0)-s2);
+ /* u+v = s*(1+...) */
+ u = s_h*t_h;
+ v = s_l*t_h+t_l*s;
+ /* 2/(3log2)*(s+...) */
+ p_h = u+v;
+ GET_FLOAT_WORD(is,p_h);
+ SET_FLOAT_WORD(p_h,is&0xfffff000);
+ p_l = v-(p_h-u);
+ z_h = cp_h*p_h; /* cp_h+cp_l = 2/(3*log2) */
+ z_l = cp_l*p_h+p_l*cp+dp_l[k];
+ /* log2(ax) = (s+..)*2/(3*log2) = n + dp_h + z_h + z_l */
+ t = (float)n;
+ t1 = (((z_h+z_l)+dp_h[k])+t);
+ GET_FLOAT_WORD(is,t1);
+ SET_FLOAT_WORD(t1,is&0xfffff000);
+ t2 = z_l-(((t1-t)-dp_h[k])-z_h);
+ }
+
+ /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
+ GET_FLOAT_WORD(is,y);
+ SET_FLOAT_WORD(y1,is&0xfffff000);
+ p_l = (y-y1)*t1+y*t2;
+ p_h = y1*t1;
+ z = p_l+p_h;
+ GET_FLOAT_WORD(j,z);
+ if (j>0x43000000) /* if z > 128 */
+ return sn*huge*huge; /* overflow */
+ else if (j==0x43000000) { /* if z == 128 */
+ if(p_l+ovt>z-p_h) return sn*huge*huge; /* overflow */
+ }
+ else if ((j&0x7fffffff)>0x43160000) /* z <= -150 */
+ return sn*tiny*tiny; /* underflow */
+ else if ((j&0xffffffff)==0xc3160000){ /* z == -150 */
+ if(p_l<=z-p_h) return sn*tiny*tiny; /* underflow */
+ }
+ /*
+ * compute 2**(p_h+p_l)
+ */
+ i = j&0x7fffffff;
+ k = (i>>23)-0x7f;
+ n = 0;
+ if(i>0x3f000000) { /* if |z| > 0.5, set n = [z+0.5] */
+ n = j+(0x00800000>>(k+1));
+ k = ((n&0x7fffffff)>>23)-0x7f; /* new k for n */
+ SET_FLOAT_WORD(t,n&~(0x007fffff>>k));
+ n = ((n&0x007fffff)|0x00800000)>>(23-k);
+ if(j<0) n = -n;
+ p_h -= t;
+ }
+ t = p_l+p_h;
+ GET_FLOAT_WORD(is,t);
+ SET_FLOAT_WORD(t,is&0xffff8000);
+ u = t*lg2_h;
+ v = (p_l-(t-p_h))*lg2+t*lg2_l;
+ z = u+v;
+ w = v-(z-u);
+ t = z*z;
+ t1 = z - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
+ r = (z*t1)/(t1-two)-(w+z*w);
+ z = one-(r-z);
+ GET_FLOAT_WORD(j,z);
+ j += (n<<23);
+ if((j>>23)<=0) z = scalbnf(z,n); /* subnormal output */
+ else SET_FLOAT_WORD(z,j);
+ return sn*z;
+}
diff --git a/mingw-w64-crt/math/softmath/softmath.c b/mingw-w64-crt/math/softmath/softmath.c
new file mode 100644
index 0000000..b537094
--- /dev/null
+++ b/mingw-w64-crt/math/softmath/softmath.c
@@ -0,0 +1,244 @@
+/*
+ This Software is provided under the Zope Public License (ZPL) Version 2.1.
+
+ Copyright (c) 2014 by the mingw-w64 project
+
+ See the AUTHORS file for the list of contributors to the mingw-w64 project.
+
+ This license has been certified as open source. It has also been designated
+ as GPL compatible by the Free Software Foundation (FSF).
+
+ Redistribution and use in source and binary forms, with or without
+ modification, are permitted provided that the following conditions are met:
+
+ 1. Redistributions in source code must retain the accompanying copyright
+ notice, this list of conditions, and the following disclaimer.
+ 2. Redistributions in binary form must reproduce the accompanying
+ copyright notice, this list of conditions, and the following disclaimer
+ in the documentation and/or other materials provided with the
+ distribution.
+ 3. Names of the copyright holders must not be used to endorse or promote
+ products derived from this software without prior written permission
+ from the copyright holders.
+ 4. The right to distribute this software or to use it for any purpose does
+ not give you the right to use Servicemarks (sm) or Trademarks (tm) of
+ the copyright holders. Use of them is covered by separate agreement
+ with the copyright holders.
+ 5. If any files are modified, you must cause the modified files to carry
+ prominent notices stating that you changed the files and the date of
+ any change.
+
+ Disclaimer
+
+ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ``AS IS'' AND ANY EXPRESSED
+ OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
+ OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO
+ EVENT SHALL THE COPYRIGHT HOLDERS BE LIABLE FOR ANY DIRECT, INDIRECT,
+ INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+ LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
+ OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+ LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+ NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
+ EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+
+#include "softmath_private.h"
+
+static inline double factorial(double number)
+{
+ if (number <= 0)
+ return 1;
+
+ return number * factorial(number - 1);
+}
+
+static inline float softmath_expf(float x)
+{
+ float result = 0.0;
+ int n;
+
+ for(n = 0; n < 8; n++)
+ {
+ result += (bsd__ieee754_powf(x, n) / factorial(n));
+ if (isnan(result) || isinf(result))
+ break;
+ }
+
+ return result;
+}
+
+static inline float softmath_logf(float x)
+{
+ int n, aprox = 8 / (x / 2);
+ float result = 0.0;
+
+ if (x == 0.0)
+ return -HUGE_VALF;
+ else if (x < 0.0001)
+ aprox = 32768;
+
+ for(n = 0; n < aprox; n++)
+ {
+ result += bsd__ieee754_powf((x - 1.0) / (x + 1.0), 2 * n + 1) * (1.0 / (2.0 * n + 1.0));
+ if (isinf(result))
+ break;
+ }
+ result *= 2;
+
+ return result;
+}
+
+float acosf(float x)
+{
+ if (x > 1.0 || x < -1.0)
+ return NAN;
+
+ return M_PI_2 - asinf(x);
+}
+
+float acoshf(float x)
+{
+ return softmath_logf(x + sqrtf(x * x * 2 - 1));
+}
+
+float asinf(float x)
+{
+ double result = 0.0;
+ int n, aprox = 11;
+
+ if (x > 1.0 || x < -1.0)
+ return NAN;
+
+ if (x > 0.95) aprox *= 8;
+ else if (x > 0.9) aprox *= 4;
+ else if (x > 0.8) aprox *= 2;
+
+ for(n = 0; n < aprox; n++)
+ {
+ double n2ff = bsd__ieee754_powf(2, n) * factorial(n);
+ result += (factorial(2 * n) / (n2ff * n2ff)) * (bsd__ieee754_powf(x, ((2 * n) + 1)) / ((2 * n) + 1));
+ }
+
+ return result;
+}
+
+float asinhf(float x)
+{
+ return softmath_logf(x + sqrtf(x * x * 2 + 1));
+}
+
+float atan2f(float y, float x)
+{
+ return atanf(y / x);
+}
+
+float atanf(float x)
+{
+ if (x < -1.0)
+ return -1 * M_PI_2 - x / (0.28 + x * x);
+ else if (x > -1.0 && x < 1.0)
+ return x / (1.0 + 0.28 * x * x);
+ else if (x > 1.0)
+ return M_PI_2 - x / (0.28 + x * x);
+ return NAN;
+}
+
+float atanhf(float x)
+{
+ if (x > 1) return NAN;
+ else if (x < -1) return -NAN;
+ return softmath_logf((1 + x) / (1 - x)) / 2;
+}
+
+float cosf(float x)
+{
+ float result = 0.0, neg = (x < 0) ? -1 : 1;
+ int n, quadrant;
+
+ x *= neg;
+ quadrant = (int)(x / M_PI_2) % 4;
+
+ x = bsd__ieee754_fmodf(x, M_PI_2);
+ if (quadrant == 1 || quadrant == 3)
+ x = M_PI_2 - x;
+
+ for(n = 0; n < 5; n++)
+ {
+ float sign = (n % 2 == 0) ? 1 : -1;
+ result += sign * (bsd__ieee754_powf(x, (2 * n)) / factorial(2 * n));
+ }
+
+ if (quadrant == 1 || quadrant == 2)
+ result *= -1;
+
+ return result;
+}
+
+float exp2f(float x)
+{
+ return bsd__ieee754_powf(2, x);
+}
+
+float expm1f(float x)
+{
+ return softmath_expf(x) - 1;
+}
+
+float fmodf(float x, float y)
+{
+ return bsd__ieee754_fmodf(x, y);
+}
+
+int ilogbf(float x)
+{
+ return floorf(log2f(x));
+}
+
+float log1pf(float x)
+{
+ return softmath_logf(x + 1);
+}
+
+float log2f(float x)
+{
+ return softmath_logf(x) / softmath_logf(2);
+}
+
+float logbf(float x)
+{
+ return floorf(log2f(x));
+}
+
+float scalbnf(float x, int exp)
+{
+ return x * exp2f(exp);
+}
+
+float sinf(float x)
+{
+ float result = 0.0, neg = (x < 0) ? -1 : 1;
+ int n, quadrant;
+
+ x *= neg;
+ quadrant = (int)(x / M_PI_2) % 4;
+
+ x = bsd__ieee754_fmodf(x, M_PI_2);
+ if (quadrant == 1 || quadrant == 3)
+ x = M_PI_2 - x;
+
+ for(n = 0; n < 5; n++)
+ {
+ float sign = (n % 2 == 0) ? 1 : -1;
+ result += sign * (bsd__ieee754_powf(x, ((2 * n) + 1)) / factorial((2 * n) + 1));
+ }
+
+ if (quadrant == 2 || quadrant == 3)
+ result *= -1;
+
+ return neg * result;
+}
+
+float tanf(float x)
+{
+ return sinf(x) / cosf(x);
+}
diff --git a/mingw-w64-crt/math/softmath/softmath_private.h b/mingw-w64-crt/math/softmath/softmath_private.h
new file mode 100644
index 0000000..2a66f59
--- /dev/null
+++ b/mingw-w64-crt/math/softmath/softmath_private.h
@@ -0,0 +1,71 @@
+/*
+* ====================================================
+* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+*
+* Developed at SunPro, a Sun Microsystems, Inc. business.
+* Permission to use, copy, modify, and distribute this
+* software is freely granted, provided that this notice
+* is preserved.
+* ====================================================
+*/
+
+#include <math.h>
+#include <inttypes.h>
+
+static const float
+ bp[] = {1.0, 1.5,},
+ dp_h[] = { 0.0, 5.84960938e-01,}, /* 0x3f15c000 */
+ dp_l[] = { 0.0, 1.56322085e-06,}, /* 0x35d1cfdc */
+ Zero[] = {0.0, -0.0,},
+ zero = 0.0,
+ one = 1.0,
+ two = 2.0,
+ two24 = 16777216.0, /* 0x4b800000 */
+ huge = 1.0e30,
+ tiny = 1.0e-30,
+ /* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */
+ L1 = 6.0000002384e-01, /* 0x3f19999a */
+ L2 = 4.2857143283e-01, /* 0x3edb6db7 */
+ L3 = 3.3333334327e-01, /* 0x3eaaaaab */
+ L4 = 2.7272811532e-01, /* 0x3e8ba305 */
+ L5 = 2.3066075146e-01, /* 0x3e6c3255 */
+ L6 = 2.0697501302e-01, /* 0x3e53f142 */
+ P1 = 1.6666667163e-01, /* 0x3e2aaaab */
+ P2 = -2.7777778450e-03, /* 0xbb360b61 */
+ P3 = 6.6137559770e-05, /* 0x388ab355 */
+ P4 = -1.6533901999e-06, /* 0xb5ddea0e */
+ P5 = 4.1381369442e-08, /* 0x3331bb4c */
+ lg2 = 6.9314718246e-01, /* 0x3f317218 */
+ lg2_h = 6.93145752e-01, /* 0x3f317200 */
+ lg2_l = 1.42860654e-06, /* 0x35bfbe8c */
+ ovt = 4.2995665694e-08, /* -(128-log2(ovfl+.5ulp)) */
+ cp = 9.6179670095e-01, /* 0x3f76384f =2/(3ln2) */
+ cp_h = 9.6191406250e-01, /* 0x3f764000 =12b cp */
+ cp_l = -1.1736857402e-04, /* 0xb8f623c6 =tail of cp_h */
+ ivln2 = 1.4426950216e+00, /* 0x3fb8aa3b =1/ln2 */
+ ivln2_h= 1.4426879883e+00, /* 0x3fb8aa00 =16b 1/ln2*/
+ ivln2_l= 7.0526075433e-06; /* 0x36eca570 =1/ln2 tail*/
+
+typedef unsigned int u_int32_t;
+
+typedef union {
+ float value;
+ u_int32_t word;
+} ieee_float_shape_type;
+
+#define GET_FLOAT_WORD(i,d) do \
+{ \
+ ieee_float_shape_type gf_u; \
+ gf_u.value = (d); \
+ (i) = gf_u.word; \
+} while(0)
+
+#define SET_FLOAT_WORD(d,i) do \
+{ \
+ ieee_float_shape_type gf_u; \
+ gf_u.word = (i); \
+ (d) = gf_u.value; \
+} while(0)
+
+extern float bsd__ieee754_fmodf(float, float);
+extern float bsd__ieee754_powf(float, float);
------------------------------------------------------------------------------
HPCC Systems Open Source Big Data Platform from LexisNexis Risk Solutions
Find What Matters Most in Your Big Data with HPCC Systems
Open Source. Fast. Scalable. Simple. Ideal for Dirty Data.
Leverages Graph Analysis for Fast Processing & Easy Data Exploration
http://p.sf.net/sfu/hpccsystems
_______________________________________________
Mingw-w64-public mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/mingw-w64-public
