-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1 Please review, i'll commit it.
btw, is it OK to add my copyright notice like: Copyright (c) 2014 by the mingw-w64 project Copyright (c) 2014 André Hentschel -----BEGIN PGP SIGNATURE----- Version: GnuPG v1 Comment: Using GnuPG with Thunderbird - http://www.enigmail.net/ iQGcBAEBAgAGBQJTvu9CAAoJEGm5GZTakYssQ38L/1sgctgWIQbQgA7A4CQLPnSm lrKXUY4ue+H+xMtPqjxwcQ1oZ52vG7C0CUR+voH/0Sx8eOUZN/qxOXBw6TyFftN0 VO0wE3Ocdzcc2/PTWlajQKfJJnssZt3k/IoGQ1QCmKSXeSr2kmoAiSEAg8gnDrZE MaU5KDnT64mheOzaS6RVs9LnuP4KMbAylvQfysSfw4FIBp8tKw5uc0a9uF/K0Iro lM0Nti1CbsGiEbgRR0i9Y7N/qy4QwZn3fDg/rWlLZIm7dGByMmGV+LrQqT/uf3fD cXsF0QiEdGOdoKkx6wL1LDzqAGEiSuXS7KZDdfOcp8gkA41zL+vBTWQG5miJgAck pk0R2UnCzVSMUkZkmhw2WFBbXlw3zTB8jQ4iqvIYkzDKgIro+MZcrzP0YKlNxQFI ur+1g7BI8gF2WR/WFGXP2lbesDYMQvifJsYw7LCu+HwPspFY2/mtiPu9m4wpZ/SW 4dcom3J8eWyVbVFU4Sq0etmh4Lnjqh7teZEUkkbNtw== =J8gG -----END PGP SIGNATURE-----
commit 9c837c7081e8efb165893a68934119c278b1517c
Author: André Hentschel <n...@dawncrow.de>
Date: Thu Jul 10 21:35:57 2014 +0200
softmath: Add acosh
diff --git a/mingw-w64-crt/Makefile.am b/mingw-w64-crt/Makefile.am
index 65e0df3..175ecdd 100644
--- a/mingw-w64-crt/Makefile.am
+++ b/mingw-w64-crt/Makefile.am
@@ -226,7 +226,7 @@ src_libmingwex=\
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/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/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 \
@@ -295,25 +295,25 @@ src_libmingwex=\
# these only go into the 64 bit version:
src_libmingwex64=\
- math/acosf.c math/acoshf.c math/asinf.c math/asinhf.c math/atan2f.c \
- math/atanf.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/acosf.c math/acosh.c math/acoshf.c math/asinf.c math/asinhf.c \
+ math/atan2f.c math/atanf.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
# these only go into the 32 bit version:
src_libmingwex32=\
- math/acosf.c math/acoshf.c math/asinf.c math/asinhf.c math/atan2f.c \
- math/atanf.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/acosf.c math/acosh.c math/acoshf.c math/asinf.c math/asinhf.c \
+ math/atan2f.c math/atanf.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
# these only go into the ARM32 version:
src_libmingwexarm32=\
math/softmath/e_fmodf.c math/softmath/e_powf.c \
- math/softmath/acosf.c math/softmath/acoshf.c math/softmath/asinf.c math/softmath/asinhf.c math/softmath/atan2f.c \
- math/softmath/atanf.c math/softmath/atanhf.c math/softmath/cosf.c math/softmath/exp2f.c math/softmath/expm1f.c \
- math/softmath/fmodf.c math/softmath/ilogbf.c math/softmath/log1pf.c math/softmath/log2f.c math/softmath/logbf.c \
- math/softmath/scalbnf.c math/softmath/sinf.c math/softmath/tanf.c
+ math/softmath/acosf.c math/softmath/acosh.c math/softmath/acoshf.c math/softmath/asinf.c math/softmath/asinhf.c \
+ math/softmath/atan2f.c math/softmath/atanf.c math/softmath/atanhf.c math/softmath/cosf.c math/softmath/exp2f.c \
+ math/softmath/expm1f.c math/softmath/fmodf.c math/softmath/ilogbf.c math/softmath/log1pf.c math/softmath/log2f.c \
+ math/softmath/logbf.c math/softmath/scalbnf.c math/softmath/sinf.c math/softmath/tanf.c
# These intrinsics are target independent:
diff --git a/mingw-w64-crt/math/softmath/acosh.c b/mingw-w64-crt/math/softmath/acosh.c
new file mode 100644
index 0000000..09bfaa6
--- /dev/null
+++ b/mingw-w64-crt/math/softmath/acosh.c
@@ -0,0 +1,50 @@
+/*
+ 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"
+
+double acosh(double x)
+{
+ return softmath_log(x + sqrt(x * x * 2 - 1));
+}
diff --git a/mingw-w64-crt/math/softmath/bsd_private.h b/mingw-w64-crt/math/softmath/bsd_private.h
new file mode 100644
index 0000000..af01999
--- /dev/null
+++ b/mingw-w64-crt/math/softmath/bsd_private.h
@@ -0,0 +1,117 @@
+/*
+* ====================================================
+* 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 double
+ bp[] = {1.0, 1.5,},
+ dp_h[] = { 0.0, 5.84962487220764160156e-01,}, /* 0x3FE2B803, 0x40000000 */
+ dp_l[] = { 0.0, 1.35003920212974897128e-08,}, /* 0x3E4CFDEB, 0x43CFD006 */
+ Zero[] = {0.0, -0.0,},
+ zero = 0.0,
+ one = 1.0,
+ two = 2.0,
+ two53 = 9007199254740992.0, /* 0x43400000, 0x00000000 */
+ huge = 1.0e300,
+ tiny = 1.0e-300,
+ /* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */
+ L1 = 5.99999999999994648725e-01, /* 0x3FE33333, 0x33333303 */
+ L2 = 4.28571428578550184252e-01, /* 0x3FDB6DB6, 0xDB6FABFF */
+ L3 = 3.33333329818377432918e-01, /* 0x3FD55555, 0x518F264D */
+ L4 = 2.72728123808534006489e-01, /* 0x3FD17460, 0xA91D4101 */
+ L5 = 2.30660745775561754067e-01, /* 0x3FCD864A, 0x93C9DB65 */
+ L6 = 2.06975017800338417784e-01, /* 0x3FCA7E28, 0x4A454EEF */
+ P1 = 1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */
+ P2 = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */
+ P3 = 6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */
+ P4 = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */
+ P5 = 4.13813679705723846039e-08, /* 0x3E663769, 0x72BEA4D0 */
+ lg2 = 6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */
+ lg2_h = 6.93147182464599609375e-01, /* 0x3FE62E43, 0x00000000 */
+ lg2_l = -1.90465429995776804525e-09, /* 0xBE205C61, 0x0CA86C39 */
+ ovt = 8.0085662595372944372e-0017, /* -(1024-log2(ovfl+.5ulp)) */
+ cp = 9.61796693925975554329e-01, /* 0x3FEEC709, 0xDC3A03FD =2/(3ln2) */
+ cp_h = 9.61796700954437255859e-01, /* 0x3FEEC709, 0xE0000000 =(float)cp */
+ cp_l = -7.02846165095275826516e-09, /* 0xBE3E2FE0, 0x145B01F5 =tail of cp_h*/
+ ivln2 = 1.44269504088896338700e+00, /* 0x3FF71547, 0x652B82FE =1/ln2 */
+ ivln2_h= 1.44269502162933349609e+00, /* 0x3FF71547, 0x60000000 =24b 1/ln2*/
+ ivln2_l= 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
+
+typedef unsigned int u_int32_t;
+
+typedef union
+{
+ double value;
+ struct
+ {
+ u_int32_t lsw;
+ u_int32_t msw;
+ } parts;
+} ieee_double_shape_type;
+
+/* Get two 32 bit ints from a double. */
+
+#define EXTRACT_WORDS(ix0,ix1,d) \
+do { \
+ ieee_double_shape_type ew_u; \
+ ew_u.value = (d); \
+ (ix0) = ew_u.parts.msw; \
+ (ix1) = ew_u.parts.lsw; \
+} while (0)
+
+/* Get the most significant 32 bit int from a double. */
+
+#define GET_HIGH_WORD(i,d) \
+do { \
+ ieee_double_shape_type gh_u; \
+ gh_u.value = (d); \
+ (i) = gh_u.parts.msw; \
+} while (0)
+
+/* Get the less significant 32 bit int from a double. */
+
+#define GET_LOW_WORD(i,d) \
+do { \
+ ieee_double_shape_type gl_u; \
+ gl_u.value = (d); \
+ (i) = gl_u.parts.lsw; \
+} while (0)
+
+/* Set a double from two 32 bit ints. */
+
+#define INSERT_WORDS(d,ix0,ix1) \
+do { \
+ ieee_double_shape_type iw_u; \
+ iw_u.parts.msw = (ix0); \
+ iw_u.parts.lsw = (ix1); \
+ (d) = iw_u.value; \
+} while (0)
+
+/* Set the more significant 32 bits of a double from an int. */
+
+#define SET_HIGH_WORD(d,v) \
+do { \
+ ieee_double_shape_type sh_u; \
+ sh_u.value = (d); \
+ sh_u.parts.msw = (v); \
+ (d) = sh_u.value; \
+} while (0)
+
+/* Set the less significant 32 bits of a double from an int. */
+
+#define SET_LOW_WORD(d,v) \
+do { \
+ ieee_double_shape_type sl_u; \
+ sl_u.value = (d); \
+ sl_u.parts.lsw = (v); \
+ (d) = sl_u.value; \
+} while (0)
diff --git a/mingw-w64-crt/math/softmath/bsd_privatef.h b/mingw-w64-crt/math/softmath/bsd_privatef.h
new file mode 100644
index 0000000..66bb415
--- /dev/null
+++ b/mingw-w64-crt/math/softmath/bsd_privatef.h
@@ -0,0 +1,68 @@
+/*
+* ====================================================
+* 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)
diff --git a/mingw-w64-crt/math/softmath/e_fmodf.c b/mingw-w64-crt/math/softmath/e_fmodf.c
index c57daf4..6a4eaf3 100644
--- a/mingw-w64-crt/math/softmath/e_fmodf.c
+++ b/mingw-w64-crt/math/softmath/e_fmodf.c
@@ -16,7 +16,7 @@
*/
#include "math.h"
-#include "softmath_private.h"
+#include "bsd_privatef.h"
float bsd__ieee754_fmodf(float x, float y)
{
diff --git a/mingw-w64-crt/math/softmath/e_pow.c b/mingw-w64-crt/math/softmath/e_pow.c
new file mode 100644
index 0000000..97fb44a
--- /dev/null
+++ b/mingw-w64-crt/math/softmath/e_pow.c
@@ -0,0 +1,223 @@
+/*
+ * ====================================================
+ * Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+
+#include "math.h"
+#include "bsd_private.h"
+
+double bsd__ieee754_pow(double x, double y)
+{
+ double z,ax,z_h,z_l,p_h,p_l;
+ double y1,t1,t2,r,s,t,u,v,w;
+ int32_t i,j,k,yisint,n;
+ int32_t hx,hy,ix,iy;
+ u_int32_t lx,ly;
+
+ EXTRACT_WORDS(hx,lx,x);
+ EXTRACT_WORDS(hy,ly,y);
+ ix = hx&0x7fffffff; iy = hy&0x7fffffff;
+
+ /* y==zero: x**0 = 1 */
+ if((iy|ly)==0) return one;
+
+ /* x==1: 1**y = 1, even if y is NaN */
+ if (hx==0x3ff00000 && lx == 0) return one;
+
+ /* y!=zero: result is NaN if either arg is NaN */
+ if(ix > 0x7ff00000 || ((ix==0x7ff00000)&&(lx!=0)) ||
+ iy > 0x7ff00000 || ((iy==0x7ff00000)&&(ly!=0)))
+ return (x+0.0)+(y+0.0);
+
+ /* 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>=0x43400000) yisint = 2; /* even integer y */
+ else if (iy>=0x3ff00000) {
+ k = (iy>>20)-0x3ff; /* exponent */
+ if(k>20) {
+ j = ly>>(52-k);
+ if((j<<(52-k))==ly) yisint = 2-(j&1);
+ } else if(ly==0) {
+ j = iy>>(20-k);
+ if((j<<(20-k))==iy) yisint = 2-(j&1);
+ }
+ }
+ }
+
+ /* special value of y */
+ if(ly==0) {
+ if (iy==0x7ff00000) { /* y is +-inf */
+ if(((ix-0x3ff00000)|lx)==0)
+ return one; /* (-1)**+-inf is NaN */
+ else if (ix >= 0x3ff00000)/* (|x|>1)**+-inf = inf,0 */
+ return (hy>=0)? y: zero;
+ else /* (|x|<1)**-,+inf = inf,0 */
+ return (hy<0)?-y: zero;
+ }
+ if(iy==0x3ff00000) { /* y is +-1 */
+ if(hy<0) return one/x; else return x;
+ }
+ if(hy==0x40000000) return x*x; /* y is 2 */
+ if(hy==0x3fe00000) { /* y is 0.5 */
+ if(hx>=0) /* x >= +0 */
+ return sqrt(x);
+ }
+ }
+
+ ax = fabs(x);
+ /* special value of x */
+ if(lx==0) {
+ if(ix==0x7ff00000||ix==0||ix==0x3ff00000){
+ z = ax; /*x is +-0,+-inf,+-1*/
+ if(hy<0) z = one/z; /* z = (1/|x|) */
+ if(hx<0) {
+ if(((ix-0x3ff00000)|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;
+ }
+ }
+
+ /* CYGNUS LOCAL + fdlibm-5.3 fix: This used to be
+ n = (hx>>31)+1;
+ but ANSI C says a right shift of a signed negative quantity is
+ implementation defined. */
+ n = ((u_int32_t)hx>>31)-1;
+
+ /* (x<0)**(non-int) is NaN */
+ if((n|yisint)==0) return (x-x)/(x-x);
+
+ s = one; /* s (sign of result -ve**odd) = -1 else = 1 */
+ if((n|(yisint-1))==0) s = -one;/* (-ve)**(odd int) */
+
+ /* |y| is huge */
+ if(iy>0x41e00000) { /* if |y| > 2**31 */
+ if(iy>0x43f00000){ /* if |y| > 2**64, must o/uflow */
+ if(ix<=0x3fefffff) return (hy<0)? huge*huge:tiny*tiny;
+ if(ix>=0x3ff00000) return (hy>0)? huge*huge:tiny*tiny;
+ }
+ /* over/underflow if x is not close to one */
+ if(ix<0x3fefffff) return (hy<0)? s*huge*huge:s*tiny*tiny;
+ if(ix>0x3ff00000) return (hy>0)? s*huge*huge:s*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-one; /* t has 20 trailing zeros */
+ w = (t*t)*(0.5-t*(0.3333333333333333333333-t*0.25));
+ u = ivln2_h*t; /* ivln2_h has 21 sig. bits */
+ v = t*ivln2_l-w*ivln2;
+ t1 = u+v;
+ SET_LOW_WORD(t1,0);
+ t2 = v-(t1-u);
+ } else {
+ double ss,s2,s_h,s_l,t_h,t_l;
+ n = 0;
+ /* take care subnormal number */
+ if(ix<0x00100000)
+ {ax *= two53; n -= 53; GET_HIGH_WORD(ix,ax); }
+ n += ((ix)>>20)-0x3ff;
+ j = ix&0x000fffff;
+ /* determine interval */
+ ix = j|0x3ff00000; /* normalize ix */
+ if(j<=0x3988E) k=0; /* |x|<sqrt(3/2) */
+ else if(j<0xBB67A) k=1; /* |x|<sqrt(3) */
+ else {k=0;n+=1;ix -= 0x00100000;}
+ SET_HIGH_WORD(ax,ix);
+
+ /* compute ss = 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]);
+ ss = u*v;
+ s_h = ss;
+ SET_LOW_WORD(s_h,0);
+ /* t_h=ax+bp[k] High */
+ t_h = zero;
+ SET_HIGH_WORD(t_h,((ix>>1)|0x20000000)+0x00080000+(k<<18));
+ t_l = ax - (t_h-bp[k]);
+ s_l = v*((u-s_h*t_h)-s_h*t_l);
+ /* compute log(ax) */
+ s2 = ss*ss;
+ r = s2*s2*(L1+s2*(L2+s2*(L3+s2*(L4+s2*(L5+s2*L6)))));
+ r += s_l*(s_h+ss);
+ s2 = s_h*s_h;
+ t_h = 3.0+s2+r;
+ SET_LOW_WORD(t_h,0);
+ t_l = r-((t_h-3.0)-s2);
+ /* u+v = ss*(1+...) */
+ u = s_h*t_h;
+ v = s_l*t_h+t_l*ss;
+ /* 2/(3log2)*(ss+...) */
+ p_h = u+v;
+ SET_LOW_WORD(p_h,0);
+ 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) = (ss+..)*2/(3*log2) = n + dp_h + z_h + z_l */
+ t = (double)n;
+ t1 = (((z_h+z_l)+dp_h[k])+t);
+ SET_LOW_WORD(t1,0);
+ t2 = z_l-(((t1-t)-dp_h[k])-z_h);
+ }
+
+ /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
+ y1 = y;
+ SET_LOW_WORD(y1,0);
+ p_l = (y-y1)*t1+y*t2;
+ p_h = y1*t1;
+ z = p_l+p_h;
+ EXTRACT_WORDS(j,i,z);
+ if (j>=0x40900000) { /* z >= 1024 */
+ if(((j-0x40900000)|i)!=0) /* if z > 1024 */
+ return s*huge*huge; /* overflow */
+ else {
+ if(p_l+ovt>z-p_h) return s*huge*huge; /* overflow */
+ }
+ } else if((j&0x7fffffff)>=0x4090cc00 ) { /* z <= -1075 */
+ if(((j-0xc090cc00)|i)!=0) /* z < -1075 */
+ return s*tiny*tiny; /* underflow */
+ else {
+ if(p_l<=z-p_h) return s*tiny*tiny; /* underflow */
+ }
+ }
+ /*
+ * compute 2**(p_h+p_l)
+ */
+ i = j&0x7fffffff;
+ k = (i>>20)-0x3ff;
+ n = 0;
+ if (i>0x3fe00000) { /* if |z| > 0.5, set n = [z+0.5] */
+ n = j+(0x00100000>>(k+1));
+ k = ((n&0x7fffffff)>>20)-0x3ff; /* new k for n */
+ t = zero;
+ SET_HIGH_WORD(t,n&~(0x000fffff>>k));
+ n = ((n&0x000fffff)|0x00100000)>>(20-k);
+ if(j<0) n = -n;
+ p_h -= t;
+ }
+ t = p_l+p_h;
+ SET_LOW_WORD(t,0);
+ 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_HIGH_WORD(j,z);
+ j += (n<<20);
+ if((j>>20)<=0) z = scalbn(z,n); /* subnormal output */
+ else SET_HIGH_WORD(z,j);
+ return s*z;
+}
diff --git a/mingw-w64-crt/math/softmath/e_powf.c b/mingw-w64-crt/math/softmath/e_powf.c
index faaf3dd..01ee105 100644
--- a/mingw-w64-crt/math/softmath/e_powf.c
+++ b/mingw-w64-crt/math/softmath/e_powf.c
@@ -10,7 +10,7 @@
*/
#include "math.h"
-#include "softmath_private.h"
+#include "bsd_privatef.h"
float bsd__ieee754_powf(float x, float y)
{
diff --git a/mingw-w64-crt/math/softmath/softmath_private.h b/mingw-w64-crt/math/softmath/softmath_private.h
index b2a2d3d..1d8476f 100644
--- a/mingw-w64-crt/math/softmath/softmath_private.h
+++ b/mingw-w64-crt/math/softmath/softmath_private.h
@@ -12,62 +12,8 @@
#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 double bsd__ieee754_pow(double, double);
extern float bsd__ieee754_powf(float, float);
static inline double softmath_fact(double number)
@@ -93,6 +39,27 @@ static inline float softmath_expf(float x)
return result;
}
+static inline double softmath_log(double x)
+{
+ int n, aprox = 8 / (x / 2);
+ double 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_pow((x - 1.0) / (x + 1.0), 2 * n + 1) * (1.0 / (2.0 * n + 1.0));
+ if (isinf(result))
+ break;
+ }
+ result *= 2;
+
+ return result;
+}
+
static inline float softmath_logf(float x)
{
int n, aprox = 8 / (x / 2);
softmath.patch.sig
Description: PGP signature
------------------------------------------------------------------------------ Open source business process management suite built on Java and Eclipse Turn processes into business applications with Bonita BPM Community Edition Quickly connect people, data, and systems into organized workflows Winner of BOSSIE, CODIE, OW2 and Gartner awards http://p.sf.net/sfu/Bonitasoft
_______________________________________________ Mingw-w64-public mailing list Mingw-w64-public@lists.sourceforge.net https://lists.sourceforge.net/lists/listinfo/mingw-w64-public
