Skip to content
GitLab
Explore
Sign in
Register
Primary navigation
Search or go to…
Project
cmucl
Manage
Activity
Members
Labels
Plan
Issues
Issue boards
Milestones
Wiki
Code
Merge requests
Repository
Branches
Commits
Tags
Repository graph
Compare revisions
Snippets
Build
Pipelines
Jobs
Pipeline schedules
Artifacts
Deploy
Releases
Package Registry
Model registry
Operate
Environments
Terraform modules
Monitor
Incidents
Service Desk
Analyze
Value stream analytics
Contributor analytics
CI/CD analytics
Repository analytics
Model experiments
Help
Help
Support
GitLab documentation
Compare GitLab plans
Community forum
Contribute to GitLab
Provide feedback
Terms and privacy
Keyboard shortcuts
?
Snippets
Groups
Projects
Show more breadcrumbs
Carl Shapiro
cmucl
Commits
92c7c5d0
Commit
92c7c5d0
authored
10 years ago
by
Raymond Toy
Browse files
Options
Downloads
Patches
Plain Diff
Add fdlibm routines e_exp and e_log, as is.
parent
99afcf7a
No related branches found
No related tags found
No related merge requests found
Changes
2
Hide whitespace changes
Inline
Side-by-side
Showing
2 changed files
src/lisp/e_exp.c
+156
-0
156 additions, 0 deletions
src/lisp/e_exp.c
src/lisp/e_log.c
+139
-0
139 additions, 0 deletions
src/lisp/e_log.c
with
295 additions
and
0 deletions
src/lisp/e_exp.c
0 → 100644
+
156
−
0
View file @
92c7c5d0
/* @(#)e_exp.c 1.6 04/04/22 */
/*
* ====================================================
* Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved.
*
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/* __ieee754_exp(x)
* Returns the exponential of x.
*
* Method
* 1. Argument reduction:
* Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658.
* Given x, find r and integer k such that
*
* x = k*ln2 + r, |r| <= 0.5*ln2.
*
* Here r will be represented as r = hi-lo for better
* accuracy.
*
* 2. Approximation of exp(r) by a special rational function on
* the interval [0,0.34658]:
* Write
* R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ...
* We use a special Remes algorithm on [0,0.34658] to generate
* a polynomial of degree 5 to approximate R. The maximum error
* of this polynomial approximation is bounded by 2**-59. In
* other words,
* R(z) ~ 2.0 + P1*z + P2*z**2 + P3*z**3 + P4*z**4 + P5*z**5
* (where z=r*r, and the values of P1 to P5 are listed below)
* and
* | 5 | -59
* | 2.0+P1*z+...+P5*z - R(z) | <= 2
* | |
* The computation of exp(r) thus becomes
* 2*r
* exp(r) = 1 + -------
* R - r
* r*R1(r)
* = 1 + r + ----------- (for better accuracy)
* 2 - R1(r)
* where
* 2 4 10
* R1(r) = r - (P1*r + P2*r + ... + P5*r ).
*
* 3. Scale back to obtain exp(x):
* From step 1, we have
* exp(x) = 2^k * exp(r)
*
* Special cases:
* exp(INF) is INF, exp(NaN) is NaN;
* exp(-INF) is 0, and
* for finite argument, only exp(0)=1 is exact.
*
* Accuracy:
* 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 exp(x) overflow
* if x < -7.45133219101941108420e+02 then exp(x) underflow
*
* 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
* to produce the hexadecimal values shown.
*/
#include
"fdlibm.h"
#ifdef __STDC__
static
const
double
#else
static
double
#endif
one
=
1
.
0
,
halF
[
2
]
=
{
0
.
5
,
-
0
.
5
,},
huge
=
1.0e+300
,
twom1000
=
9.33263618503218878990e-302
,
/* 2**-1000=0x01700000,0*/
o_threshold
=
7.09782712893383973096e+02
,
/* 0x40862E42, 0xFEFA39EF */
u_threshold
=
-
7.45133219101941108420e+02
,
/* 0xc0874910, 0xD52D3051 */
ln2HI
[
2
]
=
{
6.93147180369123816490e-01
,
/* 0x3fe62e42, 0xfee00000 */
-
6.93147180369123816490e-01
,},
/* 0xbfe62e42, 0xfee00000 */
ln2LO
[
2
]
=
{
1.90821492927058770002e-10
,
/* 0x3dea39ef, 0x35793c76 */
-
1.90821492927058770002e-10
,},
/* 0xbdea39ef, 0x35793c76 */
invln2
=
1.44269504088896338700e+00
,
/* 0x3ff71547, 0x652b82fe */
P1
=
1.66666666666666019037e-01
,
/* 0x3FC55555, 0x5555553E */
P2
=
-
2.77777777770155933842e-03
,
/* 0xBF66C16C, 0x16BEBD93 */
P3
=
6.61375632143793436117e-05
,
/* 0x3F11566A, 0xAF25DE2C */
P4
=
-
1.65339022054652515390e-06
,
/* 0xBEBBBD41, 0xC5D26BF1 */
P5
=
4.13813679705723846039e-08
;
/* 0x3E663769, 0x72BEA4D0 */
#ifdef __STDC__
double
__ieee754_exp
(
double
x
)
/* default IEEE double exp */
#else
double
__ieee754_exp
(
x
)
/* default IEEE double exp */
double
x
;
#endif
{
double
y
,
hi
,
lo
,
c
,
t
;
int
k
,
xsb
;
unsigned
hx
;
hx
=
__HI
(
x
);
/* high word of x */
xsb
=
(
hx
>>
31
)
&
1
;
/* sign bit of x */
hx
&=
0x7fffffff
;
/* high word of |x| */
/* filter out non-finite argument */
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
:
0
.
0
;
/* exp(+-inf)={inf,0} */
}
if
(
x
>
o_threshold
)
return
huge
*
huge
;
/* overflow */
if
(
x
<
u_threshold
)
return
twom1000
*
twom1000
;
/* underflow */
}
/* argument reduction */
if
(
hx
>
0x3fd62e42
)
{
/* if |x| > 0.5 ln2 */
if
(
hx
<
0x3FF0A2B2
)
{
/* and |x| < 1.5 ln2 */
hi
=
x
-
ln2HI
[
xsb
];
lo
=
ln2LO
[
xsb
];
k
=
1
-
xsb
-
xsb
;
}
else
{
k
=
(
int
)(
invln2
*
x
+
halF
[
xsb
]);
t
=
k
;
hi
=
x
-
t
*
ln2HI
[
0
];
/* t*ln2HI is exact here */
lo
=
t
*
ln2LO
[
0
];
}
x
=
hi
-
lo
;
}
else
if
(
hx
<
0x3e300000
)
{
/* when |x|<2**-28 */
if
(
huge
+
x
>
one
)
return
one
+
x
;
/* trigger inexact */
}
else
k
=
0
;
/* x is now in primary range */
t
=
x
*
x
;
c
=
x
-
t
*
(
P1
+
t
*
(
P2
+
t
*
(
P3
+
t
*
(
P4
+
t
*
P5
))));
if
(
k
==
0
)
return
one
-
((
x
*
c
)
/
(
c
-
2
.
0
)
-
x
);
else
y
=
one
-
((
lo
-
(
x
*
c
)
/
(
2
.
0
-
c
))
-
hi
);
if
(
k
>=
-
1021
)
{
__HI
(
y
)
+=
(
k
<<
20
);
/* add k to y's exponent */
return
y
;
}
else
{
__HI
(
y
)
+=
((
k
+
1000
)
<<
20
);
/* add k to y's exponent */
return
y
*
twom1000
;
}
}
This diff is collapsed.
Click to expand it.
src/lisp/e_log.c
0 → 100644
+
139
−
0
View file @
92c7c5d0
/* @(#)e_log.c 1.3 95/01/18 */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/* __ieee754_log(x)
* Return the logrithm of x
*
* Method :
* 1. Argument Reduction: find k and f such that
* x = 2^k * (1+f),
* where sqrt(2)/2 < 1+f < sqrt(2) .
*
* 2. Approximation of log(1+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) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s +Lg6*s +Lg7*s
* (the values of Lg1 to Lg7 are listed in the program)
* and
* | 2 14 | -58.45
* | Lg1*s +...+Lg7*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
* log(1+f) = f - s*(f - R) (if f is not too large)
* log(1+f) = f - (hfsq - s*(hfsq+R)). (better accuracy)
*
* 3. Finally, log(x) = k*ln2 + log(1+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:
* log(x) is NaN with signal if x < 0 (including -INF) ;
* log(+INF) is +INF; log(0) is -INF with signal;
* log(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).
*
* 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
* to produce the hexadecimal values shown.
*/
#include
"fdlibm.h"
#ifdef __STDC__
static
const
double
#else
static
double
#endif
ln2_hi
=
6.93147180369123816490e-01
,
/* 3fe62e42 fee00000 */
ln2_lo
=
1.90821492927058770002e-10
,
/* 3dea39ef 35793c76 */
two54
=
1.80143985094819840000e+16
,
/* 43500000 00000000 */
Lg1
=
6.666666666666735130e-01
,
/* 3FE55555 55555593 */
Lg2
=
3.999999999940941908e-01
,
/* 3FD99999 9997FA04 */
Lg3
=
2.857142874366239149e-01
,
/* 3FD24924 94229359 */
Lg4
=
2.222219843214978396e-01
,
/* 3FCC71C5 1D8E78AF */
Lg5
=
1.818357216161805012e-01
,
/* 3FC74664 96CB03DE */
Lg6
=
1.531383769920937332e-01
,
/* 3FC39A09 D078C69F */
Lg7
=
1.479819860511658591e-01
;
/* 3FC2F112 DF3E5244 */
static
double
zero
=
0
.
0
;
#ifdef __STDC__
double
__ieee754_log
(
double
x
)
#else
double
__ieee754_log
(
x
)
double
x
;
#endif
{
double
hfsq
,
f
,
s
,
z
,
R
,
w
,
t1
,
t2
,
dk
;
int
k
,
hx
,
i
,
j
;
unsigned
lx
;
hx
=
__HI
(
x
);
/* high word of x */
lx
=
__LO
(
x
);
/* low word of x */
k
=
0
;
if
(
hx
<
0x00100000
)
{
/* x < 2**-1022 */
if
(((
hx
&
0x7fffffff
)
|
lx
)
==
0
)
return
-
two54
/
zero
;
/* log(+-0)=-inf */
if
(
hx
<
0
)
return
(
x
-
x
)
/
zero
;
/* log(-#) = NaN */
k
-=
54
;
x
*=
two54
;
/* subnormal number, scale up x */
hx
=
__HI
(
x
);
/* high word of x */
}
if
(
hx
>=
0x7ff00000
)
return
x
+
x
;
k
+=
(
hx
>>
20
)
-
1023
;
hx
&=
0x000fffff
;
i
=
(
hx
+
0x95f64
)
&
0x100000
;
__HI
(
x
)
=
hx
|
(
i
^
0x3ff00000
);
/* normalize x or x/2 */
k
+=
(
i
>>
20
);
f
=
x
-
1
.
0
;
if
((
0x000fffff
&
(
2
+
hx
))
<
3
)
{
/* |f| < 2**-20 */
if
(
f
==
zero
)
if
(
k
==
0
)
return
zero
;
else
{
dk
=
(
double
)
k
;
return
dk
*
ln2_hi
+
dk
*
ln2_lo
;}
R
=
f
*
f
*
(
0
.
5
-
0
.
33333333333333333
*
f
);
if
(
k
==
0
)
return
f
-
R
;
else
{
dk
=
(
double
)
k
;
return
dk
*
ln2_hi
-
((
R
-
dk
*
ln2_lo
)
-
f
);}
}
s
=
f
/
(
2
.
0
+
f
);
dk
=
(
double
)
k
;
z
=
s
*
s
;
i
=
hx
-
0x6147a
;
w
=
z
*
z
;
j
=
0x6b851
-
hx
;
t1
=
w
*
(
Lg2
+
w
*
(
Lg4
+
w
*
Lg6
));
t2
=
z
*
(
Lg1
+
w
*
(
Lg3
+
w
*
(
Lg5
+
w
*
Lg7
)));
i
|=
j
;
R
=
t2
+
t1
;
if
(
i
>
0
)
{
hfsq
=
0
.
5
*
f
*
f
;
if
(
k
==
0
)
return
f
-
(
hfsq
-
s
*
(
hfsq
+
R
));
else
return
dk
*
ln2_hi
-
((
hfsq
-
(
s
*
(
hfsq
+
R
)
+
dk
*
ln2_lo
))
-
f
);
}
else
{
if
(
k
==
0
)
return
f
-
s
*
(
f
-
R
);
else
return
dk
*
ln2_hi
-
((
s
*
(
f
-
R
)
-
dk
*
ln2_lo
)
-
f
);
}
}
This diff is collapsed.
Click to expand it.
Preview
0%
Loading
Try again
or
attach a new file
.
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Save comment
Cancel
Please
register
or
sign in
to comment