[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[lmi-commits] [lmi] master 324e171e 02/22: Rectify whitespace
|
From: |
Greg Chicares |
|
Subject: |
[lmi-commits] [lmi] master 324e171e 02/22: Rectify whitespace |
|
Date: |
Fri, 20 May 2022 22:43:41 -0400 (EDT) |
branch: master
commit 324e171ef154c75b5be79fc335f196c03ba39ec4
Author: Gregory W. Chicares <gchicares@sbcglobal.net>
Commit: Gregory W. Chicares <gchicares@sbcglobal.net>
Rectify whitespace
---
fdlibm_expm1.c.txt | 378 ++++++++++++++++++++++++++---------------------------
fdlibm_log1p.c.txt | 292 ++++++++++++++++++++---------------------
2 files changed, 335 insertions(+), 335 deletions(-)
diff --git a/fdlibm_expm1.c.txt b/fdlibm_expm1.c.txt
index 6f0bd778..5b88d369 100644
--- a/fdlibm_expm1.c.txt
+++ b/fdlibm_expm1.c.txt
@@ -5,7 +5,7 @@
* Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved.
*
* Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
+ * software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
@@ -14,8 +14,8 @@
but these catch some common cases. */
#if defined(i386) || defined(i486) || \
- defined(intel) || defined(x86) || defined(i86pc) || \
- defined(__alpha) || defined(__osf__)
+ defined(intel) || defined(x86) || defined(i86pc) || \
+ defined(__alpha) || defined(__osf__)
#define __LITTLE_ENDIAN
#endif
@@ -32,9 +32,9 @@
#endif
#ifdef __STDC__
-#define __P(p) p
+#define __P(p) p
#else
-#define __P(p) ()
+#define __P(p) ()
#endif
/*
@@ -43,20 +43,20 @@
extern int signgam;
-#define MAXFLOAT ((float)3.40282346638528860e+38)
+#define MAXFLOAT ((float)3.40282346638528860e+38)
enum fdversion {fdlibm_ieee = -1, fdlibm_svid, fdlibm_xopen, fdlibm_posix};
#define _LIB_VERSION_TYPE enum fdversion
-#define _LIB_VERSION _fdlib_version
+#define _LIB_VERSION _fdlib_version
-/* if global variable _LIB_VERSION is not desirable, one may
- * change the following to be a constant by:
- * #define _LIB_VERSION_TYPE const enum version
+/* if global variable _LIB_VERSION is not desirable, one may
+ * change the following to be a constant by:
+ * #define _LIB_VERSION_TYPE const enum version
* In that case, after one initializes the value _LIB_VERSION (see
* s_lib_version.c) during compile time, it cannot be modified
* in the middle of a program
- */
+ */
extern _LIB_VERSION_TYPE _LIB_VERSION;
#define _IEEE_ fdlibm_ieee
@@ -65,28 +65,28 @@ extern _LIB_VERSION_TYPE _LIB_VERSION;
#define _POSIX_ fdlibm_posix
struct exception {
- int type;
- char *name;
- double arg1;
- double arg2;
- double retval;
+ int type;
+ char *name;
+ double arg1;
+ double arg2;
+ double retval;
};
-#define HUGE MAXFLOAT
+#define HUGE MAXFLOAT
-/*
+/*
* set X_TLOSS = pi*2**52, which is possibly defined in <values.h>
* (one may replace the following line by "#include <values.h>")
*/
-#define X_TLOSS 1.41484755040568800000e+16
+#define X_TLOSS 1.41484755040568800000e+16
-#define DOMAIN 1
-#define SING 2
-#define OVERFLOW 3
-#define UNDERFLOW 4
-#define TLOSS 5
-#define PLOSS 6
+#define DOMAIN 1
+#define SING 2
+#define OVERFLOW 3
+#define UNDERFLOW 4
+#define TLOSS 5
+#define PLOSS 6
/*
* ANSI/POSIX
@@ -173,16 +173,16 @@ extern double log1p __P((double));
#ifdef _REENTRANT
extern double gamma_r __P((double, int *));
extern double lgamma_r __P((double, int *));
-#endif /* _REENTRANT */
+#endif /* _REENTRANT */
/* ieee style elementary functions */
-extern double __ieee754_sqrt __P((double));
-extern double __ieee754_acos __P((double));
-extern double __ieee754_acosh __P((double));
-extern double __ieee754_log __P((double));
-extern double __ieee754_atanh __P((double));
-extern double __ieee754_asin __P((double));
-extern double __ieee754_atan2 __P((double,double));
+extern double __ieee754_sqrt __P((double));
+extern double __ieee754_acos __P((double));
+extern double __ieee754_acosh __P((double));
+extern double __ieee754_log __P((double));
+extern double __ieee754_atanh __P((double));
+extern double __ieee754_asin __P((double));
+extern double __ieee754_atan2 __P((double,double));
extern double __ieee754_exp __P((double));
extern double __ieee754_cosh __P((double));
extern double __ieee754_fmod __P((double,double));
@@ -209,7 +209,7 @@ extern double __ieee754_scalb __P((double,double));
#endif
/* fdlibm kernel function */
-extern double __kernel_standard __P((double,double,int));
+extern double __kernel_standard __P((double,double,int));
extern double __kernel_sin __P((double,double,int));
extern double __kernel_cos __P((double,double));
extern double __kernel_tan __P((double,double,int));
@@ -221,7 +221,7 @@ extern int __kernel_rem_pio2
__P((double*,double*,int,int,int,const int*));
* Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved.
*
* Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
+ * software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
@@ -231,92 +231,92 @@ extern int __kernel_rem_pio2
__P((double*,double*,int,int,int,const int*));
*
* Method
* 1. Argument reduction:
- * Given x, find r and integer k such that
+ * Given x, find r and integer k such that
*
- * x = k*ln2 + r, |r| <= 0.5*ln2 ~ 0.34658
+ * x = k*ln2 + r, |r| <= 0.5*ln2 ~ 0.34658
*
- * Here a correction term c will be computed to compensate
- * the error in r when rounded to a floating-point number.
+ * Here a correction term c will be computed to compensate
+ * the error in r when rounded to a floating-point number.
*
* 2. Approximating expm1(r) by a special rational function on
- * the interval [0,0.34658]:
- * Since
- * r*(exp(r)+1)/(exp(r)-1) = 2+ r^2/6 - r^4/360 + ...
- * we define R1(r*r) by
- * r*(exp(r)+1)/(exp(r)-1) = 2+ r^2/6 * R1(r*r)
- * That is,
- * R1(r**2) = 6/r *((exp(r)+1)/(exp(r)-1) - 2/r)
- * = 6/r * ( 1 + 2.0*(1/(exp(r)-1) - 1/r))
- * = 1 - r^2/60 + r^4/2520 - r^6/100800 + ...
- * We use a special Remes algorithm on [0,0.347] to generate
- * a polynomial of degree 5 in r*r to approximate R1. The
- * maximum error of this polynomial approximation is bounded
- * by 2**-61. In other words,
- * R1(z) ~ 1.0 + Q1*z + Q2*z**2 + Q3*z**3 + Q4*z**4 + Q5*z**5
- * where Q1 = -1.6666666666666567384E-2,
- * Q2 = 3.9682539681370365873E-4,
- * Q3 = -9.9206344733435987357E-6,
- * Q4 = 2.5051361420808517002E-7,
- * Q5 = -6.2843505682382617102E-9;
- * (where z=r*r, and the values of Q1 to Q5 are listed below)
- * with error bounded by
- * | 5 | -61
- * | 1.0+Q1*z+...+Q5*z - R1(z) | <= 2
- * | |
- *
- * expm1(r) = exp(r)-1 is then computed by the following
- * specific way which minimize the accumulation rounding error:
- * 2 3
- * r r [ 3 - (R1 + R1*r/2) ]
- * expm1(r) = r + --- + --- * [--------------------]
- * 2 2 [ 6 - r*(3 - R1*r/2) ]
- *
- * To compensate the error in the argument reduction, we use
- * expm1(r+c) = expm1(r) + c + expm1(r)*c
- * ~ expm1(r) + c + r*c
- * Thus c+r*c will be added in as the correction terms for
- * expm1(r+c). Now rearrange the term to avoid optimization
- * screw up:
- * ( 2 2 )
- * ({ ( r [ R1 - (3 - R1*r/2) ] ) } r )
- * expm1(r+c)~r - ({r*(--- * [--------------------]-c)-c} - --- )
- * ({ ( 2 [ 6 - r*(3 - R1*r/2) ] ) } 2 )
+ * the interval [0,0.34658]:
+ * Since
+ * r*(exp(r)+1)/(exp(r)-1) = 2+ r^2/6 - r^4/360 + ...
+ * we define R1(r*r) by
+ * r*(exp(r)+1)/(exp(r)-1) = 2+ r^2/6 * R1(r*r)
+ * That is,
+ * R1(r**2) = 6/r *((exp(r)+1)/(exp(r)-1) - 2/r)
+ * = 6/r * ( 1 + 2.0*(1/(exp(r)-1) - 1/r))
+ * = 1 - r^2/60 + r^4/2520 - r^6/100800 + ...
+ * We use a special Remes algorithm on [0,0.347] to generate
+ * a polynomial of degree 5 in r*r to approximate R1. The
+ * maximum error of this polynomial approximation is bounded
+ * by 2**-61. In other words,
+ * R1(z) ~ 1.0 + Q1*z + Q2*z**2 + Q3*z**3 + Q4*z**4 + Q5*z**5
+ * where Q1 = -1.6666666666666567384E-2,
+ * Q2 = 3.9682539681370365873E-4,
+ * Q3 = -9.9206344733435987357E-6,
+ * Q4 = 2.5051361420808517002E-7,
+ * Q5 = -6.2843505682382617102E-9;
+ * (where z=r*r, and the values of Q1 to Q5 are listed below)
+ * with error bounded by
+ * | 5 | -61
+ * | 1.0+Q1*z+...+Q5*z - R1(z) | <= 2
+ * | |
+ *
+ * expm1(r) = exp(r)-1 is then computed by the following
+ * specific way which minimize the accumulation rounding error:
+ * 2 3
+ * r r [ 3 - (R1 + R1*r/2) ]
+ * expm1(r) = r + --- + --- * [--------------------]
+ * 2 2 [ 6 - r*(3 - R1*r/2) ]
+ *
+ * To compensate the error in the argument reduction, we use
+ * expm1(r+c) = expm1(r) + c + expm1(r)*c
+ * ~ expm1(r) + c + r*c
+ * Thus c+r*c will be added in as the correction terms for
+ * expm1(r+c). Now rearrange the term to avoid optimization
+ * screw up:
+ * ( 2 2 )
+ * ({ ( r [ R1 - (3 - R1*r/2) ] ) } r )
+ * expm1(r+c)~r - ({r*(--- * [--------------------]-c)-c} - --- )
+ * ({ ( 2 [ 6 - r*(3 - R1*r/2) ] ) } 2 )
* ( )
- *
- * = r - E
+ *
+ * = r - E
* 3. Scale back to obtain expm1(x):
- * From step 1, we have
- * expm1(x) = either 2^k*[expm1(r)+1] - 1
- * = or 2^k*[expm1(r) + (1-2^-k)]
+ * From step 1, we have
+ * expm1(x) = either 2^k*[expm1(r)+1] - 1
+ * = or 2^k*[expm1(r) + (1-2^-k)]
* 4. Implementation notes:
- * (A). To save one multiplication, we scale the coefficient Qi
- * to Qi*2^i, and replace z by (x^2)/2.
- * (B). To achieve maximum accuracy, we compute expm1(x) by
- * (i) if x < -56*ln2, return -1.0, (raise inexact if x!=inf)
- * (ii) if k=0, return r-E
- * (iii) if k=-1, return 0.5*(r-E)-0.5
- * (iv) if k=1 if r < -0.25, return 2*((r+0.5)- E)
- * else return 1.0+2.0*(r-E);
- * (v) if (k<-2||k>56) return 2^k(1-(E-r)) - 1 (or exp(x)-1)
- * (vi) if k <= 20, return 2^k((1-2^-k)-(E-r)), else
- * (vii) return 2^k(1-((E+2^-k)-r))
+ * (A). To save one multiplication, we scale the coefficient Qi
+ * to Qi*2^i, and replace z by (x^2)/2.
+ * (B). To achieve maximum accuracy, we compute expm1(x) by
+ * (i) if x < -56*ln2, return -1.0, (raise inexact if x!=inf)
+ * (ii) if k=0, return r-E
+ * (iii) if k=-1, return 0.5*(r-E)-0.5
+ * (iv) if k=1 if r < -0.25, return 2*((r+0.5)- E)
+ * else return 1.0+2.0*(r-E);
+ * (v) if (k<-2||k>56) return 2^k(1-(E-r)) - 1 (or exp(x)-1)
+ * (vi) if k <= 20, return 2^k((1-2^-k)-(E-r)), else
+ * (vii) return 2^k(1-((E+2^-k)-r))
*
* Special cases:
- * expm1(INF) is INF, expm1(NaN) is NaN;
- * expm1(-INF) is -1, and
- * for finite argument, only expm1(0)=0 is exact.
+ * expm1(INF) is INF, expm1(NaN) is NaN;
+ * expm1(-INF) is -1, and
+ * for finite argument, only expm1(0)=0 is exact.
*
* Accuracy:
- * according to an error analysis, the error is always less than
- * 1 ulp (unit in the last place).
+ * according to an error analysis, the error is always less than
+ * 1 ulp (unit in the last place).
*
* Misc. info.
- * For IEEE double
- * if x > 7.09782712893383973096e+02 then expm1(x) overflow
+ * For IEEE double
+ * if x > 7.09782712893383973096e+02 then expm1(x) overflow
*
* Constants:
- * The hexadecimal values are the intended ones for the following
- * constants. The decimal values may be used, provided that the
+ * The hexadecimal values are the intended ones for the following
+ * constants. The decimal values may be used, provided that the
* compiler will convert from decimal to binary accurately enough
* to produce the hexadecimal values shown.
*/
@@ -328,14 +328,14 @@ static const double
#else
static double
#endif
-one = 1.0,
-huge = 1.0e+300,
-tiny = 1.0e-300,
-o_threshold = 7.09782712893383973096e+02,/* 0x40862E42, 0xFEFA39EF */
-ln2_hi = 6.93147180369123816490e-01,/* 0x3fe62e42, 0xfee00000 */
-ln2_lo = 1.90821492927058770002e-10,/* 0x3dea39ef, 0x35793c76 */
-invln2 = 1.44269504088896338700e+00,/* 0x3ff71547, 0x652b82fe */
- /* scaled coefficients related to expm1 */
+one = 1.0,
+huge = 1.0e+300,
+tiny = 1.0e-300,
+o_threshold = 7.09782712893383973096e+02,/* 0x40862E42, 0xFEFA39EF */
+ln2_hi = 6.93147180369123816490e-01,/* 0x3fe62e42, 0xfee00000 */
+ln2_lo = 1.90821492927058770002e-10,/* 0x3dea39ef, 0x35793c76 */
+invln2 = 1.44269504088896338700e+00,/* 0x3ff71547, 0x652b82fe */
+ /* scaled coefficients related to expm1 */
Q1 = -3.33333333333331316428e-02, /* BFA11111 111110F4 */
Q2 = 1.58730158725481460165e-03, /* 3F5A01A0 19FE5585 */
Q3 = -7.93650757867487942473e-05, /* BF14CE19 9EAADBB7 */
@@ -343,89 +343,89 @@ Q4 = 4.00821782732936239552e-06, /* 3ED0CFCA 86E65239
*/
Q5 = -2.01099218183624371326e-07; /* BE8AFDB7 6E09C32D */
#ifdef __STDC__
- double expm1(double x)
+ double expm1(double x)
#else
- double expm1(x)
- double x;
+ double expm1(x)
+ double x;
#endif
{
- double y,hi,lo,c,t,e,hxs,hfx,r1;
- int k,xsb;
- unsigned hx;
+ double y,hi,lo,c,t,e,hxs,hfx,r1;
+ int k,xsb;
+ unsigned hx;
- hx = __HI(x); /* high word of x */
- xsb = hx&0x80000000; /* sign bit of x */
- if(xsb==0) y=x; else y= -x; /* y = |x| */
- hx &= 0x7fffffff; /* high word of |x| */
+ hx = __HI(x); /* high word of x */
+ xsb = hx&0x80000000; /* sign bit of x */
+ if(xsb==0) y=x; else y= -x; /* y = |x| */
+ hx &= 0x7fffffff; /* high word of |x| */
/* filter out huge and non-finite argument */
- if(hx >= 0x4043687A) { /* if |x|>=56*ln2 */
- if(hx >= 0x40862E42) { /* if |x|>=709.78... */
- if(hx>=0x7ff00000) {
- if(((hx&0xfffff)|__LO(x))!=0)
- return x+x; /* NaN */
- else return (xsb==0)? x:-1.0;/* exp(+-inf)={inf,-1} */
- }
- if(x > o_threshold) return huge*huge; /* overflow */
- }
- if(xsb!=0) { /* x < -56*ln2, return -1.0 with inexact */
- if(x+tiny<0.0) /* raise inexact */
- return tiny-one; /* return -1 */
- }
- }
+ if(hx >= 0x4043687A) { /* if |x|>=56*ln2 */
+ if(hx >= 0x40862E42) { /* if |x|>=709.78... */
+ if(hx>=0x7ff00000) {
+ if(((hx&0xfffff)|__LO(x))!=0)
+ return x+x; /* NaN */
+ else return (xsb==0)? x:-1.0;/* exp(+-inf)={inf,-1} */
+ }
+ if(x > o_threshold) return huge*huge; /* overflow */
+ }
+ if(xsb!=0) { /* x < -56*ln2, return -1.0 with inexact */
+ if(x+tiny<0.0) /* raise inexact */
+ return tiny-one; /* return -1 */
+ }
+ }
/* argument reduction */
- if(hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */
- if(hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */
- if(xsb==0)
- {hi = x - ln2_hi; lo = ln2_lo; k = 1;}
- else
- {hi = x + ln2_hi; lo = -ln2_lo; k = -1;}
- } else {
- k = invln2*x+((xsb==0)?0.5:-0.5);
- t = k;
- hi = x - t*ln2_hi; /* t*ln2_hi is exact here */
- lo = t*ln2_lo;
- }
- x = hi - lo;
- c = (hi-x)-lo;
- }
- else if(hx < 0x3c900000) { /* when |x|<2**-54, return x */
- t = huge+x; /* return x with inexact flags when x!=0 */
- return x - (t-(huge+x));
- }
- else k = 0;
+ if(hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */
+ if(hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */
+ if(xsb==0)
+ {hi = x - ln2_hi; lo = ln2_lo; k = 1;}
+ else
+ {hi = x + ln2_hi; lo = -ln2_lo; k = -1;}
+ } else {
+ k = invln2*x+((xsb==0)?0.5:-0.5);
+ t = k;
+ hi = x - t*ln2_hi; /* t*ln2_hi is exact here */
+ lo = t*ln2_lo;
+ }
+ x = hi - lo;
+ c = (hi-x)-lo;
+ }
+ else if(hx < 0x3c900000) { /* when |x|<2**-54, return x */
+ t = huge+x; /* return x with inexact flags when x!=0 */
+ return x - (t-(huge+x));
+ }
+ else k = 0;
/* x is now in primary range */
- hfx = 0.5*x;
- hxs = x*hfx;
- r1 = one+hxs*(Q1+hxs*(Q2+hxs*(Q3+hxs*(Q4+hxs*Q5))));
- t = 3.0-r1*hfx;
- e = hxs*((r1-t)/(6.0 - x*t));
- if(k==0) return x - (x*e-hxs); /* c is 0 */
- else {
- e = (x*(e-c)-c);
- e -= hxs;
- if(k== -1) return 0.5*(x-e)-0.5;
- if(k==1)
- if(x < -0.25) return -2.0*(e-(x+0.5));
- else return one+2.0*(x-e);
- if (k <= -2 || k>56) { /* suffice to return exp(x)-1 */
- y = one-(e-x);
- __HI(y) += (k<<20); /* add k to y's exponent */
- return y-one;
- }
- t = one;
- if(k<20) {
- __HI(t) = 0x3ff00000 - (0x200000>>k); /* t=1-2^-k */
- y = t-(e-x);
- __HI(y) += (k<<20); /* add k to y's exponent */
- } else {
- __HI(t) = ((0x3ff-k)<<20); /* 2^-k */
- y = x-(e+t);
- y += one;
- __HI(y) += (k<<20); /* add k to y's exponent */
- }
- }
- return y;
+ hfx = 0.5*x;
+ hxs = x*hfx;
+ r1 = one+hxs*(Q1+hxs*(Q2+hxs*(Q3+hxs*(Q4+hxs*Q5))));
+ t = 3.0-r1*hfx;
+ e = hxs*((r1-t)/(6.0 - x*t));
+ if(k==0) return x - (x*e-hxs); /* c is 0 */
+ else {
+ e = (x*(e-c)-c);
+ e -= hxs;
+ if(k== -1) return 0.5*(x-e)-0.5;
+ if(k==1)
+ if(x < -0.25) return -2.0*(e-(x+0.5));
+ else return one+2.0*(x-e);
+ if (k <= -2 || k>56) { /* suffice to return exp(x)-1 */
+ y = one-(e-x);
+ __HI(y) += (k<<20); /* add k to y's exponent */
+ return y-one;
+ }
+ t = one;
+ if(k<20) {
+ __HI(t) = 0x3ff00000 - (0x200000>>k); /* t=1-2^-k */
+ y = t-(e-x);
+ __HI(y) += (k<<20); /* add k to y's exponent */
+ } else {
+ __HI(t) = ((0x3ff-k)<<20); /* 2^-k */
+ y = x-(e+t);
+ y += one;
+ __HI(y) += (k<<20); /* add k to y's exponent */
+ }
+ }
+ return y;
}
diff --git a/fdlibm_log1p.c.txt b/fdlibm_log1p.c.txt
index 75b3d564..48952dcf 100644
--- a/fdlibm_log1p.c.txt
+++ b/fdlibm_log1p.c.txt
@@ -5,7 +5,7 @@
* Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved.
*
* Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
+ * software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
@@ -14,8 +14,8 @@
but these catch some common cases. */
#if defined(i386) || defined(i486) || \
- defined(intel) || defined(x86) || defined(i86pc) || \
- defined(__alpha) || defined(__osf__)
+ defined(intel) || defined(x86) || defined(i86pc) || \
+ defined(__alpha) || defined(__osf__)
#define __LITTLE_ENDIAN
#endif
@@ -32,9 +32,9 @@
#endif
#ifdef __STDC__
-#define __P(p) p
+#define __P(p) p
#else
-#define __P(p) ()
+#define __P(p) ()
#endif
/*
@@ -43,20 +43,20 @@
extern int signgam;
-#define MAXFLOAT ((float)3.40282346638528860e+38)
+#define MAXFLOAT ((float)3.40282346638528860e+38)
enum fdversion {fdlibm_ieee = -1, fdlibm_svid, fdlibm_xopen, fdlibm_posix};
#define _LIB_VERSION_TYPE enum fdversion
-#define _LIB_VERSION _fdlib_version
+#define _LIB_VERSION _fdlib_version
-/* if global variable _LIB_VERSION is not desirable, one may
- * change the following to be a constant by:
- * #define _LIB_VERSION_TYPE const enum version
+/* if global variable _LIB_VERSION is not desirable, one may
+ * change the following to be a constant by:
+ * #define _LIB_VERSION_TYPE const enum version
* In that case, after one initializes the value _LIB_VERSION (see
* s_lib_version.c) during compile time, it cannot be modified
* in the middle of a program
- */
+ */
extern _LIB_VERSION_TYPE _LIB_VERSION;
#define _IEEE_ fdlibm_ieee
@@ -65,28 +65,28 @@ extern _LIB_VERSION_TYPE _LIB_VERSION;
#define _POSIX_ fdlibm_posix
struct exception {
- int type;
- char *name;
- double arg1;
- double arg2;
- double retval;
+ int type;
+ char *name;
+ double arg1;
+ double arg2;
+ double retval;
};
-#define HUGE MAXFLOAT
+#define HUGE MAXFLOAT
-/*
+/*
* set X_TLOSS = pi*2**52, which is possibly defined in <values.h>
* (one may replace the following line by "#include <values.h>")
*/
-#define X_TLOSS 1.41484755040568800000e+16
+#define X_TLOSS 1.41484755040568800000e+16
-#define DOMAIN 1
-#define SING 2
-#define OVERFLOW 3
-#define UNDERFLOW 4
-#define TLOSS 5
-#define PLOSS 6
+#define DOMAIN 1
+#define SING 2
+#define OVERFLOW 3
+#define UNDERFLOW 4
+#define TLOSS 5
+#define PLOSS 6
/*
* ANSI/POSIX
@@ -173,16 +173,16 @@ extern double log1p __P((double));
#ifdef _REENTRANT
extern double gamma_r __P((double, int *));
extern double lgamma_r __P((double, int *));
-#endif /* _REENTRANT */
+#endif /* _REENTRANT */
/* ieee style elementary functions */
-extern double __ieee754_sqrt __P((double));
-extern double __ieee754_acos __P((double));
-extern double __ieee754_acosh __P((double));
-extern double __ieee754_log __P((double));
-extern double __ieee754_atanh __P((double));
-extern double __ieee754_asin __P((double));
-extern double __ieee754_atan2 __P((double,double));
+extern double __ieee754_sqrt __P((double));
+extern double __ieee754_acos __P((double));
+extern double __ieee754_acosh __P((double));
+extern double __ieee754_log __P((double));
+extern double __ieee754_atanh __P((double));
+extern double __ieee754_asin __P((double));
+extern double __ieee754_atan2 __P((double,double));
extern double __ieee754_exp __P((double));
extern double __ieee754_cosh __P((double));
extern double __ieee754_fmod __P((double,double));
@@ -209,7 +209,7 @@ extern double __ieee754_scalb __P((double,double));
#endif
/* fdlibm kernel function */
-extern double __kernel_standard __P((double,double,int));
+extern double __kernel_standard __P((double,double,int));
extern double __kernel_sin __P((double,double,int));
extern double __kernel_cos __P((double,double));
extern double __kernel_tan __P((double,double,int));
@@ -222,74 +222,74 @@ extern int __kernel_rem_pio2
__P((double*,double*,int,int,int,const int*));
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
+ * software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/* double log1p(double x)
*
- * Method :
- * 1. Argument Reduction: find k and f such that
- * 1+x = 2^k * (1+f),
- * where sqrt(2)/2 < 1+f < sqrt(2) .
+ * Method :
+ * 1. Argument Reduction: find k and f such that
+ * 1+x = 2^k * (1+f),
+ * where sqrt(2)/2 < 1+f < sqrt(2) .
*
* Note. If k=0, then f=x is exact. However, if k!=0, then f
- * may not be representable exactly. In that case, a correction
- * term is need. Let u=1+x rounded. Let c = (1+x)-u, then
- * log(1+x) - log(u) ~ c/u. Thus, we proceed to compute log(u),
- * and add back the correction term c/u.
- * (Note: when x > 2**53, one can simply return log(x))
+ * may not be representable exactly. In that case, a correction
+ * term is need. Let u=1+x rounded. Let c = (1+x)-u, then
+ * log(1+x) - log(u) ~ c/u. Thus, we proceed to compute log(u),
+ * and add back the correction term c/u.
+ * (Note: when x > 2**53, one can simply return log(x))
*
* 2. Approximation of log1p(f).
- * Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
- * = 2s + 2/3 s**3 + 2/5 s**5 + .....,
- * = 2s + s*R
- * We use a special Reme algorithm on [0,0.1716] to generate
- * a polynomial of degree 14 to approximate R The maximum error
- * of this polynomial approximation is bounded by 2**-58.45. In
- * other words,
- * 2 4 6 8 10 12 14
- * R(z) ~ Lp1*s +Lp2*s +Lp3*s +Lp4*s +Lp5*s +Lp6*s +Lp7*s
- * (the values of Lp1 to Lp7 are listed in the program)
- * and
- * | 2 14 | -58.45
- * | Lp1*s +...+Lp7*s - R(z) | <= 2
- * | |
- * Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.
- * In order to guarantee error in log below 1ulp, we compute log
- * by
- * log1p(f) = f - (hfsq - s*(hfsq+R)).
- *
- * 3. Finally, log1p(x) = k*ln2 + log1p(f).
- * = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo)))
- * Here ln2 is split into two floating point number:
- * ln2_hi + ln2_lo,
- * where n*ln2_hi is always exact for |n| < 2000.
+ * Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
+ * = 2s + 2/3 s**3 + 2/5 s**5 + .....,
+ * = 2s + s*R
+ * We use a special Reme algorithm on [0,0.1716] to generate
+ * a polynomial of degree 14 to approximate R The maximum error
+ * of this polynomial approximation is bounded by 2**-58.45. In
+ * other words,
+ * 2 4 6 8 10 12 14
+ * R(z) ~ Lp1*s +Lp2*s +Lp3*s +Lp4*s +Lp5*s +Lp6*s +Lp7*s
+ * (the values of Lp1 to Lp7 are listed in the program)
+ * and
+ * | 2 14 | -58.45
+ * | Lp1*s +...+Lp7*s - R(z) | <= 2
+ * | |
+ * Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.
+ * In order to guarantee error in log below 1ulp, we compute log
+ * by
+ * log1p(f) = f - (hfsq - s*(hfsq+R)).
+ *
+ * 3. Finally, log1p(x) = k*ln2 + log1p(f).
+ * = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo)))
+ * Here ln2 is split into two floating point number:
+ * ln2_hi + ln2_lo,
+ * where n*ln2_hi is always exact for |n| < 2000.
*
* Special cases:
- * log1p(x) is NaN with signal if x < -1 (including -INF) ;
- * log1p(+INF) is +INF; log1p(-1) is -INF with signal;
- * log1p(NaN) is that NaN with no signal.
+ * log1p(x) is NaN with signal if x < -1 (including -INF) ;
+ * log1p(+INF) is +INF; log1p(-1) is -INF with signal;
+ * log1p(NaN) is that NaN with no signal.
*
* Accuracy:
- * according to an error analysis, the error is always less than
- * 1 ulp (unit in the last place).
+ * according to an error analysis, the error is always less than
+ * 1 ulp (unit in the last place).
*
* Constants:
- * The hexadecimal values are the intended ones for the following
- * constants. The decimal values may be used, provided that the
- * compiler will convert from decimal to binary accurately enough
+ * The hexadecimal values are the intended ones for the following
+ * constants. The decimal values may be used, provided that the
+ * compiler will convert from decimal to binary accurately enough
* to produce the hexadecimal values shown.
*
* Note: Assuming log() return accurate answer, the following
- * algorithm can be used to compute log1p(x) to within a few ULP:
- *
- * u = 1+x;
- * if(u==1.0) return x ; else
- * return log(u)*(x/(u-1.0));
+ * algorithm can be used to compute log1p(x) to within a few ULP:
+ *
+ * u = 1+x;
+ * if(u==1.0) return x ; else
+ * return log(u)*(x/(u-1.0));
*
- * See HP-15C Advanced Functions Handbook, p.193.
+ * See HP-15C Advanced Functions Handbook, p.193.
*/
#include "fdlibm.h"
@@ -299,8 +299,8 @@ static const double
#else
static double
#endif
-ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */
-ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */
+ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */
+ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */
two54 = 1.80143985094819840000e+16, /* 43500000 00000000 */
Lp1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */
Lp2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */
@@ -313,69 +313,69 @@ Lp7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
static double zero = 0.0;
#ifdef __STDC__
- double log1p(double x)
+ double log1p(double x)
#else
- double log1p(x)
- double x;
+ double log1p(x)
+ double x;
#endif
{
- double hfsq,f,c,s,z,R,u;
- int k,hx,hu,ax;
-
- hx = __HI(x); /* high word of x */
- ax = hx&0x7fffffff;
-
- k = 1;
- if (hx < 0x3FDA827A) { /* x < 0.41422 */
- if(ax>=0x3ff00000) { /* x <= -1.0 */
- if(x==-1.0) return -two54/zero; /* log1p(-1)=+inf */
- else return (x-x)/(x-x); /* log1p(x<-1)=NaN */
- }
- if(ax<0x3e200000) { /* |x| < 2**-29 */
- if(two54+x>zero /* raise inexact */
- &&ax<0x3c900000) /* |x| < 2**-54 */
- return x;
- else
- return x - x*x*0.5;
- }
- if(hx>0||hx<=((int)0xbfd2bec3)) {
- k=0;f=x;hu=1;} /* -0.2929<x<0.41422 */
- }
- if (hx >= 0x7ff00000) return x+x;
- if(k!=0) {
- if(hx<0x43400000) {
- u = 1.0+x;
- hu = __HI(u); /* high word of u */
- k = (hu>>20)-1023;
- c = (k>0)? 1.0-(u-x):x-(u-1.0);/* correction term */
- c /= u;
- } else {
- u = x;
- hu = __HI(u); /* high word of u */
- k = (hu>>20)-1023;
- c = 0;
- }
- hu &= 0x000fffff;
- if(hu<0x6a09e) {
- __HI(u) = hu|0x3ff00000; /* normalize u */
- } else {
- k += 1;
- __HI(u) = hu|0x3fe00000; /* normalize u/2 */
- hu = (0x00100000-hu)>>2;
- }
- f = u-1.0;
- }
- hfsq=0.5*f*f;
- if(hu==0) { /* |f| < 2**-20 */
- if(f==zero) if(k==0) return zero;
- else {c += k*ln2_lo; return k*ln2_hi+c;}
- R = hfsq*(1.0-0.66666666666666666*f);
- if(k==0) return f-R; else
- return k*ln2_hi-((R-(k*ln2_lo+c))-f);
- }
- s = f/(2.0+f);
- z = s*s;
- R = z*(Lp1+z*(Lp2+z*(Lp3+z*(Lp4+z*(Lp5+z*(Lp6+z*Lp7))))));
- if(k==0) return f-(hfsq-s*(hfsq+R)); else
- return k*ln2_hi-((hfsq-(s*(hfsq+R)+(k*ln2_lo+c)))-f);
+ double hfsq,f,c,s,z,R,u;
+ int k,hx,hu,ax;
+
+ hx = __HI(x); /* high word of x */
+ ax = hx&0x7fffffff;
+
+ k = 1;
+ if (hx < 0x3FDA827A) { /* x < 0.41422 */
+ if(ax>=0x3ff00000) { /* x <= -1.0 */
+ if(x==-1.0) return -two54/zero; /* log1p(-1)=+inf */
+ else return (x-x)/(x-x); /* log1p(x<-1)=NaN */
+ }
+ if(ax<0x3e200000) { /* |x| < 2**-29 */
+ if(two54+x>zero /* raise inexact */
+ &&ax<0x3c900000) /* |x| < 2**-54 */
+ return x;
+ else
+ return x - x*x*0.5;
+ }
+ if(hx>0||hx<=((int)0xbfd2bec3)) {
+ k=0;f=x;hu=1;} /* -0.2929<x<0.41422 */
+ }
+ if (hx >= 0x7ff00000) return x+x;
+ if(k!=0) {
+ if(hx<0x43400000) {
+ u = 1.0+x;
+ hu = __HI(u); /* high word of u */
+ k = (hu>>20)-1023;
+ c = (k>0)? 1.0-(u-x):x-(u-1.0);/* correction term */
+ c /= u;
+ } else {
+ u = x;
+ hu = __HI(u); /* high word of u */
+ k = (hu>>20)-1023;
+ c = 0;
+ }
+ hu &= 0x000fffff;
+ if(hu<0x6a09e) {
+ __HI(u) = hu|0x3ff00000; /* normalize u */
+ } else {
+ k += 1;
+ __HI(u) = hu|0x3fe00000; /* normalize u/2 */
+ hu = (0x00100000-hu)>>2;
+ }
+ f = u-1.0;
+ }
+ hfsq=0.5*f*f;
+ if(hu==0) { /* |f| < 2**-20 */
+ if(f==zero) if(k==0) return zero;
+ else {c += k*ln2_lo; return k*ln2_hi+c;}
+ R = hfsq*(1.0-0.66666666666666666*f);
+ if(k==0) return f-R; else
+ return k*ln2_hi-((R-(k*ln2_lo+c))-f);
+ }
+ s = f/(2.0+f);
+ z = s*s;
+ R = z*(Lp1+z*(Lp2+z*(Lp3+z*(Lp4+z*(Lp5+z*(Lp6+z*Lp7))))));
+ if(k==0) return f-(hfsq-s*(hfsq+R)); else
+ return k*ln2_hi-((hfsq-(s*(hfsq+R)+(k*ln2_lo+c)))-f);
}
- [lmi-commits] [lmi] master c4e7bd93 01/22: Import expm1() and log1p() from fdlibm, (continued)
- [lmi-commits] [lmi] master c4e7bd93 01/22: Import expm1() and log1p() from fdlibm, Greg Chicares, 2022/05/20
- [lmi-commits] [lmi] master e9f2f432 04/22: Canonicalize #include directives, Greg Chicares, 2022/05/20
- [lmi-commits] [lmi] master 5cc39664 03/22: Move a closing brace, Greg Chicares, 2022/05/20
- [lmi-commits] [lmi] master 1fd91fbf 10/22: Rename high- and low-word macros, Greg Chicares, 2022/05/20
- [lmi-commits] [lmi] master fa7bbb53 05/22: Prefer lmi style, Greg Chicares, 2022/05/20
- [lmi-commits] [lmi] master 6b073770 08/22: Expunge unreferenced fdlibm versioning, Greg Chicares, 2022/05/20
- [lmi-commits] [lmi] master cadb36d4 14/22: Don't include any fdlibm header, Greg Chicares, 2022/05/20
- [lmi-commits] [lmi] master 5b1ef5bb 15/22: Add prototypes, Greg Chicares, 2022/05/20
- [lmi-commits] [lmi] master dd3abb11 18/22: Demonstrate a defect, Greg Chicares, 2022/05/20
- [lmi-commits] [lmi] master af8d3ac0 21/22: Rationalize indentation of preprocessor directives, Greg Chicares, 2022/05/20
- [lmi-commits] [lmi] master 324e171e 02/22: Rectify whitespace,
Greg Chicares <=
- [lmi-commits] [lmi] master 15b1aa35 11/22: Add license and copyright boilerplate, Greg Chicares, 2022/05/20
- [lmi-commits] [lmi] master 4ec9da1f 12/22: Rename fdlibm sources, Greg Chicares, 2022/05/20
- [lmi-commits] [lmi] master e8a859a3 06/22: Rename prototyping macro, Greg Chicares, 2022/05/20
- [lmi-commits] [lmi] master 078e9a1c 16/22: Ignore certain diagnostics, Greg Chicares, 2022/05/20
- [lmi-commits] [lmi] master 2d514bdc 09/22: Allow some common macros, Greg Chicares, 2022/05/20
- [lmi-commits] [lmi] master 1646f4c0 13/22: Rename fdlibm entry points, Greg Chicares, 2022/05/20
- [lmi-commits] [lmi] master 75dc2686 17/22: Write unambiguous braces, Greg Chicares, 2022/05/20
- [lmi-commits] [lmi] master 9e608afe 20/22: Remove debugging output for tests that now succeed, Greg Chicares, 2022/05/20
- [lmi-commits] [lmi] master 2563b0f6 07/22: Expunge prototypes for unreferenced functions, Greg Chicares, 2022/05/20
- [lmi-commits] [lmi] master c6902b0d 19/22: Fix the defect just noted, Greg Chicares, 2022/05/20