| b_tgamma.c (92917) | b_tgamma.c (93211) |
|---|---|
| 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 38__FBSDID("$FreeBSD: head/lib/msun/bsdsrc/b_tgamma.c 92917 2002-03-21 23:54:04Z obrien $"); | 37#include <sys/cdefs.h> 38__FBSDID("$FreeBSD: head/lib/msun/bsdsrc/b_tgamma.c 93211 2002-03-26 11:59:29Z bde $"); |
| 39 40/* 41 * This code by P. McIlroy, Oct 1992; 42 * 43 * The financial support of UUNET Communications Services is greatfully 44 * acknowledged. 45 */ 46 --- 86 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 | 39 40/* 41 * This code by P. McIlroy, Oct 1992; 42 * 43 * The financial support of UUNET Communications Services is greatfully 44 * acknowledged. 45 */ 46 --- 86 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} |