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