Head

GCC Code Coverage Report

Directory:	./		Exec	Total	Coverage
File:	lib/libm/src/b_tgamma.c	Lines:	103	132	78.0 %
Date:	2017-11-13	Branches:	30	38	78.9 %


/*	$OpenBSD: b_tgamma.c,v 1.10 2016/09/12 19:47:02 guenther Exp $	*/
/*-
 * Copyright (c) 1992, 1993
 *	The Regents of the University of California.  All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 * 3. Neither the name of the University nor the names of its contributors
 *    may be used to endorse or promote products derived from this software
 *    without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
 * SUCH DAMAGE.
 */

/*
 * This code by P. McIlroy, Oct 1992;
 *
 * The financial support of UUNET Communications Services is greatfully
 * acknowledged.
 */

#include <float.h>
#include <math.h>

#include "math_private.h"

/* METHOD:
 * x < 0: Use reflection formula, G(x) = pi/(sin(pi*x)*x*G(x))
 *	At negative integers, return NaN and raise invalid.
 *
 * x < 6.5:
 *	Use argument reduction G(x+1) = xG(x) to reach the
 *	range [1.066124,2.066124].  Use a rational
 *	approximation centered at the minimum (x0+1) to
 *	ensure monotonicity.
 *
 * x >= 6.5: Use the asymptotic approximation (Stirling's formula)
 *	adjusted for equal-ripples:
 *
 *	log(G(x)) ~= (x-.5)*(log(x)-1) + .5(log(2*pi)-1) + 1/x*P(1/(x*x))
 *
 *	Keep extra precision in multiplying (x-.5)(log(x)-1), to
 *	avoid premature round-off.
 *
 * Special values:
 *	-Inf:			return NaN and raise invalid;
 *	negative integer:	return NaN and raise invalid;
 *	other x ~< -177.79:	return +-0 and raise underflow;
 *	+-0:			return +-Inf and raise divide-by-zero;
 *	finite x ~> 171.63:	return +Inf and raise overflow;
 *	+Inf:			return +Inf;
 *	NaN: 			return NaN.
 *
 * Accuracy: tgamma(x) is accurate to within
 *	x > 0:  error provably < 0.9ulp.
 *	Maximum observed in 1,000,000 trials was .87ulp.
 *	x < 0:
 *	Maximum observed error < 4ulp in 1,000,000 trials.
 */

static double neg_gam(double);
static double small_gam(double);
static double smaller_gam(double);
static struct Double large_gam(double);
static struct Double ratfun_gam(double, double);

/*
 * Rational approximation, A0 + x*x*P(x)/Q(x), on the interval
 * [1.066.., 2.066..] accurate to 4.25e-19.
 */
#define LEFT -.3955078125	/* left boundary for rat. approx */
#define x0 .461632144968362356785	/* xmin - 1 */

#define a0_hi 0.88560319441088874992
#define a0_lo -.00000000000000004996427036469019695
#define P0	 6.21389571821820863029017800727e-01
#define P1	 2.65757198651533466104979197553e-01
#define P2	 5.53859446429917461063308081748e-03
#define P3	 1.38456698304096573887145282811e-03
#define P4	 2.40659950032711365819348969808e-03
#define Q0	 1.45019531250000000000000000000e+00
#define Q1	 1.06258521948016171343454061571e+00
#define Q2	-2.07474561943859936441469926649e-01
#define Q3	-1.46734131782005422506287573015e-01
#define Q4	 3.07878176156175520361557573779e-02
#define Q5	 5.12449347980666221336054633184e-03
#define Q6	-1.76012741431666995019222898833e-03
#define Q7	 9.35021023573788935372153030556e-05
#define Q8	 6.13275507472443958924745652239e-06
/*
 * Constants for large x approximation (x in [6, Inf])
 * (Accurate to 2.8*10^-19 absolute)
 */
