82 lines · plain
1//===----------------------------------------------------------------------===//2//3// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.4// See https://llvm.org/LICENSE.txt for license information.5// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception6//7//===----------------------------------------------------------------------===//8//9// Computes natural log(x). Algorithm based on:10// Ping-Tak Peter Tang11// "Table-driven implementation of the logarithm function in IEEE floating-point12// arithmetic"13// ACM Transactions on Mathematical Software (TOMS) Volume 16, Issue 4 (December14// 1990)15//16//===----------------------------------------------------------------------===//17 18#if __CLC_FPSIZE == 6419 20#define LN0 8.33333333333317923934e-0221#define LN1 1.25000000037717509602e-0222#define LN2 2.23213998791944806202e-0323#define LN3 4.34887777707614552256e-0424 25#define LF0 8.33333333333333593622e-0226#define LF1 1.24999999978138668903e-0227#define LF2 2.23219810758559851206e-0328 29_CLC_DEF _CLC_OVERLOAD void __clc_ep_log(__CLC_GENTYPE x,30 private __CLC_INTN *xexp,31 private __CLC_GENTYPE *r1,32 private __CLC_GENTYPE *r2) {33 __CLC_LONGN near_one = x >= 0x1.e0faap-1 && x <= 0x1.1082cp+0;34 35 __CLC_ULONGN ux = __CLC_AS_ULONGN(x);36 __CLC_ULONGN uxs =37 __CLC_AS_ULONGN(__CLC_AS_GENTYPE(0x03d0000000000000UL | ux) - 0x1.0p-962);38 __CLC_LONGN c = ux < IMPBIT_DP64;39 ux = c ? uxs : ux;40 __CLC_INTN expadjust =41 __CLC_CONVERT_INTN(c ? (__CLC_LONGN)60 : (__CLC_LONGN)0);42 43 // Store the exponent of x in xexp and put f into the range [0.5,1)44 __CLC_INTN xexp1 = __CLC_CONVERT_INTN((ux >> EXPSHIFTBITS_DP64) & 0x7ff) -45 EXPBIAS_DP64 - expadjust;46 __CLC_GENTYPE f = __CLC_AS_GENTYPE(HALFEXPBITS_DP64 | (ux & MANTBITS_DP64));47 *xexp = __CLC_CONVERT_INTN(near_one) ? 0 : xexp1;48 49 __CLC_GENTYPE r = x - 1.0;50 __CLC_GENTYPE u1 = MATH_DIVIDE(r, 2.0 + r);51 __CLC_GENTYPE ru1 = -r * u1;52 u1 = u1 + u1;53 54 __CLC_INTN index = __CLC_CONVERT_INTN(ux >> 45);55 index = ((0x80 | (index & 0x7e)) >> 1) + (index & 0x1);56 57 __CLC_GENTYPE f1 = __CLC_CONVERT_GENTYPE(index) * 0x1.0p-7;58 __CLC_GENTYPE f2 = f - f1;59 __CLC_GENTYPE u2 = MATH_DIVIDE(f2, __clc_fma(0.5, f2, f1));60 61 __CLC_GENTYPE z1 = __CLC_USE_TABLE(ln_tbl_lo, (index - 64));62 __CLC_GENTYPE q = __CLC_USE_TABLE(ln_tbl_hi, (index - 64));63 64 z1 = near_one ? r : z1;65 q = near_one ? 0.0 : q;66 __CLC_GENTYPE u = near_one ? u1 : u2;67 __CLC_GENTYPE v = u * u;68 69 __CLC_GENTYPE cc = near_one ? ru1 : u2;70 71 __CLC_GENTYPE z21 =72 __clc_fma(v, __clc_fma(v, __clc_fma(v, LN3, LN2), LN1), LN0);73 __CLC_GENTYPE z22 = __clc_fma(v, __clc_fma(v, LF2, LF1), LF0);74 __CLC_GENTYPE z2 = near_one ? z21 : z22;75 z2 = __clc_fma(u * v, z2, cc) + q;76 77 *r1 = z1;78 *r2 = z2;79}80 81#endif82