xref: /illumos-gate/usr/src/lib/libm/common/C/__lgamma.c (revision 25c28e83)
1*25c28e83SPiotr Jasiukajtis /*
2*25c28e83SPiotr Jasiukajtis  * CDDL HEADER START
3*25c28e83SPiotr Jasiukajtis  *
4*25c28e83SPiotr Jasiukajtis  * The contents of this file are subject to the terms of the
5*25c28e83SPiotr Jasiukajtis  * Common Development and Distribution License (the "License").
6*25c28e83SPiotr Jasiukajtis  * You may not use this file except in compliance with the License.
7*25c28e83SPiotr Jasiukajtis  *
8*25c28e83SPiotr Jasiukajtis  * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9*25c28e83SPiotr Jasiukajtis  * or http://www.opensolaris.org/os/licensing.
10*25c28e83SPiotr Jasiukajtis  * See the License for the specific language governing permissions
11*25c28e83SPiotr Jasiukajtis  * and limitations under the License.
12*25c28e83SPiotr Jasiukajtis  *
13*25c28e83SPiotr Jasiukajtis  * When distributing Covered Code, include this CDDL HEADER in each
14*25c28e83SPiotr Jasiukajtis  * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
15*25c28e83SPiotr Jasiukajtis  * If applicable, add the following below this CDDL HEADER, with the
16*25c28e83SPiotr Jasiukajtis  * fields enclosed by brackets "[]" replaced with your own identifying
17*25c28e83SPiotr Jasiukajtis  * information: Portions Copyright [yyyy] [name of copyright owner]
18*25c28e83SPiotr Jasiukajtis  *
19*25c28e83SPiotr Jasiukajtis  * CDDL HEADER END
20*25c28e83SPiotr Jasiukajtis  */
21*25c28e83SPiotr Jasiukajtis /*
22*25c28e83SPiotr Jasiukajtis  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
23*25c28e83SPiotr Jasiukajtis  */
24*25c28e83SPiotr Jasiukajtis /*
25*25c28e83SPiotr Jasiukajtis  * Copyright 2005 Sun Microsystems, Inc.  All rights reserved.
26*25c28e83SPiotr Jasiukajtis  * Use is subject to license terms.
27*25c28e83SPiotr Jasiukajtis  */
28*25c28e83SPiotr Jasiukajtis 
29*25c28e83SPiotr Jasiukajtis /*
30*25c28e83SPiotr Jasiukajtis  * double __k_lgamma(double x, int *signgamp);
31*25c28e83SPiotr Jasiukajtis  *
32*25c28e83SPiotr Jasiukajtis  * K.C. Ng, March, 1989.
33*25c28e83SPiotr Jasiukajtis  *
34*25c28e83SPiotr Jasiukajtis  * Part of the algorithm is based on W. Cody's lgamma function.
35*25c28e83SPiotr Jasiukajtis  */
36*25c28e83SPiotr Jasiukajtis 
37*25c28e83SPiotr Jasiukajtis #include "libm.h"
38*25c28e83SPiotr Jasiukajtis 
39*25c28e83SPiotr Jasiukajtis static const double
40*25c28e83SPiotr Jasiukajtis one	= 1.0,
41*25c28e83SPiotr Jasiukajtis zero	= 0.0,
42*25c28e83SPiotr Jasiukajtis hln2pi	= 0.9189385332046727417803297,	/* log(2*pi)/2 */
43*25c28e83SPiotr Jasiukajtis pi	= 3.1415926535897932384626434,
44*25c28e83SPiotr Jasiukajtis two52	= 4503599627370496.0,		/* 43300000,00000000 (used by sin_pi) */
45*25c28e83SPiotr Jasiukajtis /*
46*25c28e83SPiotr Jasiukajtis  * Numerator and denominator coefficients for rational minimax Approximation
47*25c28e83SPiotr Jasiukajtis  * P/Q over (0.5,1.5).
48*25c28e83SPiotr Jasiukajtis  */
49*25c28e83SPiotr Jasiukajtis D1 = 	-5.772156649015328605195174e-1,
50*25c28e83SPiotr Jasiukajtis p7 =	 4.945235359296727046734888e0,
51*25c28e83SPiotr Jasiukajtis p6 =	 2.018112620856775083915565e2,
52*25c28e83SPiotr Jasiukajtis p5 =	 2.290838373831346393026739e3,
53*25c28e83SPiotr Jasiukajtis p4 =	 1.131967205903380828685045e4,
54*25c28e83SPiotr Jasiukajtis p3 =	 2.855724635671635335736389e4,
55*25c28e83SPiotr Jasiukajtis p2 =	 3.848496228443793359990269e4,
56*25c28e83SPiotr Jasiukajtis p1 =	 2.637748787624195437963534e4,
57*25c28e83SPiotr Jasiukajtis p0 =	 7.225813979700288197698961e3,
58*25c28e83SPiotr Jasiukajtis q7 =	 6.748212550303777196073036e1,
59*25c28e83SPiotr Jasiukajtis q6 =	 1.113332393857199323513008e3,
60*25c28e83SPiotr Jasiukajtis q5 =	 7.738757056935398733233834e3,
61*25c28e83SPiotr Jasiukajtis q4 =	 2.763987074403340708898585e4,
62*25c28e83SPiotr Jasiukajtis q3 =	 5.499310206226157329794414e4,
63*25c28e83SPiotr Jasiukajtis q2 =	 6.161122180066002127833352e4,
64*25c28e83SPiotr Jasiukajtis q1 =	 3.635127591501940507276287e4,
65*25c28e83SPiotr Jasiukajtis q0 =	 8.785536302431013170870835e3,
66*25c28e83SPiotr Jasiukajtis /*
67*25c28e83SPiotr Jasiukajtis  * Numerator and denominator coefficients for rational minimax Approximation
68*25c28e83SPiotr Jasiukajtis  * G/H over (1.5,4.0).
69*25c28e83SPiotr Jasiukajtis  */
70*25c28e83SPiotr Jasiukajtis D2 =	 4.227843350984671393993777e-1,
71*25c28e83SPiotr Jasiukajtis g7 =	 4.974607845568932035012064e0,
72*25c28e83SPiotr Jasiukajtis g6 =	 5.424138599891070494101986e2,
73*25c28e83SPiotr Jasiukajtis g5 =	 1.550693864978364947665077e4,
74*25c28e83SPiotr Jasiukajtis g4 =	 1.847932904445632425417223e5,
75*25c28e83SPiotr Jasiukajtis g3 =	 1.088204769468828767498470e6,
76*25c28e83SPiotr Jasiukajtis g2 =	 3.338152967987029735917223e6,
77*25c28e83SPiotr Jasiukajtis g1 =	 5.106661678927352456275255e6,
78*25c28e83SPiotr Jasiukajtis g0 =	 3.074109054850539556250927e6,
79*25c28e83SPiotr Jasiukajtis h7 =	 1.830328399370592604055942e2,
80*25c28e83SPiotr Jasiukajtis h6 =	 7.765049321445005871323047e3,
81*25c28e83SPiotr Jasiukajtis h5 =	 1.331903827966074194402448e5,
82*25c28e83SPiotr Jasiukajtis h4 =	 1.136705821321969608938755e6,
83*25c28e83SPiotr Jasiukajtis h3 =	 5.267964117437946917577538e6,
84*25c28e83SPiotr Jasiukajtis h2 =	 1.346701454311101692290052e7,
85*25c28e83SPiotr Jasiukajtis h1 =	 1.782736530353274213975932e7,
86*25c28e83SPiotr Jasiukajtis h0 =	 9.533095591844353613395747e6,
87*25c28e83SPiotr Jasiukajtis /*
88*25c28e83SPiotr Jasiukajtis  * Numerator and denominator coefficients for rational minimax Approximation
89*25c28e83SPiotr Jasiukajtis  * U/V over (4.0,12.0).
90*25c28e83SPiotr Jasiukajtis  */
91*25c28e83SPiotr Jasiukajtis D4 =	 1.791759469228055000094023e0,
92*25c28e83SPiotr Jasiukajtis u7 =	 1.474502166059939948905062e4,
93*25c28e83SPiotr Jasiukajtis u6 =	 2.426813369486704502836312e6,
94*25c28e83SPiotr Jasiukajtis u5 =	 1.214755574045093227939592e8,
95*25c28e83SPiotr Jasiukajtis u4 =	 2.663432449630976949898078e9,
96*25c28e83SPiotr Jasiukajtis u3 =	 2.940378956634553899906876e10,
97*25c28e83SPiotr Jasiukajtis u2 =	 1.702665737765398868392998e11,
98*25c28e83SPiotr Jasiukajtis u1 =	 4.926125793377430887588120e11,
99*25c28e83SPiotr Jasiukajtis u0 =	 5.606251856223951465078242e11,
100*25c28e83SPiotr Jasiukajtis v7 =	 2.690530175870899333379843e3,
101*25c28e83SPiotr Jasiukajtis v6 =	 6.393885654300092398984238e5,
102*25c28e83SPiotr Jasiukajtis v5 =	 4.135599930241388052042842e7,
103*25c28e83SPiotr Jasiukajtis v4 =	 1.120872109616147941376570e9,
104*25c28e83SPiotr Jasiukajtis v3 =	 1.488613728678813811542398e10,
105*25c28e83SPiotr Jasiukajtis v2 =	 1.016803586272438228077304e11,
106*25c28e83SPiotr Jasiukajtis v1 =	 3.417476345507377132798597e11,
107*25c28e83SPiotr Jasiukajtis v0 =	 4.463158187419713286462081e11,
108*25c28e83SPiotr Jasiukajtis /*
109*25c28e83SPiotr Jasiukajtis  * Coefficients for minimax approximation over (12, INF).
110*25c28e83SPiotr Jasiukajtis  */
111*25c28e83SPiotr Jasiukajtis c5 =	-1.910444077728e-03,
112*25c28e83SPiotr Jasiukajtis c4 =	 8.4171387781295e-04,
113*25c28e83SPiotr Jasiukajtis c3 =	-5.952379913043012e-04,
114*25c28e83SPiotr Jasiukajtis c2 =	 7.93650793500350248e-04,
115*25c28e83SPiotr Jasiukajtis c1 =	-2.777777777777681622553e-03,
116*25c28e83SPiotr Jasiukajtis c0 =	 8.333333333333333331554247e-02,
117*25c28e83SPiotr Jasiukajtis c6 =	 5.7083835261e-03;
118*25c28e83SPiotr Jasiukajtis 
119*25c28e83SPiotr Jasiukajtis /*
120*25c28e83SPiotr Jasiukajtis  * Return sin(pi*x).  We assume x is finite and negative, and if it
121*25c28e83SPiotr Jasiukajtis  * is an integer, then the sign of the zero returned doesn't matter.
122*25c28e83SPiotr Jasiukajtis  */
123*25c28e83SPiotr Jasiukajtis static double
sin_pi(double x)124*25c28e83SPiotr Jasiukajtis sin_pi(double x) {
125*25c28e83SPiotr Jasiukajtis 	double	y, z;
126*25c28e83SPiotr Jasiukajtis 	int	n;
127*25c28e83SPiotr Jasiukajtis 
128*25c28e83SPiotr Jasiukajtis 	y = -x;
129*25c28e83SPiotr Jasiukajtis 	if (y <= 0.25)
130*25c28e83SPiotr Jasiukajtis 		return (__k_sin(pi * x, 0.0));
131*25c28e83SPiotr Jasiukajtis 	if (y >= two52)
132*25c28e83SPiotr Jasiukajtis 		return (zero);
133*25c28e83SPiotr Jasiukajtis 	z = floor(y);
134*25c28e83SPiotr Jasiukajtis 	if (y == z)
135*25c28e83SPiotr Jasiukajtis 		return (zero);
136*25c28e83SPiotr Jasiukajtis 
137*25c28e83SPiotr Jasiukajtis 	/* argument reduction: set y = |x| mod 2 */
138*25c28e83SPiotr Jasiukajtis 	y *= 0.5;
139*25c28e83SPiotr Jasiukajtis 	y = 2.0 * (y - floor(y));
140*25c28e83SPiotr Jasiukajtis 
141*25c28e83SPiotr Jasiukajtis 	/* now floor(y * 4) tells which octant y is in */
142*25c28e83SPiotr Jasiukajtis 	n = (int)(y * 4.0);
143*25c28e83SPiotr Jasiukajtis 	switch (n) {
144*25c28e83SPiotr Jasiukajtis 	case 0:
145*25c28e83SPiotr Jasiukajtis 		y = __k_sin(pi * y, 0.0);
146*25c28e83SPiotr Jasiukajtis 		break;
147*25c28e83SPiotr Jasiukajtis 	case 1:
148*25c28e83SPiotr Jasiukajtis 	case 2:
149*25c28e83SPiotr Jasiukajtis 		y = __k_cos(pi * (0.5 - y), 0.0);
150*25c28e83SPiotr Jasiukajtis 		break;
151*25c28e83SPiotr Jasiukajtis 	case 3:
152*25c28e83SPiotr Jasiukajtis 	case 4:
153*25c28e83SPiotr Jasiukajtis 		y = __k_sin(pi * (1.0 - y), 0.0);
154*25c28e83SPiotr Jasiukajtis 		break;
155*25c28e83SPiotr Jasiukajtis 	case 5:
156*25c28e83SPiotr Jasiukajtis 	case 6:
157*25c28e83SPiotr Jasiukajtis 		y = -__k_cos(pi * (y - 1.5), 0.0);
158*25c28e83SPiotr Jasiukajtis 		break;
159*25c28e83SPiotr Jasiukajtis 	default:
160*25c28e83SPiotr Jasiukajtis 		y = __k_sin(pi * (y - 2.0), 0.0);
161*25c28e83SPiotr Jasiukajtis 		break;
162*25c28e83SPiotr Jasiukajtis 	}
163*25c28e83SPiotr Jasiukajtis 	return (-y);
164*25c28e83SPiotr Jasiukajtis }
165*25c28e83SPiotr Jasiukajtis 
166*25c28e83SPiotr Jasiukajtis static double
neg(double z,int * signgamp)167*25c28e83SPiotr Jasiukajtis neg(double z, int *signgamp) {
168*25c28e83SPiotr Jasiukajtis 	double	t, p;
169*25c28e83SPiotr Jasiukajtis 
170*25c28e83SPiotr Jasiukajtis 	/*
171*25c28e83SPiotr Jasiukajtis 	 * written by K.C. Ng,  Feb 2, 1989.
172*25c28e83SPiotr Jasiukajtis 	 *
173*25c28e83SPiotr Jasiukajtis 	 * Since
174*25c28e83SPiotr Jasiukajtis 	 *		-z*G(-z)*G(z) = pi/sin(pi*z),
175*25c28e83SPiotr Jasiukajtis 	 * we have
176*25c28e83SPiotr Jasiukajtis 	 * 	G(-z) = -pi/(sin(pi*z)*G(z)*z)
177*25c28e83SPiotr Jasiukajtis 	 * 	      =  pi/(sin(pi*(-z))*G(z)*z)
178*25c28e83SPiotr Jasiukajtis 	 * Algorithm
179*25c28e83SPiotr Jasiukajtis 	 *		z = |z|
180*25c28e83SPiotr Jasiukajtis 	 *		t = sin_pi(z); ...note that when z>2**52, z is an int
181*25c28e83SPiotr Jasiukajtis 	 *		and hence t=0.
182*25c28e83SPiotr Jasiukajtis 	 *
183*25c28e83SPiotr Jasiukajtis 	 *		if (t == 0.0) return 1.0/0.0;
184*25c28e83SPiotr Jasiukajtis 	 *		if (t< 0.0) *signgamp = -1; else t= -t;
185*25c28e83SPiotr Jasiukajtis 	 *		if (z+1.0 == 1.0)	...tiny z
186*25c28e83SPiotr Jasiukajtis 	 *		    return -log(z);
187*25c28e83SPiotr Jasiukajtis 	 *		else
188*25c28e83SPiotr Jasiukajtis 	 *		    return log(pi/(t*z))-__k_lgamma(z, signgamp);
189*25c28e83SPiotr Jasiukajtis 	 */
190*25c28e83SPiotr Jasiukajtis 
191*25c28e83SPiotr Jasiukajtis 	t = sin_pi(z);			/* t := sin(pi*z) */
192*25c28e83SPiotr Jasiukajtis 	if (t == zero)			/* return 1.0/0.0 = +INF */
193*25c28e83SPiotr Jasiukajtis 		return (one / fabs(t));
194*25c28e83SPiotr Jasiukajtis 	z = -z;
195*25c28e83SPiotr Jasiukajtis 	p = z + one;
196*25c28e83SPiotr Jasiukajtis 	if (p == one)
197*25c28e83SPiotr Jasiukajtis 		p = -log(z);
198*25c28e83SPiotr Jasiukajtis 	else
199*25c28e83SPiotr Jasiukajtis 		p = log(pi / (fabs(t) * z)) - __k_lgamma(z, signgamp);
200*25c28e83SPiotr Jasiukajtis 	if (t < zero)
201*25c28e83SPiotr Jasiukajtis 		*signgamp = -1;
202*25c28e83SPiotr Jasiukajtis 	return (p);
203*25c28e83SPiotr Jasiukajtis }
204*25c28e83SPiotr Jasiukajtis 
205*25c28e83SPiotr Jasiukajtis double
__k_lgamma(double x,int * signgamp)206*25c28e83SPiotr Jasiukajtis __k_lgamma(double x, int *signgamp) {
207*25c28e83SPiotr Jasiukajtis 	double	t, p, q, cr, y;
208*25c28e83SPiotr Jasiukajtis 
209*25c28e83SPiotr Jasiukajtis 	/* purge off +-inf, NaN and negative arguments */
210*25c28e83SPiotr Jasiukajtis 	if (!finite(x))
211*25c28e83SPiotr Jasiukajtis 		return (x * x);
212*25c28e83SPiotr Jasiukajtis 	*signgamp = 1;
213*25c28e83SPiotr Jasiukajtis 	if (signbit(x))
214*25c28e83SPiotr Jasiukajtis 		return (neg(x, signgamp));
215*25c28e83SPiotr Jasiukajtis 
216*25c28e83SPiotr Jasiukajtis 	/* lgamma(x) ~ log(1/x) for really tiny x */
217*25c28e83SPiotr Jasiukajtis 	t = one + x;
218*25c28e83SPiotr Jasiukajtis 	if (t == one) {
219*25c28e83SPiotr Jasiukajtis 		if (x == zero)
220*25c28e83SPiotr Jasiukajtis 			return (one / x);
221*25c28e83SPiotr Jasiukajtis 		return (-log(x));
222*25c28e83SPiotr Jasiukajtis 	}
223*25c28e83SPiotr Jasiukajtis 
224*25c28e83SPiotr Jasiukajtis 	/* for tiny < x < inf */
225*25c28e83SPiotr Jasiukajtis 	if (x <= 1.5) {
226*25c28e83SPiotr Jasiukajtis 		if (x < 0.6796875) {
227*25c28e83SPiotr Jasiukajtis 			cr = -log(x);
228*25c28e83SPiotr Jasiukajtis 			y = x;
229*25c28e83SPiotr Jasiukajtis 		} else {
230*25c28e83SPiotr Jasiukajtis 			cr = zero;
231*25c28e83SPiotr Jasiukajtis 			y = x - one;
232*25c28e83SPiotr Jasiukajtis 		}
233*25c28e83SPiotr Jasiukajtis 
234*25c28e83SPiotr Jasiukajtis 		if (x <= 0.5 || x >= 0.6796875) {
235*25c28e83SPiotr Jasiukajtis 			if (x == one)
236*25c28e83SPiotr Jasiukajtis 				return (zero);
237*25c28e83SPiotr Jasiukajtis 			p = p0+y*(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7))))));
238*25c28e83SPiotr Jasiukajtis 			q = q0+y*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*(q6+y*
239*25c28e83SPiotr Jasiukajtis 			    (q7+y)))))));
240*25c28e83SPiotr Jasiukajtis 			return (cr+y*(D1+y*(p/q)));
241*25c28e83SPiotr Jasiukajtis 		} else {
242*25c28e83SPiotr Jasiukajtis 			y = x - one;
243*25c28e83SPiotr Jasiukajtis 			p = g0+y*(g1+y*(g2+y*(g3+y*(g4+y*(g5+y*(g6+y*g7))))));
244*25c28e83SPiotr Jasiukajtis 			q = h0+y*(h1+y*(h2+y*(h3+y*(h4+y*(h5+y*(h6+y*
245*25c28e83SPiotr Jasiukajtis 			    (h7+y)))))));
246*25c28e83SPiotr Jasiukajtis 			return (cr+y*(D2+y*(p/q)));
247*25c28e83SPiotr Jasiukajtis 		}
248*25c28e83SPiotr Jasiukajtis 	} else if (x <= 4.0) {
249*25c28e83SPiotr Jasiukajtis 		if (x == 2.0)
250*25c28e83SPiotr Jasiukajtis 			return (zero);
251*25c28e83SPiotr Jasiukajtis 		y = x - 2.0;
252*25c28e83SPiotr Jasiukajtis 		p = g0+y*(g1+y*(g2+y*(g3+y*(g4+y*(g5+y*(g6+y*g7))))));
253*25c28e83SPiotr Jasiukajtis 		q = h0+y*(h1+y*(h2+y*(h3+y*(h4+y*(h5+y*(h6+y*(h7+y)))))));
254*25c28e83SPiotr Jasiukajtis 		return (y*(D2+y*(p/q)));
255*25c28e83SPiotr Jasiukajtis 	} else if (x <= 12.0) {
256*25c28e83SPiotr Jasiukajtis 		y = x - 4.0;
257*25c28e83SPiotr Jasiukajtis 		p = u0+y*(u1+y*(u2+y*(u3+y*(u4+y*(u5+y*(u6+y*u7))))));
258*25c28e83SPiotr Jasiukajtis 		q = v0+y*(v1+y*(v2+y*(v3+y*(v4+y*(v5+y*(v6+y*(v7-y)))))));
259*25c28e83SPiotr Jasiukajtis 		return (D4+y*(p/q));
260*25c28e83SPiotr Jasiukajtis 	} else if (x <= 1.0e17) {		/* x ~< 2**(prec+3) */
261*25c28e83SPiotr Jasiukajtis 		t = one / x;
262*25c28e83SPiotr Jasiukajtis 		y = t * t;
263*25c28e83SPiotr Jasiukajtis 		p = hln2pi+t*(c0+y*(c1+y*(c2+y*(c3+y*(c4+y*(c5+y*c6))))));
264*25c28e83SPiotr Jasiukajtis 		q = log(x);
265*25c28e83SPiotr Jasiukajtis 		return (x*(q-one)-(0.5*q-p));
266*25c28e83SPiotr Jasiukajtis 	} else {			/* may overflow */
267*25c28e83SPiotr Jasiukajtis 		return (x * (log(x) - 1.0));
268*25c28e83SPiotr Jasiukajtis 	}
269*25c28e83SPiotr Jasiukajtis }
270