1 /* 2 * CDDL HEADER START 3 * 4 * The contents of this file are subject to the terms of the 5 * Common Development and Distribution License (the "License"). 6 * You may not use this file except in compliance with the License. 7 * 8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE 9 * or http://www.opensolaris.org/os/licensing. 10 * See the License for the specific language governing permissions 11 * and limitations under the License. 12 * 13 * When distributing Covered Code, include this CDDL HEADER in each 14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE. 15 * If applicable, add the following below this CDDL HEADER, with the 16 * fields enclosed by brackets "[]" replaced with your own identifying 17 * information: Portions Copyright [yyyy] [name of copyright owner] 18 * 19 * CDDL HEADER END 20 */ 21 22 /* 23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved. 24 */ 25 /* 26 * Copyright 2006 Sun Microsystems, Inc. All rights reserved. 27 * Use is subject to license terms. 28 */ 29 30 #pragma weak __logl = logl 31 32 /* 33 * logl(x) 34 * Table look-up algorithm 35 * By K.C. Ng, March 6, 1989 36 * 37 * (a). For x in [31/33,33/31], using a special approximation: 38 * f = x - 1; 39 * s = f/(2.0+f); ... here |s| <= 0.03125 40 * z = s*s; 41 * return f-s*(f-z*(B1+z*(B2+z*(B3+z*(B4+...+z*B9)...)))); 42 * 43 * (b). Otherwise, normalize x = 2^n * 1.f. 44 * Use a 6-bit table look-up: find a 6 bit g that match f to 6.5 bits, 45 * then 46 * log(x) = n*ln2 + log(1.g) + log(1.f/1.g). 47 * Here the leading and trailing values of log(1.g) are obtained from 48 * a size-64 table. 49 * For log(1.f/1.g), let s = (1.f-1.g)/(1.f+1.g), then 50 * log(1.f/1.g) = log((1+s)/(1-s)) = 2s + 2/3 s^3 + 2/5 s^5 +... 51 * Note that |s|<2**-8=0.00390625. We use an odd s-polynomial 52 * approximation to compute log(1.f/1.g): 53 * s*(A1+s^2*(A2+s^2*(A3+s^2*(A4+s^2*(A5+s^2*(A6+s^2*A7)))))) 54 * (Precision is 2**-136.91 bits, absolute error) 55 * 56 * (c). The final result is computed by 57 * (n*ln2_hi+_TBL_logl_hi[j]) + 58 * ( (n*ln2_lo+_TBL_logl_lo[j]) + s*(A1+...) ) 59 * 60 * Note. 61 * For ln2_hi and _TBL_logl_hi[j], we force their last 32 bit to be zero 62 * so that n*ln2_hi + _TBL_logl_hi[j] is exact. Here 63 * _TBL_logl_hi[j] + _TBL_logl_lo[j] match log(1+j*2**-6) to 194 bits 64 * 65 * 66 * Special cases: 67 * log(x) is NaN with signal if x < 0 (including -INF) ; 68 * log(+INF) is +INF; log(0) is -INF with signal; 69 * log(NaN) is that NaN with no signal. 70 * 71 * Constants: 72 * The hexadecimal values are the intended ones for the following constants. 73 * The decimal values may be used, provided that the compiler will convert 74 * from decimal to binary accurately enough to produce the hexadecimal values 75 * shown. 76 */ 77 78 #include "libm.h" 79 80 extern const long double _TBL_logl_hi[], _TBL_logl_lo[]; 81 82 static const long double 83 zero = 0.0L, 84 one = 1.0L, 85 two = 2.0L, 86 two113 = 10384593717069655257060992658440192.0L, 87 ln2hi = 6.931471805599453094172319547495844850203e-0001L, 88 ln2lo = 1.667085920830552208890449330400379754169e-0025L, 89 A1 = 2.000000000000000000000000000000000000024e+0000L, 90 A2 = 6.666666666666666666666666666666091393804e-0001L, 91 A3 = 4.000000000000000000000000407167070220671e-0001L, 92 A4 = 2.857142857142857142730077490612903681164e-0001L, 93 A5 = 2.222222222222242577702836920812882605099e-0001L, 94 A6 = 1.818181816435493395985912667105885828356e-0001L, 95 A7 = 1.538537835211839751112067512805496931725e-0001L, 96 B1 = 6.666666666666666666666666666666961498329e-0001L, 97 B2 = 3.999999999999999999999999990037655042358e-0001L, 98 B3 = 2.857142857142857142857273426428347457918e-0001L, 99 B4 = 2.222222222222222221353229049747910109566e-0001L, 100 B5 = 1.818181818181821503532559306309070138046e-0001L, 101 B6 = 1.538461538453809210486356084587356788556e-0001L, 102 B7 = 1.333333344463358756121456892645178795480e-0001L, 103 B8 = 1.176460904783899064854645174603360383792e-0001L, 104 B9 = 1.057293869956598995326368602518056990746e-0001L; 105 106 long double 107 logl(long double x) { 108 long double f, s, z, qn, h, t; 109 int *px = (int *) &x; 110 int *pz = (int *) &z; 111 int i, j, ix, i0, i1, n; 112 113 /* get long double precision word ordering */ 114 if (*(int *) &one == 0) { 115 i0 = 3; 116 i1 = 0; 117 } else { 118 i0 = 0; 119 i1 = 3; 120 } 121 122 n = 0; 123 ix = px[i0]; 124 if (ix > 0x3ffee0f8) { /* if x > 31/33 */ 125 if (ix < 0x3fff1084) { /* if x < 33/31 */ 126 f = x - one; 127 z = f * f; 128 if (((ix - 0x3fff0000) | px[i1] | px[2] | px[1]) == 0) { 129 return (zero); /* log(1)= +0 */ 130 } 131 s = f / (two + f); /* |s|<2**-8 */ 132 z = s * s; 133 return (f - s * (f - z * (B1 + z * (B2 + z * (B3 + 134 z * (B4 + z * (B5 + z * (B6 + z * (B7 + 135 z * (B8 + z * B9)))))))))); 136 } 137 if (ix >= 0x7fff0000) 138 return (x + x); /* x is +inf or NaN */ 139 goto LARGE_N; 140 } 141 if (ix >= 0x00010000) 142 goto LARGE_N; 143 i = ix & 0x7fffffff; 144 if ((i | px[i1] | px[2] | px[1]) == 0) { 145 px[i0] |= 0x80000000; 146 return (one / x); /* log(0.0) = -inf */ 147 } 148 if (ix < 0) { 149 if ((unsigned) ix >= 0xffff0000) 150 return (x - x); /* x is -inf or NaN */ 151 return (zero / zero); /* log(x<0) is NaN */ 152 } 153 /* subnormal x */ 154 x *= two113; 155 n = -113; 156 ix = px[i0]; 157 LARGE_N: 158 n += ((ix + 0x200) >> 16) - 0x3fff; 159 ix = (ix & 0x0000ffff) | 0x3fff0000; /* scale x to [1,2] */ 160 px[i0] = ix; 161 i = ix + 0x200; 162 pz[i0] = i & 0xfffffc00; 163 pz[i1] = pz[1] = pz[2] = 0; 164 s = (x - z) / (x + z); 165 j = (i >> 10) & 0x3f; 166 z = s * s; 167 qn = (long double) n; 168 t = qn * ln2lo + _TBL_logl_lo[j]; 169 h = qn * ln2hi + _TBL_logl_hi[j]; 170 f = t + s * (A1 + z * (A2 + z * (A3 + z * (A4 + z * (A5 + 171 z * (A6 + z * A7)))))); 172 return (h + f); 173 } 174