[unix-history] / usr / src / lib / libm / common_source / log.c

/*
 * Copyright (c) 1992, 1993
 *	The Regents of the University of California.  All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 * 3. All advertising materials mentioning features or use of this software
 *    must display the following acknowledgement:
 *	This product includes software developed by the University of
 *	California, Berkeley and its contributors.
 * 4. Neither the name of the University nor the names of its contributors
 *    may be used to endorse or promote products derived from this software
 *    without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
 * SUCH DAMAGE.
 */

#ifndef lint
static char sccsid[] = "@(#)log.c	8.2 (Berkeley) 11/30/93";
#endif /* not lint */

#include <math.h>
#include <errno.h>

#include "mathimpl.h"

/* Table-driven natural logarithm.
 *
 * This code was derived, with minor modifications, from:
 *	Peter Tang, "Table-Driven Implementation of the
 *	Logarithm in IEEE Floating-Point arithmetic." ACM Trans.
 *	Math Software, vol 16. no 4, pp 378-400, Dec 1990).
 *
 * Calculates log(2^m*F*(1+f/F)), |f/j| <= 1/256,
 * where F = j/128 for j an integer in [0, 128].
 *
 * log(2^m) = log2_hi*m + log2_tail*m
 * since m is an integer, the dominant term is exact.
 * m has at most 10 digits (for subnormal numbers),
 * and log2_hi has 11 trailing zero bits.
 *
 * log(F) = logF_hi[j] + logF_lo[j] is in tabular form in log_table.h
 * logF_hi[] + 512 is exact.
 *
 * log(1+f/F) = 2*f/(2*F + f) + 1/12 * (2*f/(2*F + f))**3 + ...
 * the leading term is calculated to extra precision in two
 * parts, the larger of which adds exactly to the dominant
 * m and F terms.
 * There are two cases:
 *	1. when m, j are non-zero (m | j), use absolute
 *	   precision for the leading term.
 *	2. when m = j = 0, |1-x| < 1/256, and log(x) ~= (x-1).
 *	   In this case, use a relative precision of 24 bits.
 * (This is done differently in the original paper)
 *
 * Special cases:
 *	0	return signalling -Inf
 *	neg	return signalling NaN
 *	+Inf	return +Inf
*/

#if defined(vax) || defined(tahoe)
#define _IEEE		0
#define TRUNC(x)	x = (double) (float) (x)
#else
#define _IEEE		1
#define endian		(((*(int *) &one)) ? 1 : 0)
#define TRUNC(x)	*(((int *) &x) + endian) &= 0xf8000000
#define infnan(x)	0.0
#endif

#define N 128

/* Table of log(Fj) = logF_head[j] + logF_tail[j], for Fj = 1+j/128.
 * Used for generation of extend precision logarithms.
 * The constant 35184372088832 is 2^45, so the divide is exact.
 * It ensures correct reading of logF_head, even for inaccurate
 * decimal-to-binary conversion routines.  (Everybody gets the
 * right answer for integers less than 2^53.)
 * Values for log(F) were generated using error < 10^-57 absolute
 * with the bc -l package.
*/
static double	A1 = 	  .08333333333333178827;
static double	A2 = 	  .01250000000377174923;
static double	A3 =	 .002232139987919447809;
static double	A4 =	.0004348877777076145742;

