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

/*-
 * Copyright (c) 1992 The Regents of the University of California.
 * All rights reserved.
 *
 * %sccs.include.redist.c%
 */

#ifndef lint
static char sccsid[] = "@(#)lgamma.c	5.11 (Berkeley) %G%";
#endif /* not lint */

/*
 * Coded by Peter McIlroy, Nov 1992;
 *
 * The financial support of UUNET Communications Services is greatfully
 * acknowledged.
 */

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

#include "mathimpl.h"

/* Log gamma function.
 * Error:  x > 0 error < 1.3ulp.
 *	   x > 4, error < 1ulp.
 *	   x > 9, error < .6ulp.
 * 	   x < 0, all bets are off. (When G(x) ~ 1, log(G(x)) ~ 0)
 * Method:
 *	x > 6:
 *		Use the asymptotic expansion (Stirling's Formula)
 *	0 < x < 6:
 *		Use gamma(x+1) = x*gamma(x) for argument reduction.
 *		Use rational approximation in
 *		the range 1.2, 2.5
 *		Two approximations are used, one centered at the
 *		minimum to ensure monotonicity; one centered at 2
 *		to maintain small relative error.
 *	x < 0:
 *		Use the reflection formula,
 *		G(1-x)G(x) = PI/sin(PI*x)
 * Special values:
 *	non-positive integer	returns +Inf.
 *	NaN			returns NaN
*/
static int endian;
#if defined(vax) || defined(tahoe)
#define _IEEE		0
/* double and float have same size exponent field */
#define TRUNC(x)	x = (double) (float) (x)
#else
#define _IEEE		1
#define TRUNC(x)	*(((int *) &x) + endian) &= 0xf8000000
#define infnan(x)	0.0
#endif

extern double log1p(double);
static double small_lgam(double);
static double large_lgam(double);
static double neg_lgam(double);
static double zero = 0.0, one = 1.0;
int signgam;

#define UNDERFL (1e-1020 * 1e-1020)

#define LEFT	(1.0 - (x0 + .25))
#define RIGHT	(x0 - .218)
/*
/* Constants for approximation in [1.244,1.712]
*/
#define x0	0.461632144968362356785
#define x0_lo	-.000000000000000015522348162858676890521
#define a0_hi	-0.12148629128932952880859
#define a0_lo	.0000000007534799204229502
#define r0	-2.771227512955130520e-002
#define r1	-2.980729795228150847e-001
#define r2	-3.257411333183093394e-001
#define r3	-1.126814387531706041e-001
#define r4	-1.129130057170225562e-002
#define r5	-2.259650588213369095e-005
#define s0	 1.714457160001714442e+000
#define s1	 2.786469504618194648e+000
#define s2	 1.564546365519179805e+000
#define s3	 3.485846389981109850e-001
#define s4	 2.467759345363656348e-002
/*
 * Constants for approximation in [1.71, 2.5]
*/
#define a1_hi	4.227843350984671344505727574870e-01
#define a1_lo	4.670126436531227189e-18
#define p0	3.224670334241133695662995251041e-01
#define p1	3.569659696950364669021382724168e-01
#define p2	1.342918716072560025853732668111e-01
#define p3	1.950702176409779831089963408886e-02
#define p4	8.546740251667538090796227834289e-04
#define q0	1.000000000000000444089209850062e+00
#define q1	1.315850076960161985084596381057e+00
#define q2	6.274644311862156431658377186977e-01
#define q3	1.304706631926259297049597307705e-01
#define q4	1.102815279606722369265536798366e-02
#define q5	2.512690594856678929537585620579e-04
#define q6	-1.003597548112371003358107325598e-06
/*
 * Stirling's Formula, adjusted for equal-ripple. x in [6,Inf].
*/
#define lns2pi	.418938533204672741780329736405
#define pb0	 8.33333333333333148296162562474e-02
#define pb1	-2.77777777774548123579378966497e-03
#define pb2	 7.93650778754435631476282786423e-04
#define pb3	-5.95235082566672847950717262222e-04
#define pb4	 8.41428560346653702135821806252e-04
#define pb5	-1.89773526463879200348872089421e-03
#define pb6	 5.69394463439411649408050664078e-03
#define pb7	-1.44705562421428915453880392761e-02

