1 #include "crlibm.h"
2 #include "crlibm_private.h"
3 #include "triple-double.h"
4 #include "exp-td.h"
5
6 #define AVOID_FMA 0
7
8
9 extern void exp_td_accurate(double *polyTblh, double *polyTblm, double *polyTbll,
10 double rh, double rm, double rl,
11 double tbl1h, double tbl1m, double tbl1l,
12 double tbl2h, double tbl2m, double tbl2l);
13
14
15 /* Function exp13
16
17 Computes exp(x) with an accuracy of 113 bits as
18
19 2^exponent * (exph + expm + expl) \approx exp(x)
20
21 Unless the subnormal case for x, no special cases are
22 handled.
23
24 The triple-double exph + expm + expl is non-overlapping.
25 The domain for exph + expm + expl is 1/2..2
26 The integer exponent is in the range -1024..1024. The
27 value 2^(exponent) may therefore be non-representable
28 whereas 2^exponent * (exph + expm + expl) is.
29
30 */
31
32
exp13(int * exponent,double * exph,double * expm,double * expl,double x)33 void exp13(int *exponent, double *exph, double *expm, double *expl, double x) {
34 double rh, rm, rl, tbl1h, tbl1m, tbl1l;
35 double tbl2h, tbl2m, tbl2l;
36 double xMultLog2InvMult2L, shiftedXMult, kd;
37 double msLog2Div2LMultKh, msLog2Div2LMultKm, msLog2Div2LMultKl;
38 double t1, t2;
39 db_number shiftedXMultdb, xdb;
40 int k, M, index1, index2, xIntHi;
41
42 /* Argument reduction and filtering for special cases */
43
44 /* Compute k as a double and as an int */
45 xdb.d = x;
46 xMultLog2InvMult2L = x * log2InvMult2L;
47 shiftedXMult = xMultLog2InvMult2L + shiftConst;
48 kd = shiftedXMult - shiftConst;
49 shiftedXMultdb.d = shiftedXMult;
50
51 /* Special cases tests */
52 xIntHi = xdb.i[HI];
53 /* Test if argument is a denormal or zero */
54 if ((xIntHi & 0x7ff00000) == 0) {
55 /* We are in the RN case, return 1.0 in all cases */
56 *exph = 1.0;
57 *expm = 0.0;
58 *expl = 0.0;
59 return;
60 }
61
62 k = shiftedXMultdb.i[LO];
63 M = k >> L;
64 index1 = k & INDEXMASK1;
65 index2 = (k & INDEXMASK2) >> LHALF;
66
67 /* Table reads */
68 tbl1h = twoPowerIndex1[index1].hi;
69 tbl1m = twoPowerIndex1[index1].mi;
70 tbl2h = twoPowerIndex2[index2].hi;
71 tbl2m = twoPowerIndex2[index2].mi;
72 tbl1l = twoPowerIndex1[index1].lo;
73 tbl2l = twoPowerIndex2[index2].lo;
74
75 /* Argument reduction */
76
77 Mul133(&msLog2Div2LMultKh,&msLog2Div2LMultKm,&msLog2Div2LMultKl,kd,msLog2Div2Lh,msLog2Div2Lm,msLog2Div2Ll);
78 t1 = x + msLog2Div2LMultKh;
79 Add12Cond(rh,t2,t1,msLog2Div2LMultKm);
80 Add12Cond(rm,rl,t2,msLog2Div2LMultKl);
81
82 /* Polynomial approximation and reconstruction: factorized code */
83
84 exp_td_accurate(exph, expm, expl, rh, rm, rl, tbl1h, tbl1m, tbl1l, tbl2h, tbl2m, tbl2l);
85
86 *exponent = M;
87 }
88