static double logF_head[N+1] = {
	0.,
	.007782140442060381246,
	.015504186535963526694,
	.023167059281547608406,
	.030771658666765233647,
	.038318864302141264488,
	.045809536031242714670,
	.053244514518837604555,
	.060624621816486978786,
	.067950661908525944454,
	.075223421237524235039,
	.082443669210988446138,
	.089612158689760690322,
	.096729626458454731618,
	.103796793681567578460,
	.110814366340264314203,
	.117783035656430001836,
	.124703478501032805070,
	.131576357788617315236,
	.138402322859292326029,
	.145182009844575077295,
	.151916042025732167530,
	.158605030176659056451,
	.165249572895390883786,
	.171850256926518341060,
	.178407657472689606947,
	.184922338493834104156,
	.191394852999565046047,
	.197825743329758552135,
	.204215541428766300668,
	.210564769107350002741,
	.216873938300523150246,
	.223143551314024080056,
	.229374101064877322642,
	.235566071312860003672,
	.241719936886966024758,
	.247836163904594286577,
	.253915209980732470285,
	.259957524436686071567,
	.265963548496984003577,
	.271933715484010463114,
	.277868451003087102435,
	.283768173130738432519,
	.289633292582948342896,
	.295464212893421063199,
	.301261330578199704177,
	.307025035294827830512,
	.312755710004239517729,
	.318453731118097493890,
	.324119468654316733591,
	.329753286372579168528,
	.335355541920762334484,
	.340926586970454081892,
	.346466767346100823488,
	.351976423156884266063,
	.357455888922231679316,
	.362905493689140712376,
	.368325561158599157352,
	.373716409793814818840,
	.379078352934811846353,
	.384411698910298582632,
	.389716751140440464951,
	.394993808240542421117,
	.400243164127459749579,
	.405465108107819105498,
	.410659924985338875558,
	.415827895143593195825,
	.420969294644237379543,
	.426084395310681429691,
	.431173464818130014464,
	.436236766774527495726,
	.441274560805140936281,
	.446287102628048160113,
	.451274644139630254358,
	.456237433481874177232,
	.461175715122408291790,
	.466089729924533457960,
	.470979715219073113985,
	.475845904869856894947,
	.480688529345570714212,
	.485507815781602403149,
	.490303988045525329653,
	.495077266798034543171,
	.499827869556611403822,
	.504556010751912253908,
	.509261901790523552335,
	.513945751101346104405,
	.518607764208354637958,
	.523248143765158602036,
	.527867089620485785417,
	.532464798869114019908,
	.537041465897345915436,
	.541597282432121573947,
	.546132437597407260909,
	.550647117952394182793,
	.555141507540611200965,
	.559615787935399566777,
	.564070138285387656651,
	.568504735352689749561,
	.572919753562018740922,
	.577315365035246941260,
	.581691739635061821900,
	.586049045003164792433,
	.590387446602107957005,
	.594707107746216934174,
	.599008189645246602594,
	.603290851438941899687,
	.607555250224322662688,
	.611801541106615331955,
	.616029877215623855590,
	.620240409751204424537,
	.624433288012369303032,
	.628608659422752680256,
	.632766669570628437213,
	.636907462236194987781,
	.641031179420679109171,
	.645137961373620782978,
	.649227946625615004450,
	.653301272011958644725,
	.657358072709030238911,
	.661398482245203922502,
	.665422632544505177065,
	.669430653942981734871,
	.673422675212350441142,
	.677398823590920073911,
	.681359224807238206267,
	.685304003098281100392,
	.689233281238557538017,
	.693147180560117703862
};