#define lns2pi_hi 0.418945312500000
#define lns2pi_lo -.000006779295327258219670263595
#define Pa0	 8.33333333333333148296162562474e-02
#define Pa1	-2.77777777774548123579378966497e-03
#define Pa2	 7.93650778754435631476282786423e-04
#define Pa3	-5.95235082566672847950717262222e-04
#define Pa4	 8.41428560346653702135821806252e-04
#define Pa5	-1.89773526463879200348872089421e-03
#define Pa6	 5.69394463439411649408050664078e-03
#define Pa7	-1.44705562421428915453880392761e-02

static const double zero = 0., one = 1.0, tiny = 1e-300;

double
tgamma(double x)
{
	struct Double u;

	if (x >= 6) {
		if(x > 171.63)
			return(x/zero);
		u = large_gam(x);
		return(__exp__D(u.a, u.b));
	} else if (x >= 1.0 + LEFT + x0)
		return (small_gam(x));
	else if (x > 1.e-17)
		return (smaller_gam(x));
	else if (x > -1.e-17) {
		if (x != 0.0)
			u.a = one - tiny;	/* raise inexact */
		return (one/x);
	} else if (!isfinite(x)) {
		return (x - x);			/* x = NaN, -Inf */
	 } else
		return (neg_gam(x));
}
DEF_STD(tgamma);
LDBL_MAYBE_UNUSED_CLONE(tgamma);

/*
 * We simply call tgamma() rather than bloating the math library
 * with a float-optimized version of it.  The reason is that tgammaf()
 * is essentially useless, since the function is superexponential
 * and floats have very limited range.  -- das@freebsd.org
 */

float
tgammaf(float x)
{
	return tgamma(x);
}

/*
 * Accurate to max(ulp(1/128) absolute, 2^-66 relative) error.
 */

static struct Double
large_gam(double x)
{
	double z, p;
	struct Double t, u, v;

	z = one/(x*x);
	p = Pa0+z*(Pa1+z*(Pa2+z*(Pa3+z*(Pa4+z*(Pa5+z*(Pa6+z*Pa7))))));
	p = p/x;

	u = __log__D(x);
	u.a -= one;
	v.a = (x -= .5);
	TRUNC(v.a);
	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;
	/* return t.a + t.b + lns2pi_hi + lns2pi_lo + p */
	t.b += lns2pi_lo; t.b += p;
	u.a = lns2pi_hi + t.b; u.a += t.a;
	u.b = t.a - u.a;
	u.b += lns2pi_hi; u.b += t.b;
	return (u);
}

/*
 * Good to < 1 ulp.  (provably .90 ulp; .87 ulp on 1,000,000 runs.)
 * It also has correct monotonicity.
 */

static double
small_gam(double x)
{
	double y, ym1, t;
	struct Double yy, r;
	y = x - one;
	ym1 = y - one;
	if (y <= 1.0 + (LEFT + x0)) {
		yy = ratfun_gam(y - x0, 0);
		return (yy.a + yy.b);
	}
	r.a = y;
	TRUNC(r.a);
	yy.a = r.a - one;
	y = ym1;
	yy.b = r.b = y - yy.a;
	/* Argument reduction: G(x+1) = x*G(x) */
	for (ym1 = y-one; ym1 > LEFT + x0; y = ym1--, yy.a--) {
		t = r.a*yy.a;
		r.b = r.a*yy.b + y*r.b;
		r.a = t;
		TRUNC(r.a);
		r.b += (t - r.a);
	}
	/* Return r*tgamma(y). */
	yy = ratfun_gam(y - x0, 0);
	y = r.b*(yy.a + yy.b) + r.a*yy.b;
	y += yy.a*r.a;
	return (y);
}

/*
 * Good on (0, 1+x0+LEFT].  Accurate to 1ulp.
 */

static double
smaller_gam(double x)
{
	double t, d;
	struct Double r, xx;
	if (x < x0 + LEFT) {
		t = x;
		TRUNC(t);
		d = (t+x)*(x-t);
		t *= t;
		xx.a = (t + x);
		TRUNC(xx.a);
		xx.b = x - xx.a; xx.b += t; xx.b += d;
		t = (one-x0); t += x;
		d = (one-x0); d -= t; d += x;
		x = xx.a + xx.b;
	} else {
		xx.a =  x;
		TRUNC(xx.a);
		xx.b = x - xx.a;
		t = x - x0;
		d = (-x0 -t); d += x;
	}
	r = ratfun_gam(t, d);
	d = r.a/x;
	TRUNC(d);
	r.a -= d*xx.a; r.a -= d*xx.b; r.a += r.b;
	return (d + r.a/x);
}