double
lgamma(double x)
{
	double r;

	signgam = 1;
	endian = ((*(int *) &one)) ? 1 : 0;

	if (!finite(x))
		if (_IEEE)
			return (x+x);
		else return (infnan(EDOM));

	if (x > 6 + RIGHT) {
		r = large_lgam(x);
		return (r);
	} else if (x > 1e-16)
		return (small_lgam(x));
	else if (x > -1e-16) {
		if (x < 0)
			signgam = -1, x = -x;
		return (-log(x));
	} else
		return (neg_lgam(x));
}

static double
large_lgam(double x)
{
	double z, p, x1;
	int i;
	struct Double t, u, v;
	u = log__D(x);
	u.a -= 1.0;
	if (x > 1e15) {
		v.a = x - 0.5;
		TRUNC(v.a);
		v.b = (x - v.a) - 0.5;
		t.a = u.a*v.a;
		t.b = x*u.b + v.b*u.a;
		if (_IEEE == 0 && !finite(t.a))
			return(infnan(ERANGE));
		return(t.a + t.b);
	}
	x1 = 1./x;
	z = x1*x1;
	p = pb0+z*(pb1+z*(pb2+z*(pb3+z*(pb4+z*(pb5+z*(pb6+z*pb7))))));
					/* error in approximation = 2.8e-19 */

	p = p*x1;			/* error < 2.3e-18 absolute */
					/* 0 < p < 1/64 (at x = 5.5) */
	v.a = x = x - 0.5;
	TRUNC(v.a);			/* truncate v.a to 26 bits. */
	v.b = x - v.a;
	t.a = v.a*u.a;			/* t = (x-.5)*(log(x)-1) */
	t.b = v.b*u.a + x*u.b;
	t.b += p; t.b += lns2pi;	/* return t + lns2pi + p */
	return (t.a + t.b);
}

static double
small_lgam(double x)
{
	int x_int;
	double y, z, t, r = 0, p, q, hi, lo;
	struct Double rr;
	x_int = (x + .5);
	y = x - x_int;
	if (x_int <= 2 && y > RIGHT) {
		t = y - x0;
		y--; x_int++;
		goto CONTINUE;
	} else if (y < -LEFT) {
		t = y +(1.0-x0);
CONTINUE:
		z = t - x0_lo;
		p = r0+z*(r1+z*(r2+z*(r3+z*(r4+z*r5))));
		q = s0+z*(s1+z*(s2+z*(s3+z*s4)));
		r = t*(z*(p/q) - x0_lo);
		t = .5*t*t;
		z = 1.0;
		switch (x_int) {
		case 6:	z  = (y + 5);
		case 5:	z *= (y + 4);
		case 4:	z *= (y + 3);
		case 3:	z *= (y + 2);
			rr = log__D(z);
			rr.b += a0_lo; rr.a += a0_hi;
			return(((r+rr.b)+t+rr.a));
		case 2: return(((r+a0_lo)+t)+a0_hi);
		case 0: r -= log1p(x);
		default: rr = log__D(x);
			rr.a -= a0_hi; rr.b -= a0_lo;
			return(((r - rr.b) + t) - rr.a);
		}
	} else {
		p = p0+y*(p1+y*(p2+y*(p3+y*p4)));
		q = q0+y*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*q6)))));
		p = p*(y/q);
		t = (double)(float) y;
		z = y-t;
		hi = (double)(float) (p+a1_hi);
		lo = a1_hi - hi; lo += p; lo += a1_lo;
		r = lo*y + z*hi;	/* q + r = y*(a0+p/q) */
		q = hi*t;
		z = 1.0;
		switch (x_int) {
		case 6:	z  = (y + 5);
		case 5:	z *= (y + 4);
		case 4:	z *= (y + 3);
		case 3:	z *= (y + 2);
			rr = log__D(z);
			r += rr.b; r += q;
			return(rr.a + r);
		case 2:	return (q+ r);
		case 0: rr = log__D(x);
			r -= rr.b; r -= log1p(x);
			r += q; r-= rr.a;
			return(r);
		default: rr = log__D(x);
			r -= rr.b;
			q -= rr.a;
			return (r+q);
		}
	}
}