static double logF_tail[N+1] = {
	0.,
	-.00000000000000543229938420049,
	 .00000000000000172745674997061,
	-.00000000000001323017818229233,
	-.00000000000001154527628289872,
	-.00000000000000466529469958300,
	 .00000000000005148849572685810,
	-.00000000000002532168943117445,
	-.00000000000005213620639136504,
	-.00000000000001819506003016881,
	 .00000000000006329065958724544,
	 .00000000000008614512936087814,
	-.00000000000007355770219435028,
	 .00000000000009638067658552277,
	 .00000000000007598636597194141,
	 .00000000000002579999128306990,
	-.00000000000004654729747598444,
	-.00000000000007556920687451336,
	 .00000000000010195735223708472,
	-.00000000000017319034406422306,
	-.00000000000007718001336828098,
	 .00000000000010980754099855238,
	-.00000000000002047235780046195,
	-.00000000000008372091099235912,
	 .00000000000014088127937111135,
	 .00000000000012869017157588257,
	 .00000000000017788850778198106,
	 .00000000000006440856150696891,
	 .00000000000016132822667240822,
	-.00000000000007540916511956188,
	-.00000000000000036507188831790,
	 .00000000000009120937249914984,
	 .00000000000018567570959796010,
	-.00000000000003149265065191483,
	-.00000000000009309459495196889,
	 .00000000000017914338601329117,
	-.00000000000001302979717330866,
	 .00000000000023097385217586939,
	 .00000000000023999540484211737,
	 .00000000000015393776174455408,
	-.00000000000036870428315837678,
	 .00000000000036920375082080089,
	-.00000000000009383417223663699,
	 .00000000000009433398189512690,
	 .00000000000041481318704258568,
	-.00000000000003792316480209314,
	 .00000000000008403156304792424,
	-.00000000000034262934348285429,
	 .00000000000043712191957429145,
	-.00000000000010475750058776541,
	-.00000000000011118671389559323,
	 .00000000000037549577257259853,
	 .00000000000013912841212197565,
	 .00000000000010775743037572640,
	 .00000000000029391859187648000,
	-.00000000000042790509060060774,
	 .00000000000022774076114039555,
	 .00000000000010849569622967912,
	-.00000000000023073801945705758,
	 .00000000000015761203773969435,
	 .00000000000003345710269544082,
	-.00000000000041525158063436123,
	 .00000000000032655698896907146,
	-.00000000000044704265010452446,
	 .00000000000034527647952039772,
	-.00000000000007048962392109746,
	 .00000000000011776978751369214,
	-.00000000000010774341461609578,
	 .00000000000021863343293215910,
	 .00000000000024132639491333131,
	 .00000000000039057462209830700,
	-.00000000000026570679203560751,
	 .00000000000037135141919592021,
	-.00000000000017166921336082431,
	-.00000000000028658285157914353,
	-.00000000000023812542263446809,
	 .00000000000006576659768580062,
	-.00000000000028210143846181267,
	 .00000000000010701931762114254,
	 .00000000000018119346366441110,
	 .00000000000009840465278232627,
	-.00000000000033149150282752542,
	-.00000000000018302857356041668,
	-.00000000000016207400156744949,
	 .00000000000048303314949553201,
	-.00000000000071560553172382115,
	 .00000000000088821239518571855,
	-.00000000000030900580513238244,
	-.00000000000061076551972851496,
	 .00000000000035659969663347830,
	 .00000000000035782396591276383,
	-.00000000000046226087001544578,
	 .00000000000062279762917225156,
	 .00000000000072838947272065741,
	 .00000000000026809646615211673,
	-.00000000000010960825046059278,
	 .00000000000002311949383800537,
	-.00000000000058469058005299247,
	-.00000000000002103748251144494,
	-.00000000000023323182945587408,
	-.00000000000042333694288141916,
	-.00000000000043933937969737844,
	 .00000000000041341647073835565,
	 .00000000000006841763641591466,
	 .00000000000047585534004430641,
	 .00000000000083679678674757695,
	-.00000000000085763734646658640,
	 .00000000000021913281229340092,
	-.00000000000062242842536431148,
	-.00000000000010983594325438430,
	 .00000000000065310431377633651,
	-.00000000000047580199021710769,
	-.00000000000037854251265457040,
	 .00000000000040939233218678664,
	 .00000000000087424383914858291,
	 .00000000000025218188456842882,
	-.00000000000003608131360422557,
	-.00000000000050518555924280902,
	 .00000000000078699403323355317,
	-.00000000000067020876961949060,
	 .00000000000016108575753932458,
	 .00000000000058527188436251509,
	-.00000000000035246757297904791,
	-.00000000000018372084495629058,
	 .00000000000088606689813494916,
	 .00000000000066486268071468700,
	 .00000000000063831615170646519,
	 .00000000000025144230728376072,
	-.00000000000017239444525614834
};

double
#ifdef _ANSI_SOURCE
log(double x)
#else
log(x) double x;
#endif
{
	int m, j;
	double F, f, g, q, u, u2, v, zero = 0.0, one = 1.0;
	volatile double u1;

	/* Catch special cases */
	if (x <= 0)
		if (_IEEE && x == zero)	/* log(0) = -Inf */
			return (-one/zero);
		else if (_IEEE)		/* log(neg) = NaN */
			return (zero/zero);
		else if (x == zero)	/* NOT REACHED IF _IEEE */
			return (infnan(-ERANGE));
		else
			return (infnan(EDOM));
	else if (!finite(x))
		if (_IEEE)		/* x = NaN, Inf */
			return (x+x);
		else
			return (infnan(ERANGE));
	
	/* Argument reduction: 1 <= g < 2; x/2^m = g;	*/
	/* y = F*(1 + f/F) for |f| <= 2^-8		*/

	m = logb(x);
	g = ldexp(x, -m);
	if (_IEEE && m == -1022) {
		j = logb(g), m += j;
		g = ldexp(g, -j);
	}
	j = N*(g-1) + .5;
	F = (1.0/N) * j + 1;	/* F*128 is an integer in [128, 512] */
	f = g - F;

	/* Approximate expansion for log(1+f/F) ~= u + q */
	g = 1/(2*F+f);
	u = 2*f*g;
	v = u*u;
	q = u*v*(A1 + v*(A2 + v*(A3 + v*A4)));

    /* case 1: u1 = u rounded to 2^-43 absolute.  Since u < 2^-8,
     * 	       u1 has at most 35 bits, and F*u1 is exact, as F has < 8 bits.
     *         It also adds exactly to |m*log2_hi + log_F_head[j] | < 750
    */
	if (m | j)
		u1 = u + 513, u1 -= 513;

    /* case 2:	|1-x| < 1/256. The m- and j- dependent terms are zero;
     * 		u1 = u to 24 bits.
    */
	else
		u1 = u, TRUNC(u1);
	u2 = (2.0*(f - F*u1) - u1*f) * g;
			/* u1 + u2 = 2f/(2F+f) to extra precision.	*/

	/* log(x) = log(2^m*F*(1+f/F)) =				*/
	/* (m*log2_hi+logF_head[j]+u1) + (m*log2_lo+logF_tail[j]+q);	*/
	/* (exact) + (tiny)						*/

	u1 += m*logF_head[N] + logF_head[j];		/* exact */
	u2 = (u2 + logF_tail[j]) + q;			/* tiny */
	u2 += logF_tail[N]*m;
	return (u1 + u2);
}

