11dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project/*- 21dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * Copyright (c) 1992, 1993 31dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * The Regents of the University of California. All rights reserved. 41dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * 51dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * Redistribution and use in source and binary forms, with or without 61dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * modification, are permitted provided that the following conditions 71dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * are met: 81dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * 1. Redistributions of source code must retain the above copyright 91dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * notice, this list of conditions and the following disclaimer. 101dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * 2. Redistributions in binary form must reproduce the above copyright 111dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * notice, this list of conditions and the following disclaimer in the 121dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * documentation and/or other materials provided with the distribution. 131dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * 3. All advertising materials mentioning features or use of this software 141dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * must display the following acknowledgement: 151dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * This product includes software developed by the University of 161dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * California, Berkeley and its contributors. 171dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * 4. Neither the name of the University nor the names of its contributors 181dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * may be used to endorse or promote products derived from this software 191dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * without specific prior written permission. 201dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * 211dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND 221dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 231dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 241dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE 251dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 261dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS 271dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 281dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 291dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 301dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 311dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * SUCH DAMAGE. 321dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project */ 331dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 34a0ee07829a9ba7e99ef68e8c12551301cc797f0fElliott Hughes/* @(#)gamma.c 8.1 (Berkeley) 6/4/93 */ 351dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project#include <sys/cdefs.h> 36a0ee07829a9ba7e99ef68e8c12551301cc797f0fElliott Hughes__FBSDID("$FreeBSD$"); 371dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 381dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project/* 391dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * This code by P. McIlroy, Oct 1992; 401dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * 411dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * The financial support of UUNET Communications Services is greatfully 421dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * acknowledged. 431dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project */ 441dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 45a0ee07829a9ba7e99ef68e8c12551301cc797f0fElliott Hughes#include <math.h> 461dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project#include "mathimpl.h" 471dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 481dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project/* METHOD: 491dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * x < 0: Use reflection formula, G(x) = pi/(sin(pi*x)*x*G(x)) 50a0ee07829a9ba7e99ef68e8c12551301cc797f0fElliott Hughes * At negative integers, return NaN and raise invalid. 511dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * 521dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * x < 6.5: 531dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * Use argument reduction G(x+1) = xG(x) to reach the 541dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * range [1.066124,2.066124]. Use a rational 551dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * approximation centered at the minimum (x0+1) to 561dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * ensure monotonicity. 571dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * 581dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * x >= 6.5: Use the asymptotic approximation (Stirling's formula) 591dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * adjusted for equal-ripples: 601dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * 611dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * log(G(x)) ~= (x-.5)*(log(x)-1) + .5(log(2*pi)-1) + 1/x*P(1/(x*x)) 621dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * 631dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * Keep extra precision in multiplying (x-.5)(log(x)-1), to 641dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * avoid premature round-off. 651dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * 661dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * Special values: 67a0ee07829a9ba7e99ef68e8c12551301cc797f0fElliott Hughes * -Inf: return NaN and raise invalid; 68a0ee07829a9ba7e99ef68e8c12551301cc797f0fElliott Hughes * negative integer: return NaN and raise invalid; 69a0ee07829a9ba7e99ef68e8c12551301cc797f0fElliott Hughes * other x ~< 177.79: return +-0 and raise underflow; 70a0ee07829a9ba7e99ef68e8c12551301cc797f0fElliott Hughes * +-0: return +-Inf and raise divide-by-zero; 71a0ee07829a9ba7e99ef68e8c12551301cc797f0fElliott Hughes * finite x ~> 171.63: return +Inf and raise overflow; 72a0ee07829a9ba7e99ef68e8c12551301cc797f0fElliott Hughes * +Inf: return +Inf; 73a0ee07829a9ba7e99ef68e8c12551301cc797f0fElliott Hughes * NaN: return NaN. 741dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * 75a0ee07829a9ba7e99ef68e8c12551301cc797f0fElliott Hughes * Accuracy: tgamma(x) is accurate to within 761dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * x > 0: error provably < 0.9ulp. 771dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * Maximum observed in 1,000,000 trials was .87ulp. 781dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * x < 0: 791dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * Maximum observed error < 4ulp in 1,000,000 trials. 801dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project */ 811dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 821dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Projectstatic double neg_gam(double); 831dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Projectstatic double small_gam(double); 841dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Projectstatic double smaller_gam(double); 851dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Projectstatic struct Double large_gam(double); 861dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Projectstatic struct Double ratfun_gam(double, double); 871dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 881dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project/* 891dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * Rational approximation, A0 + x*x*P(x)/Q(x), on the interval 901dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * [1.066.., 2.066..] accurate to 4.25e-19. 911dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project */ 921dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project#define LEFT -.3955078125 /* left boundary for rat. approx */ 931dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project#define x0 .461632144968362356785 /* xmin - 1 */ 941dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 951dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project#define a0_hi 0.88560319441088874992 961dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project#define a0_lo -.00000000000000004996427036469019695 971dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project#define P0 6.21389571821820863029017800727e-01 981dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project#define P1 2.65757198651533466104979197553e-01 991dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project#define P2 5.53859446429917461063308081748e-03 1001dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project#define P3 1.38456698304096573887145282811e-03 1011dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project#define P4 2.40659950032711365819348969808e-03 1021dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project#define Q0 1.45019531250000000000000000000e+00 1031dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project#define Q1 1.06258521948016171343454061571e+00 1041dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project#define Q2 -2.07474561943859936441469926649e-01 1051dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project#define Q3 -1.46734131782005422506287573015e-01 1061dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project#define Q4 3.07878176156175520361557573779e-02 1071dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project#define Q5 5.12449347980666221336054633184e-03 1081dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project#define Q6 -1.76012741431666995019222898833e-03 1091dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project#define Q7 9.35021023573788935372153030556e-05 1101dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project#define Q8 6.13275507472443958924745652239e-06 1111dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project/* 1121dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * Constants for large x approximation (x in [6, Inf]) 1131dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * (Accurate to 2.8*10^-19 absolute) 1141dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project */ 1151dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project#define lns2pi_hi 0.418945312500000 1161dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project#define lns2pi_lo -.000006779295327258219670263595 1171dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project#define Pa0 8.33333333333333148296162562474e-02 1181dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project#define Pa1 -2.77777777774548123579378966497e-03 1191dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project#define Pa2 7.93650778754435631476282786423e-04 1201dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project#define Pa3 -5.95235082566672847950717262222e-04 1211dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project#define Pa4 8.41428560346653702135821806252e-04 1221dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project#define Pa5 -1.89773526463879200348872089421e-03 1231dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project#define Pa6 5.69394463439411649408050664078e-03 1241dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project#define Pa7 -1.44705562421428915453880392761e-02 1251dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 1261dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Projectstatic const double zero = 0., one = 1.0, tiny = 1e-300; 1271dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 1281dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Projectdouble 1291dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Projecttgamma(x) 1301dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project double x; 1311dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project{ 1321dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project struct Double u; 1331dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 1341dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if (x >= 6) { 1351dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if(x > 171.63) 136a0ee07829a9ba7e99ef68e8c12551301cc797f0fElliott Hughes return (x / zero); 1371dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project u = large_gam(x); 1381dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project return(__exp__D(u.a, u.b)); 1391dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } else if (x >= 1.0 + LEFT + x0) 1401dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project return (small_gam(x)); 1411dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project else if (x > 1.e-17) 1421dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project return (smaller_gam(x)); 1431dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project else if (x > -1.e-17) { 144a0ee07829a9ba7e99ef68e8c12551301cc797f0fElliott Hughes if (x != 0.0) 145a0ee07829a9ba7e99ef68e8c12551301cc797f0fElliott Hughes u.a = one - tiny; /* raise inexact */ 1461dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project return (one/x); 1471dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } else if (!finite(x)) 148a0ee07829a9ba7e99ef68e8c12551301cc797f0fElliott Hughes return (x - x); /* x is NaN or -Inf */ 1491dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project else 1501dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project return (neg_gam(x)); 1511dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project} 1521dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project/* 1531dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * Accurate to max(ulp(1/128) absolute, 2^-66 relative) error. 1541dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project */ 1551dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Projectstatic struct Double 1561dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Projectlarge_gam(x) 1571dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project double x; 1581dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project{ 1591dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project double z, p; 1601dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project struct Double t, u, v; 1611dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 1621dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project z = one/(x*x); 1631dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project p = Pa0+z*(Pa1+z*(Pa2+z*(Pa3+z*(Pa4+z*(Pa5+z*(Pa6+z*Pa7)))))); 1641dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project p = p/x; 1651dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 1661dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project u = __log__D(x); 1671dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project u.a -= one; 1681dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project v.a = (x -= .5); 1691dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project TRUNC(v.a); 1701dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project v.b = x - v.a; 1711dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project t.a = v.a*u.a; /* t = (x-.5)*(log(x)-1) */ 1721dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project t.b = v.b*u.a + x*u.b; 1731dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project /* return t.a + t.b + lns2pi_hi + lns2pi_lo + p */ 1741dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project t.b += lns2pi_lo; t.b += p; 1751dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project u.a = lns2pi_hi + t.b; u.a += t.a; 1761dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project u.b = t.a - u.a; 1771dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project u.b += lns2pi_hi; u.b += t.b; 1781dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project return (u); 1791dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project} 1801dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project/* 1811dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * Good to < 1 ulp. (provably .90 ulp; .87 ulp on 1,000,000 runs.) 1821dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * It also has correct monotonicity. 1831dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project */ 1841dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Projectstatic double 1851dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Projectsmall_gam(x) 1861dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project double x; 1871dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project{ 1881dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project double y, ym1, t; 1891dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project struct Double yy, r; 1901dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project y = x - one; 1911dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project ym1 = y - one; 1921dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if (y <= 1.0 + (LEFT + x0)) { 1931dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project yy = ratfun_gam(y - x0, 0); 1941dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project return (yy.a + yy.b); 1951dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } 1961dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project r.a = y; 1971dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project TRUNC(r.a); 1981dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project yy.a = r.a - one; 1991dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project y = ym1; 2001dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project yy.b = r.b = y - yy.a; 2011dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project /* Argument reduction: G(x+1) = x*G(x) */ 2021dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project for (ym1 = y-one; ym1 > LEFT + x0; y = ym1--, yy.a--) { 2031dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project t = r.a*yy.a; 2041dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project r.b = r.a*yy.b + y*r.b; 2051dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project r.a = t; 2061dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project TRUNC(r.a); 2071dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project r.b += (t - r.a); 2081dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } 2091dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project /* Return r*tgamma(y). */ 2101dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project yy = ratfun_gam(y - x0, 0); 2111dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project y = r.b*(yy.a + yy.b) + r.a*yy.b; 2121dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project y += yy.a*r.a; 2131dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project return (y); 2141dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project} 2151dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project/* 2161dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * Good on (0, 1+x0+LEFT]. Accurate to 1ulp. 2171dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project */ 2181dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Projectstatic double 2191dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Projectsmaller_gam(x) 2201dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project double x; 2211dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project{ 2221dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project double t, d; 2231dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project struct Double r, xx; 2241dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if (x < x0 + LEFT) { 2251dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project t = x, TRUNC(t); 2261dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project d = (t+x)*(x-t); 2271dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project t *= t; 2281dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project xx.a = (t + x), TRUNC(xx.a); 2291dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project xx.b = x - xx.a; xx.b += t; xx.b += d; 2301dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project t = (one-x0); t += x; 2311dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project d = (one-x0); d -= t; d += x; 2321dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project x = xx.a + xx.b; 2331dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } else { 2341dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project xx.a = x, TRUNC(xx.a); 2351dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project xx.b = x - xx.a; 2361dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project t = x - x0; 2371dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project d = (-x0 -t); d += x; 2381dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } 2391dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project r = ratfun_gam(t, d); 2401dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project d = r.a/x, TRUNC(d); 2411dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project r.a -= d*xx.a; r.a -= d*xx.b; r.a += r.b; 2421dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project return (d + r.a/x); 2431dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project} 2441dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project/* 2451dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project * returns (z+c)^2 * P(z)/Q(z) + a0 2461dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project */ 2471dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Projectstatic struct Double 2481dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Projectratfun_gam(z, c) 2491dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project double z, c; 2501dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project{ 2511dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project double p, q; 2521dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project struct Double r, t; 2531dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 2541dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project q = Q0 +z*(Q1+z*(Q2+z*(Q3+z*(Q4+z*(Q5+z*(Q6+z*(Q7+z*Q8))))))); 2551dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project p = P0 + z*(P1 + z*(P2 + z*(P3 + z*P4))); 2561dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 2571dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project /* return r.a + r.b = a0 + (z+c)^2*p/q, with r.a truncated to 26 bits. */ 2581dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project p = p/q; 2591dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project t.a = z, TRUNC(t.a); /* t ~= z + c */ 2601dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project t.b = (z - t.a) + c; 2611dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project t.b *= (t.a + z); 2621dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project q = (t.a *= t.a); /* t = (z+c)^2 */ 2631dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project TRUNC(t.a); 2641dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project t.b += (q - t.a); 2651dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project r.a = p, TRUNC(r.a); /* r = P/Q */ 2661dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project r.b = p - r.a; 2671dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project t.b = t.b*p + t.a*r.b + a0_lo; 2681dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project t.a *= r.a; /* t = (z+c)^2*(P/Q) */ 2691dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project r.a = t.a + a0_hi, TRUNC(r.a); 2701dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project r.b = ((a0_hi-r.a) + t.a) + t.b; 2711dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project return (r); /* r = a0 + t */ 2721dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project} 2731dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 2741dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Projectstatic double 2751dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Projectneg_gam(x) 2761dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project double x; 2771dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project{ 2781dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project int sgn = 1; 2791dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project struct Double lg, lsine; 2801dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project double y, z; 2811dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project 282a0ee07829a9ba7e99ef68e8c12551301cc797f0fElliott Hughes y = ceil(x); 2831dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if (y == x) /* Negative integer. */ 284a0ee07829a9ba7e99ef68e8c12551301cc797f0fElliott Hughes return ((x - x) / zero); 285a0ee07829a9ba7e99ef68e8c12551301cc797f0fElliott Hughes z = y - x; 286a0ee07829a9ba7e99ef68e8c12551301cc797f0fElliott Hughes if (z > 0.5) 287a0ee07829a9ba7e99ef68e8c12551301cc797f0fElliott Hughes z = one - z; 288a0ee07829a9ba7e99ef68e8c12551301cc797f0fElliott Hughes y = 0.5 * y; 2891dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if (y == ceil(y)) 2901dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project sgn = -1; 2911dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if (z < .25) 2921dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project z = sin(M_PI*z); 2931dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project else 2941dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project z = cos(M_PI*(0.5-z)); 2951dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project /* Special case: G(1-x) = Inf; G(x) may be nonzero. */ 2961dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if (x < -170) { 2971dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if (x < -190) 2981dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project return ((double)sgn*tiny*tiny); 2991dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project y = one - x; /* exact: 128 < |x| < 255 */ 3001dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project lg = large_gam(y); 3011dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project lsine = __log__D(M_PI/z); /* = TRUNC(log(u)) + small */ 3021dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project lg.a -= lsine.a; /* exact (opposite signs) */ 3031dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project lg.b -= lsine.b; 3041dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project y = -(lg.a + lg.b); 3051dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project z = (y + lg.a) + lg.b; 3061dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project y = __exp__D(y, z); 3071dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if (sgn < 0) y = -y; 3081dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project return (y); 3091dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project } 3101dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project y = one-x; 3111dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if (one-y == x) 3121dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project y = tgamma(y); 3131dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project else /* 1-x is inexact */ 3141dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project y = -x*tgamma(-x); 3151dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project if (sgn < 0) y = -y; 3161dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project return (M_PI / (y*z)); 3171dc9e472e19acfe6dc7f41e429236e7eef7ceda1The Android Open Source Project} 318