-
Notifications
You must be signed in to change notification settings - Fork 14
/
expm1-standalone.c
119 lines (79 loc) · 3.48 KB
/
expm1-standalone.c
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
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
#include "crlibm.h"
#include "crlibm_private.h"
#include "triple-double.h"
#include "expm1.h"
#define AVOID_FMA 0
extern void expm1_direct_td(double *expm1h, double *expm1m, double *expm1l,
double x, double xSqHalfh, double xSqHalfl, double xSqh, double xSql, int expoX);
extern void expm1_common_td(double *expm1h, double *expm1m, double *expm1l,
double rh, double rm, double rl,
double tbl1h, double tbl1m, double tbl1l,
double tbl2h, double tbl2m, double tbl2l,
int M);
/* Function expm1_13
Computes exp(x)-1 with an accuracy of 120 bits (128 for |x| <= 2^(-5)) as
(expm1h + expm1m + expm1l) \approx exp(x) - 1
There is no special case handling.
|x| is supposed to be greater than 2^(-53)
x is supposed to be in the range -37...709
The triple-double exph + expm + expl is non-overlapping.
*/
void expm1_13(double *expm1h, double *expm1m, double *expm1l, double x) {
db_number xdb, shiftedXMultdb;
int xIntHi, expoX, k, M, index1, index2;
double xSqh, xSql, xSqHalfh, xSqHalfl;
double t1, t2;
double msLog2Div2LMultKh, msLog2Div2LMultKm, msLog2Div2LMultKl;
double xMultLog2InvMult2L, shiftedXMult, kd, rh, rm, rl;
double tbl1h, tbl1m, tbl1l, tbl2h, tbl2m, tbl2l;
xdb.d = x;
xIntHi = xdb.i[HI] & 0x7fffffff;
/* Test if we have |x| <= 1/4-1/2ulp(1/4) for knowing if we use exp(x) or approximate directly */
if (xIntHi < DIRECTINTERVALBOUND) {
/* We approximate expm1 directly after a range reduction as follows
expm1(x) = (expm1(x/2) + 2) * expm1(x/2)
We perform the range reduction in such a way that finally |x| < 1/32
*/
/* Extract the exponent of |x| and add 5 (2^5 = 32) */
expoX = ((xIntHi & 0x7ff00000) >> 20) - (1023 - 5);
/* If this particularily biased exponent expoX is negative, we are already less than 1/32 */
if (expoX >= 0) {
/* If we are here, we must perform range reduction */
/* We multiply x by 2^(-expoX-1) by bit manipulation
x cannot be denormalized so there is no danger
*/
xdb.i[HI] += (-expoX-1) << 20;
/* We reassign the new x and maintain xIntHi */
xIntHi = xdb.i[HI] & 0x7fffffff;
x = xdb.d;
}
/* Here, we have always |x| < 1/32 */
/* Double precision evaluation steps and one double-double step */
Mul12(&xSqh,&xSql,x,x);
xSqHalfh = 0.5 * xSqh;
xSqHalfl = 0.5 * xSql;
expm1_direct_td(expm1h, expm1m, expm1l, x, xSqHalfh, xSqHalfl, xSqh, xSql, expoX);
} else {
/* If we are here, we can use expm1(x) = exp(x) - 1 */
/* Range reduction - exact part: compute k as double and as int */
xMultLog2InvMult2L = x * log2InvMult2L;
shiftedXMult = xMultLog2InvMult2L + shiftConst;
kd = shiftedXMult - shiftConst;
shiftedXMultdb.d = shiftedXMult;
k = shiftedXMultdb.i[LO];
M = k >> 12;
index1 = k & INDEXMASK1;
index2 = (k & INDEXMASK2) >> 6;
Mul133(&msLog2Div2LMultKh,&msLog2Div2LMultKm,&msLog2Div2LMultKl,kd,msLog2Div2Lh,msLog2Div2Lm,msLog2Div2Ll);
t1 = x + msLog2Div2LMultKh;
Add12Cond(rh,t2,t1,msLog2Div2LMultKm);
Add12Cond(rm,rl,t2,msLog2Div2LMultKl);
tbl1h = twoPowerIndex1[index1].hi;
tbl1m = twoPowerIndex1[index1].mi;
tbl2h = twoPowerIndex2[index2].hi;
tbl2m = twoPowerIndex2[index2].mi;
tbl1l = twoPowerIndex1[index1].lo;
tbl2l = twoPowerIndex2[index2].lo;
expm1_common_td(expm1h, expm1m, expm1l, rh, rm, rl, tbl1h, tbl1m, tbl1l, tbl2h, tbl2m, tbl2l, M);
}
}