/*
 * Extra precision variant, returning struct {double a, b;};
 * log(x) = a+b to 63 bits, with a is rounded to 26 bits.
 */
struct Double
#ifdef _ANSI_SOURCE
__log__D(double x)
#else
__log__D(x) double x;
#endif
{
	int m, j;
	double F, f, g, q, u, v, u2, one = 1.0;
	volatile double u1;
	struct Double r;

	/* Argument reduction: 1 <= g < 2; x/2^m = g;	*/
	/* y = F*(1 + f/F) for |f| <= 2^-8		*/

	m = logb(x);
	g = ldexp(x, -m);
	if (_IEEE && m == -1022) {
		j = logb(g), m += j;
		g = ldexp(g, -j);
	}
	j = N*(g-1) + .5;
	F = (1.0/N) * j + 1;
	f = g - F;

	g = 1/(2*F+f);
	u = 2*f*g;
	v = u*u;
	q = u*v*(A1 + v*(A2 + v*(A3 + v*A4)));
	if (m | j)
		u1 = u + 513, u1 -= 513;
	else
		u1 = u, TRUNC(u1);
	u2 = (2.0*(f - F*u1) - u1*f) * g;

	u1 += m*logF_head[N] + logF_head[j];

	u2 +=  logF_tail[j]; u2 += q;
	u2 += logF_tail[N]*m;
	r.a = u1 + u2;			/* Only difference is here */
	TRUNC(r.a);
	r.b = (u1 - r.a) + u2;
	return (r);
}
Commit	Line	Data
9b525f39	1	/*
ff0114f0 KB	2	* Copyright (c) 1992, 1993
ff0114f0 KB	3	* The Regents of the University of California. All rights reserved.
9b525f39	4	*
ad787160 C	5	* Redistribution and use in source and binary forms, with or without
	6	* modification, are permitted provided that the following conditions
	7	* are met:
	8	* 1. Redistributions of source code must retain the above copyright
	9	* notice, this list of conditions and the following disclaimer.
	10	* 2. Redistributions in binary form must reproduce the above copyright
	11	* notice, this list of conditions and the following disclaimer in the
	12	* documentation and/or other materials provided with the distribution.
	13	* 3. All advertising materials mentioning features or use of this software
	14	* must display the following acknowledgement:
	15	* This product includes software developed by the University of
	16	* California, Berkeley and its contributors.
	17	* 4. Neither the name of the University nor the names of its contributors
	18	* may be used to endorse or promote products derived from this software
	19	* without specific prior written permission.
9b525f39	20	*
ad787160 C	21	* THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
	22	* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
	23	* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
	24	* ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
	25	* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
	26	* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
	27	* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
	28	* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
	29	* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
	30	* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
	31	* SUCH DAMAGE.
59f3cb20 ZAL	32	*/
	33
	34	#ifndef lint
ed554bc5	35	static char sccsid[] = "@(#)log.c 8.2 (Berkeley) 11/30/93";
9b525f39	36	#endif /* not lint */
59f3cb20	37
3f1b1a1e KB	38	#include <math.h>
	39	#include <errno.h>
	40
