1da2e3ebdSchin #include "FEATURE/uwin"
2da2e3ebdSchin
3da2e3ebdSchin #if !_UWIN || _lib_lgamma
4da2e3ebdSchin
_STUB_lgamma()5da2e3ebdSchin void _STUB_lgamma(){}
6da2e3ebdSchin
7da2e3ebdSchin #else
8da2e3ebdSchin
9da2e3ebdSchin /*-
10da2e3ebdSchin * Copyright (c) 1992, 1993
11da2e3ebdSchin * The Regents of the University of California. All rights reserved.
12da2e3ebdSchin *
13da2e3ebdSchin * Redistribution and use in source and binary forms, with or without
14da2e3ebdSchin * modification, are permitted provided that the following conditions
15da2e3ebdSchin * are met:
16da2e3ebdSchin * 1. Redistributions of source code must retain the above copyright
17da2e3ebdSchin * notice, this list of conditions and the following disclaimer.
18da2e3ebdSchin * 2. Redistributions in binary form must reproduce the above copyright
19da2e3ebdSchin * notice, this list of conditions and the following disclaimer in the
20da2e3ebdSchin * documentation and/or other materials provided with the distribution.
21da2e3ebdSchin * 3. Neither the name of the University nor the names of its contributors
22da2e3ebdSchin * may be used to endorse or promote products derived from this software
23da2e3ebdSchin * without specific prior written permission.
24da2e3ebdSchin *
25da2e3ebdSchin * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
26da2e3ebdSchin * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27da2e3ebdSchin * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
28da2e3ebdSchin * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
29da2e3ebdSchin * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
30da2e3ebdSchin * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
31da2e3ebdSchin * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
32da2e3ebdSchin * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
33da2e3ebdSchin * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
34da2e3ebdSchin * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
35da2e3ebdSchin * SUCH DAMAGE.
36da2e3ebdSchin */
37da2e3ebdSchin
38da2e3ebdSchin #ifndef lint
39da2e3ebdSchin static char sccsid[] = "@(#)lgamma.c 8.2 (Berkeley) 11/30/93";
40da2e3ebdSchin #endif /* not lint */
41da2e3ebdSchin
42da2e3ebdSchin /*
43da2e3ebdSchin * Coded by Peter McIlroy, Nov 1992;
44da2e3ebdSchin *
45da2e3ebdSchin * The financial support of UUNET Communications Services is greatfully
46da2e3ebdSchin * acknowledged.
47da2e3ebdSchin */
48da2e3ebdSchin
49da2e3ebdSchin #define gamma ______gamma
50da2e3ebdSchin #define lgamma ______lgamma
51da2e3ebdSchin
52da2e3ebdSchin #include <math.h>
53da2e3ebdSchin #include <errno.h>
54da2e3ebdSchin #include "mathimpl.h"
55da2e3ebdSchin
56da2e3ebdSchin #undef gamma
57da2e3ebdSchin #undef lgamma
58da2e3ebdSchin
59da2e3ebdSchin /* Log gamma function.
60da2e3ebdSchin * Error: x > 0 error < 1.3ulp.
61da2e3ebdSchin * x > 4, error < 1ulp.
62da2e3ebdSchin * x > 9, error < .6ulp.
63da2e3ebdSchin * x < 0, all bets are off. (When G(x) ~ 1, log(G(x)) ~ 0)
64da2e3ebdSchin * Method:
65da2e3ebdSchin * x > 6:
66da2e3ebdSchin * Use the asymptotic expansion (Stirling's Formula)
67da2e3ebdSchin * 0 < x < 6:
68da2e3ebdSchin * Use gamma(x+1) = x*gamma(x) for argument reduction.
69da2e3ebdSchin * Use rational approximation in
70da2e3ebdSchin * the range 1.2, 2.5
71da2e3ebdSchin * Two approximations are used, one centered at the
72da2e3ebdSchin * minimum to ensure monotonicity; one centered at 2
73da2e3ebdSchin * to maintain small relative error.
74da2e3ebdSchin * x < 0:
75da2e3ebdSchin * Use the reflection formula,
76da2e3ebdSchin * G(1-x)G(x) = PI/sin(PI*x)
77da2e3ebdSchin * Special values:
78da2e3ebdSchin * non-positive integer returns +Inf.
79da2e3ebdSchin * NaN returns NaN
80da2e3ebdSchin */
81da2e3ebdSchin static int endian;
82da2e3ebdSchin #if defined(vax) || defined(tahoe)
83da2e3ebdSchin #define _IEEE 0
84da2e3ebdSchin /* double and float have same size exponent field */
85da2e3ebdSchin #define TRUNC(x) x = (double) (float) (x)
86da2e3ebdSchin #else
87da2e3ebdSchin #define _IEEE 1
88da2e3ebdSchin #define TRUNC(x) *(((int *) &x) + endian) &= 0xf8000000
89da2e3ebdSchin #define infnan(x) 0.0
90da2e3ebdSchin #endif
91da2e3ebdSchin
92da2e3ebdSchin static double small_lgam(double);
93da2e3ebdSchin static double large_lgam(double);
94da2e3ebdSchin static double neg_lgam(double);
95da2e3ebdSchin static double zero = 0.0, one = 1.0;
96da2e3ebdSchin int signgam;
97da2e3ebdSchin
98da2e3ebdSchin #define UNDERFL (1e-1020 * 1e-1020)
99da2e3ebdSchin
100da2e3ebdSchin #define LEFT (1.0 - (x0 + .25))
101da2e3ebdSchin #define RIGHT (x0 - .218)
102da2e3ebdSchin /*
103da2e3ebdSchin * Constants for approximation in [1.244,1.712]
104da2e3ebdSchin */
105da2e3ebdSchin #define x0 0.461632144968362356785
106da2e3ebdSchin #define x0_lo -.000000000000000015522348162858676890521
107da2e3ebdSchin #define a0_hi -0.12148629128932952880859
108da2e3ebdSchin #define a0_lo .0000000007534799204229502
109da2e3ebdSchin #define r0 -2.771227512955130520e-002
110da2e3ebdSchin #define r1 -2.980729795228150847e-001
111da2e3ebdSchin #define r2 -3.257411333183093394e-001
112da2e3ebdSchin #define r3 -1.126814387531706041e-001
113da2e3ebdSchin #define r4 -1.129130057170225562e-002
114da2e3ebdSchin #define r5 -2.259650588213369095e-005
115da2e3ebdSchin #define s0 1.714457160001714442e+000
116da2e3ebdSchin #define s1 2.786469504618194648e+000
117da2e3ebdSchin #define s2 1.564546365519179805e+000
118da2e3ebdSchin #define s3 3.485846389981109850e-001
119da2e3ebdSchin #define s4 2.467759345363656348e-002
120da2e3ebdSchin /*
121da2e3ebdSchin * Constants for approximation in [1.71, 2.5]
122da2e3ebdSchin */
123da2e3ebdSchin #define a1_hi 4.227843350984671344505727574870e-01
124da2e3ebdSchin #define a1_lo 4.670126436531227189e-18
125da2e3ebdSchin #define p0 3.224670334241133695662995251041e-01
126da2e3ebdSchin #define p1 3.569659696950364669021382724168e-01
127da2e3ebdSchin #define p2 1.342918716072560025853732668111e-01
128da2e3ebdSchin #define p3 1.950702176409779831089963408886e-02
129da2e3ebdSchin #define p4 8.546740251667538090796227834289e-04
130da2e3ebdSchin #define q0 1.000000000000000444089209850062e+00
131da2e3ebdSchin #define q1 1.315850076960161985084596381057e+00
132da2e3ebdSchin #define q2 6.274644311862156431658377186977e-01
133da2e3ebdSchin #define q3 1.304706631926259297049597307705e-01
134da2e3ebdSchin #define q4 1.102815279606722369265536798366e-02
135da2e3ebdSchin #define q5 2.512690594856678929537585620579e-04
136da2e3ebdSchin #define q6 -1.003597548112371003358107325598e-06
137da2e3ebdSchin /*
138da2e3ebdSchin * Stirling's Formula, adjusted for equal-ripple. x in [6,Inf].
139da2e3ebdSchin */
140da2e3ebdSchin #define lns2pi .418938533204672741780329736405
141da2e3ebdSchin #define pb0 8.33333333333333148296162562474e-02
142da2e3ebdSchin #define pb1 -2.77777777774548123579378966497e-03
143da2e3ebdSchin #define pb2 7.93650778754435631476282786423e-04
144da2e3ebdSchin #define pb3 -5.95235082566672847950717262222e-04
145da2e3ebdSchin #define pb4 8.41428560346653702135821806252e-04
146da2e3ebdSchin #define pb5 -1.89773526463879200348872089421e-03
147da2e3ebdSchin #define pb6 5.69394463439411649408050664078e-03
148da2e3ebdSchin #define pb7 -1.44705562421428915453880392761e-02
149da2e3ebdSchin
lgamma(double x)150da2e3ebdSchin extern __pure double lgamma(double x)
151da2e3ebdSchin {
152da2e3ebdSchin double r;
153da2e3ebdSchin
154da2e3ebdSchin signgam = 1;
155da2e3ebdSchin endian = ((*(int *) &one)) ? 1 : 0;
156da2e3ebdSchin
157da2e3ebdSchin if (!finite(x))
158da2e3ebdSchin if (_IEEE)
159da2e3ebdSchin return (x+x);
160da2e3ebdSchin else return (infnan(EDOM));
161da2e3ebdSchin
162da2e3ebdSchin if (x > 6 + RIGHT) {
163da2e3ebdSchin r = large_lgam(x);
164da2e3ebdSchin return (r);
165da2e3ebdSchin } else if (x > 1e-16)
166da2e3ebdSchin return (small_lgam(x));
167da2e3ebdSchin else if (x > -1e-16) {
168da2e3ebdSchin if (x < 0)
169da2e3ebdSchin signgam = -1, x = -x;
170da2e3ebdSchin return (-log(x));
171da2e3ebdSchin } else
172da2e3ebdSchin return (neg_lgam(x));
173da2e3ebdSchin }
174da2e3ebdSchin
175da2e3ebdSchin static double
large_lgam(double x)176da2e3ebdSchin large_lgam(double x)
177da2e3ebdSchin {
178da2e3ebdSchin double z, p, x1;
179da2e3ebdSchin struct Double t, u, v;
180da2e3ebdSchin u = __log__D(x);
181da2e3ebdSchin u.a -= 1.0;
182da2e3ebdSchin if (x > 1e15) {
183da2e3ebdSchin v.a = x - 0.5;
184da2e3ebdSchin TRUNC(v.a);
185da2e3ebdSchin v.b = (x - v.a) - 0.5;
186da2e3ebdSchin t.a = u.a*v.a;
187da2e3ebdSchin t.b = x*u.b + v.b*u.a;
188da2e3ebdSchin if (_IEEE == 0 && !finite(t.a))
189da2e3ebdSchin return(infnan(ERANGE));
190da2e3ebdSchin return(t.a + t.b);
191da2e3ebdSchin }
192da2e3ebdSchin x1 = 1./x;
193da2e3ebdSchin z = x1*x1;
194da2e3ebdSchin p = pb0+z*(pb1+z*(pb2+z*(pb3+z*(pb4+z*(pb5+z*(pb6+z*pb7))))));
195da2e3ebdSchin /* error in approximation = 2.8e-19 */
196da2e3ebdSchin
197da2e3ebdSchin p = p*x1; /* error < 2.3e-18 absolute */
198da2e3ebdSchin /* 0 < p < 1/64 (at x = 5.5) */
199da2e3ebdSchin v.a = x = x - 0.5;
200da2e3ebdSchin TRUNC(v.a); /* truncate v.a to 26 bits. */
201da2e3ebdSchin v.b = x - v.a;
202da2e3ebdSchin t.a = v.a*u.a; /* t = (x-.5)*(log(x)-1) */
203da2e3ebdSchin t.b = v.b*u.a + x*u.b;
204da2e3ebdSchin t.b += p; t.b += lns2pi; /* return t + lns2pi + p */
205da2e3ebdSchin return (t.a + t.b);
206da2e3ebdSchin }
207da2e3ebdSchin
208da2e3ebdSchin static double
small_lgam(double x)209da2e3ebdSchin small_lgam(double x)
210da2e3ebdSchin {
211da2e3ebdSchin int x_int;
212da2e3ebdSchin double y, z, t, r = 0, p, q, hi, lo;
213da2e3ebdSchin struct Double rr;
214da2e3ebdSchin x_int = (int)(x + .5);
215da2e3ebdSchin y = x - x_int;
216da2e3ebdSchin if (x_int <= 2 && y > RIGHT) {
217da2e3ebdSchin t = y - x0;
218da2e3ebdSchin y--; x_int++;
219da2e3ebdSchin goto CONTINUE;
220da2e3ebdSchin } else if (y < -LEFT) {
221da2e3ebdSchin t = y +(1.0-x0);
222da2e3ebdSchin CONTINUE:
223da2e3ebdSchin z = t - x0_lo;
224da2e3ebdSchin p = r0+z*(r1+z*(r2+z*(r3+z*(r4+z*r5))));
225da2e3ebdSchin q = s0+z*(s1+z*(s2+z*(s3+z*s4)));
226da2e3ebdSchin r = t*(z*(p/q) - x0_lo);
227da2e3ebdSchin t = .5*t*t;
228da2e3ebdSchin z = 1.0;
229da2e3ebdSchin switch (x_int) {
230da2e3ebdSchin case 6: z = (y + 5);
231da2e3ebdSchin case 5: z *= (y + 4);
232da2e3ebdSchin case 4: z *= (y + 3);
233da2e3ebdSchin case 3: z *= (y + 2);
234da2e3ebdSchin rr = __log__D(z);
235da2e3ebdSchin rr.b += a0_lo; rr.a += a0_hi;
236da2e3ebdSchin return(((r+rr.b)+t+rr.a));
237da2e3ebdSchin case 2: return(((r+a0_lo)+t)+a0_hi);
238da2e3ebdSchin case 0: r -= log1p(x);
239da2e3ebdSchin default: rr = __log__D(x);
240da2e3ebdSchin rr.a -= a0_hi; rr.b -= a0_lo;
241da2e3ebdSchin return(((r - rr.b) + t) - rr.a);
242da2e3ebdSchin }
243da2e3ebdSchin } else {
244da2e3ebdSchin p = p0+y*(p1+y*(p2+y*(p3+y*p4)));
245da2e3ebdSchin q = q0+y*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*q6)))));
246da2e3ebdSchin p = p*(y/q);
247da2e3ebdSchin t = (double)(float) y;
248da2e3ebdSchin z = y-t;
249da2e3ebdSchin hi = (double)(float) (p+a1_hi);
250da2e3ebdSchin lo = a1_hi - hi; lo += p; lo += a1_lo;
251da2e3ebdSchin r = lo*y + z*hi; /* q + r = y*(a0+p/q) */
252da2e3ebdSchin q = hi*t;
253da2e3ebdSchin z = 1.0;
254da2e3ebdSchin switch (x_int) {
255da2e3ebdSchin case 6: z = (y + 5);
256da2e3ebdSchin case 5: z *= (y + 4);
257da2e3ebdSchin case 4: z *= (y + 3);
258da2e3ebdSchin case 3: z *= (y + 2);
259da2e3ebdSchin rr = __log__D(z);
260da2e3ebdSchin r += rr.b; r += q;
261da2e3ebdSchin return(rr.a + r);
262da2e3ebdSchin case 2: return (q+ r);
263da2e3ebdSchin case 0: rr = __log__D(x);
264da2e3ebdSchin r -= rr.b; r -= log1p(x);
265da2e3ebdSchin r += q; r-= rr.a;
266da2e3ebdSchin return(r);
267da2e3ebdSchin default: rr = __log__D(x);
268da2e3ebdSchin r -= rr.b;
269da2e3ebdSchin q -= rr.a;
270da2e3ebdSchin return (r+q);
271da2e3ebdSchin }
272da2e3ebdSchin }
273da2e3ebdSchin }
274da2e3ebdSchin
275da2e3ebdSchin static double
neg_lgam(double x)276da2e3ebdSchin neg_lgam(double x)
277da2e3ebdSchin {
278da2e3ebdSchin int xi;
279da2e3ebdSchin double y, z, one = 1.0, zero = 0.0;
280da2e3ebdSchin extern double gamma();
281da2e3ebdSchin
282da2e3ebdSchin /* avoid destructive cancellation as much as possible */
283da2e3ebdSchin if (x > -170) {
284da2e3ebdSchin xi = (int)x;
285da2e3ebdSchin if (xi == x)
286da2e3ebdSchin if (_IEEE)
287da2e3ebdSchin return(one/zero);
288da2e3ebdSchin else
289da2e3ebdSchin return(infnan(ERANGE));
290da2e3ebdSchin y = gamma(x);
291da2e3ebdSchin if (y < 0)
292da2e3ebdSchin y = -y, signgam = -1;
293da2e3ebdSchin return (log(y));
294da2e3ebdSchin }
295da2e3ebdSchin z = floor(x + .5);
296da2e3ebdSchin if (z == x) { /* convention: G(-(integer)) -> +Inf */
297da2e3ebdSchin if (_IEEE)
298da2e3ebdSchin return (one/zero);
299da2e3ebdSchin else
300da2e3ebdSchin return (infnan(ERANGE));
301da2e3ebdSchin }
302da2e3ebdSchin y = .5*ceil(x);
303da2e3ebdSchin if (y == ceil(y))
304da2e3ebdSchin signgam = -1;
305da2e3ebdSchin x = -x;
306da2e3ebdSchin z = fabs(x + z); /* 0 < z <= .5 */
307da2e3ebdSchin if (z < .25)
308da2e3ebdSchin z = sin(M_PI*z);
309da2e3ebdSchin else
310da2e3ebdSchin z = cos(M_PI*(0.5-z));
311da2e3ebdSchin z = log(M_PI/(z*x));
312da2e3ebdSchin y = large_lgam(x);
313da2e3ebdSchin return (z - y);
314da2e3ebdSchin }
315da2e3ebdSchin
316da2e3ebdSchin #endif
317