Skip to content

Commit

Permalink
* lib/msun/Makefile:
Browse files Browse the repository at this point in the history 
  . Disconnect b_exp.c and b_log.c from the build.

* lib/msun/bsdsrc/b_exp.c:
  . Replace scalb() usage with C99's ldexp().
  . Replace finite(x) usage with C99's isfinite().
  . Whitespace changes towards style(9).
  . Remove include of "mathimpl.h".  It is no longer needed.
  . Remove #if 0 ... #endif code, which has been present since svn r93211
    (2002-03-26).
  . New minimax polynomial coefficients.
  . Add comments to explain origins of some constants.
  . Use ansi-C prototype.  Remove K&R prototype.  Add static to prototype.

* lib/msun/bsdsrc/b_log.c:
  . Remove include of "mathimpl.h".  It is no longer needed.
  . Fix comments to actually describe the code.
  . Reduce minimax polynomial from degree 4 to degree 3.
    This uses newly computed coefficients.
  . Use ansi-C prototype.  Remove K&R prototype.  Add static to prototype.
  . Remove volatile in declaration of u1.
  . Alphabetize decalaration list.
  . Whitespace changes towards style(9).
  . In argument reduction of x to g and m, replace use of logb() and
    ldexp() with a single call to frexp().  Add code to get 1 <= g < 2.
  . Remove #if 0 ... #endif code, which has been present since svn r93211
    (2002-03-26).
  . The special case m == -1022, replace logb() with ilogb().

* lib/msun/bsdsrc/b_tgamma.c:
  . Update comments.  Fix comments where needed.
  . Add float.h to get LDBL_MANT_DIG for weak reference of tgammal to tgamma.
  . Remove include of "mathimpl.h".  It is no longer needed.
  . Use "math.h" instead of <math.h>.
  . Add '#include math_private.h"
  . Add struct Double from mathimpl.h and include b_log.c and b_exp.c.
  . Remove forward declarations of neg_gam(), small_gam(), smaller_gam,
    large_gam() and ratfun_gam() by re-arranging the code to move these
    function above their first reference.
  . New minimax coefficients for polynomial in large_gam().
  . New splitting of a0 into a0hi nd a0lo, which include additional
    bits of precision.
  . Use ansi-C prototype.  Remove K&R prototype.
  . Replace the TRUNC() macro with a simple cast of a double entities
    to float before assignment (functional changes).
  . Replace sin(M_PI*z) with sinpi(z) and cos(M_PI*(0.5-z)) with cospi(0.5-z).

Submitted by:		Steve Kargl
Differential Revision:	https://reviews.freebsd.org/D33444
Reviewed by:		pfg

(cherry picked from commit 455b2cc)
  • Loading branch information
@emaste
Mark Murray authored and emaste committed Jul 26, 2023
1 parent 2b01ffb commit 4b597e5
Show file tree
Hide file tree
Showing 4 changed files with 370 additions and 396 deletions.
2 changes: 1 addition & 1 deletion lib/msun/Makefile
View file
Expand Up @@ -59,7 +59,7 @@ SHLIBDIR?= /lib
SHLIB_MAJOR= 5
WARNS?= 1
IGNORE_PRAGMA=
COMMON_SRCS= b_exp.c b_log.c b_tgamma.c \
COMMON_SRCS= b_tgamma.c \
e_acos.c e_acosf.c e_acosh.c e_acoshf.c e_asin.c e_asinf.c \
e_atan2.c e_atan2f.c e_atanh.c e_atanhf.c e_cosh.c e_coshf.c e_exp.c \
e_expf.c e_fmod.c e_fmodf.c e_gamma.c e_gamma_r.c e_gammaf.c \
Expand Down
143 changes: 48 additions & 95 deletions lib/msun/bsdsrc/b_exp.c
View file
Expand Up @@ -33,22 +33,21 @@
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");


/* EXP(X)
* RETURN THE EXPONENTIAL OF X
* DOUBLE PRECISION (IEEE 53 bits, VAX D FORMAT 56 BITS)
* CODED IN C BY K.C. NG, 1/19/85;
* REVISED BY K.C. NG on 2/6/85, 2/15/85, 3/7/85, 3/24/85, 4/16/85, 6/14/86.
*
* Required system supported functions:
* scalb(x,n)
* ldexp(x,n)
* copysign(x,y)
* finite(x)
* isfinite(x)
*
* Method:
* 1. Argument Reduction: given the input x, find r and integer k such
* that
* x = k*ln2 + r, |r| <= 0.5*ln2 .
* x = k*ln2 + r, |r| <= 0.5*ln2.
* r will be represented as r := z+c for better accuracy.
*
* 2. Compute exp(r) by
Expand All @@ -69,105 +68,59 @@ __FBSDID("$FreeBSD$");
* with 1,156,000 random arguments on a VAX, the maximum observed
* error was 0.869 ulps (units in the last place).
*/

#include "mathimpl.h"