/*
 * returns (z+c)^2 * P(z)/Q(z) + a0
 */

static struct Double
ratfun_gam(double z, double c)
{
	double p, q;
	struct Double r, t;

	q = Q0 +z*(Q1+z*(Q2+z*(Q3+z*(Q4+z*(Q5+z*(Q6+z*(Q7+z*Q8)))))));
	p = P0 + z*(P1 + z*(P2 + z*(P3 + z*P4)));

	/* return r.a + r.b = a0 + (z+c)^2*p/q, with r.a truncated to 26 bits. */
	p = p/q;
	t.a = z;
	TRUNC(t.a);			/* t ~= z + c */
	t.b = (z - t.a) + c;
	t.b *= (t.a + z);
	q = (t.a *= t.a);		/* t = (z+c)^2 */
	TRUNC(t.a);
	t.b += (q - t.a);
	r.a = p;
	TRUNC(r.a);			/* r = P/Q */
	r.b = p - r.a;
	t.b = t.b*p + t.a*r.b + a0_lo;
	t.a *= r.a;			/* t = (z+c)^2*(P/Q) */
	r.a = t.a + a0_hi;
	TRUNC(r.a);
	r.b = ((a0_hi-r.a) + t.a) + t.b;
	return (r);			/* r = a0 + t */
}

static double
neg_gam(double x)
{
	int sgn = 1;
	struct Double lg, lsine;
	double y, z;

	y = ceil(x);
	if (y == x)		/* Negative integer. */
		return ((x - x) / zero);
	z = y - x;
	if (z > 0.5)
		z = one - z;
	y = 0.5 * y;
	if (y == ceil(y))
		sgn = -1;
	if (z < .25)
		z = sin(M_PI*z);
	else
		z = cos(M_PI*(0.5-z));
	/* Special case: G(1-x) = Inf; G(x) may be nonzero. */
	if (x < -170) {
		if (x < -190)
			return ((double)sgn*tiny*tiny);
		y = one - x;		/* exact: 128 < |x| < 255 */
		lg = large_gam(y);
		lsine = __log__D(M_PI/z);	/* = TRUNC(log(u)) + small */
		lg.a -= lsine.a;		/* exact (opposite signs) */
		lg.b -= lsine.b;
		y = -(lg.a + lg.b);
		z = (y + lg.a) + lg.b;
		y = __exp__D(y, z);
		if (sgn < 0) y = -y;
		return (y);
	}
	y = one-x;
	if (one-y == x)
		y = tgamma(y);
	else		/* 1-x is inexact */
		y = -x*tgamma(-x);
	if (sgn < 0) y = -y;
	return (M_PI / (y*z));
}


Generated by: GCOVR (Version 3.3)

