1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
|
diff --git a/e_exp.cpp b/e_exp.cpp
--- a/e_exp.cpp
+++ b/e_exp.cpp
@@ -146,14 +146,17 @@ __ieee754_exp(double x) /* default IEEE double exp */
if(k >= -1021)
INSERT_WORDS(twopk,((u_int32_t)(0x3ff+k))<<20, 0);
else
INSERT_WORDS(twopk,((u_int32_t)(0x3ff+(k+1000)))<<20, 0);
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) {
- if (k==1024) return y*2.0*0x1p1023;
+ if (k==1024) {
+ double const_0x1p1023 = pow(2, 1023);
+ return y*2.0*const_0x1p1023;
+ }
return y*twopk;
} else {
return y*twopk*twom1000;
}
}
diff --git a/s_expm1.cpp b/s_expm1.cpp
--- a/s_expm1.cpp
+++ b/s_expm1.cpp
@@ -192,17 +192,20 @@ expm1(double x)
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);
- if (k == 1024) y = y*2.0*0x1p1023;
+ if (k == 1024) {
+ double const_0x1p1023 = pow(2, 1023);
+ y = y*2.0*const_0x1p1023;
+ }
else y = y*twopk;
return y-one;
}
t = one;
if(k<20) {
SET_HIGH_WORD(t,0x3ff00000 - (0x200000>>k)); /* t=1-2^-k */
y = t-(e-x);
y = y*twopk;
--- a/s_exp2.cpp 2022-12-11 21:26:36.045643147 -0500
+++ b/s_exp2.cpp_new 2022-12-11 21:28:42.058769925 -0500
@@ -381,14 +381,16 @@
INSERT_WORDS(twopk, 0x3ff00000 + k, 0);
else
INSERT_WORDS(twopkp1000, 0x3ff00000 + k + (1000 << 20), 0);
r = t + t * z * (P1 + z * (P2 + z * (P3 + z * (P4 + z * P5))));
/* Scale by 2**(k>>20). */
if(k >= -(1021 << 20)) {
- if (k == 1024 << 20)
- return (r * 2.0 * 0x1p1023);
+ if (k == 1024 << 20) {
+ double const_0x1p1023 = pow(2, 1023);
+ return (r * 2.0 * const_0x1p1023);
+ }
return (r * twopk);
} else {
return (r * twopkp1000 * twom1000);
}
}
|