a12850a3	41	#include "mathimpl.h"
3f1b1a1e KB	42
3f1b1a1e KB	43	/* Table-driven natural logarithm.
59f3cb20	44	*
3f1b1a1e KB	45	* This code was derived, with minor modifications, from:
	46	* Peter Tang, "Table-Driven Implementation of the
	47	* Logarithm in IEEE Floating-Point arithmetic." ACM Trans.
	48	* Math Software, vol 16. no 4, pp 378-400, Dec 1990).
59f3cb20	49	*
3f1b1a1e KB	50	* Calculates log(2^mF(1+f/F)), \|f/j\| <= 1/256,
3f1b1a1e KB	51	* where F = j/128 for j an integer in [0, 128].
59f3cb20	52	*
3f1b1a1e KB	53	* log(2^m) = log2_him + log2_tailm
	54	* since m is an integer, the dominant term is exact.
	55	* m has at most 10 digits (for subnormal numbers),
	56	* and log2_hi has 11 trailing zero bits.
59f3cb20	57	*
3f1b1a1e KB	58	* log(F) = logF_hi[j] + logF_lo[j] is in tabular form in log_table.h
3f1b1a1e KB	59	* logF_hi[] + 512 is exact.
59f3cb20	60	*
3f1b1a1e KB	61	* log(1+f/F) = 2f/(2F + f) + 1/12 * (2f/(2F + f))**3 + ...
	62	* the leading term is calculated to extra precision in two
	63	* parts, the larger of which adds exactly to the dominant
	64	* m and F terms.
	65	* There are two cases:
	66	* 1. when m, j are non-zero (m \| j), use absolute
	67	* precision for the leading term.
	68	* 2. when m = j = 0, \|1-x\| < 1/256, and log(x) ~= (x-1).
	69	* In this case, use a relative precision of 24 bits.
	70	* (This is done differently in the original paper)
59f3cb20 ZAL	71	*
59f3cb20 ZAL	72	* Special cases:
3f1b1a1e KB	73	* 0 return signalling -Inf
	74	* neg return signalling NaN
	75	* +Inf return +Inf
	76	*/
	77
	78	#if defined(vax) \|\| defined(tahoe)
df4f4b2e PM	79	#define _IEEE 0
df4f4b2e PM	80	#define TRUNC(x) x = (double) (float) (x)
3f1b1a1e	81	#else
75148696	82	#define _IEEE 1
df4f4b2e PM	83	#define endian ((((int ) &one)) ? 1 : 0)
	84	#define TRUNC(x) (((int ) &x) + endian) &= 0xf8000000
	85	#define infnan(x) 0.0
9eda3584 KB	86	#endif
9eda3584 KB	87
75148696 KB	88	#define N 128
	89
	90	/* Table of log(Fj) = logF_head[j] + logF_tail[j], for Fj = 1+j/128.
	91	* Used for generation of extend precision logarithms.
	92	* The constant 35184372088832 is 2^45, so the divide is exact.
	93	* It ensures correct reading of logF_head, even for inaccurate
	94	* decimal-to-binary conversion routines. (Everybody gets the
	95	* right answer for integers less than 2^53.)
	96	* Values for log(F) were generated using error < 10^-57 absolute
	97	* with the bc -l package.
	98	*/
7af77d91 KB	99	static double A1 = .08333333333333178827;
	100	static double A2 = .01250000000377174923;
	101	static double A3 = .002232139987919447809;
	102	static double A4 = .0004348877777076145742;
