/*
* CDDL HEADER START
*
* The contents of this file are subject to the terms of the
* Common Development and Distribution License (the "License").
* You may not use this file except in compliance with the License.
*
* You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
* See the License for the specific language governing permissions
* and limitations under the License.
*
* When distributing Covered Code, include this CDDL HEADER in each
* file and include the License file at usr/src/OPENSOLARIS.LICENSE.
* If applicable, add the following below this CDDL HEADER, with the
* fields enclosed by brackets "[]" replaced with your own identifying
* information: Portions Copyright [yyyy] [name of copyright owner]
*
* CDDL HEADER END
*/
/*
* Copyright 2011 Nexenta Systems, Inc. All rights reserved.
*/
/*
* Copyright 2006 Sun Microsystems, Inc. All rights reserved.
* Use is subject to license terms.
*/
/*
* logl(x)
* Table look-up algorithm
* By K.C. Ng, March 6, 1989
*
* (a). For x in [31/33,33/31], using a special approximation:
* f = x - 1;
* s = f/(2.0+f); ... here |s| <= 0.03125
* z = s*s;
* return f-s*(f-z*(B1+z*(B2+z*(B3+z*(B4+...+z*B9)...))));
*
* (b). Otherwise, normalize x = 2^n * 1.f.
* Use a 6-bit table look-up: find a 6 bit g that match f to 6.5 bits,
* then
* log(x) = n*ln2 + log(1.g) + log(1.f/1.g).
* Here the leading and trailing values of log(1.g) are obtained from
* a size-64 table.
* For log(1.f/1.g), let s = (1.f-1.g)/(1.f+1.g), then
* log(1.f/1.g) = log((1+s)/(1-s)) = 2s + 2/3 s^3 + 2/5 s^5 +...
* Note that |s|<2**-8=0.00390625. We use an odd s-polynomial
* approximation to compute log(1.f/1.g):
* s*(A1+s^2*(A2+s^2*(A3+s^2*(A4+s^2*(A5+s^2*(A6+s^2*A7))))))
* (Precision is 2**-136.91 bits, absolute error)
*
* (c). The final result is computed by
* (n*ln2_hi+_TBL_logl_hi[j]) +
* ( (n*ln2_lo+_TBL_logl_lo[j]) + s*(A1+...) )
*
* Note.
* For ln2_hi and _TBL_logl_hi[j], we force their last 32 bit to be zero
* so that n*ln2_hi + _TBL_logl_hi[j] is exact. Here
* _TBL_logl_hi[j] + _TBL_logl_lo[j] match log(1+j*2**-6) to 194 bits
*
*
* Special cases:
* log(x) is NaN with signal if x < 0 (including -INF) ;
* log(+INF) is +INF; log(0) is -INF with signal;
* log(NaN) is that NaN with no signal.
*
* Constants:
* The hexadecimal values are the intended ones for the following constants.
* The decimal values may be used, provided that the compiler will convert
* from decimal to binary accurately enough to produce the hexadecimal values
* shown.
*/
#include "libm.h"
extern const long double _TBL_logl_hi[], _TBL_logl_lo[];
static const long double
long double
logl(long double x) {
long double f, s, z, qn, h, t;
int *px = (int *) &x;
int *pz = (int *) &z;
/* get long double precision word ordering */
if (*(int *) &one == 0) {
i0 = 3;
i1 = 0;
} else {
i0 = 0;
i1 = 3;
}
n = 0;
f = x - one;
z = f * f;
return (zero); /* log(1)= +0 */
}
s = f / (two + f); /* |s|<2**-8 */
z = s * s;
}
if (ix >= 0x7fff0000)
return (x + x); /* x is +inf or NaN */
goto LARGE_N;
}
if (ix >= 0x00010000)
goto LARGE_N;
i = ix & 0x7fffffff;
return (one / x); /* log(0.0) = -inf */
}
if (ix < 0) {
if ((unsigned) ix >= 0xffff0000)
return (x - x); /* x is -inf or NaN */
}
/* subnormal x */
x *= two113;
n = -113;
i = ix + 0x200;
s = (x - z) / (x + z);
j = (i >> 10) & 0x3f;
z = s * s;
qn = (long double) n;
return (h + f);
}