|
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