75148696	103
7af77d91	104	static double logF_head[N+1] = {
75148696 KB	105	0.,
	106	.007782140442060381246,
	107	.015504186535963526694,
	108	.023167059281547608406,
	109	.030771658666765233647,
	110	.038318864302141264488,
	111	.045809536031242714670,
	112	.053244514518837604555,
	113	.060624621816486978786,
	114	.067950661908525944454,
	115	.075223421237524235039,
	116	.082443669210988446138,
	117	.089612158689760690322,
	118	.096729626458454731618,
	119	.103796793681567578460,
	120	.110814366340264314203,
	121	.117783035656430001836,
	122	.124703478501032805070,
	123	.131576357788617315236,
	124	.138402322859292326029,
	125	.145182009844575077295,
	126	.151916042025732167530,
	127	.158605030176659056451,
	128	.165249572895390883786,
	129	.171850256926518341060,
	130	.178407657472689606947,
	131	.184922338493834104156,
	132	.191394852999565046047,
	133	.197825743329758552135,
	134	.204215541428766300668,
	135	.210564769107350002741,
	136	.216873938300523150246,
	137	.223143551314024080056,
	138	.229374101064877322642,
	139	.235566071312860003672,
	140	.241719936886966024758,
	141	.247836163904594286577,
	142	.253915209980732470285,
	143	.259957524436686071567,
	144	.265963548496984003577,
	145	.271933715484010463114,
	146	.277868451003087102435,
	147	.283768173130738432519,
	148	.289633292582948342896,
	149	.295464212893421063199,
	150	.301261330578199704177,
	151	.307025035294827830512,
	152	.312755710004239517729,
	153	.318453731118097493890,
	154	.324119468654316733591,
	155	.329753286372579168528,
	156	.335355541920762334484,
	157	.340926586970454081892,
	158	.346466767346100823488,
	159	.351976423156884266063,
	160	.357455888922231679316,
	161	.362905493689140712376,
	162	.368325561158599157352,
	163	.373716409793814818840,
	164	.379078352934811846353,
	165	.384411698910298582632,
	166	.389716751140440464951,
	167	.394993808240542421117,
	168	.400243164127459749579,
169	.405465108107819105498,
170	.410659924985338875558,
171	.415827895143593195825,
172	.420969294644237379543,
173	.426084395310681429691,
174	.431173464818130014464,
175	.436236766774527495726,
176	.441274560805140936281,
177	.446287102628048160113,
178	.451274644139630254358,
179	.456237433481874177232,
180	.461175715122408291790,
181	.466089729924533457960,
182	.470979715219073113985,
183	.475845904869856894947,
184	.480688529345570714212,
185	.485507815781602403149,
186	.490303988045525329653,
187	.495077266798034543171,
188	.499827869556611403822,
189	.504556010751912253908,
190	.509261901790523552335,
191	.513945751101346104405,
192	.518607764208354637958,
193	.523248143765158602036,
194	.527867089620485785417,
195	.532464798869114019908,
196	.537041465897345915436,
197	.541597282432121573947,
198	.546132437597407260909,
199	.550647117952394182793,
200	.555141507540611200965,
201	.559615787935399566777,
202	.564070138285387656651,
203	.568504735352689749561,
204	.572919753562018740922,
205	.577315365035246941260,
206	.581691739635061821900,
207	.586049045003164792433,
208	.590387446602107957005,
209	.594707107746216934174,
210	.599008189645246602594,
211	.603290851438941899687,
212	.607555250224322662688,
213	.611801541106615331955,
214	.616029877215623855590,
215	.620240409751204424537,
216	.624433288012369303032,
217	.628608659422752680256,
218	.632766669570628437213,
219	.636907462236194987781,
220	.641031179420679109171,
221	.645137961373620782978,
222	.649227946625615004450,
223	.653301272011958644725,
224	.657358072709030238911,
225	.661398482245203922502,
226	.665422632544505177065,
227	.669430653942981734871,
228	.673422675212350441142,
229	.677398823590920073911,
230	.681359224807238206267,
231	.685304003098281100392,
232	.689233281238557538017,
233	.693147180560117703862
234	};
235
7af77d91	236	static double logF_tail[N+1] = {
75148696 KB	237	0.,
	238	-.00000000000000543229938420049,
	239	.00000000000000172745674997061,
	240	-.00000000000001323017818229233,
	241	-.00000000000001154527628289872,
	242	-.00000000000000466529469958300,
	243	.00000000000005148849572685810,
	244	-.00000000000002532168943117445,
	245	-.00000000000005213620639136504,
	246	-.00000000000001819506003016881,
	247	.00000000000006329065958724544,
	248	.00000000000008614512936087814,
	249	-.00000000000007355770219435028,
	250	.00000000000009638067658552277,
	251	.00000000000007598636597194141,
	252	.00000000000002579999128306990,
	253	-.00000000000004654729747598444,
	254	-.00000000000007556920687451336,
	255	.00000000000010195735223708472,
	256	-.00000000000017319034406422306,
	257	-.00000000000007718001336828098,
	258	.00000000000010980754099855238,
	259	-.00000000000002047235780046195,
	260	-.00000000000008372091099235912,
	261	.00000000000014088127937111135,
	262	.00000000000012869017157588257,
	263	.00000000000017788850778198106,
	264	.00000000000006440856150696891,
	265	.00000000000016132822667240822,
	266	-.00000000000007540916511956188,
	267	-.00000000000000036507188831790,
	268	.00000000000009120937249914984,
	269	.00000000000018567570959796010,
	270	-.00000000000003149265065191483,
	271	-.00000000000009309459495196889,
	272	.00000000000017914338601329117,
	273	-.00000000000001302979717330866,
	274	.00000000000023097385217586939,
	275	.00000000000023999540484211737,
	276	.00000000000015393776174455408,
	277	-.00000000000036870428315837678,
	278	.00000000000036920375082080089,
	279	-.00000000000009383417223663699,
	280	.00000000000009433398189512690,
	281	.00000000000041481318704258568,
	282	-.00000000000003792316480209314,
	283	.00000000000008403156304792424,
	284	-.00000000000034262934348285429,
	285	.00000000000043712191957429145,
	286	-.00000000000010475750058776541,
	287	-.00000000000011118671389559323,
	288	.00000000000037549577257259853,
	289	.00000000000013912841212197565,
	290	.00000000000010775743037572640,
	291	.00000000000029391859187648000,
	292	-.00000000000042790509060060774,
	293	.00000000000022774076114039555,
	294	.00000000000010849569622967912,
	295	-.00000000000023073801945705758,
	296	.00000000000015761203773969435,
	297	.00000000000003345710269544082,
	298	-.00000000000041525158063436123,
	299	.00000000000032655698896907146,
	300	-.00000000000044704265010452446,
301	.00000000000034527647952039772,
302	-.00000000000007048962392109746,
303	.00000000000011776978751369214,
304	-.00000000000010774341461609578,
305	.00000000000021863343293215910,
306	.00000000000024132639491333131,
307	.00000000000039057462209830700,
308	-.00000000000026570679203560751,
309	.00000000000037135141919592021,
310	-.00000000000017166921336082431,
311	-.00000000000028658285157914353,
312	-.00000000000023812542263446809,
313	.00000000000006576659768580062,
314	-.00000000000028210143846181267,
315	.00000000000010701931762114254,
316	.00000000000018119346366441110,
317	.00000000000009840465278232627,
318	-.00000000000033149150282752542,
319	-.00000000000018302857356041668,
320	-.00000000000016207400156744949,
321	.00000000000048303314949553201,
322	-.00000000000071560553172382115,
323	.00000000000088821239518571855,
324	-.00000000000030900580513238244,
325	-.00000000000061076551972851496,
326	.00000000000035659969663347830,
327	.00000000000035782396591276383,
328	-.00000000000046226087001544578,
329	.00000000000062279762917225156,
330	.00000000000072838947272065741,
331	.00000000000026809646615211673,
332	-.00000000000010960825046059278,
333	.00000000000002311949383800537,
334	-.00000000000058469058005299247,
335	-.00000000000002103748251144494,
336	-.00000000000023323182945587408,
337	-.00000000000042333694288141916,
338	-.00000000000043933937969737844,
339	.00000000000041341647073835565,
340	.00000000000006841763641591466,
341	.00000000000047585534004430641,
342	.00000000000083679678674757695,
343	-.00000000000085763734646658640,
344	.00000000000021913281229340092,
345	-.00000000000062242842536431148,
346	-.00000000000010983594325438430,
347	.00000000000065310431377633651,
348	-.00000000000047580199021710769,
349	-.00000000000037854251265457040,
350	.00000000000040939233218678664,
351	.00000000000087424383914858291,
352	.00000000000025218188456842882,
353	-.00000000000003608131360422557,
354	-.00000000000050518555924280902,
355	.00000000000078699403323355317,
356	-.00000000000067020876961949060,
357	.00000000000016108575753932458,
358	.00000000000058527188436251509,
359	-.00000000000035246757297904791,
360	-.00000000000018372084495629058,
361	.00000000000088606689813494916,
362	.00000000000066486268071468700,
363	.00000000000063831615170646519,
364	.00000000000025144230728376072,
365	-.00000000000017239444525614834
366	};
367
3f1b1a1e KB	368	double
	369	#ifdef _ANSI_SOURCE
	370	log(double x)
	371	#else
	372	log(x) double x;
	373	#endif
