-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1 Please review this initial implementation for the remaining math functions I'm unsure about the folder name and the license for taylor_private.h, it mostly contains defines/typedefs/consts from FreeBSD... -----BEGIN PGP SIGNATURE----- Version: GnuPG v1 Comment: Using GnuPG with Thunderbird - http://www.enigmail.net/
iQGcBAEBAgAGBQJTn0KlAAoJEGm5GZTakYssggYL/iNc6+WuWR8Nx+k06Zjy4yxl CCB3/c9AL/uqSZdUaqHr0HGTZsy6ny8iLXUy082XnY5wE/87QInWPZ7hZjJ/5WKc 3pay4H1FufDKfrHJRBkcWD1z6PDiYNWboTKkHVOrNLcA6rQV8iOjFocIdepQjiF7 nIJDx+/i+YQ2VDNt3KgHhYnuR3smZFC2vhKSOR61XVx72v+CZ78IBE072ONvJyfM UxKa/FNEzKx5yF+ebC60nX15KldS+xBGcnsSsq8rdQpRN8dog0Z/vfQepV3yrzLV mx3s7nR3wU3b4j6rnxwnAkIAqfjkZk6JbPVB6I+rx+zKekkg966ZkNZX/TQeh2BM WyX+jlnRpwolLRPXN77Ywt/FQ6UJ8S8AoCaIWh7rpPtBADUh6z/mZxhH1F+h0Zs/ rTSE6w+dGWj3Szh1K/boqSdM6RdBg8gqrDr7SzNFy2XKgvuRjp2O/nUT3SB/nJHh nNoDvEt4LerCYnplOa/rfFDBd9kLAQiSacVw1iCPow== =4G6y -----END PGP SIGNATURE-----
commit 6d9072f9de43493e177248e2cf080899567f2be4 Author: André Hentschel <[email protected]> Date: Mon Jun 16 20:59:39 2014 +0200 sinf cosf draft diff --git a/mingw-w64-crt/Makefile.am b/mingw-w64-crt/Makefile.am index b530a7b..c7ac221 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 \ \ @@ -224,12 +225,12 @@ src_libmingwex=\ 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/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/scalbnf.S math/scalbnl.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/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/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/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 \ @@ -245,7 +246,7 @@ src_libmingwex=\ 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/sin.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/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 \ @@ -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/cosf.c math/sinl_internal.S math/sinf.c + +# these only go into the 32 bit version: +src_libmingwex32=math/cosf.c math/sinl_internal.S math/sinf.c + +# these only go into the ARM32 version: +src_libmingwexarm32=math/taylor/e_fmodf.c math/taylor/e_powf.c math/taylor/taylor.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/taylor/e_fmodf.c b/mingw-w64-crt/math/taylor/e_fmodf.c new file mode 100644 index 0000000..ae097b0 --- /dev/null +++ b/mingw-w64-crt/math/taylor/e_fmodf.c @@ -0,0 +1,102 @@ +/* e_fmodf.c -- float version of e_fmod.c. + * Conversion to float by Ian Lance Taylor, Cygnus Support, [email protected]. + */ + +/* + * ==================================================== + * 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 "taylor_private.h" + + + + +float +__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/taylor/e_powf.c b/mingw-w64-crt/math/taylor/e_powf.c new file mode 100644 index 0000000..e5e1e01 --- /dev/null +++ b/mingw-w64-crt/math/taylor/e_powf.c @@ -0,0 +1,216 @@ +/* e_powf.c -- float version of e_pow.c. + * Conversion to float by Ian Lance Taylor, Cygnus Support, [email protected]. + */ + +/* + * ==================================================== + * 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 "taylor_private.h" + + + + +float +__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/taylor/taylor.c b/mingw-w64-crt/math/taylor/taylor.c new file mode 100644 index 0000000..971b28a --- /dev/null +++ b/mingw-w64-crt/math/taylor/taylor.c @@ -0,0 +1,93 @@ +/* + 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 "taylor_private.h" + +static inline float factorial(float number) +{ + if (number == 0) + return 1; + + return number * factorial(number - 1); +} + +float cosf(float x) +{ + float result = 0.0; + int n, quadrant = (int)(x / M_PI_2) % 4; + + x = __ieee754_fmodf(x, M_PI_2); + if (quadrant == 1 || quadrant == 3) + x = M_PI_2 - x; + + for(n = 1; n < 6; n++) + { + float sign = ((n + 1) % 2 == 0) ? 1 : -1; + result += sign * (__ieee754_powf(x, ((2 * n) - 2)) / factorial((2 * n) - 2)); + } + + if (quadrant == 1 || quadrant == 2) + return result * -1; + return result; +} + +float sinf(float x) +{ + float result = 0.0; + int n, quadrant = (int)(x / M_PI_2) % 4; + + x = __ieee754_fmodf(x, M_PI_2); + if (quadrant == 1 || quadrant == 3) + x = M_PI_2 - x; + + for(n = 1; n < 6; n++) + { + float sign = ((n + 1) % 2 == 0) ? 1 : -1; + result += sign * (__ieee754_powf(x, ((2 * n) - 1)) / factorial((2 * n) - 1)); + } + + if (quadrant == 2 || quadrant == 3) + return result * -1; + return result; +} diff --git a/mingw-w64-crt/math/taylor/taylor_private.h b/mingw-w64-crt/math/taylor/taylor_private.h new file mode 100644 index 0000000..80deabd --- /dev/null +++ b/mingw-w64-crt/math/taylor/taylor_private.h @@ -0,0 +1,60 @@ +#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 __ieee754_fmodf(float, float); +extern float __ieee754_powf(float, float);
tmp.diff.sig
Description: PGP signature
------------------------------------------------------------------------------ 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
