2009-08-23 Ozkan Sezer <sezeroz@gmail.com>
* math/*.c: Many whitespace/indentation tidy-ups.
Replaced K&R declarations with ANSI ones. Fixed
double semicolons.
git-svn-id: svn+ssh://svn.code.sf.net/p/mingw-w64/code/trunk@1280 4407c894-4637-0410-b4f5-ada5f102cad1
diff --git a/mingw-w64-crt/math/cbrt.c b/mingw-w64-crt/math/cbrt.c
index 290bf62..56e9e46 100644
--- a/mingw-w64-crt/math/cbrt.c
+++ b/mingw-w64-crt/math/cbrt.c
@@ -12,104 +12,93 @@
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;
+double cbrt(double x)
{
-int e, rem, sign;
-double z;
+ int e, rem, sign;
+ double z;
#ifdef __MINGW32__
-if (!isfinite (x) || x == 0 )
- return x;
+ if (!isfinite (x) || x == 0)
+ return x;
#else
-
#ifdef NANS
-if( isnan(x) )
- return x;
+ if (isnan(x))
+ return x;
#endif
#ifdef INFINITIES
-if( !isfinite(x) )
- return x;
+ if (!isfinite(x))
+ return x;
#endif
-if( x == 0 )
- return( x );
-
+ if (x == 0)
+ return (x);
#endif /* __MINGW32__ */
-if( x > 0 )
- sign = 1;
-else
+ if (x > 0)
+ sign = 1;
+ else
{
- sign = -1;
- x = -x;
+ sign = -1;
+ x = -x;
}
-z = x;
-/* extract power of 2, leaving
- * mantissa between 0.5 and 1
- */
-x = frexp( x, &e );
+ 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;
+ /* 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 )
+ /* 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;
+ 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);
-/* 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;
+ /* Newton iteration */
+ x -= ( x - (z/(x*x)) )*0.33333333333333333333;
#ifdef DEC
-x -= ( x - (z/(x*x)) )/3.0;
+ x -= ( x - (z/(x*x)) )/3.0;
#else
-x -= ( x - (z/(x*x)) )*0.33333333333333333333;
+ x -= ( x - (z/(x*x)) )*0.33333333333333333333;
#endif
-if( sign < 0 )
- x = -x;
-return(x);
+ if (sign < 0)
+ x = -x;
+ return (x);
}
diff --git a/mingw-w64-crt/math/cbrtf.c b/mingw-w64-crt/math/cbrtf.c
index cb3d59b..03c285e 100644
--- a/mingw-w64-crt/math/cbrtf.c
+++ b/mingw-w64-crt/math/cbrtf.c
@@ -11,68 +11,65 @@
float cbrtf (float x)
{
-int e, rem, sign;
-float z;
-if (!isfinite (x) || x == 0.0F )
- return x;
-if( x > 0 )
- sign = 1;
-else
+ int e, rem, sign;
+ float z;
+ if (!isfinite (x) || x == 0.0F)
+ return x;
+ if (x > 0)
+ sign = 1;
+ else
{
- sign = -1;
- x = -x;
+ sign = -1;
+ x = -x;
}
-z = x;
-/* extract power of 2, leaving
- * mantissa between 0.5 and 1
- */
-x = frexpf( x, &e );
+ 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;
+ /* 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 )
+ /* 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;
+ rem = e;
+ e /= 3;
+ rem -= 3*e;
+ if (rem == 1)
+ x *= CBRT2;
+ else if (rem == 2)
+ x *= CBRT4;
}
-
-
/* argument less than 1 */
-
-else
+ else
{
- e = -e;
- rem = e;
- e /= 3;
- rem -= 3*e;
- if( rem == 1 )
- x /= CBRT2;
- else if( rem == 2 )
- x /= CBRT4;
- e = -e;
+ 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 );
+ /* multiply by power of 2 */
+ x = ldexpf(x, e);
-/* Newton iteration */
-x -= ( x - (z/(x*x)) ) * 0.333333333333;
+ /* Newton iteration */
+ x -= ( x - (z/(x*x)) ) * 0.333333333333;
-if( sign < 0 )
- x = -x;
-return(x);
+ if (sign < 0)
+ x = -x;
+ return (x);
}
diff --git a/mingw-w64-crt/math/cbrtl.c b/mingw-w64-crt/math/cbrtl.c
index dc7bae4..3841032 100644
--- a/mingw-w64-crt/math/cbrtl.c
+++ b/mingw-w64-crt/math/cbrtl.c
@@ -12,72 +12,71 @@
extern long double ldexpl(long double,int);
-long double cbrtl(x)
-long double x;
+long double cbrtl(long double x)
{
-int e, rem, sign;
-long double z;
+ int e, rem, sign;
+ long double z;
-if (!isfinite (x) || x == 0.0L)
- return(x);
+ if (!isfinite (x) || x == 0.0L)
+ return (x);
-if( x > 0 )
- sign = 1;
-else
+ if (x > 0)
+ sign = 1;
+ else
{
- sign = -1;
- x = -x;
+ sign = -1;
+ x = -x;
}
-z = x;
-/* extract power of 2, leaving
- * mantissa between 0.5 and 1
- */
-x = frexpl( x, &e );
+ 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;
+ /* 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 )
+ /* 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;
+ rem = e;
+ e /= 3;
+ rem -= 3*e;
+ if (rem == 1)
+ x *= CBRT2;
+ else if (rem == 2)
+ x *= CBRT4;
}
-else
+ 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;
+ 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 );
+ /* multiply by power of 2 */
+ x = ldexpl(x, e);
-/* Newton iteration */
+ /* Newton iteration */
-x -= ( x - (z/(x*x)) )*0.3333333333333333333333L;
-x -= ( x - (z/(x*x)) )*0.3333333333333333333333L;
+ x -= ( x - (z/(x*x)) )*0.3333333333333333333333L;
+ x -= ( x - (z/(x*x)) )*0.3333333333333333333333L;
-if( sign < 0 )
- x = -x;
-return(x);
+ if (sign < 0)
+ x = -x;
+ return (x);
}
diff --git a/mingw-w64-crt/math/coshf.c b/mingw-w64-crt/math/coshf.c
index 32e3534..14ebf0b 100644
--- a/mingw-w64-crt/math/coshf.c
+++ b/mingw-w64-crt/math/coshf.c
@@ -5,4 +5,6 @@
*/
#include <math.h>
float coshf (float x)
- {return (float) cosh (x);}
+{
+ return (float) cosh (x);
+}
diff --git a/mingw-w64-crt/math/coshl.c b/mingw-w64-crt/math/coshl.c
index c3fb112..fee84e5 100644
--- a/mingw-w64-crt/math/coshl.c
+++ b/mingw-w64-crt/math/coshl.c
@@ -9,37 +9,36 @@
#define _SET_ERRNO(x)
#endif
-long double coshl(x)
-long double x;
+long double coshl(long double x)
{
-long double y;
+ long double y;
#ifdef NANS
-if( isnanl(x) )
+ if (isnanl(x))
{
- _SET_ERRNO(EDOM);
- return(x);
+ _SET_ERRNO(EDOM);
+ return (x);
}
#endif
-if( x < 0 )
- x = -x;
-if( x > (MAXLOGL + LOGE2L) )
+ if (x < 0)
+ x = -x;
+ if (x > (MAXLOGL + LOGE2L))
{
- mtherr( "coshl", OVERFLOW );
- _SET_ERRNO(ERANGE);
+ mtherr( "coshl", OVERFLOW );
+ _SET_ERRNO(ERANGE);
#ifdef INFINITIES
- return( INFINITYL );
+ return (INFINITYL);
#else
- return( MAXNUML );
+ 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 );
+ 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/mingw-w64-crt/math/erfl.c b/mingw-w64-crt/math/erfl.c
index 165062f..61d39c8 100644
--- a/mingw-w64-crt/math/erfl.c
+++ b/mingw-w64-crt/math/erfl.c
@@ -44,7 +44,7 @@
* IEEE 0,1 50000 2.0e-19 5.7e-20
*
*/
-
+
/* erfcl.c
*
* Complementary error function
@@ -95,7 +95,7 @@
*
*
*/
-
+
/*
Modified from file ndtrl.c
@@ -215,89 +215,89 @@
static long double expx2l (long double x)
{
- long double u, u1, m, f;
+ long double u, u1, m, f;
- x = fabsl (x);
- /* Represent x as an exact multiple of M plus a residual.
- M is a power of 2 chosen so that exp(m * m) does not overflow
- or underflow and so that |x - m| is small. */
- m = MINV * floorl(M * x + 0.5L);
- f = x - m;
+ x = fabsl (x);
+ /* Represent x as an exact multiple of M plus a residual.
+ M is a power of 2 chosen so that exp(m * m) does not overflow
+ or underflow and so that |x - m| is small. */
+ m = MINV * floorl(M * x + 0.5L);
+ f = x - m;
- /* x^2 = m^2 + 2mf + f^2 */
- u = m * m;
- u1 = 2 * m * f + f * f;
+ /* x^2 = m^2 + 2mf + f^2 */
+ u = m * m;
+ u1 = 2 * m * f + f * f;
- if ((u+u1) > MAXLOGL)
- return (INFINITYL);
+ if ((u + u1) > MAXLOGL)
+ return (INFINITYL);
- /* u is exact, u1 is small. */
- u = expl(u) * expl(u1);
- return(u);
+ /* u is exact, u1 is small. */
+ u = expl(u) * expl(u1);
+ return (u);
}
long double erfcl(long double a)
{
-long double p,q,x,y,z;
+ long double p, q, x, y, z;
-if (isinf (a))
- return (signbit (a) ? 2.0 : 0.0);
+ if (isinf (a))
+ return (signbit(a) ? 2.0 : 0.0);
-x = fabsl (a);
+ x = fabsl (a);
-if (x < 1.0L)
- return (1.0L - erfl(a));
+ if (x < 1.0L)
+ return (1.0L - erfl(a));
-z = a * a;
+ z = a * a;
-if( z > MAXLOGL )
+ if (z > MAXLOGL)
{
under:
- mtherr( "erfcl", UNDERFLOW );
- errno = ERANGE;
- return (signbit (a) ? 2.0 : 0.0);
+ mtherr("erfcl", UNDERFLOW);
+ errno = ERANGE;
+ return (signbit(a) ? 2.0 : 0.0);
}
-/* Compute z = expl(a * a). */
-z = expx2l (a);
-y = 1.0L/x;
+ /* Compute z = expl(a * a). */
+ z = expx2l(a);
+ y = 1.0L/x;
-if (x < 8.0L)
+ if (x < 8.0L)
{
- p = polevll (y, P, 9);
- q = p1evll (y, Q, 10);
+ p = polevll(y, P, 9);
+ q = p1evll(y, Q, 10);
}
-else
+ else
{
- q = y * y;
- p = y * polevll (q, R, 4);
- q = p1evll (q, S, 5);
+ q = y * y;
+ p = y * polevll(q, R, 4);
+ q = p1evll(q, S, 5);
}
-y = p/(q * z);
+ y = p/(q * z);
-if (a < 0.0L)
- y = 2.0L - y;
+ if (a < 0.0L)
+ y = 2.0L - y;
-if (y == 0.0L)
- goto under;
+ if (y == 0.0L)
+ goto under;
-return (y);
+ return (y);
}
long double erfl(long double x)
{
-long double y, z;
+ long double y, z;
-if( x == 0.0L )
- return (x);
+ if (x == 0.0L)
+ return (x);
-if (isinf (x))
- return (signbit (x) ? -1.0L : 1.0L);
+ if (isinf (x))
+ return (signbit(x) ? -1.0L : 1.0L);
-if (fabsl(x) > 1.0L)
- return (1.0L - erfcl (x));
+ if (fabsl(x) > 1.0L)
+ return (1.0L - erfcl(x));
-z = x * x;
-y = x * polevll( z, T, 6 ) / p1evll( z, U, 6 );
-return( y );
+ z = x * x;
+ y = x * polevll(z, T, 6) / p1evll(z, U, 6);
+ return (y);
}
diff --git a/mingw-w64-crt/math/expf.c b/mingw-w64-crt/math/expf.c
index 593050a..926dae8 100644
--- a/mingw-w64-crt/math/expf.c
+++ b/mingw-w64-crt/math/expf.c
@@ -5,4 +5,6 @@
*/
#include <math.h>
float expf (float x)
- {return (float) exp (x);}
+{
+ return (float) exp (x);
+}
diff --git a/mingw-w64-crt/math/fp_consts.c b/mingw-w64-crt/math/fp_consts.c
index eea647b..bae4167 100644
--- a/mingw-w64-crt/math/fp_consts.c
+++ b/mingw-w64-crt/math/fp_consts.c
@@ -6,7 +6,7 @@
#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 __INF = { __DOUBLE_INF_REP };
const union _ieee_rep __DENORM = { __DOUBLE_DENORM_REP };
/* ISO C99 */
@@ -14,6 +14,7 @@
/* FIXME */
double nan (const char *);
double nan (const char * tagp __attribute__((unused)) )
- { return __QNAN.double_val; }
-
+{
+ return __QNAN.double_val;
+}
diff --git a/mingw-w64-crt/math/fp_consts.h b/mingw-w64-crt/math/fp_consts.h
index fa8305e..5fe2cc7 100644
--- a/mingw-w64-crt/math/fp_consts.h
+++ b/mingw-w64-crt/math/fp_consts.h
@@ -1,53 +1,54 @@
-/**
- * This file has no copyright assigned and is placed in the Public Domain.
- * This file is part of the w64 mingw-runtime package.
- * No warranty is given; refer to the file DISCLAIMER within this package.
- */
-#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
+/**
+ * This file has no copyright assigned and is placed in the Public Domain.
+ * This file is part of the w64 mingw-runtime package.
+ * No warranty is given; refer to the file DISCLAIMER within this package.
+ */
+#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 /* _FP_CONSTS_H */
+
diff --git a/mingw-w64-crt/math/fp_constsf.c b/mingw-w64-crt/math/fp_constsf.c
index 7f1f847..1f88b47 100644
--- a/mingw-w64-crt/math/fp_constsf.c
+++ b/mingw-w64-crt/math/fp_constsf.c
@@ -5,10 +5,10 @@
*/
#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 };
+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
@@ -16,4 +16,7 @@
float nanf(const char *);
float nanf(const char * tagp __attribute__((unused)) )
- { return __QNANF.float_val;}
+{
+ return __QNANF.float_val;
+}
+
diff --git a/mingw-w64-crt/math/fp_constsl.c b/mingw-w64-crt/math/fp_constsl.c
index 06d23ef..63fbbab 100644
--- a/mingw-w64-crt/math/fp_constsl.c
+++ b/mingw-w64-crt/math/fp_constsl.c
@@ -6,13 +6,15 @@
#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 __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 *);
long double nanl (const char * tagp __attribute__((unused)) )
- { return __QNANL.ldouble_val; }
+{
+ return __QNANL.ldouble_val;
+}
+
diff --git a/mingw-w64-crt/math/frexpf.c b/mingw-w64-crt/math/frexpf.c
index b3ce1f8..4c4f189 100644
--- a/mingw-w64-crt/math/frexpf.c
+++ b/mingw-w64-crt/math/frexpf.c
@@ -4,7 +4,10 @@
* No warranty is given; refer to the file DISCLAIMER within this package.
*/
extern double __cdecl frexp(double _X,int *_Y);
-float frexpf (float, int *);
+float frexpf (float, int *);
float frexpf (float x, int *expn)
- {return (float)frexp(x, expn);}
+{
+ return (float)frexp(x, expn);
+}
+
diff --git a/mingw-w64-crt/math/hypotf.c b/mingw-w64-crt/math/hypotf.c
index 9bbd2a7..a5e7f37 100644
--- a/mingw-w64-crt/math/hypotf.c
+++ b/mingw-w64-crt/math/hypotf.c
@@ -6,4 +6,7 @@
#include <math.h>
float hypotf (float x, float y)
- { return (float) _hypot (x, y);}
+{
+ return (float) _hypot (x, y);
+}
+
diff --git a/mingw-w64-crt/math/ldexpf.c b/mingw-w64-crt/math/ldexpf.c
index 67ecbe6..0510fd5 100644
--- a/mingw-w64-crt/math/ldexpf.c
+++ b/mingw-w64-crt/math/ldexpf.c
@@ -4,7 +4,10 @@
* No warranty is given; refer to the file DISCLAIMER within this package.
*/
extern double __cdecl ldexp(double _X,int _Y);
-float ldexpf (float x, int expn);
+float ldexpf (float x, int expn);
float ldexpf (float x, int expn)
- {return (float) ldexp (x, expn);}
+{
+ return (float) ldexp (x, expn);
+}
+
diff --git a/mingw-w64-crt/math/lgamma.c b/mingw-w64-crt/math/lgamma.c
index c8f4cf8..04e8097 100644
--- a/mingw-w64-crt/math/lgamma.c
+++ b/mingw-w64-crt/math/lgamma.c
@@ -10,11 +10,11 @@
*/
#ifdef UNK
static uD A[] = {
- { { 8.11614167470508450300E-4 } },
+ { { 8.11614167470508450300E-4 } },
{ { -5.95061904284301438324E-4 } },
- { { 7.93650340457716943945E-4 } },
+ { { 7.93650340457716943945E-4 } },
{ { -2.77777777730099687205E-3 } },
- { { 8.33333333333331927722E-2 } }
+ { { 8.33333333333331927722E-2 } }
};
static uD B[] = {
{ { -1.37825152569120859100E3 } },
@@ -160,118 +160,118 @@
double __lgamma_r(double x, int* sgngam)
{
-double p, q, u, w, z;
-int i;
+ double p, q, u, w, z;
+ int i;
-*sgngam = 1;
+ *sgngam = 1;
#ifdef NANS
-if( isnan(x) )
- return(x);
+ if (isnan(x))
+ return (x);
#endif
#ifdef INFINITIES
-if( !isfinite(x) )
- return(INFINITY);
+ if (!isfinite(x))
+ return (INFINITY);
#endif
-if( x < -34.0 )
+ if (x < -34.0)
{
- q = -x;
- w = __lgamma_r(q, sgngam); /* note this modifies sgngam! */
- p = floor(q);
- if( p == q )
+ q = -x;
+ w = __lgamma_r(q, sgngam); /* note this modifies sgngam! */
+ p = floor(q);
+ if (p == q)
{
lgsing:
- _SET_ERRNO(EDOM);
- mtherr( "lgam", SING );
+ _SET_ERRNO(EDOM);
+ mtherr( "lgam", SING );
#ifdef INFINITIES
- return (INFINITY);
+ return (INFINITY);
#else
- return (MAXNUM);
+ return (MAXNUM);
#endif
}
- i = p;
- if( (i & 1) == 0 )
- *sgngam = -1;
- else
- *sgngam = 1;
- z = q - p;
- if( z > 0.5 )
+ i = p;
+ if ((i & 1) == 0)
+ *sgngam = -1;
+ else
+ *sgngam = 1;
+ z = q - p;
+ if (z > 0.5)
{
- p += 1.0;
- z = p - q;
+ p += 1.0;
+ z = p - q;
}
- z = q * sin( PI * z );
- if( z == 0.0 )
- goto lgsing;
-/* z = log(PI) - log( z ) - w;*/
- z = LOGPI - log( z ) - w;
- return( z );
- }
-
-if( x < 13.0 )
- {
- z = 1.0;
- p = 0.0;
- u = x;
- while( u >= 3.0 )
- {
- p -= 1.0;
- u = x + p;
- z *= u;
- }
- while( u < 2.0 )
- {
- if( u == 0.0 )
+ z = q * sin( PI * z );
+ if (z == 0.0)
goto lgsing;
- z /= u;
- p += 1.0;
- u = x + p;
- }
- if( z < 0.0 )
- {
- *sgngam = -1;
- z = -z;
- }
- else
- *sgngam = 1;
- if( u == 2.0 )
- return( log(z) );
- p -= 2.0;
- x = x + p;
- p = x * polevl( x, B, 5 ) / p1evl( x, C, 6);
- return( log(z) + p );
+ /* z = log(PI) - log( z ) - w;*/
+ z = LOGPI - log( z ) - w;
+ return (z);
}
-if( x > MAXLGM )
+ if (x < 13.0)
{
- _SET_ERRNO(ERANGE);
- mtherr( "lgamma", OVERFLOW );
+ z = 1.0;
+ p = 0.0;
+ u = x;
+ while (u >= 3.0)
+ {
+ p -= 1.0;
+ u = x + p;
+ z *= u;
+ }
+ while (u < 2.0)
+ {
+ if (u == 0.0)
+ goto lgsing;
+ z /= u;
+ p += 1.0;
+ u = x + p;
+ }
+ if (z < 0.0)
+ {
+ *sgngam = -1;
+ z = -z;
+ }
+ else
+ *sgngam = 1;
+ if (u == 2.0)
+ return ( log(z) );
+ p -= 2.0;
+ x = x + p;
+ p = x * polevl(x, B, 5) / p1evl(x, C, 6);
+ return ( log(z) + p );
+ }
+
+ if (x > MAXLGM)
+ {
+ _SET_ERRNO(ERANGE);
+ mtherr("lgamma", OVERFLOW);
#ifdef INFINITIES
- return( *sgngam * INFINITY );
+ return (*sgngam * INFINITY);
#else
- return( *sgngam * MAXNUM );
+ return (*sgngam * MAXNUM);
#endif
}
-q = ( x - 0.5 ) * log(x) - x + LS2PI;
-if( x > 1.0e8 )
- return( q );
+ q = (x - 0.5) * log(x) - x + LS2PI;
+ if (x > 1.0e8)
+ return (q);
-p = 1.0/(x*x);
-if( x >= 1000.0 )
- q += (( 7.9365079365079365079365e-4 * p
- - 2.7777777777777777777778e-3) *p
- + 0.0833333333333333333333) / x;
-else
- q += polevl( p, A, 4 ) / x;
-return( q );
+ p = 1.0/(x*x);
+ if (x >= 1000.0)
+ q += (( 7.9365079365079365079365e-4 * p
+ - 2.7777777777777777777778e-3) *p
+ + 0.0833333333333333333333) / x;
+ else
+ q += polevl( p, A, 4 ) / x;
+ return (q);
}
/* This is the C99 version */
-
double lgamma(double x)
{
- int local_sgngam=0;
- return (__lgamma_r(x, &local_sgngam));
+ int local_sgngam = 0;
+ return (__lgamma_r(x, &local_sgngam));
}
+
diff --git a/mingw-w64-crt/math/lgammaf.c b/mingw-w64-crt/math/lgammaf.c
index a3b187f..042f8d1 100644
--- a/mingw-w64-crt/math/lgammaf.c
+++ b/mingw-w64-crt/math/lgammaf.c
@@ -5,26 +5,26 @@
*/
/* log gamma(x+2), -.5 < x < .5 */
static const float B[] = {
- 6.055172732649237E-004,
--1.311620815545743E-003,
- 2.863437556468661E-003,
--7.366775108654962E-003,
- 2.058355474821512E-002,
--6.735323259371034E-002,
- 3.224669577325661E-001,
- 4.227843421859038E-001
+ 6.055172732649237E-004,
+ -1.311620815545743E-003,
+ 2.863437556468661E-003,
+ -7.366775108654962E-003,
+ 2.058355474821512E-002,
+ -6.735323259371034E-002,
+ 3.224669577325661E-001,
+ 4.227843421859038E-001
};
/* log gamma(x+1), -.25 < x < .25 */
static const float C[] = {
- 1.369488127325832E-001,
--1.590086327657347E-001,
- 1.692415923504637E-001,
--2.067882815621965E-001,
- 2.705806208275915E-001,
--4.006931650563372E-001,
- 8.224670749082976E-001,
--5.772156501719101E-001
+ 1.369488127325832E-001,
+ -1.590086327657347E-001,
+ 1.692415923504637E-001,
+ -2.067882815621965E-001,
+ 2.705806208275915E-001,
+ -4.006931650563372E-001,
+ 8.224670749082976E-001,
+ -5.772156501719101E-001
};
/* log( sqrt( 2*pi ) ) */
@@ -36,154 +36,153 @@
/* Reentrant version */
/* Logarithm of gamma function */
-float __lgammaf_r( float x, int* sgngamf );
+float __lgammaf_r(float x, int* sgngamf);
-float __lgammaf_r( float x, int* sgngamf )
+float __lgammaf_r(float x, int* sgngamf)
{
-float p, q, w, z;
-float nx, tx;
-int i, direction;
+ float p, q, w, z;
+ float nx, tx;
+ int i, direction;
-*sgngamf = 1;
+ *sgngamf = 1;
#ifdef NANS
-if( isnan(x) )
- return(x);
+ if (isnan(x))
+ return (x);
#endif
#ifdef INFINITIES
-if( !isfinite(x) )
- return(x);
+ if (!isfinite(x))
+ return (x);
#endif
-
-if( x < 0.0 )
+ if (x < 0.0)
{
- q = -x;
- w = __lgammaf_r(q, sgngamf); /* note this modifies sgngam! */
- p = floorf(q);
- if( p == q )
+ q = -x;
+ w = __lgammaf_r(q, sgngamf); /* note this modifies sgngam! */
+ p = floorf(q);
+ if (p == q)
{
lgsing:
- _SET_ERRNO(EDOM);
- mtherr( "lgamf", SING );
+ _SET_ERRNO(EDOM);
+ mtherr("lgamf", SING);
#ifdef INFINITIES
- return (INFINITYF);
+ return (INFINITYF);
#else
- return( *sgngamf * MAXNUMF );
+ return( *sgngamf * MAXNUMF );
#endif
}
- i = p;
- if( (i & 1) == 0 )
- *sgngamf = -1;
- else
- *sgngamf = 1;
- z = q - p;
- if( z > 0.5 )
+ i = p;
+ if ((i & 1) == 0)
+ *sgngamf = -1;
+ else
+ *sgngamf = 1;
+ z = q - p;
+ if (z > 0.5)
{
- p += 1.0;
- z = p - q;
+ p += 1.0;
+ z = p - q;
}
- z = q * sinf( PIF * z );
- if( z == 0.0 )
- goto lgsing;
- z = -logf( PIINV*z ) - w;
- return( z );
+ z = q * sinf(PIF * z);
+ if (z == 0.0)
+ goto lgsing;
+ z = -logf(PIINV * z) - w;
+ return (z);
}
-if( x < 6.5 )
+ if (x < 6.5)
{
- direction = 0;
- z = 1.0;
- tx = x;
- nx = 0.0;
- if( x >= 1.5 )
+ direction = 0;
+ z = 1.0;
+ tx = x;
+ nx = 0.0;
+ if (x >= 1.5)
{
- while( tx > 2.5 )
+ while (tx > 2.5)
{
- nx -= 1.0;
- tx = x + nx;
- z *=tx;
+ nx -= 1.0;
+ tx = x + nx;
+ z *=tx;
}
- x += nx - 2.0;
+ x += nx - 2.0;
iv1r5:
- p = x * polevlf( x, B, 7 );
- goto cont;
+ p = x * polevlf( x, B, 7 );
+ goto cont;
}
- if( x >= 1.25 )
+ if (x >= 1.25)
{
- z *= x;
- x -= 1.0; /* x + 1 - 2 */
+ z *= x;
+ x -= 1.0; /* x + 1 - 2 */
+ direction = 1;
+ goto iv1r5;
+ }
+ if (x >= 0.75)
+ {
+ x -= 1.0;
+ p = x * polevlf( x, C, 7 );
+ q = 0.0;
+ goto contz;
+ }
+ while (tx < 1.5)
+ {
+ if (tx == 0.0)
+ goto lgsing;
+ z *=tx;
+ nx += 1.0;
+ tx = x + nx;
+ }
direction = 1;
- goto iv1r5;
- }
- if( x >= 0.75 )
- {
- x -= 1.0;
- p = x * polevlf( x, C, 7 );
- q = 0.0;
- goto contz;
- }
- while( tx < 1.5 )
- {
- if( tx == 0.0 )
- goto lgsing;
- z *=tx;
- nx += 1.0;
- tx = x + nx;
- }
- direction = 1;
- x += nx - 2.0;
- p = x * polevlf( x, B, 7 );
+ x += nx - 2.0;
+ p = x * polevlf(x, B, 7);
cont:
- if( z < 0.0 )
+ if (z < 0.0)
{
- *sgngamf = -1;
- z = -z;
+ *sgngamf = -1;
+ z = -z;
}
- else
+ else
{
- *sgngamf = 1;
+ *sgngamf = 1;
}
- q = logf(z);
- if( direction )
- q = -q;
+ q = logf(z);
+ if (direction)
+ q = -q;
contz:
- return( p + q );
+ return( p + q );
}
-if( x > MAXLGM )
+ if (x > MAXLGM)
{
- _SET_ERRNO(ERANGE);
- mtherr( "lgamf", OVERFLOW );
+ _SET_ERRNO(ERANGE);
+ mtherr("lgamf", OVERFLOW);
#ifdef INFINITIES
- return( *sgngamf * INFINITYF );
+ return (*sgngamf * INFINITYF);
#else
- return( *sgngamf * MAXNUMF );
+ return (*sgngamf * MAXNUMF);
#endif
}
/* Note, though an asymptotic formula could be used for x >= 3,
* there is cancellation error in the following if x < 6.5. */
-q = LS2PI - x;
-q += ( x - 0.5 ) * logf(x);
+ q = LS2PI - x;
+ q += (x - 0.5) * logf(x);
-if( x <= 1.0e4 )
+ if (x <= 1.0e4)
{
- z = 1.0/x;
- p = z * z;
- q += (( 6.789774945028216E-004 * p
- - 2.769887652139868E-003 ) * p
- + 8.333316229807355E-002 ) * z;
+ z = 1.0/x;
+ p = z * z;
+ q += (( 6.789774945028216E-004 * p
+ - 2.769887652139868E-003 ) * p
+ + 8.333316229807355E-002 ) * z;
}
-return( q );
+ return (q);
}
/* This is the C99 version */
-
float lgammaf(float x)
{
- int local_sgngamf=0;
- return (__lgammaf_r(x, &local_sgngamf));
+ int local_sgngamf = 0;
+ return (__lgammaf_r(x, &local_sgngamf));
}
+
diff --git a/mingw-w64-crt/math/lgammal.c b/mingw-w64-crt/math/lgammal.c
index da2affe..7fdf6ee 100644
--- a/mingw-w64-crt/math/lgammal.c
+++ b/mingw-w64-crt/math/lgammal.c
@@ -8,14 +8,14 @@
#if UNK
static uLD S[9] = {
{ { -1.193945051381510095614E-3L } },
- { { 7.220599478036909672331E-3L } },
+ { { 7.220599478036909672331E-3L } },
{ { -9.622023360406271645744E-3L } },
{ { -4.219773360705915470089E-2L } },
- { { 1.665386113720805206758E-1L } },
+ { { 1.665386113720805206758E-1L } },
{ { -4.200263503403344054473E-2L } },
{ { -6.558780715202540684668E-1L } },
- { { 5.772156649015328608253E-1L } },
- { { 1.000000000000000000000E0L } }
+ { { 5.772156649015328608253E-1L } },
+ { { 1.000000000000000000000E0L } }
};
#endif
#if IBMPC
@@ -47,14 +47,14 @@
#if UNK
static uLD SN[9] = {
- { { 1.133374167243894382010E-3L } },
- { { 7.220837261893170325704E-3L } },
- { { 9.621911155035976733706E-3L } },
+ { { 1.133374167243894382010E-3L } },
+ { { 7.220837261893170325704E-3L } },
+ { { 9.621911155035976733706E-3L } },
{ { -4.219773343731191721664E-2L } },
{ { -1.665386113944413519335E-1L } },
{ { -4.200263503402112910504E-2L } },
- { { 6.558780715202536547116E-1L } },
- { { 5.772156649015328608727E-1L } },
+ { { 6.558780715202536547116E-1L } },
+ { { 5.772156649015328608727E-1L } },
{ { -1.000000000000000000000E0L } }
};
#endif
@@ -98,13 +98,13 @@
*/
#if UNK
static uLD A[7] = {
- { { 4.885026142432270781165E-3L } },
+ { { 4.885026142432270781165E-3L } },
{ { -1.880801938119376907179E-3L } },
- { { 8.412723297322498080632E-4L } },
+ { { 8.412723297322498080632E-4L } },
{ { -5.952345851765688514613E-4L } },
- { { 7.936507795855070755671E-4L } },
+ { { 7.936507795855070755671E-4L } },
{ { -2.777777777750349603440E-3L } },
- { { 8.333333333333331447505E-2L } }
+ { { 8.333333333333331447505E-2L } }
};
#endif
#if IBMPC
@@ -207,130 +207,129 @@
long double __lgammal_r(long double x, int* sgngaml)
{
- long double p, q, w, z, f, nx;
- int i;
+ long double p, q, w, z, f, nx;
+ int i;
- *sgngaml = 1;
+ *sgngaml = 1;
#ifdef NANS
- if( isnanl(x) )
- return(NANL);
+ if (isnanl(x))
+ return(NANL);
#endif
#ifdef INFINITIES
- if( !isfinitel(x) )
- return(INFINITYL);
+ if (!isfinitel(x))
+ return (INFINITYL);
#endif
- if( x < -34.0L )
- {
- q = -x;
- w = __lgammal_r(q, sgngaml); /* note this modifies sgngam! */
- p = floorl(q);
- if( p == q )
- {
+ if (x < -34.0L)
+ {
+ q = -x;
+ w = __lgammal_r(q, sgngaml); /* note this modifies sgngam! */
+ p = floorl(q);
+ if (p == q)
+ {
lgsing:
- _SET_ERRNO(EDOM);
- mtherr( "lgammal", SING );
+ _SET_ERRNO(EDOM);
+ mtherr( "lgammal", SING );
#ifdef INFINITIES
- return (INFINITYL);
+ return (INFINITYL);
#else
- return (MAXNUML);
+ return (MAXNUML);
#endif
- }
- i = p;
- if( (i & 1) == 0 )
- *sgngaml = -1;
- else
- *sgngaml = 1;
- z = q - p;
- if( z > 0.5L )
- {
- p += 1.0L;
- z = p - q;
- }
- z = q * sinl( PIL * z );
- if( z == 0.0L )
- goto lgsing;
- /* z = LOGPI - logl( z ) - w; */
- z = logl( PIL/z ) - w;
- return( z );
- }
+ }
+ i = p;
+ if ((i & 1) == 0)
+ *sgngaml = -1;
+ else
+ *sgngaml = 1;
+ z = q - p;
+ if (z > 0.5L)
+ {
+ p += 1.0L;
+ z = p - q;
+ }
+ z = q * sinl(PIL * z);
+ if (z == 0.0L)
+ goto lgsing;
+ /* z = LOGPI - logl( z ) - w; */
+ z = logl(PIL/z) - w;
+ return (z);
+ }
- if( x < 13.0L )
- {
- z = 1.0L;
- nx = floorl( x + 0.5L );
- f = x - nx;
- while( x >= 3.0L )
- {
- nx -= 1.0L;
- x = nx + f;
- z *= x;
- }
- while( x < 2.0L )
- {
- if( fabsl(x) <= 0.03125 )
- goto lsmall;
- z /= nx + f;
- nx += 1.0L;
- x = nx + f;
- }
- if( z < 0.0L )
- {
- *sgngaml = -1;
- z = -z;
- }
- else
- *sgngaml = 1;
- if( x == 2.0L )
- return( logl(z) );
- x = (nx - 2.0L) + f;
- p = x * polevll( x, B, 6 ) / p1evll( x, C, 7);
- return( logl(z) + p );
- }
+ if (x < 13.0L)
+ {
+ z = 1.0L;
+ nx = floorl(x + 0.5L);
+ f = x - nx;
+ while (x >= 3.0L)
+ {
+ nx -= 1.0L;
+ x = nx + f;
+ z *= x;
+ }
+ while (x < 2.0L)
+ {
+ if (fabsl(x) <= 0.03125)
+ goto lsmall;
+ z /= nx + f;
+ nx += 1.0L;
+ x = nx + f;
+ }
+ if (z < 0.0L)
+ {
+ *sgngaml = -1;
+ z = -z;
+ }
+ else
+ *sgngaml = 1;
+ if (x == 2.0L)
+ return ( logl(z) );
+ x = (nx - 2.0L) + f;
+ p = x * polevll(x, B, 6) / p1evll(x, C, 7);
+ return ( logl(z) + p );
+ }
- if( x > MAXLGM )
- {
- _SET_ERRNO(ERANGE);
- mtherr( "lgammal", OVERFLOW );
+ if (x > MAXLGM)
+ {
+ _SET_ERRNO(ERANGE);
+ mtherr("lgammal", OVERFLOW);
#ifdef INFINITIES
- return( *sgngaml * INFINITYL );
+ return (*sgngaml * INFINITYL);
#else
- return( *sgngaml * MAXNUML );
+ return (*sgngaml * MAXNUML);
#endif
- }
+ }
- q = ( x - 0.5L ) * logl(x) - x + LS2PI;
- if( x > 1.0e10L )
- return(q);
- p = 1.0L/(x*x);
- q += polevll( p, A, 6 ) / x;
- return( q );
-
+ q = (x - 0.5L) * logl(x) - x + LS2PI;
+ if (x > 1.0e10L)
+ return(q);
+ p = 1.0L/(x*x);
+ q += polevll(p, A, 6) / x;
+ return (q);
lsmall:
- if( x == 0.0L )
- goto lgsing;
- if( x < 0.0L )
- {
- x = -x;
- q = z / (x * polevll( x, SN, 8 ));
- }
- else
- q = z / (x * polevll( x, S, 8 ));
- if( q < 0.0L )
- {
- *sgngaml = -1;
- q = -q;
- }
- else
- *sgngaml = 1;
- q = logl( q );
- return(q);
+ if (x == 0.0L)
+ goto lgsing;
+ if (x < 0.0L)
+ {
+ x = -x;
+ q = z / (x * polevll(x, SN, 8));
+ }
+ else
+ q = z / (x * polevll(x, S, 8));
+ if (q < 0.0L)
+ {
+ *sgngaml = -1;
+ q = -q;
+ }
+ else
+ *sgngaml = 1;
+ q = logl(q);
+ return (q);
}
/* This is the C99 version */
-
long double lgammal(long double x)
{
- int local_sgngaml=0;
- return (__lgammal_r(x, &local_sgngaml));
+ int local_sgngaml = 0;
+ return (__lgammal_r(x, &local_sgngaml));
}
+
diff --git a/mingw-w64-crt/math/llround.c b/mingw-w64-crt/math/llround.c
index 769cbb9..99efaa2 100644
--- a/mingw-w64-crt/math/llround.c
+++ b/mingw-w64-crt/math/llround.c
@@ -11,27 +11,28 @@
llround (double x)
{
double res;
-
+
if (x >= 0.0)
{
res = ceil (x);
if (res - x > 0.5)
- res -= 1.0;
+ res -= 1.0;
}
else
{
res = ceil (-x);
if (res + x > 0.5)
- res -= 1.0;
- res = -res;;
+ res -= 1.0;
+ res = -res;
}
if (!isfinite (res)
|| res > (double) LONG_LONG_MAX
|| res < (double) LONG_LONG_MIN)
- {
+ {
errno = ERANGE;
- /* Undefined behaviour, so we could return anything. */
+ /* Undefined behaviour, so we could return anything. */
/* return res > 0.0 ? LONG_LONG_MAX : LONG_LONG_MIN; */
}
return (long long) res;
}
+
diff --git a/mingw-w64-crt/math/llroundf.c b/mingw-w64-crt/math/llroundf.c
index 8277c18..3c449e6 100644
--- a/mingw-w64-crt/math/llroundf.c
+++ b/mingw-w64-crt/math/llroundf.c
@@ -11,7 +11,7 @@
llroundf (float x)
{
float res;
-
+
if (x >= 0.0F)
{
res = ceilf (x);
@@ -34,4 +34,5 @@
/* return res > 0.0F ? LONG_LONG_MAX : LONG_LONG_MIN; */
}
return (long long) res;
-}
+}
+
diff --git a/mingw-w64-crt/math/llroundl.c b/mingw-w64-crt/math/llroundl.c
index b39de26..b8aad25 100644
--- a/mingw-w64-crt/math/llroundl.c
+++ b/mingw-w64-crt/math/llroundl.c
@@ -11,7 +11,7 @@
llroundl (long double x)
{
long double res;
-
+
if (x >= 0.0L)
{
res = ceill (x);
@@ -35,3 +35,4 @@
}
return (long long) res;
}
+
diff --git a/mingw-w64-crt/math/lroundf.c b/mingw-w64-crt/math/lroundf.c
index d6475be..d03b352 100644
--- a/mingw-w64-crt/math/lroundf.c
+++ b/mingw-w64-crt/math/lroundf.c
@@ -34,4 +34,4 @@
/* return res > 0.0F ? LONG_MAX : LONG_MIN; */
}
return (long) res;
-}
+}
diff --git a/mingw-w64-crt/math/lroundl.c b/mingw-w64-crt/math/lroundl.c
index 64303b1..633d382 100644
--- a/mingw-w64-crt/math/lroundl.c
+++ b/mingw-w64-crt/math/lroundl.c
@@ -23,7 +23,7 @@
res = ceill (-x);
if (res + x > 0.5L)
res -= 1.0L;
- res = -res;;
+ res = -res;
}
if (!isfinite (res)
|| res > (long double)LONG_MAX
diff --git a/mingw-w64-crt/math/nextafterl.c b/mingw-w64-crt/math/nextafterl.c
index 1a4bbf7..3d320c9 100644
--- a/mingw-w64-crt/math/nextafterl.c
+++ b/mingw-w64-crt/math/nextafterl.c
@@ -59,7 +59,6 @@
/* If we have updated the expn of a normal number,
or moved from denormal to normal, [re]set the normal bit. */
-
if (u.parts.expn & 0x7fff)
u.parts.mantissa |= normal_bit;
@@ -69,3 +68,4 @@
/* nexttowardl is the same function with a different name. */
long double
nexttowardl (long double, long double) __attribute__ ((alias("nextafterl")));
+
diff --git a/mingw-w64-crt/math/pow.c b/mingw-w64-crt/math/pow.c
index 4a60550..92874e1 100644
--- a/mingw-w64-crt/math/pow.c
+++ b/mingw-w64-crt/math/pow.c
@@ -276,7 +276,6 @@
#else /* __MINGW32__ */
-#ifdef ANSIPROT
extern double floor ( double );
extern double fabs ( double );
extern double frexp ( double, int * );
@@ -288,12 +287,6 @@
extern int isnan ( double );
extern int isfinite ( double );
static double reduc ( double );
-#else
-double floor(), fabs(), frexp(), ldexp();
-double polevl(), p1evl(), __powi();
-int signbit(), isnan(), isfinite();
-static double reduc();
-#endif
extern double MAXNUM;
#ifdef INFINITIES
extern double INFINITY;
@@ -307,397 +300,390 @@
#endif /* __MINGW32__ */
-double pow( x, y )
-double x, y;
+double pow(double x, double y)
{
-double w, z, W, Wa, Wb, ya, yb, u;
-/* double F, Fa, Fb, G, Ga, Gb, H, Ha, Hb */
-double aw, ay, wy;
-int e, i, nflg, iyflg, yoddint;
+ double w, z, W, Wa, Wb, ya, yb, u;
+/* double F, Fa, Fb, G, Ga, Gb, H, Ha, Hb */
+ double aw, ay, wy;
+ int e, i, nflg, iyflg, yoddint;
-if( y == 0.0 )
- return( 1.0 );
+ if (y == 0.0)
+ return (1.0);
#ifdef NANS
-if( isnan(x) || isnan(y) )
+ if (isnan(x) || isnan(y))
{
- _SET_ERRNO (EDOM);
- return( x + y );
+ _SET_ERRNO (EDOM);
+ return (x + y);
}
#endif
-if( y == 1.0 )
- return( x );
-
+ if (y == 1.0)
+ return (x);
#ifdef INFINITIES
-if( !isfinite(y) && (x == 1.0 || x == -1.0) )
+ if (!isfinite(y) && (x == 1.0 || x == -1.0))
{
- mtherr( "pow", DOMAIN );
+ mtherr( "pow", DOMAIN );
#ifdef NANS
- return( NAN );
+ return( NAN );
#else
- return( INFINITY );
-#endif
- }
-#endif
-
-if( x == 1.0 )
- return( 1.0 );
-
-if( y >= MAXNUM )
- {
- _SET_ERRNO (ERANGE);
-#ifdef INFINITIES
- if( x > 1.0 )
return( INFINITY );
-#else
- if( x > 1.0 )
- return( MAXNUM );
#endif
- if( x > 0.0 && x < 1.0 )
- return( 0.0);
- if( x < -1.0 )
+ }
+#endif
+
+ if (x == 1.0)
+ return (1.0);
+
+ if (y >= MAXNUM)
+ {
+ _SET_ERRNO (ERANGE);
+#ifdef INFINITIES
+ if (x > 1.0)
+ return (INFINITY);
+#else
+ if (x > 1.0)
+ return (MAXNUM);
+#endif
+ if (x > 0.0 && x < 1.0)
+ return (0.0);
+ if (x < -1.0)
{
#ifdef INFINITIES
- return( INFINITY );
+ return (INFINITY);
#else
- return( MAXNUM );
+ return (MAXNUM);
#endif
}
- if( x > -1.0 && x < 0.0 )
- return( 0.0 );
+ if (x > -1.0 && x < 0.0)
+ return (0.0);
}
-if( y <= -MAXNUM )
+ if (y <= -MAXNUM)
{
- _SET_ERRNO (ERANGE);
- if( x > 1.0 )
- return( 0.0 );
+ _SET_ERRNO (ERANGE);
+ if (x > 1.0)
+ return (0.0);
#ifdef INFINITIES
- if( x > 0.0 && x < 1.0 )
- return( INFINITY );
+ if (x > 0.0 && x < 1.0)
+ return (INFINITY);
#else
- if( x > 0.0 && x < 1.0 )
- return( MAXNUM );
+ if (x > 0.0 && x < 1.0)
+ return (MAXNUM);
#endif
- if( x < -1.0 )
- return( 0.0 );
+ if (x < -1.0)
+ return (0.0);
#ifdef INFINITIES
- if( x > -1.0 && x < 0.0 )
- return( INFINITY );
+ if (x > -1.0 && x < 0.0)
+ return (INFINITY);
#else
- if( x > -1.0 && x < 0.0 )
- return( MAXNUM );
+ if (x > -1.0 && x < 0.0)
+ return (MAXNUM);
#endif
}
-if( x >= MAXNUM )
+ if (x >= MAXNUM)
{
#if INFINITIES
- if( y > 0.0 )
- return( INFINITY );
+ if (y > 0.0)
+ return (INFINITY);
#else
- if( y > 0.0 )
- return( MAXNUM );
+ if (y > 0.0)
+ return (MAXNUM);
#endif
- return(0.0);
+ return (0.0);
}
-/* Set iyflg to 1 if y is an integer. */
-iyflg = 0;
-w = floor(y);
-if( w == y )
- iyflg = 1;
+ /* Set iyflg to 1 if y is an integer. */
+ iyflg = 0;
+ w = floor(y);
+ if (w == y)
+ iyflg = 1;
/* Test for odd integer y. */
-yoddint = 0;
-if( iyflg )
+ yoddint = 0;
+ if (iyflg)
{
- ya = fabs(y);
- ya = floor(0.5 * ya);
- yb = 0.5 * fabs(w);
- if( ya != yb )
- yoddint = 1;
+ ya = fabs(y);
+ ya = floor(0.5 * ya);
+ yb = 0.5 * fabs(w);
+ if (ya != yb)
+ yoddint = 1;
}
-if( x <= -MAXNUM )
+ if (x <= -MAXNUM)
{
- if( y > 0.0 )
+ if (y > 0.0)
{
#ifdef INFINITIES
- if( yoddint )
- return( -INFINITY );
- return( INFINITY );
+ if (yoddint)
+ return (-INFINITY);
+ return (INFINITY);
#else
- if( yoddint )
- return( -MAXNUM );
- return( MAXNUM );
+ if (yoddint)
+ return (-MAXNUM);
+ return (MAXNUM);
#endif
}
- if( y < 0.0 )
+ if (y < 0.0)
{
#ifdef MINUSZERO
- if( yoddint )
- return( NEGZERO );
+ if (yoddint)
+ return (NEGZERO);
#endif
- return( 0.0 );
+ return (0.0);
}
}
-nflg = 0; /* flag = 1 if x<0 raised to integer power */
-if( x <= 0.0 )
+ nflg = 0; /* flag = 1 if x<0 raised to integer power */
+ if (x <= 0.0)
{
- if( x == 0.0 )
+ if (x == 0.0)
{
- if( y < 0.0 )
+ if (y < 0.0)
{
#ifdef MINUSZERO
- if( signbit(x) && yoddint )
- return( -INFINITY );
+ if (signbit(x) && yoddint)
+ return (-INFINITY);
#endif
#ifdef INFINITIES
- return( INFINITY );
+ return (INFINITY);
#else
- return( MAXNUM );
+ return (MAXNUM);
#endif
}
- if( y > 0.0 )
+ if (y > 0.0)
{
#ifdef MINUSZERO
- if( signbit(x) && yoddint )
- return( NEGZERO );
+ if (signbit(x) && yoddint)
+ return (NEGZERO);
#endif
- return( 0.0 );
+ return (0.0);
}
- return( 1.0 );
+ return (1.0);
}
- else
+ else
{
- if( iyflg == 0 )
+ if (iyflg == 0)
{ /* noninteger power of negative number */
- mtherr( fname, DOMAIN );
- _SET_ERRNO (EDOM);
+ mtherr(fname, DOMAIN);
+ _SET_ERRNO (EDOM);
#ifdef NANS
- return(NAN);
+ return (NAN);
#else
- return(0.0L);
+ return (0.0L);
#endif
}
- nflg = 1;
+ nflg = 1;
}
}
-/* Integer power of an integer. */
+ /* Integer power of an integer. */
-if( iyflg )
+ if (iyflg)
{
- i = w;
- w = floor(x);
- if( (w == x) && (fabs(y) < 32768.0) )
+ i = w;
+ w = floor(x);
+ if( (w == x) && (fabs(y) < 32768.0) )
{
- w = __powi( x, (int) y );
- return( w );
+ w = __powi(x, (int) y);
+ return (w);
}
}
-if( nflg )
- x = fabs(x);
+ if (nflg)
+ x = fabs(x);
-/* For results close to 1, use a series expansion. */
-w = x - 1.0;
-aw = fabs(w);
-ay = fabs(y);
-wy = w * y;
-ya = fabs(wy);
-if((aw <= 1.0e-3 && ay <= 1.0)
- || (ya <= 1.0e-3 && ay >= 1.0))
+ /* For results close to 1, use a series expansion. */
+ w = x - 1.0;
+ aw = fabs(w);
+ ay = fabs(y);
+ wy = w * y;
+ ya = fabs(wy);
+ if ((aw <= 1.0e-3 && ay <= 1.0)
+ || (ya <= 1.0e-3 && ay >= 1.0))
{
- z = (((((w*(y-5.)/720. + 1./120.)*w*(y-4.) + 1./24.)*w*(y-3.)
- + 1./6.)*w*(y-2.) + 0.5)*w*(y-1.) )*wy + wy + 1.;
- goto done;
+ z = (((((w*(y-5.)/720. + 1./120.)*w*(y-4.) + 1./24.)*w*(y-3.)
+ + 1./6.)*w*(y-2.) + 0.5)*w*(y-1.) )*wy + wy + 1.;
+ goto done;
}
-/* These are probably too much trouble. */
+ /* These are probably too much trouble. */
#if 0
-w = y * log(x);
-if (aw > 1.0e-3 && fabs(w) < 1.0e-3)
- {
- z = ((((((
- w/7. + 1.)*w/6. + 1.)*w/5. + 1.)*w/4. + 1.)*w/3. + 1.)*w/2. + 1.)*w + 1.;
- goto done;
- }
+ w = y * log(x);
+ if (aw > 1.0e-3 && fabs(w) < 1.0e-3)
+ {
+ z = ((((((w/7. + 1.)*w/6. + 1.)*w/5. + 1.)*w/4. + 1.)*w/3. + 1.)*w/2. + 1.)*w + 1.;
+ goto done;
+ }
-if(ya <= 1.0e-3 && aw <= 1.0e-4)
- {
- z = (((((
- wy*1./720.
- + (-w*1./48. + 1./120.) )*wy
- + ((w*17./144. - 1./12.)*w + 1./24.) )*wy
- + (((-w*5./16. + 7./24.)*w - 1./4.)*w + 1./6.) )*wy
- + ((((w*137./360. - 5./12.)*w + 11./24.)*w - 1./2.)*w + 1./2.) )*wy
- + (((((-w*1./6. + 1./5.)*w - 1./4)*w + 1./3.)*w -1./2.)*w ) )*wy
- + wy + 1.0;
- goto done;
- }
+ if (ya <= 1.0e-3 && aw <= 1.0e-4)
+ {
+ z = (((((wy*1./720.
+ + (-w*1./48. + 1./120.) )*wy
+ + ((w*17./144. - 1./12.)*w + 1./24.) )*wy
+ + (((-w*5./16. + 7./24.)*w - 1./4.)*w + 1./6.) )*wy
+ + ((((w*137./360. - 5./12.)*w + 11./24.)*w - 1./2.)*w + 1./2.) )*wy
+ + (((((-w*1./6. + 1./5.)*w - 1./4)*w + 1./3.)*w -1./2.)*w ) )*wy
+ + wy + 1.0;
+ goto done;
+ }
#endif
-/* separate significand from exponent */
-x = frexp( x, &e );
+ /* separate significand from exponent */
+ x = frexp(x, &e);
#if 0
-/* For debugging, check for gross overflow. */
-if( (e * y) > (MEXP + 1024) )
- goto overflow;
+ /* For debugging, check for gross overflow. */
+ if ((e * y) > (MEXP + 1024))
+ goto overflow;
#endif
-/* Find significand of x in antilog table A[]. */
-i = 1;
-if( x <= douba(9) )
- i = 9;
-if( x <= douba(i+4) )
- i += 4;
-if( x <= douba(i+2) )
- i += 2;
-if( x >= douba(1) )
- i = -1;
-i += 1;
+ /* Find significand of x in antilog table A[]. */
+ i = 1;
+ if (x <= douba(9))
+ i = 9;
+ 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);
-/* 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 * polevl( x, P, 3 ) / p1evl( x, Q, 4 ) );
+ w = w - ldexp( z, -1 ); /* w - 0.5 * z */
+ /* Convert to base 2 logarithm:
+ * multiply by log2(e)
+ */
+ w = w + LOG2EA * w;
+ /* Note x was not yet added in
+ * to above rational approximation,
+ * so do it now, while multiplying
+ * by log2(e).
+ */
+ z = w + LOG2EA * x;
+ z = z + x;
-/* 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 * polevl( x, P, 3 ) / p1evl( x, Q, 4 ) );
-w = w - ldexp( z, -1 ); /* w - 0.5 * z */
+ /* Compute exponent term of the base 2 logarithm. */
+ w = -i;
+ w = ldexp( w, -4 ); /* divide by 16 */
+ w += e;
+ /* Now base 2 log of x is w + z. */
-/* Convert to base 2 logarithm:
- * multiply by log2(e)
- */
-w = w + LOG2EA * w;
-/* Note x was not yet added in
- * to above rational approximation,
- * so do it now, while multiplying
- * by log2(e).
- */
-z = w + LOG2EA * x;
-z = z + x;
+ /* Multiply base 2 log by y, in extended precision. */
-/* Compute exponent term of the base 2 logarithm. */
-w = -i;
-w = ldexp( w, -4 ); /* divide by 16 */
-w += e;
-/* Now base 2 log of x is w + z. */
+ /* separate y into large part ya
+ * and small part yb less than 1/16
+ */
+ ya = reduc(y);
+ yb = y - ya;
-/* Multiply base 2 log by y, in extended precision. */
+ F = z * y + w * yb;
+ Fa = reduc(F);
+ Fb = F - Fa;
-/* separate y into large part ya
- * and small part yb less than 1/16
- */
-ya = reduc(y);
-yb = y - ya;
+ G = Fa + w * ya;
+ Ga = reduc(G);
+ Gb = G - Ga;
+ H = Fb + Gb;
+ Ha = reduc(H);
+ w = ldexp(Ga + Ha, 4);
-F = z * y + w * yb;
-Fa = reduc(F);
-Fb = F - Fa;
-
-G = Fa + w * ya;
-Ga = reduc(G);
-Gb = G - Ga;
-
-H = Fb + Gb;
-Ha = reduc(H);
-w = ldexp( Ga+Ha, 4 );
-
-/* Test the power of 2 for overflow */
-if( w > MEXP )
+ /* Test the power of 2 for overflow */
+ if (w > MEXP)
{
#ifndef INFINITIES
- mtherr( fname, OVERFLOW );
+ mtherr(fname, OVERFLOW);
#endif
#ifdef INFINITIES
- if( nflg && yoddint )
- return( -INFINITY );
- return( INFINITY );
+ if (nflg && yoddint)
+ return (-INFINITY);
+ return (INFINITY);
#else
- if( nflg && yoddint )
- return( -MAXNUM );
- return( MAXNUM );
+ if (nflg && yoddint)
+ return (-MAXNUM);
+ return (MAXNUM);
#endif
}
-if( w < (MNEXP - 1) )
+ if (w < (MNEXP - 1))
{
#ifndef DENORMAL
- mtherr( fname, UNDERFLOW );
+ mtherr(fname, UNDERFLOW);
#endif
#ifdef MINUSZERO
- if( nflg && yoddint )
- return( NEGZERO );
+ if (nflg && yoddint)
+ return (NEGZERO);
#endif
- return( 0.0 );
+ return (0.0);
}
-e = w;
-Hb = H - Ha;
+ e = w;
+ Hb = H - Ha;
-if( Hb > 0.0 )
+ if (Hb > 0.0)
{
- e += 1;
- Hb -= 0.0625;
+ e += 1;
+ Hb -= 0.0625;
}
-/* Now the product y * log2(x) = Hb + e/16.0.
- *
- * Compute base 2 exponential of Hb,
- * where -0.0625 <= Hb <= 0.
- */
-z = Hb * polevl( Hb, R, 6 ); /* z = 2**Hb - 1 */
+ /* Now the product y * log2(x) = Hb + e/16.0.
+ *
+ * Compute base 2 exponential of Hb,
+ * where -0.0625 <= Hb <= 0.
+ */
+ z = Hb * polevl(Hb, R, 6); /* z = 2**Hb - 1 */
-/* Express e/16 as an integer plus a negative number of 16ths.
- * Find lookup table entry for the fractional power of 2.
- */
-if( e < 0 )
- i = 0;
-else
- i = 1;
-i = e/16 + i;
-e = 16*i - e;
-w = douba( e );
-z = w + w * z; /* 2**-e * ( 1 + (2**Hb-1) ) */
-z = ldexp( z, i ); /* multiply by integer power of 2 */
+ /* Express e/16 as an integer plus a negative number of 16ths.
+ * Find lookup table entry for the fractional power of 2.
+ */
+ if (e < 0)
+ i = 0;
+ else
+ i = 1;
+ i = e/16 + i;
+ e = 16*i - e;
+ w = douba(e);
+ z = w + w * z; /* 2**-e * ( 1 + (2**Hb-1) ) */
+ z = ldexp(z, i); /* multiply by integer power of 2 */
done:
-
-/* Negate if odd integer power of negative number */
-if( nflg && yoddint )
+ /* Negate if odd integer power of negative number */
+ if (nflg && yoddint)
{
#ifdef MINUSZERO
- if( z == 0.0 )
- z = NEGZERO;
- else
+ if (z == 0.0)
+ z = NEGZERO;
+ else
#endif
- z = -z;
+ z = -z;
}
-return( z );
+ return (z);
}
/* Find a multiple of 1/16 that is within 1/16 of x. */
static double reduc(double x)
{
-double t;
+ double t;
-t = ldexp( x, 4 );
-t = floor( t );
-t = ldexp( t, -4 );
-return(t);
+ t = ldexp(x, 4);
+ t = floor(t);
+ t = ldexp(t, -4);
+ return (t);
}
+
diff --git a/mingw-w64-crt/math/powi.c b/mingw-w64-crt/math/powi.c
index d8c521c..ee2dce8 100644
--- a/mingw-w64-crt/math/powi.c
+++ b/mingw-w64-crt/math/powi.c
@@ -12,130 +12,128 @@
double __powi (double x,int nn)
{
-int n, e, sign, asign, lx;
-double w, y, s;
+ int n, e, sign, asign, lx;
+ double w, y, s;
-/* See pow.c for these tests. */
-if( x == 0.0 )
+ /* See pow.c for these tests. */
+ if (x == 0.0)
{
- if( nn == 0 )
- return( 1.0 );
- else if( nn < 0 )
- return( INFINITY );
+ if (nn == 0)
+ return (1.0);
+ else if (nn < 0)
+ return (INFINITY);
+ else
+ {
+ if (nn & 1)
+ return (x);
+ else
+ return (0.0);
+ }
+ }
+
+ if (nn == 0)
+ return (1.0);
+
+ if (nn == -1)
+ return (1.0/x);
+
+ if (x < 0.0)
+ {
+ asign = -1;
+ x = -x;
+ }
else
- {
- if( nn & 1 )
- return( x );
- else
- return( 0.0 );
- }
+ asign = 0;
+
+ if (nn < 0)
+ {
+ sign = -1;
+ n = -nn;
+ }
+ else
+ {
+ sign = 1;
+ n = nn;
}
-if( nn == 0 )
- return( 1.0 );
+ /* Even power will be positive. */
+ if ((n & 1) == 0)
+ asign = 0;
-if( nn == -1 )
- return( 1.0/x );
+ /* Overflow detection */
-if( x < 0.0 )
+ /* Calculate approximate logarithm of answer */
+ s = frexp(x, &lx);
+ e = (lx - 1)*n;
+ if ((e == 0) || (e > 64) || (e < -64))
{
- asign = -1;
- x = -x;
+ s = (s - 7.0710678118654752e-1) / (s + 7.0710678118654752e-1);
+ s = (2.9142135623730950 * s - 0.5 + lx) * nn * LOGE2;
}
-else
- asign = 0;
-
-
-if( nn < 0 )
+ else
{
- sign = -1;
- n = -nn;
- }
-else
- {
- sign = 1;
- n = nn;
+ s = LOGE2 * e;
}
-/* Even power will be positive. */
-if( (n & 1) == 0 )
- asign = 0;
-
-/* Overflow detection */
-
-/* Calculate approximate logarithm of answer */
-s = frexp( x, &lx );
-e = (lx - 1)*n;
-if( (e == 0) || (e > 64) || (e < -64) )
+ if (s > MAXLOG)
{
- s = (s - 7.0710678118654752e-1) / (s + 7.0710678118654752e-1);
- s = (2.9142135623730950 * s - 0.5 + lx) * nn * LOGE2;
- }
-else
- {
- s = LOGE2 * e;
- }
-
-if( s > MAXLOG )
- {
- mtherr( "powi", OVERFLOW );
- _SET_ERRNO(ERANGE);
- y = INFINITY;
- goto done;
+ mtherr("powi", OVERFLOW);
+ _SET_ERRNO(ERANGE);
+ y = INFINITY;
+ goto done;
}
#if DENORMAL
-if( s < MINLOG )
+ if (s < MINLOG)
{
- y = 0.0;
- goto done;
+ y = 0.0;
+ goto done;
}
-/* 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 < (-MAXLOG+2.0)) && (sign < 0) )
+ /* Handle tiny denormal answer, but with less accuracy
+ * since roundoff error in 1.0/x will be amplified.
+ * The precise demarcation should be the gradual underflow threshold.
+ */
+ if ((s < (-MAXLOG+2.0)) && (sign < 0))
{
- x = 1.0/x;
- sign = -sign;
+ x = 1.0/x;
+ sign = -sign;
}
#else
-/* do not produce denormal answer */
-if( s < -MAXLOG )
- return(0.0);
+ /* do not produce denormal answer */
+ if (s < -MAXLOG)
+ return (0.0);
#endif
-
/* First bit of the power */
-if( n & 1 )
- y = x;
-
-else
- y = 1.0;
+ if (n & 1)
+ y = x;
+ else
+ y = 1.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;
+ 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;
}
-if( sign < 0 )
- y = 1.0/y;
+ if (sign < 0)
+ y = 1.0/y;
done:
-if( asign )
+ if (asign)
{
- /* odd power of negative number */
- if( y == 0.0 )
- y = NEGZERO;
- else
- y = -y;
+ /* odd power of negative number */
+ if (y == 0.0)
+ y = NEGZERO;
+ else
+ y = -y;
}
-return(y);
+ return (y);
}
+
diff --git a/mingw-w64-crt/math/powif.c b/mingw-w64-crt/math/powif.c
index a0fb1eb..4a72acb 100644
--- a/mingw-w64-crt/math/powif.c
+++ b/mingw-w64-crt/math/powif.c
@@ -12,130 +12,128 @@
float __powif (float x, int nn)
{
-int n, e, sign, asign, lx;
-float w, y, s;
+ int n, e, sign, asign, lx;
+ float w, y, s;
-/* See pow.c for these tests. */
-if( x == 0.0F )
+ /* See pow.c for these tests. */
+ if (x == 0.0F)
{
- if( nn == 0 )
- return( 1.0F );
- else if( nn < 0 )
- return( INFINITYF );
+ if (nn == 0)
+ return (1.0F );
+ else if (nn < 0)
+ return (INFINITYF);
+ else
+ {
+ if (nn & 1)
+ return (x);
+ else
+ return (0.0);
+ }
+ }
+
+ if (nn == 0)
+ return (1.0);
+
+ if (nn == -1)
+ return (1.0/x);
+
+ if (x < 0.0)
+ {
+ asign = -1;
+ x = -x;
+ }
else
- {
- if( nn & 1 )
- return( x );
- else
- return( 0.0 );
- }
+ asign = 0;
+
+ if (nn < 0)
+ {
+ sign = -1;
+ n = -nn;
+ }
+ else
+ {
+ sign = 1;
+ n = nn;
}
-if( nn == 0 )
- return( 1.0 );
+ /* Even power will be positive. */
+ if ((n & 1) == 0)
+ asign = 0;
-if( nn == -1 )
- return( 1.0/x );
+ /* Overflow detection */
-if( x < 0.0 )
+ /* Calculate approximate logarithm of answer */
+ s = frexpf(x, &lx);
+ e = (lx - 1)*n;
+ if ((e == 0) || (e > 64) || (e < -64))
{
- asign = -1;
- x = -x;
+ s = (s - 7.0710678118654752e-1) / (s + 7.0710678118654752e-1);
+ s = (2.9142135623730950 * s - 0.5 + lx) * nn * LOGE2F;
}
-else
- asign = 0;
-
-
-if( nn < 0 )
+ else
{
- sign = -1;
- n = -nn;
- }
-else
- {
- sign = 1;
- n = nn;
+ s = LOGE2F * e;
}
-/* Even power will be positive. */
-if( (n & 1) == 0 )
- asign = 0;
-
-/* Overflow detection */
-
-/* Calculate approximate logarithm of answer */
-s = frexpf( x, &lx );
-e = (lx - 1)*n;
-if( (e == 0) || (e > 64) || (e < -64) )
+ if (s > MAXLOGF)
{
- s = (s - 7.0710678118654752e-1) / (s + 7.0710678118654752e-1);
- s = (2.9142135623730950 * s - 0.5 + lx) * nn * LOGE2F;
- }
-else
- {
- s = LOGE2F * e;
- }
-
-if( s > MAXLOGF )
- {
- mtherr( "__powif", OVERFLOW );
- _SET_ERRNO(ERANGE);
- y = INFINITYF;
- goto done;
+ mtherr("__powif", OVERFLOW);
+ _SET_ERRNO(ERANGE);
+ y = INFINITYF;
+ goto done;
}
#if DENORMAL
-if( s < MINLOGF )
+ if (s < MINLOGF)
{
- y = 0.0;
- goto done;
+ y = 0.0;
+ goto done;
}
-/* Handle tiny denormal answer, but with less accuracy
- * since roundoff error in 1.0/x will be amplified.
- * The precise demarcation should be the gradual underflow threshold.
- */
-if( (s < (-MAXLOGF+2.0)) && (sign < 0) )
+ /* Handle tiny denormal answer, but with less accuracy
+ * since roundoff error in 1.0/x will be amplified.
+ * The precise demarcation should be the gradual underflow threshold.
+ */
+ if ((s < (-MAXLOGF+2.0)) && (sign < 0))
{
- x = 1.0/x;
- sign = -sign;
+ x = 1.0/x;
+ sign = -sign;
}
#else
-/* do not produce denormal answer */
-if( s < -MAXLOGF )
- return(0.0);
+ /* do not produce denormal answer */
+ if (s < -MAXLOGF)
+ return (0.0);
#endif
+ /* First bit of the power */
+ if (n & 1)
+ y = x;
+ else
+ y = 1.0;
-/* First bit of the power */
-if( n & 1 )
- y = x;
-
-else
- y = 1.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;
+ 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;
}
-if( sign < 0 )
- y = 1.0/y;
+ if (sign < 0)
+ y = 1.0/y;
done:
-if( asign )
+ if (asign)
{
- /* odd power of negative number */
- if( y == 0.0 )
- y = NEGZEROF;
- else
- y = -y;
+ /* odd power of negative number */
+ if (y == 0.0)
+ y = NEGZEROF;
+ else
+ y = -y;
}
-return(y);
+ return (y);
}
+
diff --git a/mingw-w64-crt/math/powil.c b/mingw-w64-crt/math/powil.c
index 790e629..1e1e510 100644
--- a/mingw-w64-crt/math/powil.c
+++ b/mingw-w64-crt/math/powil.c
@@ -12,110 +12,107 @@
long double __powil (long double x,int nn)
{
-long double w, y;
-long double s;
-int n, e, sign, asign, lx;
+ long double w, y;
+ long double s;
+ int n, e, sign, asign, lx;
-if( x == 0.0L )
+ if (x == 0.0L)
{
- if( nn == 0 )
- return( 1.0L );
- else if( nn < 0 )
- return( INFINITYL );
+ 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
- return( 0.0L );
+ asign = 0;
+
+ if (nn < 0)
+ {
+ sign = -1;
+ n = -nn;
+ }
+ else
+ {
+ sign = 1;
+ n = nn;
}
-if( nn == 0 )
- return( 1.0L );
+ /* Overflow detection */
-
-if( x < 0.0L )
+ /* Calculate approximate logarithm of answer */
+ s = x;
+ s = frexpl(s, &lx);
+ e = (lx - 1)*n;
+ if ((e == 0) || (e > 64) || (e < -64))
{
- asign = -1;
- x = -x;
+ s = (s - 7.0710678118654752e-1L) / (s + 7.0710678118654752e-1L);
+ s = (2.9142135623730950L * s - 0.5L + lx) * nn * LOGE2L;
}
-else
- asign = 0;
-
-
-if( nn < 0 )
+ else
{
- sign = -1;
- n = -nn;
- }
-else
- {
- sign = 1;
- n = nn;
+ s = LOGE2L * e;
}
-/* Overflow detection */
-
-/* Calculate approximate logarithm of answer */
-s = x;
-s = frexpl( s, &lx );
-e = (lx - 1)*n;
-if( (e == 0) || (e > 64) || (e < -64) )
+ if (s > MAXLOGL)
{
- s = (s - 7.0710678118654752e-1L) / (s + 7.0710678118654752e-1L);
- s = (2.9142135623730950L * s - 0.5L + lx) * nn * LOGE2L;
- }
-else
- {
- s = LOGE2L * e;
+ mtherr("__powil", OVERFLOW);
+ _SET_ERRNO(ERANGE);
+ y = INFINITYL;
+ goto done;
}
-if( s > MAXLOGL )
+ if (s < MINLOGL)
{
- mtherr( "__powil", OVERFLOW );
- _SET_ERRNO(ERANGE);
- y = INFINITYL;
- goto done;
+ 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;
}
-if( s < MINLOGL )
+ /* First bit of the power */
+ if (n & 1)
+ y = x;
+ else
{
- 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;
+ y = 1.0L;
+ asign = 0;
}
-/* 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;
+ 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);
+ if (asign)
+ y = -y; /* odd power of negative number */
+ if (sign < 0)
+ y = 1.0L/y;
+ return (y);
}
+
diff --git a/mingw-w64-crt/math/powl.c b/mingw-w64-crt/math/powl.c
index 6030465..09a5780 100644
--- a/mingw-w64-crt/math/powl.c
+++ b/mingw-w64-crt/math/powl.c
@@ -314,9 +314,9 @@
static VOLATILE long double z;
static long double w, W, Wa, Wb, ya, yb, u;
-static __inline__ long double reducl( long double );
-extern long double __powil ( long double, int );
-extern long double powl ( long double x, long double y);
+static __inline__ long double reducl(long double);
+extern long double __powil(long double, int);
+extern long double powl(long double, long double);
/* No error checking. We handle Infs and zeros ourselves. */
static __inline__ long double
@@ -331,358 +331,355 @@
#define ldexpl __fast_ldexpl
-long double powl( x, y )
-long double x, y;
+long double powl(long double x, long doubley)
{
-/* double F, Fa, Fb, G, Ga, Gb, H, Ha, Hb */
-int i, nflg, iyflg, yoddint;
-long e;
+/* double F, Fa, Fb, G, Ga, Gb, H, Ha, Hb */
+ int i, nflg, iyflg, yoddint;
+ long e;
-if( y == 0.0L )
- return( 1.0L );
+ if (y == 0.0L)
+ return (1.0L);
#ifdef NANS
-if( isnanl(x) )
+ if (isnanl(x))
{
- _SET_ERRNO (EDOM);
- return( x );
+ _SET_ERRNO (EDOM);
+ return (x);
}
-if( isnanl(y) )
+ if (isnanl(y))
{
- _SET_ERRNO (EDOM);
- return( y );
+ _SET_ERRNO (EDOM);
+ return (y);
}
#endif
-if( y == 1.0L )
- return( x );
+ if (y == 1.0L)
+ return (x);
-if( isinfl(y) && (x == -1.0L || x == 1.0L) )
- return( y );
+ if (isinfl(y) && (x == -1.0L || x == 1.0L))
+ return (y);
-if( x == 1.0L )
- return( 1.0L );
+ if (x == 1.0L)
+ return (1.0L);
-if( y >= MAXNUML )
+ if (y >= MAXNUML)
{
- _SET_ERRNO (ERANGE);
+ _SET_ERRNO (ERANGE);
#ifdef INFINITIES
- if( x > 1.0L )
- return( INFINITYL );
+ if (x > 1.0L)
+ return (INFINITYL);
#else
- if( x > 1.0L )
- return( MAXNUML );
+ if (x > 1.0L)
+ return (MAXNUML);
#endif
- if( x > 0.0L && x < 1.0L )
- return( 0.0L );
+ if (x > 0.0L && x < 1.0L)
+ return (0.0L);
#ifdef INFINITIES
- if( x < -1.0L )
- return( INFINITYL );
+ if (x < -1.0L)
+ return (INFINITYL);
#else
- if( x < -1.0L )
- return( MAXNUML );
+ if (x < -1.0L)
+ return (MAXNUML);
#endif
- if( x > -1.0L && x < 0.0L )
- return( 0.0L );
+ if (x > -1.0L && x < 0.0L)
+ return (0.0L);
}
-if( y <= -MAXNUML )
+ if (y <= -MAXNUML)
{
- _SET_ERRNO (ERANGE);
- if( x > 1.0L )
- return( 0.0L );
+ _SET_ERRNO (ERANGE);
+ if (x > 1.0L)
+ return (0.0L);
#ifdef INFINITIES
- if( x > 0.0L && x < 1.0L )
- return( INFINITYL );
+ if (x > 0.0L && x < 1.0L)
+ return (INFINITYL);
#else
- if( x > 0.0L && x < 1.0L )
- return( MAXNUML );
+ if (x > 0.0L && x < 1.0L)
+ return (MAXNUML);
#endif
- if( x < -1.0L )
- return( 0.0L );
+ if (x < -1.0L)
+ return (0.0L);
#ifdef INFINITIES
- if( x > -1.0L && x < 0.0L )
- return( INFINITYL );
+ if (x > -1.0L && x < 0.0L)
+ return (INFINITYL);
#else
- if( x > -1.0L && x < 0.0L )
- return( MAXNUML );
+ if (x > -1.0L && x < 0.0L)
+ return (MAXNUML);
#endif
}
-if( x >= MAXNUML )
+ if (x >= MAXNUML)
{
#if INFINITIES
- if( y > 0.0L )
- return( INFINITYL );
+ if (y > 0.0L)
+ return (INFINITYL);
#else
- if( y > 0.0L )
- return( MAXNUML );
+ if (y > 0.0L)
+ return (MAXNUML);
#endif
- return( 0.0L );
+ return (0.0L);
}
-w = floorl(y);
-/* Set iyflg to 1 if y is an integer. */
-iyflg = 0;
-if( w == y )
- iyflg = 1;
+ 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 )
+ /* 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;
+ ya = fabsl(y);
+ ya = floorl(0.5L * ya);
+ yb = 0.5L * fabsl(w);
+ if (ya != yb)
+ yoddint = 1;
}
-if( x <= -MAXNUML )
+ if (x <= -MAXNUML)
{
- if( y > 0.0L )
+ if (y > 0.0L)
{
#ifdef INFINITIES
- if( yoddint )
- return( -INFINITYL );
- return( INFINITYL );
+ if (yoddint)
+ return (-INFINITYL);
+ return (INFINITYL);
#else
- if( yoddint )
- return( -MAXNUML );
- return( MAXNUML );
+ if (yoddint)
+ return (-MAXNUML);
+ return (MAXNUML);
#endif
}
- if( y < 0.0L )
+ if (y < 0.0L)
{
#ifdef MINUSZERO
- if( yoddint )
- return( NEGZEROL );
+ if (yoddint)
+ return (NEGZEROL);
#endif
- return( 0.0 );
+ return (0.0);
}
- }
+ }
-nflg = 0; /* flag = 1 if x<0 raised to integer power */
-if( x <= 0.0L )
+ nflg = 0; /* flag = 1 if x<0 raised to integer power */
+ if (x <= 0.0L)
{
- if( x == 0.0L )
+ if (x == 0.0L)
{
- if( y < 0.0 )
+ if (y < 0.0)
{
#ifdef MINUSZERO
- if( signbitl(x) && yoddint )
- return( -INFINITYL );
+ if (signbitl(x) && yoddint)
+ return (-INFINITYL);
#endif
#ifdef INFINITIES
- return( INFINITYL );
+ return (INFINITYL);
#else
- return( MAXNUML );
+ return (MAXNUML);
#endif
}
- if( y > 0.0 )
+ if (y > 0.0)
{
#ifdef MINUSZERO
- if( signbitl(x) && yoddint )
- return( NEGZEROL );
+ if (signbitl(x) && yoddint)
+ return (NEGZEROL);
#endif
- return( 0.0 );
+ return (0.0);
}
- if( y == 0.0L )
- return( 1.0L ); /* 0**0 */
- else
- return( 0.0L ); /* 0**y */
+ if (y == 0.0L)
+ return (1.0L); /* 0**0 */
+ else
+ return (0.0L); /* 0**y */
}
- else
+ else
{
- if( iyflg == 0 )
+ if (iyflg == 0)
{ /* noninteger power of negative number */
- mtherr( fname, DOMAIN );
- _SET_ERRNO (EDOM);
+ mtherr(fname, DOMAIN);
+ _SET_ERRNO (EDOM);
#ifdef NANS
- return(NANL);
+ return (NANL);
#else
- return(0.0L);
+ return (0.0L);
#endif
}
- nflg = 1;
+ nflg = 1;
}
}
/* Integer power of an integer. */
-if( iyflg )
+ if (iyflg)
{
- i = w;
- w = floorl(x);
- if( (w == x) && (fabsl(y) < 32768.0) )
+ i = w;
+ w = floorl(x);
+ if ((w == x) && (fabsl(y) < 32768.0))
{
- w = __powil( x, (int) y );
- return( w );
+ w = __powil(x, (int) y);
+ return (w);
}
}
+ if (nflg)
+ x = fabsl(x);
-if( nflg )
- x = fabsl(x);
+ /* separate significand from exponent */
+ x = frexpl( x, &i );
+ e = i;
-/* 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 )
- {
- _SET_ERRNO (ERANGE);
- mtherr( fname, OVERFLOW );
- return( MAXNUML );
- }
-
-if( w < MNEXP )
- {
- _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
+ /* find significand in antilog table A[] */
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 (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;
-if( nflg )
+ /* 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)
{
-/* 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 */
+ _SET_ERRNO (ERANGE);
+ mtherr(fname, OVERFLOW);
+ return (MAXNUML);
}
-return( z );
+ if (w < MNEXP)
+ {
+ _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);
}
-static __inline__ long double
+static __inline__ long double
__convert_inf_to_maxnum(long double x)
{
- if (isinf(x))
- return (x > 0.0L ? MAXNUML : -MAXNUML);
- else
- return x;
+ if (isinf(x))
+ return (x > 0.0L ? MAXNUML : -MAXNUML);
+ else
+ return x;
}
-
/* Find a multiple of 1/NXT that is within 1/NXT of x. */
static long double reducl(long double x)
{
-long double t;
+ long double t;
-/* If the call to ldexpl overflows, set it to MAXNUML.
- This avoids Inf - Inf = Nan result when calculating the 'small'
- part of a reduction. Instead, the small part becomes Inf,
- causing under/overflow when adding it to the 'large' part.
- There must be a cleaner way of doing this. */
-t = __convert_inf_to_maxnum (ldexpl( x, LNXT ));
-t = floorl( t );
-t = ldexpl( t, -LNXT );
-return(t);
+ /* If the call to ldexpl overflows, set it to MAXNUML.
+ This avoids Inf - Inf = Nan result when calculating the 'small'
+ part of a reduction. Instead, the small part becomes Inf,
+ causing under/overflow when adding it to the 'large' part.
+ There must be a cleaner way of doing this. */
+ t = __convert_inf_to_maxnum (ldexpl( x, LNXT ));
+ t = floorl(t);
+ t = ldexpl(t, -LNXT);
+ return (t);
}
+
diff --git a/mingw-w64-crt/math/rint.c b/mingw-w64-crt/math/rint.c
index 6134e2c..c797506 100644
--- a/mingw-w64-crt/math/rint.c
+++ b/mingw-w64-crt/math/rint.c
@@ -4,7 +4,7 @@
* No warranty is given; refer to the file DISCLAIMER within this package.
*/
#include <math.h>
-double rint (double x){
+double rint (double x) {
double retval = 0.0;
__asm__ ("frndint;" : "=t" (retval) : "0" (x));
return retval;
diff --git a/mingw-w64-crt/math/rintf.c b/mingw-w64-crt/math/rintf.c
index 74da29f..9db9d33 100644
--- a/mingw-w64-crt/math/rintf.c
+++ b/mingw-w64-crt/math/rintf.c
@@ -5,7 +5,7 @@
*/
#include <math.h>
-float rintf (float x){
+float rintf (float x) {
float retval = 0.0F;
__asm__ ("frndint;": "=t" (retval) : "0" (x));
return retval;
diff --git a/mingw-w64-crt/math/rintl.c b/mingw-w64-crt/math/rintl.c
index 2a3aa52..05dd5cc 100644
--- a/mingw-w64-crt/math/rintl.c
+++ b/mingw-w64-crt/math/rintl.c
@@ -5,7 +5,7 @@
*/
#include <math.h>
-long double rintl (long double x){
+long double rintl (long double x) {
long double retval = 0.0L;
__asm__ ("frndint;": "=t" (retval) : "0" (x));
return retval;
diff --git a/mingw-w64-crt/math/roundl.c b/mingw-w64-crt/math/roundl.c
index aaf3273..df7b843 100644
--- a/mingw-w64-crt/math/roundl.c
+++ b/mingw-w64-crt/math/roundl.c
@@ -20,7 +20,7 @@
res = ceill (-x);
if (res + x > 0.5L)
res -= 1.0L;
- res = -res;;
+ res = -res;
}
return res;
}
diff --git a/mingw-w64-crt/math/s_erf.c b/mingw-w64-crt/math/s_erf.c
index e19da83..384f898 100644
--- a/mingw-w64-crt/math/s_erf.c
+++ b/mingw-w64-crt/math/s_erf.c
@@ -37,7 +37,6 @@
* Q is an odd poly of degree 10.
* -57.90
* | R - (erf(x)-x)/x | <= 2
- *
*
* Remark. The formula is derived by noting
* erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....)
@@ -46,14 +45,14 @@
* is close to one. The interval is chosen because the fix
* point of erf(x) is near 0.6174 (i.e., erf(x)=x when x is
* near 0.6174), and by some experiment, 0.84375 is chosen to
- * guarantee the error is less than one ulp for erf.
+ * guarantee the error is less than one ulp for erf.
*
- * 2. For |x| in [0.84375,1.25], let s = |x| - 1, and
+ * 2. For |x| in [0.84375,1.25], let s = |x| - 1, and
* c = 0.84506291151 rounded to single (24 bits)
- * erf(x) = sign(x) * (c + P1(s)/Q1(s))
- * erfc(x) = (1-c) - P1(s)/Q1(s) if x > 0
+ * erf(x) = sign(x) * (c + P1(s)/Q1(s))
+ * erfc(x) = (1-c) - P1(s)/Q1(s) if x > 0
* 1+(c+P1(s)/Q1(s)) if x < 0
- * |P1/Q1 - (erf(|x|)-c)| <= 2**-59.06
+ * |P1/Q1 - (erf(|x|)-c)| <= 2**-59.06
* Remark: here we use the taylor series expansion at x=1.
* erf(1+s) = erf(1) + s*Poly(s)
* = 0.845.. + P1(s)/Q1(s)
@@ -64,49 +63,49 @@
* P1(s) = degree 6 poly in s
* Q1(s) = degree 6 poly in s
*
- * 3. For x in [1.25,1/0.35(~2.857143)],
- * erfc(x) = (1/x)*exp(-x*x-0.5625+R1/S1)
- * erf(x) = 1 - erfc(x)
+ * 3. For x in [1.25,1/0.35(~2.857143)],
+ * erfc(x) = (1/x)*exp(-x*x-0.5625+R1/S1)
+ * erf(x) = 1 - erfc(x)
* where
* R1(z) = degree 7 poly in z, (z=1/x^2)
* S1(z) = degree 8 poly in z
*
- * 4. For x in [1/0.35,28]
- * erfc(x) = (1/x)*exp(-x*x-0.5625+R2/S2) if x > 0
+ * 4. For x in [1/0.35,28]
+ * erfc(x) = (1/x)*exp(-x*x-0.5625+R2/S2) if x > 0
* = 2.0 - (1/x)*exp(-x*x-0.5625+R2/S2) if -6<x<0
* = 2.0 - tiny (if x <= -6)
- * erf(x) = sign(x)*(1.0 - erfc(x)) if x < 6, else
- * erf(x) = sign(x)*(1.0 - tiny)
+ * erf(x) = sign(x)*(1.0 - erfc(x)) if x < 6, else
+ * erf(x) = sign(x)*(1.0 - tiny)
* where
* R2(z) = degree 6 poly in z, (z=1/x^2)
* S2(z) = degree 7 poly in z
*
- * Note1:
+ * Note1:
* To compute exp(-x*x-0.5625+R/S), let s be a single
* precision number and s := x; then
* -x*x = -s*s + (s-x)*(s+x)
- * exp(-x*x-0.5626+R/S) =
+ * exp(-x*x-0.5626+R/S) =
* exp(-s*s-0.5625)*exp((s-x)*(s+x)+R/S);
- * Note2:
+ * Note2:
* Here 4 and 5 make use of the asymptotic series
* exp(-x*x)
* erfc(x) ~ ---------- * ( 1 + Poly(1/x^2) )
* x*sqrt(pi)
* We use rational approximation to approximate
- * g(s)=f(1/x^2) = log(erfc(x)*x) - x*x + 0.5625
+ * g(s)=f(1/x^2) = log(erfc(x)*x) - x*x + 0.5625
* Here is the error bound for R1/S1 and R2/S2
- * |R1/S1 - f(x)| < 2**(-62.57)
- * |R2/S2 - f(x)| < 2**(-61.52)
+ * |R1/S1 - f(x)| < 2**(-62.57)
+ * |R2/S2 - f(x)| < 2**(-61.52)
*
- * 5. For inf > x >= 28
- * erf(x) = sign(x) *(1 - tiny) (raise inexact)
- * erfc(x) = tiny*tiny (raise underflow) if x > 0
+ * 5. For inf > x >= 28
+ * erf(x) = sign(x) *(1 - tiny) (raise inexact)
+ * erfc(x) = tiny*tiny (raise underflow) if x > 0
* = 2 - tiny if x<0
*
- * 7. Special case:
- * erf(0) = 0, erf(inf) = 1, erf(-inf) = -1,
- * erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2,
- * erfc/erf(NaN) is NaN
+ * 7. Special case:
+ * erf(0) = 0, erf(inf) = 1, erf(-inf) = -1,
+ * erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2,
+ * erfc/erf(NaN) is NaN
*/
@@ -118,10 +117,10 @@
#define __ieee754_exp exp
-typedef union
+typedef union
{
double value;
- struct
+ struct
{
uint32_t lsw;
uint32_t msw;
@@ -145,103 +144,95 @@
}
-#ifdef __STDC__
static const double
-#else
-static double
-#endif
-tiny = 1e-300,
-half= 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
-one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
-two = 2.00000000000000000000e+00, /* 0x40000000, 0x00000000 */
- /* c = (float)0.84506291151 */
-erx = 8.45062911510467529297e-01, /* 0x3FEB0AC1, 0x60000000 */
-/*
- * Coefficients for approximation to erf on [0,0.84375]
- */
-efx = 1.28379167095512586316e-01, /* 0x3FC06EBA, 0x8214DB69 */
-efx8= 1.02703333676410069053e+00, /* 0x3FF06EBA, 0x8214DB69 */
-pp0 = 1.28379167095512558561e-01, /* 0x3FC06EBA, 0x8214DB68 */
-pp1 = -3.25042107247001499370e-01, /* 0xBFD4CD7D, 0x691CB913 */
-pp2 = -2.84817495755985104766e-02, /* 0xBF9D2A51, 0xDBD7194F */
-pp3 = -5.77027029648944159157e-03, /* 0xBF77A291, 0x236668E4 */
-pp4 = -2.37630166566501626084e-05, /* 0xBEF8EAD6, 0x120016AC */
-qq1 = 3.97917223959155352819e-01, /* 0x3FD97779, 0xCDDADC09 */
-qq2 = 6.50222499887672944485e-02, /* 0x3FB0A54C, 0x5536CEBA */
-qq3 = 5.08130628187576562776e-03, /* 0x3F74D022, 0xC4D36B0F */
-qq4 = 1.32494738004321644526e-04, /* 0x3F215DC9, 0x221C1A10 */
-qq5 = -3.96022827877536812320e-06, /* 0xBED09C43, 0x42A26120 */
-/*
- * Coefficients for approximation to erf in [0.84375,1.25]
- */
-pa0 = -2.36211856075265944077e-03, /* 0xBF6359B8, 0xBEF77538 */
-pa1 = 4.14856118683748331666e-01, /* 0x3FDA8D00, 0xAD92B34D */
-pa2 = -3.72207876035701323847e-01, /* 0xBFD7D240, 0xFBB8C3F1 */
-pa3 = 3.18346619901161753674e-01, /* 0x3FD45FCA, 0x805120E4 */
-pa4 = -1.10894694282396677476e-01, /* 0xBFBC6398, 0x3D3E28EC */
-pa5 = 3.54783043256182359371e-02, /* 0x3FA22A36, 0x599795EB */
-pa6 = -2.16637559486879084300e-03, /* 0xBF61BF38, 0x0A96073F */
-qa1 = 1.06420880400844228286e-01, /* 0x3FBB3E66, 0x18EEE323 */
-qa2 = 5.40397917702171048937e-01, /* 0x3FE14AF0, 0x92EB6F33 */
-qa3 = 7.18286544141962662868e-02, /* 0x3FB2635C, 0xD99FE9A7 */
-qa4 = 1.26171219808761642112e-01, /* 0x3FC02660, 0xE763351F */
-qa5 = 1.36370839120290507362e-02, /* 0x3F8BEDC2, 0x6B51DD1C */
-qa6 = 1.19844998467991074170e-02, /* 0x3F888B54, 0x5735151D */
-/*
- * Coefficients for approximation to erfc in [1.25,1/0.35]
- */
-ra0 = -9.86494403484714822705e-03, /* 0xBF843412, 0x600D6435 */
-ra1 = -6.93858572707181764372e-01, /* 0xBFE63416, 0xE4BA7360 */
-ra2 = -1.05586262253232909814e+01, /* 0xC0251E04, 0x41B0E726 */
-ra3 = -6.23753324503260060396e+01, /* 0xC04F300A, 0xE4CBA38D */
-ra4 = -1.62396669462573470355e+02, /* 0xC0644CB1, 0x84282266 */
-ra5 = -1.84605092906711035994e+02, /* 0xC067135C, 0xEBCCABB2 */
-ra6 = -8.12874355063065934246e+01, /* 0xC0545265, 0x57E4D2F2 */
-ra7 = -9.81432934416914548592e+00, /* 0xC023A0EF, 0xC69AC25C */
-sa1 = 1.96512716674392571292e+01, /* 0x4033A6B9, 0xBD707687 */
-sa2 = 1.37657754143519042600e+02, /* 0x4061350C, 0x526AE721 */
-sa3 = 4.34565877475229228821e+02, /* 0x407B290D, 0xD58A1A71 */
-sa4 = 6.45387271733267880336e+02, /* 0x40842B19, 0x21EC2868 */
-sa5 = 4.29008140027567833386e+02, /* 0x407AD021, 0x57700314 */
-sa6 = 1.08635005541779435134e+02, /* 0x405B28A3, 0xEE48AE2C */
-sa7 = 6.57024977031928170135e+00, /* 0x401A47EF, 0x8E484A93 */
-sa8 = -6.04244152148580987438e-02, /* 0xBFAEEFF2, 0xEE749A62 */
-/*
- * Coefficients for approximation to erfc in [1/.35,28]
- */
-rb0 = -9.86494292470009928597e-03, /* 0xBF843412, 0x39E86F4A */
-rb1 = -7.99283237680523006574e-01, /* 0xBFE993BA, 0x70C285DE */
-rb2 = -1.77579549177547519889e+01, /* 0xC031C209, 0x555F995A */
-rb3 = -1.60636384855821916062e+02, /* 0xC064145D, 0x43C5ED98 */
-rb4 = -6.37566443368389627722e+02, /* 0xC083EC88, 0x1375F228 */
-rb5 = -1.02509513161107724954e+03, /* 0xC0900461, 0x6A2E5992 */
-rb6 = -4.83519191608651397019e+02, /* 0xC07E384E, 0x9BDC383F */
-sb1 = 3.03380607434824582924e+01, /* 0x403E568B, 0x261D5190 */
-sb2 = 3.25792512996573918826e+02, /* 0x40745CAE, 0x221B9F0A */
-sb3 = 1.53672958608443695994e+03, /* 0x409802EB, 0x189D5118 */
-sb4 = 3.19985821950859553908e+03, /* 0x40A8FFB7, 0x688C246A */
-sb5 = 2.55305040643316442583e+03, /* 0x40A3F219, 0xCEDF3BE6 */
-sb6 = 4.74528541206955367215e+02, /* 0x407DA874, 0xE79FE763 */
-sb7 = -2.24409524465858183362e+01; /* 0xC03670E2, 0x42712D62 */
+ tiny= 1e-300,
+ half= 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
+ one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
+ two = 2.00000000000000000000e+00, /* 0x40000000, 0x00000000 */
+ /* c = (float)0.84506291151 */
+ erx = 8.45062911510467529297e-01, /* 0x3FEB0AC1, 0x60000000 */
+ /*
+ * Coefficients for approximation to erf on [0,0.84375]
+ */
+ efx = 1.28379167095512586316e-01, /* 0x3FC06EBA, 0x8214DB69 */
+ efx8= 1.02703333676410069053e+00, /* 0x3FF06EBA, 0x8214DB69 */
+ pp0 = 1.28379167095512558561e-01, /* 0x3FC06EBA, 0x8214DB68 */
+ pp1 = -3.25042107247001499370e-01, /* 0xBFD4CD7D, 0x691CB913 */
+ pp2 = -2.84817495755985104766e-02, /* 0xBF9D2A51, 0xDBD7194F */
+ pp3 = -5.77027029648944159157e-03, /* 0xBF77A291, 0x236668E4 */
+ pp4 = -2.37630166566501626084e-05, /* 0xBEF8EAD6, 0x120016AC */
+ qq1 = 3.97917223959155352819e-01, /* 0x3FD97779, 0xCDDADC09 */
+ qq2 = 6.50222499887672944485e-02, /* 0x3FB0A54C, 0x5536CEBA */
+ qq3 = 5.08130628187576562776e-03, /* 0x3F74D022, 0xC4D36B0F */
+ qq4 = 1.32494738004321644526e-04, /* 0x3F215DC9, 0x221C1A10 */
+ qq5 = -3.96022827877536812320e-06, /* 0xBED09C43, 0x42A26120 */
+ /*
+ * Coefficients for approximation to erf in [0.84375,1.25]
+ */
+ pa0 = -2.36211856075265944077e-03, /* 0xBF6359B8, 0xBEF77538 */
+ pa1 = 4.14856118683748331666e-01, /* 0x3FDA8D00, 0xAD92B34D */
+ pa2 = -3.72207876035701323847e-01, /* 0xBFD7D240, 0xFBB8C3F1 */
+ pa3 = 3.18346619901161753674e-01, /* 0x3FD45FCA, 0x805120E4 */
+ pa4 = -1.10894694282396677476e-01, /* 0xBFBC6398, 0x3D3E28EC */
+ pa5 = 3.54783043256182359371e-02, /* 0x3FA22A36, 0x599795EB */
+ pa6 = -2.16637559486879084300e-03, /* 0xBF61BF38, 0x0A96073F */
+ qa1 = 1.06420880400844228286e-01, /* 0x3FBB3E66, 0x18EEE323 */
+ qa2 = 5.40397917702171048937e-01, /* 0x3FE14AF0, 0x92EB6F33 */
+ qa3 = 7.18286544141962662868e-02, /* 0x3FB2635C, 0xD99FE9A7 */
+ qa4 = 1.26171219808761642112e-01, /* 0x3FC02660, 0xE763351F */
+ qa5 = 1.36370839120290507362e-02, /* 0x3F8BEDC2, 0x6B51DD1C */
+ qa6 = 1.19844998467991074170e-02, /* 0x3F888B54, 0x5735151D */
+ /*
+ * Coefficients for approximation to erfc in [1.25,1/0.35]
+ */
+ ra0 = -9.86494403484714822705e-03, /* 0xBF843412, 0x600D6435 */
+ ra1 = -6.93858572707181764372e-01, /* 0xBFE63416, 0xE4BA7360 */
+ ra2 = -1.05586262253232909814e+01, /* 0xC0251E04, 0x41B0E726 */
+ ra3 = -6.23753324503260060396e+01, /* 0xC04F300A, 0xE4CBA38D */
+ ra4 = -1.62396669462573470355e+02, /* 0xC0644CB1, 0x84282266 */
+ ra5 = -1.84605092906711035994e+02, /* 0xC067135C, 0xEBCCABB2 */
+ ra6 = -8.12874355063065934246e+01, /* 0xC0545265, 0x57E4D2F2 */
+ ra7 = -9.81432934416914548592e+00, /* 0xC023A0EF, 0xC69AC25C */
+ sa1 = 1.96512716674392571292e+01, /* 0x4033A6B9, 0xBD707687 */
+ sa2 = 1.37657754143519042600e+02, /* 0x4061350C, 0x526AE721 */
+ sa3 = 4.34565877475229228821e+02, /* 0x407B290D, 0xD58A1A71 */
+ sa4 = 6.45387271733267880336e+02, /* 0x40842B19, 0x21EC2868 */
+ sa5 = 4.29008140027567833386e+02, /* 0x407AD021, 0x57700314 */
+ sa6 = 1.08635005541779435134e+02, /* 0x405B28A3, 0xEE48AE2C */
+ sa7 = 6.57024977031928170135e+00, /* 0x401A47EF, 0x8E484A93 */
+ sa8 = -6.04244152148580987438e-02, /* 0xBFAEEFF2, 0xEE749A62 */
+ /*
+ * Coefficients for approximation to erfc in [1/.35,28]
+ */
+ rb0 = -9.86494292470009928597e-03, /* 0xBF843412, 0x39E86F4A */
+ rb1 = -7.99283237680523006574e-01, /* 0xBFE993BA, 0x70C285DE */
+ rb2 = -1.77579549177547519889e+01, /* 0xC031C209, 0x555F995A */
+ rb3 = -1.60636384855821916062e+02, /* 0xC064145D, 0x43C5ED98 */
+ rb4 = -6.37566443368389627722e+02, /* 0xC083EC88, 0x1375F228 */
+ rb5 = -1.02509513161107724954e+03, /* 0xC0900461, 0x6A2E5992 */
+ rb6 = -4.83519191608651397019e+02, /* 0xC07E384E, 0x9BDC383F */
+ sb1 = 3.03380607434824582924e+01, /* 0x403E568B, 0x261D5190 */
+ sb2 = 3.25792512996573918826e+02, /* 0x40745CAE, 0x221B9F0A */
+ sb3 = 1.53672958608443695994e+03, /* 0x409802EB, 0x189D5118 */
+ sb4 = 3.19985821950859553908e+03, /* 0x40A8FFB7, 0x688C246A */
+ sb5 = 2.55305040643316442583e+03, /* 0x40A3F219, 0xCEDF3BE6 */
+ sb6 = 4.74528541206955367215e+02, /* 0x407DA874, 0xE79FE763 */
+ sb7 = -2.24409524465858183362e+01; /* 0xC03670E2, 0x42712D62 */
-#ifdef __STDC__
- double erf(double x)
-#else
- double erf(x)
- double x;
-#endif
+
+double erf(double x)
{
- int hx,ix,i;
- double R,S,P,Q,s,y,z,r;
+ int hx, ix, i;
+ double R, S, P, Q, s, y, z, r;
hx = __get_hi_word(x);
- ix = hx&0x7fffffff;
- if(ix>=0x7ff00000) { /* erf(nan)=nan */
+ ix = hx & 0x7fffffff;
+ if (ix >= 0x7ff00000) { /* erf(nan)=nan */
i = ((unsigned)hx>>31)<<1;
return (double)(1-i)+one/x; /* erf(+-inf)=+-1 */
}
- if(ix < 0x3feb0000) { /* |x|<0.84375 */
- if(ix < 0x3e300000) { /* |x|<2**-28 */
+ if (ix < 0x3feb0000) { /* |x|<0.84375 */
+ if (ix < 0x3e300000) { /* |x|<2**-28 */
if (ix < 0x00800000)
return 0.125*(8.0*x+efx8*x); /*avoid underflow */
return x + efx*x;
@@ -252,58 +243,62 @@
y = r/s;
return x + x*y;
}
- if(ix < 0x3ff40000) { /* 0.84375 <= |x| < 1.25 */
+ if (ix < 0x3ff40000) { /* 0.84375 <= |x| < 1.25 */
s = fabs(x)-one;
P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
- if(hx>=0) return erx + P/Q; else return -erx - P/Q;
+ if (hx >= 0)
+ return erx + P/Q;
+ else
+ return -erx - P/Q;
}
if (ix >= 0x40180000) { /* inf>|x|>=6 */
- if(hx>=0) return one-tiny; else return tiny-one;
+ if (hx >= 0)
+ return one-tiny;
+ else
+ return tiny-one;
}
x = fabs(x);
s = one/(x*x);
- if(ix< 0x4006DB6E) { /* |x| < 1/0.35 */
- R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
+ if (ix < 0x4006DB6E) { /* |x| < 1/0.35 */
+ R = ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
ra5+s*(ra6+s*ra7))))));
- S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
+ S = one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
sa5+s*(sa6+s*(sa7+s*sa8)))))));
} else { /* |x| >= 1/0.35 */
- R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
+ R = rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
rb5+s*rb6)))));
- S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
+ S = one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
sb5+s*(sb6+s*sb7))))));
}
- z = x;
+ z = x;
__trunc_lo_word(&z);
- r = __ieee754_exp(-z*z-0.5625)*__ieee754_exp((z-x)*(z+x)+R/S);
- if(hx>=0) return one-r/x; else return r/x-one;
+ r = __ieee754_exp(-z*z-0.5625)*__ieee754_exp((z-x)*(z+x)+R/S);
+ if (hx >= 0)
+ return one-r/x;
+ else
+ return r/x-one;
}
-#ifdef __STDC__
- double erfc(double x)
-#else
- double erfc(x)
- double x;
-#endif
+double erfc(double x)
{
int hx,ix;
double R,S,P,Q,s,y,z,r;
hx = __get_hi_word(x);
ix = hx&0x7fffffff;
- if(ix>=0x7ff00000) { /* erfc(nan)=nan */
+ if (ix >= 0x7ff00000) { /* erfc(nan)=nan */
/* erfc(+-inf)=0,2 */
return (double)(((unsigned)hx>>31)<<1)+one/x;
}
- if(ix < 0x3feb0000) { /* |x|<0.84375 */
- if(ix < 0x3c700000) /* |x|<2**-56 */
+ if (ix < 0x3feb0000) { /* |x|<0.84375 */
+ if (ix < 0x3c700000) /* |x|<2**-56 */
return one-x;
z = x*x;
r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
y = r/s;
- if(hx < 0x3fd00000) { /* x<1/4 */
+ if (hx < 0x3fd00000) { /* x<1/4 */
return one-(x+x*y);
} else {
r = x*y;
@@ -311,11 +306,11 @@
return half - r ;
}
}
- if(ix < 0x3ff40000) { /* 0.84375 <= |x| < 1.25 */
+ if (ix < 0x3ff40000) { /* 0.84375 <= |x| < 1.25 */
s = fabs(x)-one;
P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
- if(hx>=0) {
+ if (hx >= 0) {
z = one-erx; return z - P/Q;
} else {
z = erx+P/Q; return one+z;
@@ -324,26 +319,34 @@
if (ix < 0x403c0000) { /* |x|<28 */
x = fabs(x);
s = one/(x*x);
- if(ix< 0x4006DB6D) { /* |x| < 1/.35 ~ 2.857143*/
- R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
+ if (ix < 0x4006DB6D) { /* |x| < 1/.35 ~ 2.857143*/
+ R = ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
ra5+s*(ra6+s*ra7))))));
- S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
+ S = one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
sa5+s*(sa6+s*(sa7+s*sa8)))))));
} else { /* |x| >= 1/.35 ~ 2.857143 */
- if(hx<0&&ix>=0x40180000) return two-tiny;/* x < -6 */
- R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
+ if (hx < 0 && ix >= 0x40180000)
+ return two-tiny; /* x < -6 */
+ R = rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
rb5+s*rb6)))));
- S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
+ S = one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
sb5+s*(sb6+s*sb7))))));
}
- z = x;
+ z = x;
__trunc_lo_word(&z);
- r = __ieee754_exp(-z*z-0.5625)*
+ r = __ieee754_exp(-z*z-0.5625)*
__ieee754_exp((z-x)*(z+x)+R/S);
- if(hx>0) return r/x; else return two-r/x;
+ if (hx > 0)
+ return r/x;
+ else
+ return two-r/x;
} else {
/* set range error */
errno = ERANGE;
- if(hx>0) return tiny*tiny; else return two-tiny;
+ if (hx > 0)
+ return tiny*tiny;
+ else
+ return two-tiny;
}
}
+
diff --git a/mingw-w64-crt/math/sf_erf.c b/mingw-w64-crt/math/sf_erf.c
index 5b3ff40..f0e8f03 100644
--- a/mingw-w64-crt/math/sf_erf.c
+++ b/mingw-w64-crt/math/sf_erf.c
@@ -27,8 +27,6 @@
#define __ieee754_expf expf
-
-
typedef union
{
float value;
@@ -47,17 +45,17 @@
/* Set a float from a 32 bit int. */
-#define SET_FLOAT_WORD(d,i) \
-do { \
- ieee_float_shape_type sf_u; \
- sf_u.word = (i); \
- (d) = sf_u.value; \
+#define SET_FLOAT_WORD(d,i) \
+do { \
+ ieee_float_shape_type sf_u; \
+ sf_u.word = (i); \
+ (d) = sf_u.value; \
} while (0)
static inline void __trunc_float_word(float * x)
{
ieee_float_shape_type u;
- u.value = * x;
+ u.value = * x;
u.word &= 0xfffff000;
}
@@ -65,104 +63,95 @@
#define const
#endif
-#ifdef __STDC__
static const float
-#else
-static float
-#endif
-tiny = 1e-30,
-half= 5.0000000000e-01, /* 0x3F000000 */
-one = 1.0000000000e+00, /* 0x3F800000 */
-two = 2.0000000000e+00, /* 0x40000000 */
+ tiny= 1e-30,
+ half= 5.0000000000e-01, /* 0x3F000000 */
+ one = 1.0000000000e+00, /* 0x3F800000 */
+ two = 2.0000000000e+00, /* 0x40000000 */
/* c = (subfloat)0.84506291151 */
-erx = 8.4506291151e-01, /* 0x3f58560b */
-/*
- * Coefficients for approximation to erf on [0,0.84375]
- */
-efx = 1.2837916613e-01, /* 0x3e0375d4 */
-efx8= 1.0270333290e+00, /* 0x3f8375d4 */
-pp0 = 1.2837916613e-01, /* 0x3e0375d4 */
-pp1 = -3.2504209876e-01, /* 0xbea66beb */
-pp2 = -2.8481749818e-02, /* 0xbce9528f */
-pp3 = -5.7702702470e-03, /* 0xbbbd1489 */
-pp4 = -2.3763017452e-05, /* 0xb7c756b1 */
-qq1 = 3.9791721106e-01, /* 0x3ecbbbce */
-qq2 = 6.5022252500e-02, /* 0x3d852a63 */
-qq3 = 5.0813062117e-03, /* 0x3ba68116 */
-qq4 = 1.3249473704e-04, /* 0x390aee49 */
-qq5 = -3.9602282413e-06, /* 0xb684e21a */
-/*
- * Coefficients for approximation to erf in [0.84375,1.25]
- */
-pa0 = -2.3621185683e-03, /* 0xbb1acdc6 */
-pa1 = 4.1485610604e-01, /* 0x3ed46805 */
-pa2 = -3.7220788002e-01, /* 0xbebe9208 */
-pa3 = 3.1834661961e-01, /* 0x3ea2fe54 */
-pa4 = -1.1089469492e-01, /* 0xbde31cc2 */
-pa5 = 3.5478305072e-02, /* 0x3d1151b3 */
-pa6 = -2.1663755178e-03, /* 0xbb0df9c0 */
-qa1 = 1.0642088205e-01, /* 0x3dd9f331 */
-qa2 = 5.4039794207e-01, /* 0x3f0a5785 */
-qa3 = 7.1828655899e-02, /* 0x3d931ae7 */
-qa4 = 1.2617121637e-01, /* 0x3e013307 */
-qa5 = 1.3637083583e-02, /* 0x3c5f6e13 */
-qa6 = 1.1984500103e-02, /* 0x3c445aa3 */
-/*
- * Coefficients for approximation to erfc in [1.25,1/0.35]
- */
-ra0 = -9.8649440333e-03, /* 0xbc21a093 */
-ra1 = -6.9385856390e-01, /* 0xbf31a0b7 */
-ra2 = -1.0558626175e+01, /* 0xc128f022 */
-ra3 = -6.2375331879e+01, /* 0xc2798057 */
-ra4 = -1.6239666748e+02, /* 0xc322658c */
-ra5 = -1.8460508728e+02, /* 0xc3389ae7 */
-ra6 = -8.1287437439e+01, /* 0xc2a2932b */
-ra7 = -9.8143291473e+00, /* 0xc11d077e */
-sa1 = 1.9651271820e+01, /* 0x419d35ce */
-sa2 = 1.3765776062e+02, /* 0x4309a863 */
-sa3 = 4.3456588745e+02, /* 0x43d9486f */
-sa4 = 6.4538726807e+02, /* 0x442158c9 */
-sa5 = 4.2900814819e+02, /* 0x43d6810b */
-sa6 = 1.0863500214e+02, /* 0x42d9451f */
-sa7 = 6.5702495575e+00, /* 0x40d23f7c */
-sa8 = -6.0424413532e-02, /* 0xbd777f97 */
-/*
- * Coefficients for approximation to erfc in [1/.35,28]
- */
-rb0 = -9.8649431020e-03, /* 0xbc21a092 */
-rb1 = -7.9928326607e-01, /* 0xbf4c9dd4 */
-rb2 = -1.7757955551e+01, /* 0xc18e104b */
-rb3 = -1.6063638306e+02, /* 0xc320a2ea */
-rb4 = -6.3756646729e+02, /* 0xc41f6441 */
-rb5 = -1.0250950928e+03, /* 0xc480230b */
-rb6 = -4.8351919556e+02, /* 0xc3f1c275 */
-sb1 = 3.0338060379e+01, /* 0x41f2b459 */
-sb2 = 3.2579251099e+02, /* 0x43a2e571 */
-sb3 = 1.5367296143e+03, /* 0x44c01759 */
-sb4 = 3.1998581543e+03, /* 0x4547fdbb */
-sb5 = 2.5530502930e+03, /* 0x451f90ce */
-sb6 = 4.7452853394e+02, /* 0x43ed43a7 */
-sb7 = -2.2440952301e+01; /* 0xc1b38712 */
+ erx = 8.4506291151e-01, /* 0x3f58560b */
+ /*
+ * Coefficients for approximation to erf on [0,0.84375]
+ */
+ efx = 1.2837916613e-01, /* 0x3e0375d4 */
+ efx8= 1.0270333290e+00, /* 0x3f8375d4 */
+ pp0 = 1.2837916613e-01, /* 0x3e0375d4 */
+ pp1 = -3.2504209876e-01, /* 0xbea66beb */
+ pp2 = -2.8481749818e-02, /* 0xbce9528f */
+ pp3 = -5.7702702470e-03, /* 0xbbbd1489 */
+ pp4 = -2.3763017452e-05, /* 0xb7c756b1 */
+ qq1 = 3.9791721106e-01, /* 0x3ecbbbce */
+ qq2 = 6.5022252500e-02, /* 0x3d852a63 */
+ qq3 = 5.0813062117e-03, /* 0x3ba68116 */
+ qq4 = 1.3249473704e-04, /* 0x390aee49 */
+ qq5 = -3.9602282413e-06, /* 0xb684e21a */
+ /*
+ * Coefficients for approximation to erf in [0.84375,1.25]
+ */
+ pa0 = -2.3621185683e-03, /* 0xbb1acdc6 */
+ pa1 = 4.1485610604e-01, /* 0x3ed46805 */
+ pa2 = -3.7220788002e-01, /* 0xbebe9208 */
+ pa3 = 3.1834661961e-01, /* 0x3ea2fe54 */
+ pa4 = -1.1089469492e-01, /* 0xbde31cc2 */
+ pa5 = 3.5478305072e-02, /* 0x3d1151b3 */
+ pa6 = -2.1663755178e-03, /* 0xbb0df9c0 */
+ qa1 = 1.0642088205e-01, /* 0x3dd9f331 */
+ qa2 = 5.4039794207e-01, /* 0x3f0a5785 */
+ qa3 = 7.1828655899e-02, /* 0x3d931ae7 */
+ qa4 = 1.2617121637e-01, /* 0x3e013307 */
+ qa5 = 1.3637083583e-02, /* 0x3c5f6e13 */
+ qa6 = 1.1984500103e-02, /* 0x3c445aa3 */
+ /*
+ * Coefficients for approximation to erfc in [1.25,1/0.35]
+ */
+ ra0 = -9.8649440333e-03, /* 0xbc21a093 */
+ ra1 = -6.9385856390e-01, /* 0xbf31a0b7 */
+ ra2 = -1.0558626175e+01, /* 0xc128f022 */
+ ra3 = -6.2375331879e+01, /* 0xc2798057 */
+ ra4 = -1.6239666748e+02, /* 0xc322658c */
+ ra5 = -1.8460508728e+02, /* 0xc3389ae7 */
+ ra6 = -8.1287437439e+01, /* 0xc2a2932b */
+ ra7 = -9.8143291473e+00, /* 0xc11d077e */
+ sa1 = 1.9651271820e+01, /* 0x419d35ce */
+ sa2 = 1.3765776062e+02, /* 0x4309a863 */
+ sa3 = 4.3456588745e+02, /* 0x43d9486f */
+ sa4 = 6.4538726807e+02, /* 0x442158c9 */
+ sa5 = 4.2900814819e+02, /* 0x43d6810b */
+ sa6 = 1.0863500214e+02, /* 0x42d9451f */
+ sa7 = 6.5702495575e+00, /* 0x40d23f7c */
+ sa8 = -6.0424413532e-02, /* 0xbd777f97 */
+ /*
+ * Coefficients for approximation to erfc in [1/.35,28]
+ */
+ rb0 = -9.8649431020e-03, /* 0xbc21a092 */
+ rb1 = -7.9928326607e-01, /* 0xbf4c9dd4 */
+ rb2 = -1.7757955551e+01, /* 0xc18e104b */
+ rb3 = -1.6063638306e+02, /* 0xc320a2ea */
+ rb4 = -6.3756646729e+02, /* 0xc41f6441 */
+ rb5 = -1.0250950928e+03, /* 0xc480230b */
+ rb6 = -4.8351919556e+02, /* 0xc3f1c275 */
+ sb1 = 3.0338060379e+01, /* 0x41f2b459 */
+ sb2 = 3.2579251099e+02, /* 0x43a2e571 */
+ sb3 = 1.5367296143e+03, /* 0x44c01759 */
+ sb4 = 3.1998581543e+03, /* 0x4547fdbb */
+ sb5 = 2.5530502930e+03, /* 0x451f90ce */
+ sb6 = 4.7452853394e+02, /* 0x43ed43a7 */
+ sb7 = -2.2440952301e+01; /* 0xc1b38712 */
-#ifdef __STDC__
- float erff(float x)
-#else
- float erff(x)
- float x;
-#endif
+float erff(float x)
{
- int32_t hx,ix,i;
- float R,S,P,Q,s,y,z,r;
+ int32_t hx, ix, i;
+ float R, S, P, Q, s, y, z, r;
hx = __get_float_word(x);
- ix = hx&0x7fffffff;
- if(!(ix<0x7f800000L)) { /* erf(nan)=nan */
+ ix = hx & 0x7fffffff;
+ if (!(ix<0x7f800000L)) { /* erf(nan)=nan */
i = ((uint32_t)hx>>31)<<1;
return (float)(1-i)+one/x; /* erf(+-inf)=+-1 */
}
- if(ix < 0x3f580000) { /* |x|<0.84375 */
- if(ix < 0x31800000) { /* |x|<2**-28 */
- if (ix < 0x04000000)
+ if (ix < 0x3f580000) { /* |x|<0.84375 */
+ if (ix < 0x31800000) { /* |x|<2**-28 */
+ if (ix < 0x04000000)
/*avoid underflow */
return (float)0.125*((float)8.0*x+efx8*x);
return x + efx*x;
@@ -173,71 +162,75 @@
y = r/s;
return x + x*y;
}
- if(ix < 0x3fa00000) { /* 0.84375 <= |x| < 1.25 */
+ if (ix < 0x3fa00000) { /* 0.84375 <= |x| < 1.25 */
s = fabsf(x)-one;
P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
- if(hx>=0) return erx + P/Q; else return -erx - P/Q;
+ if (hx >= 0)
+ return erx + P/Q;
+ else
+ return -erx - P/Q;
}
if (ix >= 0x40c00000) { /* inf>|x|>=6 */
- if(hx>=0) return one-tiny; else return tiny-one;
+ if (hx >= 0)
+ return one-tiny;
+ else
+ return tiny-one;
}
x = fabsf(x);
s = one/(x*x);
- if(ix< 0x4036DB6E) { /* |x| < 1/0.35 */
- R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
+ if (ix< 0x4036DB6E) { /* |x| < 1/0.35 */
+ R = ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
ra5+s*(ra6+s*ra7))))));
- S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
+ S = one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
sa5+s*(sa6+s*(sa7+s*sa8)))))));
} else { /* |x| >= 1/0.35 */
- R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
+ R = rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
rb5+s*rb6)))));
- S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
+ S = one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
sb5+s*(sb6+s*sb7))))));
}
z = x;
__trunc_float_word (&z);
- r = __ieee754_expf(-z*z-(float)0.5625)*__ieee754_expf((z-x)*(z+x)+R/S);
- if(hx>=0) return one-r/x; else return r/x-one;
+ r = __ieee754_expf(-z*z-(float)0.5625)*__ieee754_expf((z-x)*(z+x)+R/S);
+ if (hx >= 0)
+ return one-r/x;
+ else
+ return r/x-one;
}
-#ifdef __STDC__
- float erfcf(float x)
-#else
- float erfcf(x)
- float x;
-#endif
+float erfcf(float x)
{
- int32_t hx,ix;
- float R,S,P,Q,s,y,z,r;
+ int32_t hx, ix;
+ float R, S, P, Q, s, y, z, r;
hx = __get_float_word(x);
- ix = hx&0x7fffffff;
- if(!(ix<0x7f800000L)) { /* erfc(nan)=nan */
+ ix = hx & 0x7fffffff;
+ if (!(ix < 0x7f800000L)) { /* erfc(nan)=nan */
/* erfc(+-inf)=0,2 */
return (float)(((uint32_t)hx>>31)<<1)+one/x;
}
- if(ix < 0x3f580000) { /* |x|<0.84375 */
- if(ix < 0x23800000) /* |x|<2**-56 */
+ if (ix < 0x3f580000) { /* |x|<0.84375 */
+ if (ix < 0x23800000) /* |x|<2**-56 */
return one-x;
z = x*x;
r = pp0+z*(pp1+z*(pp2+z*(pp3+z*pp4)));
s = one+z*(qq1+z*(qq2+z*(qq3+z*(qq4+z*qq5))));
y = r/s;
- if(hx < 0x3e800000) { /* x<1/4 */
+ if (hx < 0x3e800000) { /* x<1/4 */
return one-(x+x*y);
} else {
r = x*y;
r += (x-half);
- return half - r ;
+ return half - r ;
}
}
- if(ix < 0x3fa00000) { /* 0.84375 <= |x| < 1.25 */
+ if (ix < 0x3fa00000) { /* 0.84375 <= |x| < 1.25 */
s = fabsf(x)-one;
P = pa0+s*(pa1+s*(pa2+s*(pa3+s*(pa4+s*(pa5+s*pa6)))));
Q = one+s*(qa1+s*(qa2+s*(qa3+s*(qa4+s*(qa5+s*qa6)))));
- if(hx>=0) {
+ if (hx >= 0) {
z = one-erx; return z - P/Q;
} else {
z = erx+P/Q; return one+z;
@@ -247,27 +240,34 @@
if (ix < 0x41e00000) { /* |x|<28 */
x = fabsf(x);
s = one/(x*x);
- if(ix< 0x4036DB6D) { /* |x| < 1/.35 ~ 2.857143*/
- R=ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
+ if (ix< 0x4036DB6D) { /* |x| < 1/.35 ~ 2.857143*/
+ R = ra0+s*(ra1+s*(ra2+s*(ra3+s*(ra4+s*(
ra5+s*(ra6+s*ra7))))));
- S=one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
+ S = one+s*(sa1+s*(sa2+s*(sa3+s*(sa4+s*(
sa5+s*(sa6+s*(sa7+s*sa8)))))));
} else { /* |x| >= 1/.35 ~ 2.857143 */
- if(hx<0&&ix>=0x40c00000) return two-tiny;/* x < -6 */
- R=rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
+ if (hx < 0 && ix >= 0x40c00000)
+ return two-tiny;/* x < -6 */
+ R = rb0+s*(rb1+s*(rb2+s*(rb3+s*(rb4+s*(
rb5+s*rb6)))));
- S=one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
+ S = one+s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(
sb5+s*(sb6+s*sb7))))));
}
z = x;
__trunc_float_word (&z);
- r = __ieee754_expf(-z*z-(float)0.5625)*
+ r = __ieee754_expf(-z*z-(float)0.5625)*
__ieee754_expf((z-x)*(z+x)+R/S);
- if(hx>0) return r/x; else return two-r/x;
+ if (hx > 0)
+ return r/x;
+ else
+ return two-r/x;
} else {
/* set range error */
errno = ERANGE;
- if(hx>0) return tiny*tiny; else return two-tiny;
+ if (hx > 0)
+ return tiny*tiny;
+ else
+ return two-tiny;
}
}
diff --git a/mingw-w64-crt/math/sinhf.c b/mingw-w64-crt/math/sinhf.c
index 4fde4e6..81afd18 100644
--- a/mingw-w64-crt/math/sinhf.c
+++ b/mingw-w64-crt/math/sinhf.c
@@ -5,4 +5,6 @@
*/
#include <math.h>
float sinhf (float x)
- {return (float) sinh (x);}
+{
+ return (float) sinh (x);
+}
diff --git a/mingw-w64-crt/math/sinhl.c b/mingw-w64-crt/math/sinhl.c
index 7d51fb2..a3da29c 100644
--- a/mingw-w64-crt/math/sinhl.c
+++ b/mingw-w64-crt/math/sinhl.c
@@ -62,8 +62,8 @@
long double a;
#ifdef MINUSZERO
- if( x == 0.0 )
- return(x);
+ if (x == 0.0)
+ return (x);
#endif
#ifdef NANS
if (isnanl(x))
@@ -72,39 +72,40 @@
}
#endif
a = fabsl(x);
- if( (x > (MAXLOGL + LOGE2L)) || (x > -(MINLOGL-LOGE2L) ) )
+ if ((x > (MAXLOGL + LOGE2L)) || (x > -(MINLOGL-LOGE2L)))
{
- mtherr( "sinhl", DOMAIN );
+ mtherr("sinhl", DOMAIN);
_SET_ERRNO(ERANGE);
#ifdef INFINITIES
- if( x > 0.0L )
- return( INFINITYL );
+ if (x > 0.0L)
+ return (INFINITYL);
else
- return( -INFINITYL );
+ return (-INFINITYL);
#else
- if( x > 0.0L )
- return( MAXNUML );
+ if (x > 0.0L)
+ return (MAXNUML);
else
- return( -MAXNUML );
+ return (-MAXNUML);
#endif
}
- if( a > 1.0L )
+ if (a > 1.0L)
{
- if( a >= (MAXLOGL - LOGE2L) )
+ if (a >= (MAXLOGL - LOGE2L))
{
a = expl(0.5L*a);
a = (0.5L * a) * a;
- if( x < 0.0L )
+ if (x < 0.0L)
a = -a;
- return(a);
+ return (a);
}
a = expl(a);
a = 0.5L*a - (0.5L/a);
- if( x < 0.0L )
+ if (x < 0.0L)
a = -a;
- return(a);
+ return (a);
}
a *= a;
- return( x + x * a * (polevll(a,P,3)/polevll(a,Q,4)) );
+ return (x + x * a * (polevll(a,P,3)/polevll(a,Q,4)));
}
+
diff --git a/mingw-w64-crt/math/tanhf.c b/mingw-w64-crt/math/tanhf.c
index 738ed12..6a5f9d8 100644
--- a/mingw-w64-crt/math/tanhf.c
+++ b/mingw-w64-crt/math/tanhf.c
@@ -5,4 +5,6 @@
*/
#include <math.h>
float tanhf (float x)
- {return (float) tanh (x);}
+{
+ return (float) tanh (x);
+}
diff --git a/mingw-w64-crt/math/tanhl.c b/mingw-w64-crt/math/tanhl.c
index 2b6da46..d22a3b7 100644
--- a/mingw-w64-crt/math/tanhl.c
+++ b/mingw-w64-crt/math/tanhl.c
@@ -55,8 +55,8 @@
long double s, z;
#ifdef MINUSZERO
- if( x == 0.0L )
- return(x);
+ if (x == 0.0L)
+ return (x);
#endif
if (isnanl(x))
{
@@ -65,19 +65,19 @@
}
z = fabsl(x);
- if( z > 0.5L * MAXLOGL )
+ if (z > 0.5L * MAXLOGL)
{
_SET_ERRNO (ERANGE);
- if( x > 0 )
- return( 1.0L );
+ if (x > 0)
+ return (1.0L);
else
- return( -1.0L );
+ return (-1.0L);
}
- if( z >= 0.625L )
+ if (z >= 0.625L)
{
s = expl(2.0*z);
z = 1.0L - 2.0/(s + 1.0L);
- if( x < 0 )
+ if (x < 0)
z = -z;
}
else
@@ -87,5 +87,6 @@
z = x * s * z;
z = x + z;
}
- return( z );
+ return (z);
}
+
diff --git a/mingw-w64-crt/math/tgamma.c b/mingw-w64-crt/math/tgamma.c
index c93fbed..55c33dc 100644
--- a/mingw-w64-crt/math/tgamma.c
+++ b/mingw-w64-crt/math/tgamma.c
@@ -119,29 +119,29 @@
#define SQTPI SQT.d
#endif
-static double stirf ( double );
+static double stirf (double);
/* Gamma function computed by Stirling's formula.
* The polynomial STIR is valid for 33 <= x <= 172.
*/
static double stirf(double x)
{
- double y, w, v;
+ double y, w, v;
- w = 1.0/x;
- w = 1.0 + w * polevl( w, STIR, 4 );
- y = exp(x);
- if( x > MAXSTIR )
- { /* Avoid overflow in pow() */
- v = pow( x, 0.5 * x - 0.25 );
- y = v * (v / y);
- }
- else
- {
- y = pow( x, x - 0.5 ) / y;
- }
- y = SQTPI * y * w;
- return( y );
+ w = 1.0/x;
+ w = 1.0 + w * polevl(w, STIR, 4);
+ y = exp(x);
+ if (x > MAXSTIR)
+ { /* Avoid overflow in pow() */
+ v = pow(x, 0.5 * x - 0.25);
+ y = v * (v / y);
+ }
+ else
+ {
+ y = pow(x, x - 0.5) / y;
+ }
+ y = SQTPI * y * w;
+ return (y);
}
@@ -149,117 +149,117 @@
double __tgamma_r(double x, int *sgngam)
{
- double p, q, z;
- int i;
+ double p, q, z;
+ int i;
- *sgngam = 1;
+ *sgngam = 1;
#ifdef NANS
- if( isnan(x) )
- return(x);
+ if (isnan(x))
+ return (x);
#endif
#ifdef INFINITIES
#ifdef NANS
- if( x == INFINITY )
- return(x);
- if( x == -INFINITY )
- return(NAN);
+ if (x == INFINITY)
+ return (x);
+ if (x == -INFINITY)
+ return (NAN);
#else
- if( !isfinite(x) )
- return(x);
+ if (!isfinite(x))
+ return (x);
#endif
#endif
- q = fabs(x);
+ q = fabs(x);
- if( q > 33.0 )
- {
- if( x < 0.0 )
- {
- p = floor(q);
- if( p == q )
- {
+ if (q > 33.0)
+ {
+ if (x < 0.0)
+ {
+ p = floor(q);
+ if (p == q)
+ {
gsing:
- _SET_ERRNO(EDOM);
- mtherr( "tgamma", SING );
+ _SET_ERRNO(EDOM);
+ mtherr("tgamma", SING);
#ifdef INFINITIES
- return (INFINITY);
+ return (INFINITY);
#else
- return (MAXNUM);
+ return (MAXNUM);
#endif
- }
- i = p;
- if( (i & 1) == 0 )
- *sgngam = -1;
- z = q - p;
- if( z > 0.5 )
- {
- p += 1.0;
- z = q - p;
- }
- z = q * sin( PI * z );
- if( z == 0.0 )
- {
- _SET_ERRNO(ERANGE);
- mtherr( "tgamma", OVERFLOW );
+ }
+ i = p;
+ if ((i & 1) == 0)
+ *sgngam = -1;
+ z = q - p;
+ if (z > 0.5)
+ {
+ p += 1.0;
+ z = q - p;
+ }
+ z = q * sin(PI * z);
+ if (z == 0.0)
+ {
+ _SET_ERRNO(ERANGE);
+ mtherr("tgamma", OVERFLOW);
#ifdef INFINITIES
- return( *sgngam * INFINITY);
+ return (*sgngam * INFINITY);
#else
- return( *sgngam * MAXNUM);
+ return (*sgngam * MAXNUM);
#endif
- }
- z = fabs(z);
- z = PI/(z * stirf(q) );
- }
- else
- {
- z = stirf(x);
- }
- return( *sgngam * z );
- }
+ }
+ z = fabs(z);
+ z = PI/(z * stirf(q));
+ }
+ else
+ {
+ z = stirf(x);
+ }
+ return (*sgngam * z);
+ }
- z = 1.0;
- while( x >= 3.0 )
- {
- x -= 1.0;
- z *= x;
- }
+ z = 1.0;
+ while (x >= 3.0)
+ {
+ x -= 1.0;
+ z *= x;
+ }
- while( x < 0.0 )
- {
- if( x > -1.E-9 )
- goto Small;
- z /= x;
- x += 1.0;
- }
+ while (x < 0.0)
+ {
+ if (x > -1.E-9)
+ goto Small;
+ z /= x;
+ x += 1.0;
+ }
- while( x < 2.0 )
- {
- if( x < 1.e-9 )
- goto Small;
- z /= x;
- x += 1.0;
- }
+ while (x < 2.0)
+ {
+ if (x < 1.e-9)
+ goto Small;
+ z /= x;
+ x += 1.0;
+ }
- if( x == 2.0 )
- return(z);
+ if (x == 2.0)
+ return (z);
- x -= 2.0;
- p = polevl( x, P, 6 );
- q = polevl( x, Q, 7 );
- return( z * p / q );
+ x -= 2.0;
+ p = polevl( x, P, 6 );
+ q = polevl( x, Q, 7 );
+ return (z * p / q);
Small:
- if( x == 0.0 )
- {
- goto gsing;
- }
- else
- return( z/((1.0 + 0.5772156649015329 * x) * x) );
+ if (x == 0.0)
+ {
+ goto gsing;
+ }
+ else
+ return (z/((1.0 + 0.5772156649015329 * x) * x));
}
/* This is the C99 version */
-
double tgamma(double x)
{
- int local_sgngam=0;
- return (__tgamma_r(x, &local_sgngam));
+ int local_sgngam = 0;
+ return (__tgamma_r(x, &local_sgngam));
}
+
diff --git a/mingw-w64-crt/math/tgammaf.c b/mingw-w64-crt/math/tgammaf.c
index d199cf6..ff335b2 100644
--- a/mingw-w64-crt/math/tgammaf.c
+++ b/mingw-w64-crt/math/tgammaf.c
@@ -13,9 +13,9 @@
* relative error < 1.9e-11
*/
static const float STIR[] = {
--2.705194986674176E-003,
- 3.473255786154910E-003,
- 8.333331788340907E-002,
+ -2.705194986674176E-003,
+ 3.473255786154910E-003,
+ 8.333331788340907E-002,
};
static const float MAXSTIR = 26.77;
static const float SQTPIF = 2.50662827463100050242; /* sqrt( 2 pi ) */
@@ -28,168 +28,167 @@
*/
static float stirf( float x )
{
-float y, w, v;
+ float y, w, v;
-w = 1.0/x;
-w = 1.0 + w * polevlf( w, STIR, 2 );
-y = expf( -x );
-if( x > MAXSTIR )
+ w = 1.0/x;
+ w = 1.0 + w * polevlf(w, STIR, 2);
+ y = expf(-x);
+ if (x > MAXSTIR)
{ /* Avoid overflow in pow() */
- v = powf( x, 0.5 * x - 0.25 );
- y *= v;
- y *= v;
+ v = powf(x, 0.5 * x - 0.25);
+ y *= v;
+ y *= v;
}
-else
+ else
{
- y = powf( x, x - 0.5 ) * y;
+ y = powf(x, x - 0.5) * y;
}
-y = SQTPIF * y * w;
-return( y );
+ y = SQTPIF * y * w;
+ return (y);
}
/* gamma(x+2), 0 < x < 1 */
static const float P[] = {
- 1.536830450601906E-003,
- 5.397581592950993E-003,
- 4.130370201859976E-003,
- 7.232307985516519E-002,
- 8.203960091619193E-002,
- 4.117857447645796E-001,
- 4.227867745131584E-001,
- 9.999999822945073E-001,
+ 1.536830450601906E-003,
+ 5.397581592950993E-003,
+ 4.130370201859976E-003,
+ 7.232307985516519E-002,
+ 8.203960091619193E-002,
+ 4.117857447645796E-001,
+ 4.227867745131584E-001,
+ 9.999999822945073E-001,
};
float __tgammaf_r( float x, int* sgngamf);
float __tgammaf_r( float x, int* sgngamf)
{
-float p, q, z, nz;
-int i, direction, negative;
+ float p, q, z, nz;
+ int i, direction, negative;
#ifdef NANS
-if( isnan(x) )
- return(x);
+ if (isnan(x))
+ return (x);
#endif
#ifdef INFINITIES
#ifdef NANS
-if( x == INFINITYF )
- return(x);
-if( x == -INFINITYF )
- return(NANF);
+ if (x == INFINITYF)
+ return (x);
+ if (x == -INFINITYF)
+ return (NANF);
#else
-if( !isfinite(x) )
- return(x);
+ if (!isfinite(x))
+ return (x);
#endif
#endif
-*sgngamf = 1;
-negative = 0;
-nz = 0.0;
-if( x < 0.0 )
+ *sgngamf = 1;
+ negative = 0;
+ nz = 0.0;
+ if (x < 0.0)
{
- negative = 1;
- q = -x;
- p = floorf(q);
- if( p == q )
+ negative = 1;
+ q = -x;
+ p = floorf(q);
+ if (p == q)
{
gsing:
- _SET_ERRNO(EDOM);
- mtherr( "tgammaf", SING );
+ _SET_ERRNO(EDOM);
+ mtherr("tgammaf", SING);
#ifdef INFINITIES
- return (INFINITYF);
+ return (INFINITYF);
#else
- return (MAXNUMF);
+ return (MAXNUMF);
#endif
- }
- i = p;
- if( (i & 1) == 0 )
- *sgngamf = -1;
- nz = q - p;
- if( nz > 0.5 )
- {
- p += 1.0;
+ }
+ i = p;
+ if ((i & 1) == 0)
+ *sgngamf = -1;
nz = q - p;
- }
- nz = q * sinf( PIF * nz );
- if( nz == 0.0 )
+ if (nz > 0.5)
{
- _SET_ERRNO(ERANGE);
- mtherr( "tgamma", OVERFLOW );
+ p += 1.0;
+ nz = q - p;
+ }
+ nz = q * sinf(PIF * nz);
+ if (nz == 0.0)
+ {
+ _SET_ERRNO(ERANGE);
+ mtherr("tgamma", OVERFLOW);
#ifdef INFINITIES
- return( *sgngamf * INFINITYF);
+ return(*sgngamf * INFINITYF);
#else
- return( *sgngamf * MAXNUMF);
+ return(*sgngamf * MAXNUMF);
#endif
}
- if( nz < 0 )
- nz = -nz;
- x = q;
+ if (nz < 0)
+ nz = -nz;
+ x = q;
}
-if( x >= 10.0 )
+ if (x >= 10.0)
{
- z = stirf(x);
+ z = stirf(x);
}
-if( x < 2.0 )
- direction = 1;
-else
- direction = 0;
-z = 1.0;
-while( x >= 3.0 )
+ if (x < 2.0)
+ direction = 1;
+ else
+ direction = 0;
+ z = 1.0;
+ while (x >= 3.0)
{
- x -= 1.0;
- z *= x;
+ x -= 1.0;
+ z *= x;
}
-/*
-while( x < 0.0 )
+ /*
+ while (x < 0.0)
{
- if( x > -1.E-4 )
- goto Small;
- z *=x;
- x += 1.0;
+ if (x > -1.E-4)
+ goto Small;
+ z *=x;
+ x += 1.0;
}
-*/
-while( x < 2.0 )
+ */
+ while (x < 2.0)
{
- if( x < 1.e-4 )
- goto Small;
- z *=x;
- x += 1.0;
+ if (x < 1.e-4)
+ goto Small;
+ z *=x;
+ x += 1.0;
}
-if( direction )
- z = 1.0/z;
+ if (direction)
+ z = 1.0/z;
-if( x == 2.0 )
- return(z);
+ if (x == 2.0)
+ return (z);
-x -= 2.0;
-p = z * polevlf( x, P, 7 );
+ x -= 2.0;
+ p = z * polevlf(x, P, 7);
gdone:
-
-if( negative )
+ if (negative)
{
- p = *sgngamf * PIF/(nz * p );
+ p = *sgngamf * PIF/(nz * p );
}
-return(p);
+ return (p);
Small:
-if( x == 0.0 )
+ if (x == 0.0)
{
- goto gsing;
+ goto gsing;
}
-else
+ else
{
- p = z / ((1.0 + 0.5772156649015329 * x) * x);
- goto gdone;
+ p = z / ((1.0 + 0.5772156649015329 * x) * x);
+ goto gdone;
}
}
/* This is the C99 version */
-
float tgammaf(float x)
{
- int local_sgngamf=0;
- return (__tgammaf_r(x, &local_sgngamf));
+ int local_sgngamf = 0;
+ return (__tgammaf_r(x, &local_sgngamf));
}
+
diff --git a/mingw-w64-crt/math/tgammal.c b/mingw-w64-crt/math/tgammal.c
index 3aa9129..f9ca66b 100644
--- a/mingw-w64-crt/math/tgammal.c
+++ b/mingw-w64-crt/math/tgammal.c
@@ -228,168 +228,166 @@
};
#endif
-static long double stirf ( long double);
+static long double stirf (long double);
/* Gamma function computed by Stirling's formula. */
static long double stirf(long double x)
{
- long double y, w, v;
+ long double y, w, v;
- w = 1.0L/x;
- /* For large x, use rational coefficients from the analytical expansion. */
- if( x > 1024.0L )
- w = (((((6.97281375836585777429E-5L * w
- + 7.84039221720066627474E-4L) * w
- - 2.29472093621399176955E-4L) * w
- - 2.68132716049382716049E-3L) * w
- + 3.47222222222222222222E-3L) * w
- + 8.33333333333333333333E-2L) * w
- + 1.0L;
- else
- w = 1.0L + w * polevll( w, STIR, 8 );
- y = expl(x);
- if( x > MAXSTIR )
- { /* Avoid overflow in pow() */
- v = powl( x, 0.5L * x - 0.25L );
- y = v * (v / y);
- }
- else
- {
- y = powl( x, x - 0.5L ) / y;
- }
- y = SQTPI * y * w;
- return( y );
+ w = 1.0L/x;
+ /* For large x, use rational coefficients from the analytical expansion. */
+ if (x > 1024.0L)
+ w = (((((6.97281375836585777429E-5L * w
+ + 7.84039221720066627474E-4L) * w
+ - 2.29472093621399176955E-4L) * w
+ - 2.68132716049382716049E-3L) * w
+ + 3.47222222222222222222E-3L) * w
+ + 8.33333333333333333333E-2L) * w
+ + 1.0L;
+ else
+ w = 1.0L + w * polevll( w, STIR, 8 );
+ y = expl(x);
+ if (x > MAXSTIR)
+ { /* Avoid overflow in pow() */
+ v = powl(x, 0.5L * x - 0.25L);
+ y = v * (v / y);
+ }
+ else
+ {
+ y = powl(x, x - 0.5L) / y;
+ }
+ y = SQTPI * y * w;
+ return (y);
}
-long double __tgammal_r(long double x, int* sgngaml);
+long double __tgammal_r(long double, int *);
long double __tgammal_r(long double x, int* sgngaml)
{
-long double p, q, z;
-int i;
+ long double p, q, z;
+ int i;
-*sgngaml = 1;
+ *sgngaml = 1;
#ifdef NANS
-if( isnanl(x) )
- return(NANL);
+ if (isnanl(x))
+ return (NANL);
#endif
#ifdef INFINITIES
#ifdef NANS
-if( x == INFINITYL )
- return(x);
-if( x == -INFINITYL )
- return(NANL);
+ if (x == INFINITYL)
+ return (x);
+ if (x == -INFINITYL)
+ return (NANL);
#else
-if( !isfinite(x) )
- return(x);
+ if (!isfinite(x))
+ return (x);
#endif
#endif
-q = fabsl(x);
+ q = fabsl(x);
-if( q > 13.0L )
+ if (q > 13.0L)
{
- if( q > MAXGAML )
- goto goverf;
- if( x < 0.0L )
+ if (q > MAXGAML)
+ goto goverf;
+ if (x < 0.0L)
{
- p = floorl(q);
- if( p == q )
+ p = floorl(q);
+ if (p == q)
{
gsing:
- _SET_ERRNO(EDOM);
- mtherr( "tgammal", SING );
+ _SET_ERRNO(EDOM);
+ mtherr("tgammal", SING);
#ifdef INFINITIES
- return (INFINITYL);
+ return (INFINITYL);
#else
- return( *sgngaml * MAXNUML);
+ return (*sgngaml * MAXNUML);
#endif
}
- i = p;
- if( (i & 1) == 0 )
- *sgngaml = -1;
- z = q - p;
- if( z > 0.5L )
- {
- p += 1.0L;
+ i = p;
+ if ((i & 1) == 0)
+ *sgngaml = -1;
z = q - p;
+ if (z > 0.5L)
+ {
+ p += 1.0L;
+ z = q - p;
}
- z = q * sinl( PIL * z );
- z = fabsl(z) * stirf(q);
- if( z <= PIL/MAXNUML )
+ z = q * sinl(PIL * z);
+ z = fabsl(z) * stirf(q);
+ if (z <= PIL/MAXNUML)
{
goverf:
- _SET_ERRNO(ERANGE);
- mtherr( "tgammal", OVERFLOW );
+ _SET_ERRNO(ERANGE);
+ mtherr("tgammal", OVERFLOW);
#ifdef INFINITIES
- return( *sgngaml * INFINITYL);
+ return(*sgngaml * INFINITYL);
#else
- return( *sgngaml * MAXNUML);
+ return(*sgngaml * MAXNUML);
#endif
}
- z = PIL/z;
+ z = PIL/z;
}
- else
+ else
{
- z = stirf(x);
+ z = stirf(x);
}
- return( *sgngaml * z );
+ return (*sgngaml * z);
}
-z = 1.0L;
-while( x >= 3.0L )
+ z = 1.0L;
+ while (x >= 3.0L)
{
- x -= 1.0L;
- z *= x;
+ x -= 1.0L;
+ z *= x;
}
-while( x < -0.03125L )
+ while (x < -0.03125L)
{
- z /= x;
- x += 1.0L;
+ z /= x;
+ x += 1.0L;
}
-if( x <= 0.03125L )
- goto Small;
+ if (x <= 0.03125L)
+ goto Small;
-while( x < 2.0L )
+ while (x < 2.0L)
{
- z /= x;
- x += 1.0L;
+ z /= x;
+ x += 1.0L;
}
-if( x == 2.0L )
- return(z);
+ if (x == 2.0L)
+ return (z);
-x -= 2.0L;
-p = polevll( x, P, 7 );
-q = polevll( x, Q, 8 );
-return( z * p / q );
+ x -= 2.0L;
+ p = polevll( x, P, 7 );
+ q = polevll( x, Q, 8 );
+ return (z * p / q);
Small:
-if( x == 0.0L )
+ if (x == 0.0L)
{
- goto gsing;
+ goto gsing;
}
-else
- {
- if( x < 0.0L )
- {
- x = -x;
- q = z / (x * polevll( x, SN, 8 ));
- }
else
- q = z / (x * polevll( x, S, 8 ));
+ {
+ if (x < 0.0L)
+ {
+ x = -x;
+ q = z / (x * polevll(x, SN, 8));
+ }
+ else
+ q = z / (x * polevll(x, S, 8));
}
-return q;
+ return q;
}
-
/* This is the C99 version. */
-
long double tgammal(long double x)
{
- int local_sgngaml=0;
- return (__tgammal_r(x, &local_sgngaml));
+ int local_sgngaml = 0;
+ return (__tgammal_r(x, &local_sgngaml));
}