59f3cb20	374	{
3f1b1a1e KB	375	int m, j;
3f1b1a1e KB	376	double F, f, g, q, u, u2, v, zero = 0.0, one = 1.0;
3f1b1a1e KB	377	volatile double u1;
	378
	379	/* Catch special cases */
	380	if (x <= 0)
	381	if (_IEEE && x == zero) /* log(0) = -Inf */
	382	return (-one/zero);
	383	else if (_IEEE) /* log(neg) = NaN */
	384	return (zero/zero);
	385	else if (x == zero) /* NOT REACHED IF _IEEE */
	386	return (infnan(-ERANGE));
59f3cb20	387	else
3f1b1a1e KB	388	return (infnan(EDOM));
	389	else if (!finite(x))
	390	if (_IEEE) /* x = NaN, Inf */
	391	return (x+x);
	392	else
	393	return (infnan(ERANGE));
	394
	395	/* Argument reduction: 1 <= g < 2; x/2^m = g; */
	396	/* y = F(1 + f/F) for \|f\| <= 2^-8 /
	397
	398	m = logb(x);
	399	g = ldexp(x, -m);
	400	if (_IEEE && m == -1022) {
	401	j = logb(g), m += j;
	402	g = ldexp(g, -j);
59f3cb20	403	}
3f1b1a1e KB	404	j = N*(g-1) + .5;
	405	F = (1.0/N) * j + 1; /* F128 is an integer in [128, 512] /
	406	f = g - F;
	407
	408	/* Approximate expansion for log(1+f/F) ~= u + q */
	409	g = 1/(2*F+f);
	410	u = 2fg;
	411	v = u*u;
	412	q = uv(A1 + v(A2 + v(A3 + v*A4)));
	413
	414	/* case 1: u1 = u rounded to 2^-43 absolute. Since u < 2^-8,
	415	* u1 has at most 35 bits, and F*u1 is exact, as F has < 8 bits.
	416	* It also adds exactly to \|m*log2_hi + log_F_head[j] \| < 750
	417	*/
	418	if (m \| j)
	419	u1 = u + 513, u1 -= 513;
	420
	421	/* case 2: \|1-x\| < 1/256. The m- and j- dependent terms are zero;
	422	* u1 = u to 24 bits.
	423	*/
	424	else
	425	u1 = u, TRUNC(u1);
	426	u2 = (2.0(f - Fu1) - u1f) g;
	427	/* u1 + u2 = 2f/(2F+f) to extra precision. */
	428
	429	/* log(x) = log(2^mF(1+f/F)) = */
	430	/* (mlog2_hi+logF_head[j]+u1) + (mlog2_lo+logF_tail[j]+q); */
	431	/* (exact) + (tiny) */
	432
	433	u1 += mlogF_head[N] + logF_head[j]; / exact */
	434	u2 = (u2 + logF_tail[j]) + q; /* tiny */
	435	u2 += logF_tail[N]*m;
	436	return (u1 + u2);
	437	}