Line	Branch	Exec	Source
1			/* $OpenBSD: b_tgamma.c,v 1.10 2016/09/12 19:47:02 guenther Exp $ */
2			/*-
3			* Copyright (c) 1992, 1993
4			* The Regents of the University of California. All rights reserved.
5			*
6			* Redistribution and use in source and binary forms, with or without
7			* modification, are permitted provided that the following conditions
8			* are met:
9			* 1. Redistributions of source code must retain the above copyright
10			* notice, this list of conditions and the following disclaimer.
11			* 2. Redistributions in binary form must reproduce the above copyright
12			* notice, this list of conditions and the following disclaimer in the
13			* documentation and/or other materials provided with the distribution.
14			* 3. Neither the name of the University nor the names of its contributors
15			* may be used to endorse or promote products derived from this software
16			* without specific prior written permission.
17			*
18			* THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
19			* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
20			* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
21			* ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
22			* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
23			* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
24			* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
25			* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
26			* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
27			* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
28			* SUCH DAMAGE.
29			*/
30
31			/*
32			* This code by P. McIlroy, Oct 1992;
33			*
34			* The financial support of UUNET Communications Services is greatfully
35			* acknowledged.
36			*/
37
38			#include <float.h>
39			#include <math.h>
40
41			#include "math_private.h"
42
43			/* METHOD:
44			* x < 0: Use reflection formula, G(x) = pi/(sin(pix)x*G(x))
45			* At negative integers, return NaN and raise invalid.
46			*
47			* x < 6.5:
48			* Use argument reduction G(x+1) = xG(x) to reach the
49			* range [1.066124,2.066124]. Use a rational
50			* approximation centered at the minimum (x0+1) to
51			* ensure monotonicity.
52			*
53			* x >= 6.5: Use the asymptotic approximation (Stirling's formula)
54			* adjusted for equal-ripples:
55			*
56			* log(G(x)) ~= (x-.5)(log(x)-1) + .5(log(2pi)-1) + 1/xP(1/(xx))
57			*
58			* Keep extra precision in multiplying (x-.5)(log(x)-1), to
59			* avoid premature round-off.
60			*
61			* Special values:
62			* -Inf: return NaN and raise invalid;
63			* negative integer: return NaN and raise invalid;
64			* other x ~< -177.79: return +-0 and raise underflow;
65			* +-0: return +-Inf and raise divide-by-zero;
66			* finite x ~> 171.63: return +Inf and raise overflow;
67			* +Inf: return +Inf;
68			* NaN: return NaN.
69			*
70			* Accuracy: tgamma(x) is accurate to within
71			* x > 0: error provably < 0.9ulp.
72			* Maximum observed in 1,000,000 trials was .87ulp.
73			* x < 0:
74			* Maximum observed error < 4ulp in 1,000,000 trials.
75			*/
76
77			static double neg_gam(double);
78			static double small_gam(double);
79			static double smaller_gam(double);
80			static struct Double large_gam(double);
81			static struct Double ratfun_gam(double, double);
82
83			/*
84			* Rational approximation, A0 + xxP(x)/Q(x), on the interval
85			* [1.066.., 2.066..] accurate to 4.25e-19.
86			*/
87			#define LEFT -.3955078125 /* left boundary for rat. approx */
88			#define x0 .461632144968362356785 /* xmin - 1 */
89
90			#define a0_hi 0.88560319441088874992
91			#define a0_lo -.00000000000000004996427036469019695
92			#define P0 6.21389571821820863029017800727e-01
93			#define P1 2.65757198651533466104979197553e-01
94			#define P2 5.53859446429917461063308081748e-03
95			#define P3 1.38456698304096573887145282811e-03
96			#define P4 2.40659950032711365819348969808e-03
97			#define Q0 1.45019531250000000000000000000e+00
98			#define Q1 1.06258521948016171343454061571e+00
99			#define Q2 -2.07474561943859936441469926649e-01
100			#define Q3 -1.46734131782005422506287573015e-01
101			#define Q4 3.07878176156175520361557573779e-02
102			#define Q5 5.12449347980666221336054633184e-03
103			#define Q6 -1.76012741431666995019222898833e-03
104			#define Q7 9.35021023573788935372153030556e-05
105			#define Q8 6.13275507472443958924745652239e-06
106			/*
107			* Constants for large x approximation (x in [6, Inf])
108			* (Accurate to 2.8*10^-19 absolute)
109			*/
110			#define lns2pi_hi 0.418945312500000
111			#define lns2pi_lo -.000006779295327258219670263595
112			#define Pa0 8.33333333333333148296162562474e-02
113			#define Pa1 -2.77777777774548123579378966497e-03
114			#define Pa2 7.93650778754435631476282786423e-04
115			#define Pa3 -5.95235082566672847950717262222e-04
116			#define Pa4 8.41428560346653702135821806252e-04
117			#define Pa5 -1.89773526463879200348872089421e-03
118			#define Pa6 5.69394463439411649408050664078e-03
119			#define Pa7 -1.44705562421428915453880392761e-02
120
121			static const double zero = 0., one = 1.0, tiny = 1e-300;
122
123			double
124			tgamma(double x)
125			{
126			struct Double u;
127
128	✓✓	72	if (x >= 6) {
129	✓✓	9	if(x > 171.63)
130		6	return(x/zero);
131		3	u = large_gam(x);
132		3	return(__exp__D(u.a, u.b));
133	✓✓	27	} else if (x >= 1.0 + LEFT + x0)
134		6	return (small_gam(x));
135	✗✓	21	else if (x > 1.e-17)
136			return (smaller_gam(x));
137	✓✓	21	else if (x > -1.e-17) {
138	✗✓	6	if (x != 0.0)
139			u.a = one - tiny; /* raise inexact */
140		6	return (one/x);
141	✓✓	15	} else if (!isfinite(x)) {
142		6	return (x - x); /* x = NaN, -Inf */
143			} else
144		9	return (neg_gam(x));
145		36	}
146			DEF_STD(tgamma);
147			LDBL_MAYBE_UNUSED_CLONE(tgamma);
148
149			/*
150			* We simply call tgamma() rather than bloating the math library
151			* with a float-optimized version of it. The reason is that tgammaf()
152			* is essentially useless, since the function is superexponential
153			* and floats have very limited range. -- das@freebsd.org
154			*/
155
156			float
157			tgammaf(float x)
158			{
159			return tgamma(x);
160			}
161
162			/*
163			* Accurate to max(ulp(1/128) absolute, 2^-66 relative) error.
164			*/
165
166			static struct Double
167			large_gam(double x)
168			{
169			double z, p;
170		12	struct Double t, u, v;
171
172		6	z = one/(x*x);
173		6	p = Pa0+z(Pa1+z(Pa2+z(Pa3+z(Pa4+z(Pa5+z(Pa6+z*Pa7))))));
174		6	p = p/x;
175
176		6	u = __log__D(x);
177		6	u.a -= one;
178		6	v.a = (x -= .5);
179		6	TRUNC(v.a);
180		6	v.b = x - v.a;
181		6	t.a = v.au.a; / t = (x-.5)(log(x)-1) /
182		6	t.b = v.bu.a + xu.b;
183			/* return t.a + t.b + lns2pi_hi + lns2pi_lo + p */
184		6	t.b += lns2pi_lo; t.b += p;
185		6	u.a = lns2pi_hi + t.b; u.a += t.a;
186		6	u.b = t.a - u.a;
187		6	u.b += lns2pi_hi; u.b += t.b;
188			return (u);
189		6	}
190
191			/*
192			* Good to < 1 ulp. (provably .90 ulp; .87 ulp on 1,000,000 runs.)
193			* It also has correct monotonicity.
194			*/
195
196			static double
197			small_gam(double x)
198			{
199			double y, ym1, t;
200		12	struct Double yy, r;
201		6	y = x - one;
202		6	ym1 = y - one;
203	✓✓	6	if (y <= 1.0 + (LEFT + x0)) {
204		3	yy = ratfun_gam(y - x0, 0);
205		3	return (yy.a + yy.b);
206			}
207		3	r.a = y;
208		3	TRUNC(r.a);
209		3	yy.a = r.a - one;
210			y = ym1;
211		3	yy.b = r.b = y - yy.a;
212			/* Argument reduction: G(x+1) = xG(x) /
213	✓✓	12	for (ym1 = y-one; ym1 > LEFT + x0; y = ym1--, yy.a--) {
214		3	t = r.a*yy.a;
215		3	r.b = r.ayy.b + yr.b;
216		3	r.a = t;
217		3	TRUNC(r.a);
218		3	r.b += (t - r.a);
219			}
220			/* Return rtgamma(y). /
221		3	yy = ratfun_gam(y - x0, 0);
222		3	y = r.b(yy.a + yy.b) + r.ayy.b;
223		3	y += yy.a*r.a;
224		3	return (y);
225		6	}
226
227			/*
228			* Good on (0, 1+x0+LEFT]. Accurate to 1ulp.
229			*/
230
231			static double
232			smaller_gam(double x)
233			{
234			double t, d;
235			struct Double r, xx;
236			if (x < x0 + LEFT) {
237			t = x;
238			TRUNC(t);
239			d = (t+x)*(x-t);
240			t *= t;
241			xx.a = (t + x);
242			TRUNC(xx.a);
243			xx.b = x - xx.a; xx.b += t; xx.b += d;
244			t = (one-x0); t += x;
245			d = (one-x0); d -= t; d += x;
246			x = xx.a + xx.b;
247			} else {
248			xx.a = x;
249			TRUNC(xx.a);
250			xx.b = x - xx.a;
251			t = x - x0;
252			d = (-x0 -t); d += x;
253			}
254			r = ratfun_gam(t, d);
255			d = r.a/x;
256			TRUNC(d);
257			r.a -= dxx.a; r.a -= dxx.b; r.a += r.b;
258			return (d + r.a/x);
259			}
260
261			/*
262			* returns (z+c)^2 * P(z)/Q(z) + a0
263			*/
264
265			static struct Double
266			ratfun_gam(double z, double c)
267			{
268			double p, q;
269		12	struct Double r, t;
270
271		6	q = Q0 +z(Q1+z(Q2+z(Q3+z(Q4+z(Q5+z(Q6+z(Q7+zQ8)))))));
272		6	p = P0 + z(P1 + z(P2 + z(P3 + zP4)));
273
274			/* return r.a + r.b = a0 + (z+c)^2p/q, with r.a truncated to 26 bits. /
275		6	p = p/q;
276		6	t.a = z;
277		6	TRUNC(t.a); /* t ~= z + c */
278		6	t.b = (z - t.a) + c;
279		6	t.b *= (t.a + z);
280		6	q = (t.a = t.a); / t = (z+c)^2 */
281		6	TRUNC(t.a);
282		6	t.b += (q - t.a);
283		6	r.a = p;
284		6	TRUNC(r.a); /* r = P/Q */
285		6	r.b = p - r.a;
286		6	t.b = t.bp + t.ar.b + a0_lo;
287		6	t.a = r.a; / t = (z+c)^2(P/Q) /
288		6	r.a = t.a + a0_hi;
289		6	TRUNC(r.a);
290		6	r.b = ((a0_hi-r.a) + t.a) + t.b;
291		6	return (r); /* r = a0 + t */
292		6	}
293
294			static double
295			neg_gam(double x)
296			{
297			int sgn = 1;
298			struct Double lg, lsine;
299			double y, z;
300
301		18	y = ceil(x);
302	✓✓	9	if (y == x) /* Negative integer. */
303		3	return ((x - x) / zero);
304		6	z = y - x;
305	✓✓	6	if (z > 0.5)
306		3	z = one - z;
307		6	y = 0.5 * y;
308	✓✓	6	if (y == ceil(y))
309		3	sgn = -1;
310	✓✓	6	if (z < .25)
311		3	z = sin(M_PI*z);
312			else
313		3	z = cos(M_PI*(0.5-z));
314			/* Special case: G(1-x) = Inf; G(x) may be nonzero. */
315	✓✓	6	if (x < -170) {
316	✗✓	3	if (x < -190)
317			return ((double)sgntinytiny);
318		3	y = one - x; /* exact: 128 < \|x\| < 255 */
319		3	lg = large_gam(y);
320		3	lsine = __log__D(M_PI/z); /* = TRUNC(log(u)) + small */
321		3	lg.a -= lsine.a; /* exact (opposite signs) */
322		3	lg.b -= lsine.b;
323		3	y = -(lg.a + lg.b);
324		3	z = (y + lg.a) + lg.b;
325		3	y = __exp__D(y, z);
326	✗✓	3	if (sgn < 0) y = -y;
327		3	return (y);
328			}
329		3	y = one-x;
330	✓✗	3	if (one-y == x)
331		3	y = tgamma(y);
332			else /* 1-x is inexact */
333			y = -x*tgamma(-x);
334	✓✗	6	if (sgn < 0) y = -y;
335		3	return (M_PI / (y*z));
336		9	}