static const double p1 = 0x1.555555555553ep-3;
static const double p2 = -0x1.6c16c16bebd93p-9;
static const double p3 = 0x1.1566aaf25de2cp-14;
static const double p4 = -0x1.bbd41c5d26bf1p-20;
static const double p5 = 0x1.6376972bea4d0p-25;
static const double ln2hi = 0x1.62e42fee00000p-1;
static const double ln2lo = 0x1.a39ef35793c76p-33;
static const double lnhuge = 0x1.6602b15b7ecf2p9;
static const double lntiny = -0x1.77af8ebeae354p9;
static const double invln2 = 0x1.71547652b82fep0;

#if 0
double exp(x)
double x;
{
double z,hi,lo,c;
int k;

#if !defined(vax)&&!defined(tahoe)
if(x!=x) return(x); /* x is NaN */
#endif /* !defined(vax)&&!defined(tahoe) */
if( x <= lnhuge ) {
if( x >= lntiny ) {

/* argument reduction : x --> x - k*ln2 */

k=invln2*x+copysign(0.5,x); /* k=NINT(x/ln2) */

/* express x-k*ln2 as hi-lo and let x=hi-lo rounded */

hi=x-k*ln2hi;
x=hi-(lo=k*ln2lo);

/* return 2^k*[1+x+x*c/(2+c)] */
z=x*x;
c= x - z*(p1+z*(p2+z*(p3+z*(p4+z*p5))));
return scalb(1.0+(hi-(lo-(x*c)/(2.0-c))),k);

}
/* end of x > lntiny */

else
/* exp(-big#) underflows to zero */
if(finite(x)) return(scalb(1.0,-5000));

/* exp(-INF) is zero */
else return(0.0);
}
/* end of x < lnhuge */

else
/* exp(INF) is INF, exp(+big#) overflows to INF */
return( finite(x) ? scalb(1.0,5000) : x);
}
#endif
static const double
p1 = 1.6666666666666660e-01, /* 0x3fc55555, 0x55555553 */
p2 = -2.7777777777564776e-03, /* 0xbf66c16c, 0x16c0ac3c */
p3 = 6.6137564717940088e-05, /* 0x3f11566a, 0xb5c2ba0d */
p4 = -1.6534060280704225e-06, /* 0xbebbbd53, 0x273e8fb7 */
p5 = 4.1437773411069054e-08; /* 0x3e663f2a, 0x09c94b6c */

static const double
ln2hi = 0x1.62e42fee00000p-1, /* High 32 bits round-down. */
ln2lo = 0x1.a39ef35793c76p-33; /* Next 53 bits round-to-nearst. */

static const double
lnhuge = 0x1.6602b15b7ecf2p9, /* (DBL_MAX_EXP + 9) * log(2.) */
lntiny = -0x1.77af8ebeae354p9, /* (DBL_MIN_EXP - 53 - 10) * log(2.) */
invln2 = 0x1.71547652b82fep0; /* 1 / log(2.) */

/* returns exp(r = x + c) for |c| < |x| with no overlap. */

double __exp__D(x, c)
double x, c;
static double
__exp__D(double x, double c)
{
double z,hi,lo;
double hi, lo, z;
int k;

if (x != x) /* x is NaN */
if (x != x) /* x is NaN. */
return(x);
if ( x <= lnhuge ) {
if ( x >= lntiny ) {

/* argument reduction : x --> x - k*ln2 */
z = invln2*x;
k = z + copysign(.5, x);

/* express (x+c)-k*ln2 as hi-lo and let x=hi-lo rounded */

hi=(x-k*ln2hi); /* Exact. */
x= hi - (lo = k*ln2lo-c);
/* return 2^k*[1+x+x*c/(2+c)] */
z=x*x;
c= x - z*(p1+z*(p2+z*(p3+z*(p4+z*p5))));
c = (x*c)/(2.0-c);

return scalb(1.+(hi-(lo - c)), k);
if (x <= lnhuge) {
if (x >= lntiny) {
/* argument reduction: x --> x - k*ln2 */
z = invln2 * x;
k = z + copysign(0.5, x);

/*
* Express (x + c) - k * ln2 as hi - lo.
* Let x = hi - lo rounded.
*/
hi = x - k * ln2hi; /* Exact. */
lo = k * ln2lo - c;
x = hi - lo;

/* Return 2^k*[1+x+x*c/(2+c)] */
z = x * x;
c = x - z * (p1 + z * (p2 + z * (p3 + z * (p4 +
z * p5))));
c = (x * c) / (2 - c);

return (ldexp(1 + (hi - (lo - c)), k));
} else {
/* exp(-INF) is 0. exp(-big) underflows to 0. */
return (isfinite(x) ? ldexp(1., -5000) : 0);
}
/* end of x > lntiny */

else
/* exp(-big#) underflows to zero */
if(finite(x)) return(scalb(1.0,-5000));

/* exp(-INF) is zero */
else return(0.0);
}
/* end of x < lnhuge */

else
} else
/* exp(INF) is INF, exp(+big#) overflows to INF */
return( finite(x) ? scalb(1.0,5000) : x);
return (isfinite(x) ? ldexp(1., 5000) : x);
}

0 comments on commit 4b597e5

Please sign in to comment.