static double
neg_lgam(double x)
{
	int xi;
	double y, z, one = 1.0, zero = 0.0;
	extern double gamma();

	/* avoid destructive cancellation as much as possible */
	if (x > -170) {
		xi = x;
		if (xi == x)
			if (_IEEE)
				return(one/zero);
			else
				return(infnan(ERANGE));
		y = gamma(x);
		if (y < 0)
			y = -y, signgam = -1;
		return (log(y));
	}
	z = floor(x + .5);
	if (z == x) {		/* convention: G(-(integer)) -> +Inf */
		if (_IEEE)
			return (one/zero);
		else
			return (infnan(ERANGE));
	}
	y = .5*ceil(x);
	if (y == ceil(y))
		signgam = -1;
	x = -x;
	z = fabs(x + z);	/* 0 < z <= .5 */
	if (z < .25)
		z = sin(M_PI*z);
	else
		z = cos(M_PI*(0.5-z));
	z = log(M_PI/(z*x));
	y = large_lgam(x);
	return (z - y);
}
Commit	Line	Data
830b4dec	1	/*-
a7784218	2	* Copyright (c) 1992 The Regents of the University of California.
830b4dec KB	3	* All rights reserved.
830b4dec KB	4	*
a7784218	5	* %sccs.include.redist.c%
b6be2d90 KB	6	*/
b6be2d90 KB	7
0ebe4989	8	#ifndef lint
9ab6b304	9	static char sccsid[] = "@(#)lgamma.c 5.11 (Berkeley) %G%";
b6be2d90	10	#endif /* not lint */
0ebe4989	11
19d48448 KB	12	/*
	13	* Coded by Peter McIlroy, Nov 1992;
	14	*
	15	* The financial support of UUNET Communications Services is greatfully
	16	* acknowledged.
	17	*/
	18
05909853 KB	19	#include <math.h>
05909853 KB	20	#include <errno.h>
0ebe4989	21
05909853	22	#include "mathimpl.h"
0ebe4989	23
0662f309 PM	24	/* Log gamma function.
	25	* Error: x > 0 error < 1.3ulp.
	26	* x > 4, error < 1ulp.
	27	* x > 9, error < .6ulp.
19d48448	28	* x < 0, all bets are off. (When G(x) ~ 1, log(G(x)) ~ 0)
0662f309 PM	29	* Method:
	30	* x > 6:
	31	* Use the asymptotic expansion (Stirling's Formula)
	32	* 0 < x < 6:
19d48448	33	* Use gamma(x+1) = x*gamma(x) for argument reduction.
0662f309 PM	34	* Use rational approximation in
0662f309 PM	35	* the range 1.2, 2.5
19d48448 KB	36	* Two approximations are used, one centered at the
	37	* minimum to ensure monotonicity; one centered at 2
	38	* to maintain small relative error.
0662f309 PM	39	* x < 0:
	40	* Use the reflection formula,
	41	* G(1-x)G(x) = PI/sin(PI*x)
	42	* Special values:
	43	* non-positive integer returns +Inf.
	44	* NaN returns NaN
0ebe4989	45	*/
19d48448	46	static int endian;
05909853	47	#if defined(vax) \|\| defined(tahoe)
19d48448	48	#define _IEEE 0
0662f309	49	/* double and float have same size exponent field */
19d48448	50	#define TRUNC(x) x = (double) (float) (x)
05909853	51	#else
19d48448 KB	52	#define _IEEE 1
	53	#define TRUNC(x) (((int ) &x) + endian) &= 0xf8000000
	54	#define infnan(x) 0.0
05909853 KB	55	#endif
05909853 KB	56
0662f309 PM	57	extern double log1p(double);
	58	static double small_lgam(double);
	59	static double large_lgam(double);
	60	static double neg_lgam(double);
	61	static double zero = 0.0, one = 1.0;
