|
38 | 38 | #define LF1 1.24999999978138668903e-02 |
39 | 39 | #define LF2 2.23219810758559851206e-03 |
40 | 40 |
|
41 | | -_CLC_DEF void __clc_ep_log(double x, int *xexp, double *r1, double *r2) |
42 | | -{ |
43 | | - // Computes natural log(x). Algorithm based on: |
44 | | - // Ping-Tak Peter Tang |
45 | | - // "Table-driven implementation of the logarithm function in IEEE |
46 | | - // floating-point arithmetic" |
47 | | - // ACM Transactions on Mathematical Software (TOMS) |
48 | | - // Volume 16, Issue 4 (December 1990) |
49 | | - int near_one = x >= 0x1.e0faap-1 & x <= 0x1.1082cp+0; |
50 | | - |
51 | | - ulong ux = as_ulong(x); |
52 | | - ulong uxs = as_ulong(as_double(0x03d0000000000000UL | ux) - 0x1.0p-962); |
53 | | - int c = ux < IMPBIT_DP64; |
54 | | - ux = c ? uxs : ux; |
55 | | - int expadjust = c ? 60 : 0; |
56 | | - |
57 | | - // Store the exponent of x in xexp and put f into the range [0.5,1) |
58 | | - int xexp1 = ((as_int2(ux).hi >> 20) & 0x7ff) - EXPBIAS_DP64 - expadjust; |
59 | | - double f = as_double(HALFEXPBITS_DP64 | (ux & MANTBITS_DP64)); |
60 | | - *xexp = near_one ? 0 : xexp1; |
61 | | - |
62 | | - double r = x - 1.0; |
63 | | - double u1 = MATH_DIVIDE(r, 2.0 + r); |
64 | | - double ru1 = -r * u1; |
65 | | - u1 = u1 + u1; |
66 | | - |
67 | | - int index = as_int2(ux).hi >> 13; |
68 | | - index = ((0x80 | (index & 0x7e)) >> 1) + (index & 0x1); |
69 | | - |
70 | | - double f1 = index * 0x1.0p-7; |
71 | | - double f2 = f - f1; |
72 | | - double u2 = MATH_DIVIDE(f2, fma(0.5, f2, f1)); |
73 | | - |
74 | | - double2 tv = USE_TABLE(ln_tbl, (index - 64)); |
75 | | - double z1 = tv.s0; |
76 | | - double q = tv.s1; |
77 | | - |
78 | | - z1 = near_one ? r : z1; |
79 | | - q = near_one ? 0.0 : q; |
80 | | - double u = near_one ? u1 : u2; |
81 | | - double v = u*u; |
82 | | - |
83 | | - double cc = near_one ? ru1 : u2; |
84 | | - |
85 | | - double z21 = fma(v, fma(v, fma(v, LN3, LN2), LN1), LN0); |
86 | | - double z22 = fma(v, fma(v, LF2, LF1), LF0); |
87 | | - double z2 = near_one ? z21 : z22; |
88 | | - z2 = fma(u*v, z2, cc) + q; |
89 | | - |
90 | | - *r1 = z1; |
91 | | - *r2 = z2; |
| 41 | +_CLC_DEF void __clc_ep_log(double x, int *xexp, double *r1, double *r2) { |
| 42 | + // Computes natural log(x). Algorithm based on: |
| 43 | + // Ping-Tak Peter Tang |
| 44 | + // "Table-driven implementation of the logarithm function in IEEE |
| 45 | + // floating-point arithmetic" |
| 46 | + // ACM Transactions on Mathematical Software (TOMS) |
| 47 | + // Volume 16, Issue 4 (December 1990) |
| 48 | + int near_one = x >= 0x1.e0faap-1 & x <= 0x1.1082cp+0; |
| 49 | + |
| 50 | + ulong ux = as_ulong(x); |
| 51 | + ulong uxs = as_ulong(as_double(0x03d0000000000000UL | ux) - 0x1.0p-962); |
| 52 | + int c = ux < IMPBIT_DP64; |
| 53 | + ux = c ? uxs : ux; |
| 54 | + int expadjust = c ? 60 : 0; |
| 55 | + |
| 56 | + // Store the exponent of x in xexp and put f into the range [0.5,1) |
| 57 | + int xexp1 = ((as_int2(ux).hi >> 20) & 0x7ff) - EXPBIAS_DP64 - expadjust; |
| 58 | + double f = as_double(HALFEXPBITS_DP64 | (ux & MANTBITS_DP64)); |
| 59 | + *xexp = near_one ? 0 : xexp1; |
| 60 | + |
| 61 | + double r = x - 1.0; |
| 62 | + double u1 = MATH_DIVIDE(r, 2.0 + r); |
| 63 | + double ru1 = -r * u1; |
| 64 | + u1 = u1 + u1; |
| 65 | + |
| 66 | + int index = as_int2(ux).hi >> 13; |
| 67 | + index = ((0x80 | (index & 0x7e)) >> 1) + (index & 0x1); |
| 68 | + |
| 69 | + double f1 = index * 0x1.0p-7; |
| 70 | + double f2 = f - f1; |
| 71 | + double u2 = MATH_DIVIDE(f2, fma(0.5, f2, f1)); |
| 72 | + |
| 73 | + double2 tv = USE_TABLE(ln_tbl, (index - 64)); |
| 74 | + double z1 = tv.s0; |
| 75 | + double q = tv.s1; |
| 76 | + |
| 77 | + z1 = near_one ? r : z1; |
| 78 | + q = near_one ? 0.0 : q; |
| 79 | + double u = near_one ? u1 : u2; |
| 80 | + double v = u * u; |
| 81 | + |
| 82 | + double cc = near_one ? ru1 : u2; |
| 83 | + |
| 84 | + double z21 = fma(v, fma(v, fma(v, LN3, LN2), LN1), LN0); |
| 85 | + double z22 = fma(v, fma(v, LF2, LF1), LF0); |
| 86 | + double z2 = near_one ? z21 : z22; |
| 87 | + z2 = fma(u * v, z2, cc) + q; |
| 88 | + |
| 89 | + *r1 = z1; |
| 90 | + *r2 = z2; |
92 | 91 | } |
93 | 92 |
|
94 | 93 | #endif |
0 commit comments