xref: /illumos-gate/usr/src/lib/libm/common/Q/logl.c (revision ddc0e0b5)
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 
30*ddc0e0b5SRichard Lowe #pragma weak __logl = logl
3125c28e83SPiotr Jasiukajtis 
3225c28e83SPiotr Jasiukajtis /*
3325c28e83SPiotr Jasiukajtis  * logl(x)
3425c28e83SPiotr Jasiukajtis  * Table look-up algorithm
3525c28e83SPiotr Jasiukajtis  * By K.C. Ng, March 6, 1989
3625c28e83SPiotr Jasiukajtis  *
3725c28e83SPiotr Jasiukajtis  * (a). For x in [31/33,33/31], using a special approximation:
3825c28e83SPiotr Jasiukajtis  *	f = x - 1;
3925c28e83SPiotr Jasiukajtis  *	s = f/(2.0+f);	... here |s| <= 0.03125
4025c28e83SPiotr Jasiukajtis  *	z = s*s;
4125c28e83SPiotr Jasiukajtis  *	return f-s*(f-z*(B1+z*(B2+z*(B3+z*(B4+...+z*B9)...))));
4225c28e83SPiotr Jasiukajtis  *
4325c28e83SPiotr Jasiukajtis  * (b). Otherwise, normalize x = 2^n * 1.f.
4425c28e83SPiotr Jasiukajtis  *	Use a 6-bit table look-up: find a 6 bit g that match f to 6.5 bits,
4525c28e83SPiotr Jasiukajtis  *	then
4625c28e83SPiotr Jasiukajtis  *	    log(x) = n*ln2 + log(1.g) + log(1.f/1.g).
4725c28e83SPiotr Jasiukajtis  *	Here the leading and trailing values of log(1.g) are obtained from
4825c28e83SPiotr Jasiukajtis  *	a size-64 table.
4925c28e83SPiotr Jasiukajtis  *	For log(1.f/1.g), let s = (1.f-1.g)/(1.f+1.g), then
5025c28e83SPiotr Jasiukajtis  *	    log(1.f/1.g) = log((1+s)/(1-s)) = 2s + 2/3 s^3 + 2/5 s^5 +...
5125c28e83SPiotr Jasiukajtis  *	Note that |s|<2**-8=0.00390625. We use an odd s-polynomial
5225c28e83SPiotr Jasiukajtis  *	approximation to compute log(1.f/1.g):
5325c28e83SPiotr Jasiukajtis  *		s*(A1+s^2*(A2+s^2*(A3+s^2*(A4+s^2*(A5+s^2*(A6+s^2*A7))))))
5425c28e83SPiotr Jasiukajtis  *	(Precision is 2**-136.91 bits, absolute error)
5525c28e83SPiotr Jasiukajtis  *
5625c28e83SPiotr Jasiukajtis  * (c). The final result is computed by
5725c28e83SPiotr Jasiukajtis  *		(n*ln2_hi+_TBL_logl_hi[j]) +
5825c28e83SPiotr Jasiukajtis  *			( (n*ln2_lo+_TBL_logl_lo[j]) + s*(A1+...) )
5925c28e83SPiotr Jasiukajtis  *
6025c28e83SPiotr Jasiukajtis  * Note.
6125c28e83SPiotr Jasiukajtis  *	For ln2_hi and _TBL_logl_hi[j], we force their last 32 bit to be zero
6225c28e83SPiotr Jasiukajtis  *	so that n*ln2_hi + _TBL_logl_hi[j] is exact. Here
6325c28e83SPiotr Jasiukajtis  *	_TBL_logl_hi[j] + _TBL_logl_lo[j] match log(1+j*2**-6) to 194 bits
6425c28e83SPiotr Jasiukajtis  *
6525c28e83SPiotr Jasiukajtis  *
6625c28e83SPiotr Jasiukajtis  * Special cases:
6725c28e83SPiotr Jasiukajtis  *	log(x) is NaN with signal if x < 0 (including -INF) ;
6825c28e83SPiotr Jasiukajtis  *	log(+INF) is +INF; log(0) is -INF with signal;
6925c28e83SPiotr Jasiukajtis  *	log(NaN) is that NaN with no signal.
7025c28e83SPiotr Jasiukajtis  *
7125c28e83SPiotr Jasiukajtis  * Constants:
7225c28e83SPiotr Jasiukajtis  * The hexadecimal values are the intended ones for the following constants.
7325c28e83SPiotr Jasiukajtis  * The decimal values may be used, provided that the compiler will convert
7425c28e83SPiotr Jasiukajtis  * from decimal to binary accurately enough to produce the hexadecimal values
7525c28e83SPiotr Jasiukajtis  * shown.
7625c28e83SPiotr Jasiukajtis  */
7725c28e83SPiotr Jasiukajtis 
7825c28e83SPiotr Jasiukajtis #include "libm.h"
7925c28e83SPiotr Jasiukajtis 
8025c28e83SPiotr Jasiukajtis extern const long double _TBL_logl_hi[], _TBL_logl_lo[];
8125c28e83SPiotr Jasiukajtis 
8225c28e83SPiotr Jasiukajtis static const long double
8325c28e83SPiotr Jasiukajtis 	zero	=   0.0L,
8425c28e83SPiotr Jasiukajtis 	one	=   1.0L,
8525c28e83SPiotr Jasiukajtis 	two	=   2.0L,
8625c28e83SPiotr Jasiukajtis 	two113  =   10384593717069655257060992658440192.0L,
8725c28e83SPiotr Jasiukajtis 	ln2hi	=   6.931471805599453094172319547495844850203e-0001L,
8825c28e83SPiotr Jasiukajtis 	ln2lo	=   1.667085920830552208890449330400379754169e-0025L,
8925c28e83SPiotr Jasiukajtis 	A1	=   2.000000000000000000000000000000000000024e+0000L,
9025c28e83SPiotr Jasiukajtis 	A2	=   6.666666666666666666666666666666091393804e-0001L,
9125c28e83SPiotr Jasiukajtis 	A3	=   4.000000000000000000000000407167070220671e-0001L,
9225c28e83SPiotr Jasiukajtis 	A4	=   2.857142857142857142730077490612903681164e-0001L,
9325c28e83SPiotr Jasiukajtis 	A5	=   2.222222222222242577702836920812882605099e-0001L,
9425c28e83SPiotr Jasiukajtis 	A6	=   1.818181816435493395985912667105885828356e-0001L,
9525c28e83SPiotr Jasiukajtis 	A7	=   1.538537835211839751112067512805496931725e-0001L,
9625c28e83SPiotr Jasiukajtis 	B1	=   6.666666666666666666666666666666961498329e-0001L,
9725c28e83SPiotr Jasiukajtis 	B2	=   3.999999999999999999999999990037655042358e-0001L,
9825c28e83SPiotr Jasiukajtis 	B3	=   2.857142857142857142857273426428347457918e-0001L,
9925c28e83SPiotr Jasiukajtis 	B4	=   2.222222222222222221353229049747910109566e-0001L,
10025c28e83SPiotr Jasiukajtis 	B5	=   1.818181818181821503532559306309070138046e-0001L,
10125c28e83SPiotr Jasiukajtis 	B6	=   1.538461538453809210486356084587356788556e-0001L,
10225c28e83SPiotr Jasiukajtis 	B7	=   1.333333344463358756121456892645178795480e-0001L,
10325c28e83SPiotr Jasiukajtis 	B8	=   1.176460904783899064854645174603360383792e-0001L,
10425c28e83SPiotr Jasiukajtis 	B9	=   1.057293869956598995326368602518056990746e-0001L;
10525c28e83SPiotr Jasiukajtis 
10625c28e83SPiotr Jasiukajtis long double
logl(long double x)10725c28e83SPiotr Jasiukajtis logl(long double x) {
10825c28e83SPiotr Jasiukajtis 	long double f, s, z, qn, h, t;
10925c28e83SPiotr Jasiukajtis 	int *px = (int *) &x;
11025c28e83SPiotr Jasiukajtis 	int *pz = (int *) &z;
11125c28e83SPiotr Jasiukajtis 	int i, j, ix, i0, i1, n;
11225c28e83SPiotr Jasiukajtis 
11325c28e83SPiotr Jasiukajtis 	/* get long double precision word ordering */
11425c28e83SPiotr Jasiukajtis 	if (*(int *) &one == 0) {
11525c28e83SPiotr Jasiukajtis 		i0 = 3;
11625c28e83SPiotr Jasiukajtis 		i1 = 0;
11725c28e83SPiotr Jasiukajtis 	} else {
11825c28e83SPiotr Jasiukajtis 		i0 = 0;
11925c28e83SPiotr Jasiukajtis 		i1 = 3;
12025c28e83SPiotr Jasiukajtis 	}
12125c28e83SPiotr Jasiukajtis 
12225c28e83SPiotr Jasiukajtis 	n = 0;
12325c28e83SPiotr Jasiukajtis 	ix = px[i0];
12425c28e83SPiotr Jasiukajtis 	if (ix > 0x3ffee0f8) {	/* if x >  31/33 */
12525c28e83SPiotr Jasiukajtis 		if (ix < 0x3fff1084) {	/* if x < 33/31 */
12625c28e83SPiotr Jasiukajtis 			f = x - one;
12725c28e83SPiotr Jasiukajtis 			z = f * f;
12825c28e83SPiotr Jasiukajtis 			if (((ix - 0x3fff0000) | px[i1] | px[2] | px[1]) == 0) {
12925c28e83SPiotr Jasiukajtis 				return (zero);	/* log(1)= +0 */
13025c28e83SPiotr Jasiukajtis 			}
13125c28e83SPiotr Jasiukajtis 			s = f / (two + f);	/* |s|<2**-8 */
13225c28e83SPiotr Jasiukajtis 			z = s * s;
13325c28e83SPiotr Jasiukajtis 			return (f - s * (f - z * (B1 + z * (B2 + z * (B3 +
13425c28e83SPiotr Jasiukajtis 				z * (B4 + z * (B5 + z * (B6 + z * (B7 +
13525c28e83SPiotr Jasiukajtis 				z * (B8 + z * B9))))))))));
13625c28e83SPiotr Jasiukajtis 		}
13725c28e83SPiotr Jasiukajtis 		if (ix >= 0x7fff0000)
13825c28e83SPiotr Jasiukajtis 			return (x + x);	/* x is +inf or NaN */
13925c28e83SPiotr Jasiukajtis 		goto LARGE_N;
14025c28e83SPiotr Jasiukajtis 	}
14125c28e83SPiotr Jasiukajtis 	if (ix >= 0x00010000)
14225c28e83SPiotr Jasiukajtis 		goto LARGE_N;
14325c28e83SPiotr Jasiukajtis 	i = ix & 0x7fffffff;
14425c28e83SPiotr Jasiukajtis 	if ((i | px[i1] | px[2] | px[1]) == 0) {
14525c28e83SPiotr Jasiukajtis 		px[i0] |= 0x80000000;
14625c28e83SPiotr Jasiukajtis 		return (one / x);	/* log(0.0) = -inf */
14725c28e83SPiotr Jasiukajtis 	}
14825c28e83SPiotr Jasiukajtis 	if (ix < 0) {
14925c28e83SPiotr Jasiukajtis 		if ((unsigned) ix >= 0xffff0000)
15025c28e83SPiotr Jasiukajtis 			return (x - x);	/* x is -inf or NaN */
15125c28e83SPiotr Jasiukajtis 		return (zero / zero);	/* log(x<0) is NaN  */
15225c28e83SPiotr Jasiukajtis 	}
15325c28e83SPiotr Jasiukajtis 	/* subnormal x */
15425c28e83SPiotr Jasiukajtis 	x *= two113;
15525c28e83SPiotr Jasiukajtis 	n = -113;
15625c28e83SPiotr Jasiukajtis 	ix = px[i0];
15725c28e83SPiotr Jasiukajtis LARGE_N:
15825c28e83SPiotr Jasiukajtis 	n += ((ix + 0x200) >> 16) - 0x3fff;
15925c28e83SPiotr Jasiukajtis 	ix = (ix & 0x0000ffff) | 0x3fff0000;	/* scale x to [1,2] */
16025c28e83SPiotr Jasiukajtis 	px[i0] = ix;
16125c28e83SPiotr Jasiukajtis 	i = ix + 0x200;
16225c28e83SPiotr Jasiukajtis 	pz[i0] = i & 0xfffffc00;
16325c28e83SPiotr Jasiukajtis 	pz[i1] = pz[1] = pz[2] = 0;
16425c28e83SPiotr Jasiukajtis 	s = (x - z) / (x + z);
16525c28e83SPiotr Jasiukajtis 	j = (i >> 10) & 0x3f;
16625c28e83SPiotr Jasiukajtis 	z = s * s;
16725c28e83SPiotr Jasiukajtis 	qn = (long double) n;
16825c28e83SPiotr Jasiukajtis 	t = qn * ln2lo + _TBL_logl_lo[j];
16925c28e83SPiotr Jasiukajtis 	h = qn * ln2hi + _TBL_logl_hi[j];
17025c28e83SPiotr Jasiukajtis 	f = t + s * (A1 + z * (A2 + z * (A3 + z * (A4 + z * (A5 +
17125c28e83SPiotr Jasiukajtis 		z * (A6 + z * A7))))));
17225c28e83SPiotr Jasiukajtis 	return (h + f);
17325c28e83SPiotr Jasiukajtis }
174