b_tgamma.c (69160b1eb704d1931b536a65fcf50757ddbfc469) b_tgamma.c (46d7c2979ef66b28bbe1d5e3db922f3b47784ccd)
1/*-
2 * Copyright (c) 1992, 1993
3 * The Regents of the University of California. All rights reserved.
4 *
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

--- 20 unchanged lines hidden (view full) ---

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.
32 */
33
34#ifndef lint
35static char sccsid[] = "@(#)gamma.c 8.1 (Berkeley) 6/4/93";
36#endif /* not lint */
1/*-
2 * Copyright (c) 1992, 1993
3 * The Regents of the University of California. All rights reserved.
4 *
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

--- 20 unchanged lines hidden (view full) ---

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.
32 */
33
34#ifndef lint
35static char sccsid[] = "@(#)gamma.c 8.1 (Berkeley) 6/4/93";
36#endif /* not lint */
37include <sys/cdefs.h>
37#include <sys/cdefs.h>
38__FBSDID("$FreeBSD$");
39
40/*
41 * This code by P. McIlroy, Oct 1992;
42 *
43 * The financial support of UUNET Communications Services is greatfully
44 * acknowledged.
45 */

--- 87 unchanged lines hidden (view full) ---

133#define TRUNC(x) x = (double) (float) (x)
134#else
135#define _IEEE 1
136#define TRUNC(x) *(((int *) &x) + endian) &= 0xf8000000
137#define infnan(x) 0.0
138#endif
139
140double
38__FBSDID("$FreeBSD$");
39
40/*
41 * This code by P. McIlroy, Oct 1992;
42 *
43 * The financial support of UUNET Communications Services is greatfully
44 * acknowledged.
45 */

--- 87 unchanged lines hidden (view full) ---

133#define TRUNC(x) x = (double) (float) (x)
134#else
135#define _IEEE 1
136#define TRUNC(x) *(((int *) &x) + endian) &= 0xf8000000
137#define infnan(x) 0.0
138#endif
139
140double
141gamma(x)
141tgamma(x)
142 double x;
143{
144 struct Double u;
145 endian = (*(int *) &one) ? 1 : 0;
146
147 if (x >= 6) {
148 if(x > 171.63)
149 return(one/zero);

--- 70 unchanged lines hidden (view full) ---

220 /* Argument reduction: G(x+1) = x*G(x) */
221 for (ym1 = y-one; ym1 > LEFT + x0; y = ym1--, yy.a--) {
222 t = r.a*yy.a;
223 r.b = r.a*yy.b + y*r.b;
224 r.a = t;
225 TRUNC(r.a);
226 r.b += (t - r.a);
227 }
142 double x;
143{
144 struct Double u;
145 endian = (*(int *) &one) ? 1 : 0;
146
147 if (x >= 6) {
148 if(x > 171.63)
149 return(one/zero);

--- 70 unchanged lines hidden (view full) ---

220 /* Argument reduction: G(x+1) = x*G(x) */
221 for (ym1 = y-one; ym1 > LEFT + x0; y = ym1--, yy.a--) {
222 t = r.a*yy.a;
223 r.b = r.a*yy.b + y*r.b;
224 r.a = t;
225 TRUNC(r.a);
226 r.b += (t - r.a);
227 }
228 /* Return r*gamma(y). */
228 /* Return r*tgamma(y). */
229 yy = ratfun_gam(y - x0, 0);
230 y = r.b*(yy.a + yy.b) + r.a*yy.b;
231 y += yy.a*r.a;
232 return (y);
233}
234/*
235 * Good on (0, 1+x0+LEFT]. Accurate to 1ulp.
236 */

--- 88 unchanged lines hidden (view full) ---

325 y = -(lg.a + lg.b);
326 z = (y + lg.a) + lg.b;
327 y = __exp__D(y, z);
328 if (sgn < 0) y = -y;
329 return (y);
330 }
331 y = one-x;
332 if (one-y == x)
229 yy = ratfun_gam(y - x0, 0);
230 y = r.b*(yy.a + yy.b) + r.a*yy.b;
231 y += yy.a*r.a;
232 return (y);
233}
234/*
235 * Good on (0, 1+x0+LEFT]. Accurate to 1ulp.
236 */

--- 88 unchanged lines hidden (view full) ---

325 y = -(lg.a + lg.b);
326 z = (y + lg.a) + lg.b;
327 y = __exp__D(y, z);
328 if (sgn < 0) y = -y;
329 return (y);
330 }
331 y = one-x;
332 if (one-y == x)
333 y = gamma(y);
333 y = tgamma(y);
334 else /* 1-x is inexact */
334 else /* 1-x is inexact */
335 y = -x*gamma(-x);
335 y = -x*tgamma(-x);
336 if (sgn < 0) y = -y;
337 return (M_PI / (y*z));
338}
336 if (sgn < 0) y = -y;
337 return (M_PI / (y*z));
338}