125c28e83SPiotr Jasiukajtis /* 225c28e83SPiotr Jasiukajtis * CDDL HEADER START 325c28e83SPiotr Jasiukajtis * 425c28e83SPiotr Jasiukajtis * The contents of this file are subject to the terms of the 525c28e83SPiotr Jasiukajtis * Common Development and Distribution License (the "License"). 625c28e83SPiotr Jasiukajtis * You may not use this file except in compliance with the License. 725c28e83SPiotr Jasiukajtis * 825c28e83SPiotr Jasiukajtis * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE 925c28e83SPiotr Jasiukajtis * or http://www.opensolaris.org/os/licensing. 1025c28e83SPiotr Jasiukajtis * See the License for the specific language governing permissions 1125c28e83SPiotr Jasiukajtis * and limitations under the License. 1225c28e83SPiotr Jasiukajtis * 1325c28e83SPiotr Jasiukajtis * When distributing Covered Code, include this CDDL HEADER in each 1425c28e83SPiotr Jasiukajtis * file and include the License file at usr/src/OPENSOLARIS.LICENSE. 1525c28e83SPiotr Jasiukajtis * If applicable, add the following below this CDDL HEADER, with the 1625c28e83SPiotr Jasiukajtis * fields enclosed by brackets "[]" replaced with your own identifying 1725c28e83SPiotr Jasiukajtis * information: Portions Copyright [yyyy] [name of copyright owner] 1825c28e83SPiotr Jasiukajtis * 1925c28e83SPiotr Jasiukajtis * CDDL HEADER END 2025c28e83SPiotr Jasiukajtis */ 2125c28e83SPiotr Jasiukajtis 2225c28e83SPiotr Jasiukajtis /* 2325c28e83SPiotr Jasiukajtis * Copyright 2011 Nexenta Systems, Inc. All rights reserved. 2425c28e83SPiotr Jasiukajtis */ 2525c28e83SPiotr Jasiukajtis /* 2625c28e83SPiotr Jasiukajtis * Copyright 2006 Sun Microsystems, Inc. All rights reserved. 2725c28e83SPiotr Jasiukajtis * Use is subject to license terms. 2825c28e83SPiotr Jasiukajtis */ 2925c28e83SPiotr Jasiukajtis 3025c28e83SPiotr Jasiukajtis #ifdef __LITTLE_ENDIAN 3125c28e83SPiotr Jasiukajtis #define H0(x) *(3 + (int *) &x) 3225c28e83SPiotr Jasiukajtis #define H1(x) *(2 + (int *) &x) 3325c28e83SPiotr Jasiukajtis #define H2(x) *(1 + (int *) &x) 3425c28e83SPiotr Jasiukajtis #define H3(x) *(int *) &x 3525c28e83SPiotr Jasiukajtis #else 3625c28e83SPiotr Jasiukajtis #define H0(x) *(int *) &x 3725c28e83SPiotr Jasiukajtis #define H1(x) *(1 + (int *) &x) 3825c28e83SPiotr Jasiukajtis #define H2(x) *(2 + (int *) &x) 3925c28e83SPiotr Jasiukajtis #define H3(x) *(3 + (int *) &x) 4025c28e83SPiotr Jasiukajtis #endif 4125c28e83SPiotr Jasiukajtis 4225c28e83SPiotr Jasiukajtis /* 4325c28e83SPiotr Jasiukajtis * log1pl(x) 4425c28e83SPiotr Jasiukajtis * Table look-up algorithm by modifying logl.c 4525c28e83SPiotr Jasiukajtis * By K.C. Ng, July 6, 1995 4625c28e83SPiotr Jasiukajtis * 4725c28e83SPiotr Jasiukajtis * (a). For 1+x in [31/33,33/31], using a special approximation: 4825c28e83SPiotr Jasiukajtis * s = x/(2.0+x); ... here |s| <= 0.03125 4925c28e83SPiotr Jasiukajtis * z = s*s; 5025c28e83SPiotr Jasiukajtis * return x-s*(x-z*(B1+z*(B2+z*(B3+z*(B4+...+z*B9)...)))); 5125c28e83SPiotr Jasiukajtis * (i.e., x is in [-2/33,2/31]) 5225c28e83SPiotr Jasiukajtis * 5325c28e83SPiotr Jasiukajtis * (b). Otherwise, normalize 1+x = 2^n * 1.f. 5425c28e83SPiotr Jasiukajtis * Here we may need a correction term for 1+x rounded. 5525c28e83SPiotr Jasiukajtis * Use a 6-bit table look-up: find a 6 bit g that match f to 6.5 bits, 5625c28e83SPiotr Jasiukajtis * then 5725c28e83SPiotr Jasiukajtis * log(1+x) = n*ln2 + log(1.g) + log(1.f/1.g). 5825c28e83SPiotr Jasiukajtis * Here the leading and trailing values of log(1.g) are obtained from 5925c28e83SPiotr Jasiukajtis * a size-64 table. 6025c28e83SPiotr Jasiukajtis * For log(1.f/1.g), let s = (1.f-1.g)/(1.f+1.g). Note that 6125c28e83SPiotr Jasiukajtis * 1.f = 2^-n(1+x) 6225c28e83SPiotr Jasiukajtis * 6325c28e83SPiotr Jasiukajtis * then 6425c28e83SPiotr Jasiukajtis * log(1.f/1.g) = log((1+s)/(1-s)) = 2s + 2/3 s^3 + 2/5 s^5 +... 6525c28e83SPiotr Jasiukajtis * Note that |s|<2**-8=0.00390625. We use an odd s-polynomial 6625c28e83SPiotr Jasiukajtis * approximation to compute log(1.f/1.g): 6725c28e83SPiotr Jasiukajtis * s*(A1+s^2*(A2+s^2*(A3+s^2*(A4+s^2*(A5+s^2*(A6+s^2*A7)))))) 6825c28e83SPiotr Jasiukajtis * (Precision is 2**-136.91 bits, absolute error) 6925c28e83SPiotr Jasiukajtis * 7025c28e83SPiotr Jasiukajtis * CAUTION: 7125c28e83SPiotr Jasiukajtis * For x>=1, compute 1+x will lost one bit (OK). 7225c28e83SPiotr Jasiukajtis * For x in [-0.5,-1), 1+x is exact. 7325c28e83SPiotr Jasiukajtis * For x in (-0.5,-2/33]U[2/31,1), up to 4 last bits of x will be lost 7425c28e83SPiotr Jasiukajtis * in 1+x. Therefore, to recover the lost bits, one need to compute 7525c28e83SPiotr Jasiukajtis * 1.f-1.g accurately. 7625c28e83SPiotr Jasiukajtis * 7725c28e83SPiotr Jasiukajtis * Let hx = HI(x), m = (hx>>16)-0x3fff (=ilogbl(x)), note that 7825c28e83SPiotr Jasiukajtis * -2/33 = -0.0606...= 2^-5 * 1.939..., 7925c28e83SPiotr Jasiukajtis * 2/31 = 0.09375 = 2^-4 * 1.500..., 8025c28e83SPiotr Jasiukajtis * so for x in (-0.5,-2/33], -5<=m<=-2, n= -1, 1+f=2*(1+x) 8125c28e83SPiotr Jasiukajtis * for x in [2/33,1), -4<=m<=-1, n= 0, f=x 8225c28e83SPiotr Jasiukajtis * 8325c28e83SPiotr Jasiukajtis * In short: 8425c28e83SPiotr Jasiukajtis * if x>0, let g: hg= ((hx + (0x200<<(-m)))>>(10-m))<<(10-m) 8525c28e83SPiotr Jasiukajtis * then 1.f-1.g = x-g 8625c28e83SPiotr Jasiukajtis * if x<0, let g': hg' =((ix-(0x200)<<(-m-1))>>(9-m))<<(9-m) 8725c28e83SPiotr Jasiukajtis * (ix=hx&0x7fffffff) 8825c28e83SPiotr Jasiukajtis * then 1.f-1.g = 2*(g'+x), 8925c28e83SPiotr Jasiukajtis * 9025c28e83SPiotr Jasiukajtis * (c). The final result is computed by 9125c28e83SPiotr Jasiukajtis * (n*ln2_hi+_TBL_logl_hi[j]) + 9225c28e83SPiotr Jasiukajtis * ( (n*ln2_lo+_TBL_logl_lo[j]) + s*(A1+...) ) 9325c28e83SPiotr Jasiukajtis * 9425c28e83SPiotr Jasiukajtis * Note. 9525c28e83SPiotr Jasiukajtis * For ln2_hi and _TBL_logl_hi[j], we force their last 32 bit to be zero 9625c28e83SPiotr Jasiukajtis * so that n*ln2_hi + _TBL_logl_hi[j] is exact. Here 9725c28e83SPiotr Jasiukajtis * _TBL_logl_hi[j] + _TBL_logl_lo[j] match log(1+j*2**-6) to 194 bits 9825c28e83SPiotr Jasiukajtis * 9925c28e83SPiotr Jasiukajtis * 10025c28e83SPiotr Jasiukajtis * Special cases: 10125c28e83SPiotr Jasiukajtis * log(x) is NaN with signal if x < 0 (including -INF) ; 10225c28e83SPiotr Jasiukajtis * log(+INF) is +INF; log(0) is -INF with signal; 10325c28e83SPiotr Jasiukajtis * log(NaN) is that NaN with no signal. 10425c28e83SPiotr Jasiukajtis * 10525c28e83SPiotr Jasiukajtis * Constants: 10625c28e83SPiotr Jasiukajtis * The hexadecimal values are the intended ones for the following constants. 10725c28e83SPiotr Jasiukajtis * The decimal values may be used, provided that the compiler will convert 10825c28e83SPiotr Jasiukajtis * from decimal to binary accurately enough to produce the hexadecimal values 10925c28e83SPiotr Jasiukajtis * shown. 11025c28e83SPiotr Jasiukajtis */ 11125c28e83SPiotr Jasiukajtis 112*ddc0e0b5SRichard Lowe #pragma weak __log1pl = log1pl 11325c28e83SPiotr Jasiukajtis 11425c28e83SPiotr Jasiukajtis #include "libm.h" 11525c28e83SPiotr Jasiukajtis 11625c28e83SPiotr Jasiukajtis extern const long double _TBL_logl_hi[], _TBL_logl_lo[]; 11725c28e83SPiotr Jasiukajtis 11825c28e83SPiotr Jasiukajtis static const long double 11925c28e83SPiotr Jasiukajtis zero = 0.0L, 12025c28e83SPiotr Jasiukajtis one = 1.0L, 12125c28e83SPiotr Jasiukajtis two = 2.0L, 12225c28e83SPiotr Jasiukajtis ln2hi = 6.931471805599453094172319547495844850203e-0001L, 12325c28e83SPiotr Jasiukajtis ln2lo = 1.667085920830552208890449330400379754169e-0025L, 12425c28e83SPiotr Jasiukajtis A1 = 2.000000000000000000000000000000000000024e+0000L, 12525c28e83SPiotr Jasiukajtis A2 = 6.666666666666666666666666666666091393804e-0001L, 12625c28e83SPiotr Jasiukajtis A3 = 4.000000000000000000000000407167070220671e-0001L, 12725c28e83SPiotr Jasiukajtis A4 = 2.857142857142857142730077490612903681164e-0001L, 12825c28e83SPiotr Jasiukajtis A5 = 2.222222222222242577702836920812882605099e-0001L, 12925c28e83SPiotr Jasiukajtis A6 = 1.818181816435493395985912667105885828356e-0001L, 13025c28e83SPiotr Jasiukajtis A7 = 1.538537835211839751112067512805496931725e-0001L, 13125c28e83SPiotr Jasiukajtis B1 = 6.666666666666666666666666666666961498329e-0001L, 13225c28e83SPiotr Jasiukajtis B2 = 3.999999999999999999999999990037655042358e-0001L, 13325c28e83SPiotr Jasiukajtis B3 = 2.857142857142857142857273426428347457918e-0001L, 13425c28e83SPiotr Jasiukajtis B4 = 2.222222222222222221353229049747910109566e-0001L, 13525c28e83SPiotr Jasiukajtis B5 = 1.818181818181821503532559306309070138046e-0001L, 13625c28e83SPiotr Jasiukajtis B6 = 1.538461538453809210486356084587356788556e-0001L, 13725c28e83SPiotr Jasiukajtis B7 = 1.333333344463358756121456892645178795480e-0001L, 13825c28e83SPiotr Jasiukajtis B8 = 1.176460904783899064854645174603360383792e-0001L, 13925c28e83SPiotr Jasiukajtis B9 = 1.057293869956598995326368602518056990746e-0001L; 14025c28e83SPiotr Jasiukajtis 14125c28e83SPiotr Jasiukajtis long double 14225c28e83SPiotr Jasiukajtis log1pl(long double x) { 14325c28e83SPiotr Jasiukajtis long double f, s, z, qn, h, t, y, g; 14425c28e83SPiotr Jasiukajtis int i, j, ix, iy, n, hx, m; 14525c28e83SPiotr Jasiukajtis 14625c28e83SPiotr Jasiukajtis hx = H0(x); 14725c28e83SPiotr Jasiukajtis ix = hx & 0x7fffffff; 14825c28e83SPiotr Jasiukajtis if (ix < 0x3ffaf07c) { /* |x|<2/33 */ 14925c28e83SPiotr Jasiukajtis if (ix <= 0x3f8d0000) { /* x <= 2**-114, return x */ 15025c28e83SPiotr Jasiukajtis if ((int) x == 0) 15125c28e83SPiotr Jasiukajtis return (x); 15225c28e83SPiotr Jasiukajtis } 15325c28e83SPiotr Jasiukajtis s = x / (two + x); /* |s|<2**-8 */ 15425c28e83SPiotr Jasiukajtis z = s * s; 15525c28e83SPiotr Jasiukajtis return (x - s * (x - z * (B1 + z * (B2 + z * (B3 + z * (B4 + 15625c28e83SPiotr Jasiukajtis z * (B5 + z * (B6 + z * (B7 + z * (B8 + z * B9)))))))))); 15725c28e83SPiotr Jasiukajtis } 15825c28e83SPiotr Jasiukajtis if (ix >= 0x7fff0000) { /* x is +inf or NaN */ 15925c28e83SPiotr Jasiukajtis return (x + fabsl(x)); 16025c28e83SPiotr Jasiukajtis } 16125c28e83SPiotr Jasiukajtis if (hx < 0 && ix >= 0x3fff0000) { 16225c28e83SPiotr Jasiukajtis if (ix > 0x3fff0000 || (H1(x) | H2(x) | H3(x)) != 0) 16325c28e83SPiotr Jasiukajtis x = zero; 16425c28e83SPiotr Jasiukajtis return (x / zero); /* log1p(x) is NaN if x<-1 */ 16525c28e83SPiotr Jasiukajtis /* log1p(-1) is -inf */ 16625c28e83SPiotr Jasiukajtis } 16725c28e83SPiotr Jasiukajtis if (ix >= 0x7ffeffff) 16825c28e83SPiotr Jasiukajtis y = x; /* avoid spurious overflow */ 16925c28e83SPiotr Jasiukajtis else 17025c28e83SPiotr Jasiukajtis y = one + x; 17125c28e83SPiotr Jasiukajtis iy = H0(y); 17225c28e83SPiotr Jasiukajtis n = ((iy + 0x200) >> 16) - 0x3fff; 17325c28e83SPiotr Jasiukajtis iy = (iy & 0x0000ffff) | 0x3fff0000; /* scale 1+x to [1,2] */ 17425c28e83SPiotr Jasiukajtis H0(y) = iy; 17525c28e83SPiotr Jasiukajtis z = zero; 17625c28e83SPiotr Jasiukajtis m = (ix >> 16) - 0x3fff; 17725c28e83SPiotr Jasiukajtis /* HI(1+x) = (((hx&0xffff)|0x10000)>>(-m))|0x3fff0000 */ 17825c28e83SPiotr Jasiukajtis if (n == 0) { /* x in [2/33,1) */ 17925c28e83SPiotr Jasiukajtis g = zero; 18025c28e83SPiotr Jasiukajtis H0(g) = ((hx + (0x200 << (-m))) >> (10 - m)) << (10 - m); 18125c28e83SPiotr Jasiukajtis t = x - g; 18225c28e83SPiotr Jasiukajtis i = (((((hx & 0xffff) | 0x10000) >> (-m)) | 0x3fff0000) + 18325c28e83SPiotr Jasiukajtis 0x200) >> 10; 18425c28e83SPiotr Jasiukajtis H0(z) = i << 10; 18525c28e83SPiotr Jasiukajtis 18625c28e83SPiotr Jasiukajtis } else if ((1 + n) == 0 && (ix < 0x3ffe0000)) { /* x in (-0.5,-2/33] */ 18725c28e83SPiotr Jasiukajtis g = zero; 18825c28e83SPiotr Jasiukajtis H0(g) = ((ix + (0x200 << (-m - 1))) >> (9 - m)) << (9 - m); 18925c28e83SPiotr Jasiukajtis t = g + x; 19025c28e83SPiotr Jasiukajtis t = t + t; 19125c28e83SPiotr Jasiukajtis /* 19225c28e83SPiotr Jasiukajtis * HI(2*(1+x)) = 19325c28e83SPiotr Jasiukajtis * ((0x10000-(((hx&0xffff)|0x10000)>>(-m)))<<1)|0x3fff0000 19425c28e83SPiotr Jasiukajtis */ 19525c28e83SPiotr Jasiukajtis /* 19625c28e83SPiotr Jasiukajtis * i = 19725c28e83SPiotr Jasiukajtis * ((((0x10000-(((hx&0xffff)|0x10000)>>(-m)))<<1)|0x3fff0000)+ 19825c28e83SPiotr Jasiukajtis * 0x200)>>10; H0(z)=i<<10; 19925c28e83SPiotr Jasiukajtis */ 20025c28e83SPiotr Jasiukajtis z = two * (one - g); 20125c28e83SPiotr Jasiukajtis i = H0(z) >> 10; 20225c28e83SPiotr Jasiukajtis } else { 20325c28e83SPiotr Jasiukajtis i = (iy + 0x200) >> 10; 20425c28e83SPiotr Jasiukajtis H0(z) = i << 10; 20525c28e83SPiotr Jasiukajtis t = y - z; 20625c28e83SPiotr Jasiukajtis } 20725c28e83SPiotr Jasiukajtis 20825c28e83SPiotr Jasiukajtis s = t / (y + z); 20925c28e83SPiotr Jasiukajtis j = i & 0x3f; 21025c28e83SPiotr Jasiukajtis z = s * s; 21125c28e83SPiotr Jasiukajtis qn = (long double) n; 21225c28e83SPiotr Jasiukajtis t = qn * ln2lo + _TBL_logl_lo[j]; 21325c28e83SPiotr Jasiukajtis h = qn * ln2hi + _TBL_logl_hi[j]; 21425c28e83SPiotr Jasiukajtis f = t + s * (A1 + z * (A2 + z * (A3 + z * (A4 + z * (A5 + z * (A6 + 21525c28e83SPiotr Jasiukajtis z * A7)))))); 21625c28e83SPiotr Jasiukajtis return (h + f); 21725c28e83SPiotr Jasiukajtis } 218