59f3cb20	438
75148696 KB	439	/*
	440	* Extra precision variant, returning struct {double a, b;};
	441	* log(x) = a+b to 63 bits, with a is rounded to 26 bits.
3f1b1a1e KB	442	*/
	443	struct Double
	444	#ifdef _ANSI_SOURCE
5acba3ee	445	__log__D(double x)
3f1b1a1e	446	#else
5acba3ee	447	__log__D(x) double x;
3f1b1a1e KB	448	#endif
	449	{
	450	int m, j;
75148696	451	double F, f, g, q, u, v, u2, one = 1.0;
3f1b1a1e KB	452	volatile double u1;
	453	struct Double r;
	454
	455	/* Argument reduction: 1 <= g < 2; x/2^m = g; */
	456	/* y = F(1 + f/F) for \|f\| <= 2^-8 /
	457
	458	m = logb(x);
	459	g = ldexp(x, -m);
	460	if (_IEEE && m == -1022) {
	461	j = logb(g), m += j;
	462	g = ldexp(g, -j);
	463	}
	464	j = N*(g-1) + .5;
	465	F = (1.0/N) * j + 1;
	466	f = g - F;
	467
	468	g = 1/(2*F+f);
	469	u = 2fg;
	470	v = u*u;
	471	q = uv(A1 + v(A2 + v(A3 + v*A4)));
	472	if (m \| j)
	473	u1 = u + 513, u1 -= 513;
	474	else
	475	u1 = u, TRUNC(u1);
	476	u2 = (2.0(f - Fu1) - u1f) g;
	477
	478	u1 += m*logF_head[N] + logF_head[j];
	479
	480	u2 += logF_tail[j]; u2 += q;
	481	u2 += logF_tail[N]*m;
	482	r.a = u1 + u2; /* Only difference is here */
	483	TRUNC(r.a);
	484	r.b = (u1 - r.a) + u2;
	485	return (r);
59f3cb20	486	}