05909853 KB	62	int signgam;
05909853 KB	63
0662f309	64	#define UNDERFL (1e-1020 * 1e-1020)
0ebe4989	65
0662f309 PM	66	#define LEFT (1.0 - (x0 + .25))
	67	#define RIGHT (x0 - .218)
	68	/*
	69	/* Constants for approximation in [1.244,1.712]
	70	*/
	71	#define x0 0.461632144968362356785
	72	#define x0_lo -.000000000000000015522348162858676890521
	73	#define a0_hi -0.12148629128932952880859
	74	#define a0_lo .0000000007534799204229502
	75	#define r0 -2.771227512955130520e-002
	76	#define r1 -2.980729795228150847e-001
	77	#define r2 -3.257411333183093394e-001
	78	#define r3 -1.126814387531706041e-001
	79	#define r4 -1.129130057170225562e-002
	80	#define r5 -2.259650588213369095e-005
	81	#define s0 1.714457160001714442e+000
	82	#define s1 2.786469504618194648e+000
	83	#define s2 1.564546365519179805e+000
	84	#define s3 3.485846389981109850e-001
	85	#define s4 2.467759345363656348e-002
	86	/*
	87	* Constants for approximation in [1.71, 2.5]
	88	*/
	89	#define a1_hi 4.227843350984671344505727574870e-01
	90	#define a1_lo 4.670126436531227189e-18
	91	#define p0 3.224670334241133695662995251041e-01
	92	#define p1 3.569659696950364669021382724168e-01
	93	#define p2 1.342918716072560025853732668111e-01
	94	#define p3 1.950702176409779831089963408886e-02
	95	#define p4 8.546740251667538090796227834289e-04
	96	#define q0 1.000000000000000444089209850062e+00
	97	#define q1 1.315850076960161985084596381057e+00
	98	#define q2 6.274644311862156431658377186977e-01
	99	#define q3 1.304706631926259297049597307705e-01
	100	#define q4 1.102815279606722369265536798366e-02
	101	#define q5 2.512690594856678929537585620579e-04
	102	#define q6 -1.003597548112371003358107325598e-06
	103	/*
	104	* Stirling's Formula, adjusted for equal-ripple. x in [6,Inf].
	105	*/
19d48448 KB	106	#define lns2pi .418938533204672741780329736405
	107	#define pb0 8.33333333333333148296162562474e-02
	108	#define pb1 -2.77777777774548123579378966497e-03
	109	#define pb2 7.93650778754435631476282786423e-04
	110	#define pb3 -5.95235082566672847950717262222e-04
	111	#define pb4 8.41428560346653702135821806252e-04
	112	#define pb5 -1.89773526463879200348872089421e-03
	113	#define pb6 5.69394463439411649408050664078e-03
	114	#define pb7 -1.44705562421428915453880392761e-02
9eda3584	115
0ebe4989	116	double
0662f309	117	lgamma(double x)
0ebe4989	118	{
05909853	119	double r;
19d48448	120
05909853	121	signgam = 1;
19d48448 KB	122	endian = (((int ) &one)) ? 1 : 0;
19d48448 KB	123
0662f309 PM	124	if (!finite(x))
	125	if (_IEEE)
	126	return (x+x);
	127	else return (infnan(EDOM));
	128
05909853	129	if (x > 6 + RIGHT) {
05909853 KB	130	r = large_lgam(x);
05909853 KB	131	return (r);
0662f309	132	} else if (x > 1e-16)
05909853	133	return (small_lgam(x));
0662f309	134	else if (x > -1e-16) {
05909853 KB	135	if (x < 0)
	136	signgam = -1, x = -x;
	137	return (-log(x));
	138	} else
	139	return (neg_lgam(x));
0ebe4989 ZAL	140	}
	141
	142	static double
0662f309	143	large_lgam(double x)
0ebe4989	144	{
05909853	145	double z, p, x1;
0ebe4989	146	int i;
05909853	147	struct Double t, u, v;
0662f309 PM	148	u = log__D(x);
	149	u.a -= 1.0;
	150	if (x > 1e15) {
	151	v.a = x - 0.5;
	152	TRUNC(v.a);
	153	v.b = (x - v.a) - 0.5;
	154	t.a = u.a*v.a;
	155	t.b = xu.b + v.bu.a;
	156	if (_IEEE == 0 && !finite(t.a))
	157	return(infnan(ERANGE));
	158	return(t.a + t.b);
	159	}
	160	x1 = 1./x;
	161	z = x1*x1;
	162	p = pb0+z(pb1+z(pb2+z(pb3+z(pb4+z(pb5+z(pb6+z*pb7))))));
	163	/* error in approximation = 2.8e-19 */
05909853	164
0662f309 PM	165	p = px1; / error < 2.3e-18 absolute */
0662f309 PM	166	/* 0 < p < 1/64 (at x = 5.5) */
b26726d3	167	v.a = x = x - 0.5;
0662f309	168	TRUNC(v.a); /* truncate v.a to 26 bits. */
05909853 KB	169	v.b = x - v.a;
	170	t.a = v.au.a; / t = (x-.5)(log(x)-1) /
	171	t.b = v.bu.a + xu.b;
0662f309 PM	172	t.b += p; t.b += lns2pi; /* return t + lns2pi + p */
0662f309 PM	173	return (t.a + t.b);
0ebe4989	174	}
0662f309	175
0ebe4989	176	static double
0662f309	177	small_lgam(double x)
0ebe4989	178	{
0662f309 PM	179	int x_int;
0662f309 PM	180	double y, z, t, r = 0, p, q, hi, lo;
05909853	181	struct Double rr;
0662f309 PM	182	x_int = (x + .5);
	183	y = x - x_int;
	184	if (x_int <= 2 && y > RIGHT) {
	185	t = y - x0;
	186	y--; x_int++;
	187	goto CONTINUE;
	188	} else if (y < -LEFT) {
	189	t = y +(1.0-x0);
	190	CONTINUE:
05909853 KB	191	z = t - x0_lo;
	192	p = r0+z(r1+z(r2+z(r3+z(r4+z*r5))));
	193	q = s0+z(s1+z(s2+z(s3+zs4)));
0662f309 PM	194	r = t(z(p/q) - x0_lo);
	195	t = .5tt;
	196	z = 1.0;
	197	switch (x_int) {
	198	case 6: z = (y + 5);
	199	case 5: z *= (y + 4);
	200	case 4: z *= (y + 3);
	201	case 3: z *= (y + 2);
	202	rr = log__D(z);
	203	rr.b += a0_lo; rr.a += a0_hi;
	204	return(((r+rr.b)+t+rr.a));
	205	case 2: return(((r+a0_lo)+t)+a0_hi);
	206	case 0: r -= log1p(x);
	207	default: rr = log__D(x);
	208	rr.a -= a0_hi; rr.b -= a0_lo;
	209	return(((r - rr.b) + t) - rr.a);
	210	}
05909853	211	} else {
0662f309	212	p = p0+y(p1+y(p2+y(p3+yp4)));
05909853	213	q = q0+y(q1+y(q2+y(q3+y(q4+y(q5+yq6)))));
0662f309 PM	214	p = p*(y/q);
	215	t = (double)(float) y;
	216	z = y-t;
	217	hi = (double)(float) (p+a1_hi);
	218	lo = a1_hi - hi; lo += p; lo += a1_lo;
	219	r = loy + zhi; /* q + r = y(a0+p/q) /
	220	q = hi*t;
	221	z = 1.0;
	222	switch (x_int) {
	223	case 6: z = (y + 5);
	224	case 5: z *= (y + 4);
	225	case 4: z *= (y + 3);
	226	case 3: z *= (y + 2);
	227	rr = log__D(z);
	228	r += rr.b; r += q;
	229	return(rr.a + r);
	230	case 2: return (q+ r);
	231	case 0: rr = log__D(x);
	232	r -= rr.b; r -= log1p(x);
	233	r += q; r-= rr.a;
	234	return(r);
	235	default: rr = log__D(x);
	236	r -= rr.b;
	237	q -= rr.a;
	238	return (r+q);
	239	}
0ebe4989	240	}
0ebe4989 ZAL	241	}
	242
	243	static double
0662f309	244	neg_lgam(double x)
0ebe4989	245	{
19d48448	246	int xi;
05909853	247	double y, z, one = 1.0, zero = 0.0;
19d48448	248	extern double gamma();
0ebe4989	249
19d48448 KB	250	/* avoid destructive cancellation as much as possible */
	251	if (x > -170) {
	252	xi = x;
	253	if (xi == x)
	254	if (_IEEE)
	255	return(one/zero);
	256	else
	257	return(infnan(ERANGE));
	258	y = gamma(x);
	259	if (y < 0)
	260	y = -y, signgam = -1;
	261	return (log(y));
	262	}
05909853	263	z = floor(x + .5);
0662f309 PM	264	if (z == x) { /* convention: G(-(integer)) -> +Inf */
	265	if (_IEEE)
	266	return (one/zero);
	267	else
	268	return (infnan(ERANGE));
0ebe4989	269	}
0662f309 PM	270	y = .5*ceil(x);
0662f309 PM	271	if (y == ceil(y))
05909853	272	signgam = -1;
05909853 KB	273	x = -x;
	274	z = fabs(x + z); /* 0 < z <= .5 */
	275	if (z < .25)
	276	z = sin(M_PI*z);
	277	else
0662f309	278	z = cos(M_PI*(0.5-z));
19d48448	279	z = log(M_PI/(z*x));
9ab6b304	280	y = large_lgam(x);
19d48448	281	return (z - y);
0ebe4989	282	}