1*25c28e83SPiotr Jasiukajtis /*
2*25c28e83SPiotr Jasiukajtis * CDDL HEADER START
3*25c28e83SPiotr Jasiukajtis *
4*25c28e83SPiotr Jasiukajtis * The contents of this file are subject to the terms of the
5*25c28e83SPiotr Jasiukajtis * Common Development and Distribution License (the "License").
6*25c28e83SPiotr Jasiukajtis * You may not use this file except in compliance with the License.
7*25c28e83SPiotr Jasiukajtis *
8*25c28e83SPiotr Jasiukajtis * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9*25c28e83SPiotr Jasiukajtis * or http://www.opensolaris.org/os/licensing.
10*25c28e83SPiotr Jasiukajtis * See the License for the specific language governing permissions
11*25c28e83SPiotr Jasiukajtis * and limitations under the License.
12*25c28e83SPiotr Jasiukajtis *
13*25c28e83SPiotr Jasiukajtis * When distributing Covered Code, include this CDDL HEADER in each
14*25c28e83SPiotr Jasiukajtis * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
15*25c28e83SPiotr Jasiukajtis * If applicable, add the following below this CDDL HEADER, with the
16*25c28e83SPiotr Jasiukajtis * fields enclosed by brackets "[]" replaced with your own identifying
17*25c28e83SPiotr Jasiukajtis * information: Portions Copyright [yyyy] [name of copyright owner]
18*25c28e83SPiotr Jasiukajtis *
19*25c28e83SPiotr Jasiukajtis * CDDL HEADER END
20*25c28e83SPiotr Jasiukajtis */
21*25c28e83SPiotr Jasiukajtis
22*25c28e83SPiotr Jasiukajtis /*
23*25c28e83SPiotr Jasiukajtis * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
24*25c28e83SPiotr Jasiukajtis */
25*25c28e83SPiotr Jasiukajtis /*
26*25c28e83SPiotr Jasiukajtis * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
27*25c28e83SPiotr Jasiukajtis * Use is subject to license terms.
28*25c28e83SPiotr Jasiukajtis */
29*25c28e83SPiotr Jasiukajtis
30*25c28e83SPiotr Jasiukajtis /*
31*25c28e83SPiotr Jasiukajtis * long double __k_lgammal(long double x, int *signgamlp);
32*25c28e83SPiotr Jasiukajtis * K.C. Ng, August, 1989.
33*25c28e83SPiotr Jasiukajtis *
34*25c28e83SPiotr Jasiukajtis * We choose [1.5,2.5] to be the primary interval. Our algorithms
35*25c28e83SPiotr Jasiukajtis * are mainly derived from
36*25c28e83SPiotr Jasiukajtis *
37*25c28e83SPiotr Jasiukajtis *
38*25c28e83SPiotr Jasiukajtis * zeta(2)-1 2 zeta(3)-1 3
39*25c28e83SPiotr Jasiukajtis * lgamma(2+s) = s*(1-euler) + --------- * s - --------- * s + ...
40*25c28e83SPiotr Jasiukajtis * 2 3
41*25c28e83SPiotr Jasiukajtis *
42*25c28e83SPiotr Jasiukajtis *
43*25c28e83SPiotr Jasiukajtis * Note 1. Since gamma(1+s)=s*gamma(s), hence
44*25c28e83SPiotr Jasiukajtis * lgamma(1+s) = log(s) + lgamma(s), or
45*25c28e83SPiotr Jasiukajtis * lgamma(s) = lgamma(1+s) - log(s).
46*25c28e83SPiotr Jasiukajtis * When s is really tiny (like roundoff), lgamma(1+s) ~ s(1-enler)
47*25c28e83SPiotr Jasiukajtis * Hence lgamma(s) ~ -log(s) for tiny s
48*25c28e83SPiotr Jasiukajtis *
49*25c28e83SPiotr Jasiukajtis */
50*25c28e83SPiotr Jasiukajtis
51*25c28e83SPiotr Jasiukajtis #include "libm.h"
52*25c28e83SPiotr Jasiukajtis #include "longdouble.h"
53*25c28e83SPiotr Jasiukajtis
54*25c28e83SPiotr Jasiukajtis static long double neg(long double, int *);
55*25c28e83SPiotr Jasiukajtis static long double poly(long double, const long double *, int);
56*25c28e83SPiotr Jasiukajtis static long double polytail(long double);
57*25c28e83SPiotr Jasiukajtis static long double primary(long double);
58*25c28e83SPiotr Jasiukajtis
59*25c28e83SPiotr Jasiukajtis static const long double
60*25c28e83SPiotr Jasiukajtis c0 = 0.0L,
61*25c28e83SPiotr Jasiukajtis ch = 0.5L,
62*25c28e83SPiotr Jasiukajtis c1 = 1.0L,
63*25c28e83SPiotr Jasiukajtis c2 = 2.0L,
64*25c28e83SPiotr Jasiukajtis c3 = 3.0L,
65*25c28e83SPiotr Jasiukajtis c4 = 4.0L,
66*25c28e83SPiotr Jasiukajtis c5 = 5.0L,
67*25c28e83SPiotr Jasiukajtis c6 = 6.0L,
68*25c28e83SPiotr Jasiukajtis pi = 3.1415926535897932384626433832795028841971L,
69*25c28e83SPiotr Jasiukajtis tiny = 1.0e-40L;
70*25c28e83SPiotr Jasiukajtis
71*25c28e83SPiotr Jasiukajtis long double
__k_lgammal(long double x,int * signgamlp)72*25c28e83SPiotr Jasiukajtis __k_lgammal(long double x, int *signgamlp) {
73*25c28e83SPiotr Jasiukajtis long double t, y;
74*25c28e83SPiotr Jasiukajtis int i;
75*25c28e83SPiotr Jasiukajtis
76*25c28e83SPiotr Jasiukajtis /* purge off +-inf, NaN and negative arguments */
77*25c28e83SPiotr Jasiukajtis if (!finitel(x))
78*25c28e83SPiotr Jasiukajtis return (x*x);
79*25c28e83SPiotr Jasiukajtis *signgamlp = 1;
80*25c28e83SPiotr Jasiukajtis if (signbitl(x))
81*25c28e83SPiotr Jasiukajtis return (neg(x, signgamlp));
82*25c28e83SPiotr Jasiukajtis
83*25c28e83SPiotr Jasiukajtis /* for x < 8.0 */
84*25c28e83SPiotr Jasiukajtis if (x < 8.0L) {
85*25c28e83SPiotr Jasiukajtis y = anintl(x);
86*25c28e83SPiotr Jasiukajtis i = (int) y;
87*25c28e83SPiotr Jasiukajtis switch (i) {
88*25c28e83SPiotr Jasiukajtis case 0:
89*25c28e83SPiotr Jasiukajtis if (x < 1.0e-40L)
90*25c28e83SPiotr Jasiukajtis return (-logl(x));
91*25c28e83SPiotr Jasiukajtis else
92*25c28e83SPiotr Jasiukajtis return (primary(x)-log1pl(x))-logl(x);
93*25c28e83SPiotr Jasiukajtis case 1:
94*25c28e83SPiotr Jasiukajtis return (primary(x-y)-logl(x));
95*25c28e83SPiotr Jasiukajtis case 2:
96*25c28e83SPiotr Jasiukajtis return (primary(x-y));
97*25c28e83SPiotr Jasiukajtis case 3:
98*25c28e83SPiotr Jasiukajtis return (primary(x-y)+logl(x-c1));
99*25c28e83SPiotr Jasiukajtis case 4:
100*25c28e83SPiotr Jasiukajtis return (primary(x-y)+logl((x-c1)*(x-c2)));
101*25c28e83SPiotr Jasiukajtis case 5:
102*25c28e83SPiotr Jasiukajtis return (primary(x-y)+logl((x-c1)*(x-c2)*(x-c3)));
103*25c28e83SPiotr Jasiukajtis case 6:
104*25c28e83SPiotr Jasiukajtis return (primary(x-y)+logl((x-c1)*(x-c2)*(x-c3)*(x-c4)));
105*25c28e83SPiotr Jasiukajtis case 7:
106*25c28e83SPiotr Jasiukajtis return (primary(x-y)+logl((x-c1)*(x-c2)*(x-c3)*(x-c4)*(x-c5)));
107*25c28e83SPiotr Jasiukajtis case 8:
108*25c28e83SPiotr Jasiukajtis return primary(x-y)+
109*25c28e83SPiotr Jasiukajtis logl((x-c1)*(x-c2)*(x-c3)*(x-c4)*(x-c5)*(x-c6));
110*25c28e83SPiotr Jasiukajtis }
111*25c28e83SPiotr Jasiukajtis }
112*25c28e83SPiotr Jasiukajtis
113*25c28e83SPiotr Jasiukajtis /* 8.0 <= x < 1.0e40 */
114*25c28e83SPiotr Jasiukajtis if (x < 1.0e40L) {
115*25c28e83SPiotr Jasiukajtis t = logl(x);
116*25c28e83SPiotr Jasiukajtis return (x*(t-c1)-(ch*t-polytail(c1/x)));
117*25c28e83SPiotr Jasiukajtis }
118*25c28e83SPiotr Jasiukajtis
119*25c28e83SPiotr Jasiukajtis /* 1.0e40 <= x <= inf */
120*25c28e83SPiotr Jasiukajtis return (x*(logl(x)-c1));
121*25c28e83SPiotr Jasiukajtis }
122*25c28e83SPiotr Jasiukajtis
123*25c28e83SPiotr Jasiukajtis static const long double an1[] = { /* 20 terms */
124*25c28e83SPiotr Jasiukajtis -0.0772156649015328606065120900824024309741L,
125*25c28e83SPiotr Jasiukajtis 3.224670334241132182362075833230130289059e-0001L,
126*25c28e83SPiotr Jasiukajtis -6.735230105319809513324605383668929964120e-0002L,
127*25c28e83SPiotr Jasiukajtis 2.058080842778454787900092432928910226297e-0002L,
128*25c28e83SPiotr Jasiukajtis -7.385551028673985266273054086081102125704e-0003L,
129*25c28e83SPiotr Jasiukajtis 2.890510330741523285758867304409628648727e-0003L,
130*25c28e83SPiotr Jasiukajtis -1.192753911703260976581414338096267498555e-0003L,
131*25c28e83SPiotr Jasiukajtis 5.096695247430424562831956662855697824035e-0004L,
132*25c28e83SPiotr Jasiukajtis -2.231547584535777978926798502084300123638e-0004L,
133*25c28e83SPiotr Jasiukajtis 9.945751278186384670278268034322157947635e-0005L,
134*25c28e83SPiotr Jasiukajtis -4.492623673665547726647838474125147631082e-0005L,
135*25c28e83SPiotr Jasiukajtis 2.050721280617796810096993154281561168706e-0005L,
136*25c28e83SPiotr Jasiukajtis -9.439487785617396552092393234044767313568e-0006L,
137*25c28e83SPiotr Jasiukajtis 4.374872903516051510689234173139793159340e-0006L,
138*25c28e83SPiotr Jasiukajtis -2.039156676413643091040459825776029327487e-0006L,
139*25c28e83SPiotr Jasiukajtis 9.555777181318621470466563543806211523634e-0007L,
140*25c28e83SPiotr Jasiukajtis -4.468344919709630637558538313482398989638e-0007L,
141*25c28e83SPiotr Jasiukajtis 2.216738086090045781773004477831059444178e-0007L,
142*25c28e83SPiotr Jasiukajtis -7.472783403418388455860445842543843485916e-0008L,
143*25c28e83SPiotr Jasiukajtis 8.777317930927149922056782132706238921648e-0008L,
144*25c28e83SPiotr Jasiukajtis };
145*25c28e83SPiotr Jasiukajtis
146*25c28e83SPiotr Jasiukajtis static const long double an2[] = { /* 20 terms */
147*25c28e83SPiotr Jasiukajtis -.0772156649015328606062692723698127607018L,
148*25c28e83SPiotr Jasiukajtis 3.224670334241132182635552349060279118047e-0001L,
149*25c28e83SPiotr Jasiukajtis -6.735230105319809367555642883133994818325e-0002L,
150*25c28e83SPiotr Jasiukajtis 2.058080842778459676880822202762143671813e-0002L,
151*25c28e83SPiotr Jasiukajtis -7.385551028672828216011343150077846918930e-0003L,
152*25c28e83SPiotr Jasiukajtis 2.890510330762060607399561536905727853178e-0003L,
153*25c28e83SPiotr Jasiukajtis -1.192753911419623262328187532759756368041e-0003L,
154*25c28e83SPiotr Jasiukajtis 5.096695278636456678258091134532258618614e-0004L,
155*25c28e83SPiotr Jasiukajtis -2.231547306817535743052975194022893369135e-0004L,
156*25c28e83SPiotr Jasiukajtis 9.945771461633313282744264853986643877087e-0005L,
157*25c28e83SPiotr Jasiukajtis -4.492503279458972037926876061257489481619e-0005L,
158*25c28e83SPiotr Jasiukajtis 2.051311416812082875492678651369394595613e-0005L,
159*25c28e83SPiotr Jasiukajtis -9.415778282365955203915850761537462941165e-0006L,
160*25c28e83SPiotr Jasiukajtis 4.452428829045147098722932981088650055919e-0006L,
161*25c28e83SPiotr Jasiukajtis -1.835024727987632579886951760650722695781e-0006L,
162*25c28e83SPiotr Jasiukajtis 1.379783080658545009579060714946381462565e-0006L,
163*25c28e83SPiotr Jasiukajtis 2.282637532109775156769736768748402175238e-0007L,
164*25c28e83SPiotr Jasiukajtis 1.002577375515900191362119718128149880168e-0006L,
165*25c28e83SPiotr Jasiukajtis 5.177028794262638311939991106423220002463e-0007L,
166*25c28e83SPiotr Jasiukajtis 3.127947245174847104122426445937830555755e-0007L,
167*25c28e83SPiotr Jasiukajtis };
168*25c28e83SPiotr Jasiukajtis
169*25c28e83SPiotr Jasiukajtis static const long double an3[] = { /* 20 terms */
170*25c28e83SPiotr Jasiukajtis -.0772156649015328227870646417729220690875L,
171*25c28e83SPiotr Jasiukajtis 3.224670334241156699881788955959915250365e-0001L,
172*25c28e83SPiotr Jasiukajtis -6.735230105312273571375431059744975563170e-0002L,
173*25c28e83SPiotr Jasiukajtis 2.058080842924464587662846071337083809005e-0002L,
174*25c28e83SPiotr Jasiukajtis -7.385551008677271654723604653956131791619e-0003L,
175*25c28e83SPiotr Jasiukajtis 2.890510536479782086197110272583833176602e-0003L,
176*25c28e83SPiotr Jasiukajtis -1.192752262076857692740571567808259138697e-0003L,
177*25c28e83SPiotr Jasiukajtis 5.096800771149805289371135155128380707889e-0004L,
178*25c28e83SPiotr Jasiukajtis -2.231000836682831335505058492409860123647e-0004L,
179*25c28e83SPiotr Jasiukajtis 9.968912171073936803871803966360595275047e-0005L,
180*25c28e83SPiotr Jasiukajtis -4.412020779327746243544387946167256187258e-0005L,
181*25c28e83SPiotr Jasiukajtis 2.281374113541454151067016632998630209049e-0005L,
182*25c28e83SPiotr Jasiukajtis -4.028361291428629491824694655287954266830e-0006L,
183*25c28e83SPiotr Jasiukajtis 1.470694920619518924598956849226530750139e-0005L,
184*25c28e83SPiotr Jasiukajtis 1.381686137617987197975289545582377713772e-0005L,
185*25c28e83SPiotr Jasiukajtis 2.012493539265777728944759982054970441601e-0005L,
186*25c28e83SPiotr Jasiukajtis 1.723917864208965490251560644681933675799e-0005L,
187*25c28e83SPiotr Jasiukajtis 1.202954035243788300138608765425123713395e-0005L,
188*25c28e83SPiotr Jasiukajtis 5.079851887558623092776296577030850938146e-0006L,
189*25c28e83SPiotr Jasiukajtis 1.220657945824153751555138592006604026282e-0006L,
190*25c28e83SPiotr Jasiukajtis };
191*25c28e83SPiotr Jasiukajtis
192*25c28e83SPiotr Jasiukajtis static const long double an4[] = { /* 21 terms */
193*25c28e83SPiotr Jasiukajtis -.0772156649015732285350261816697540392371L,
194*25c28e83SPiotr Jasiukajtis 3.224670334221752060691751340365212226097e-0001L,
195*25c28e83SPiotr Jasiukajtis -6.735230109744009693977755991488196368279e-0002L,
196*25c28e83SPiotr Jasiukajtis 2.058080778913037626909954141611580783216e-0002L,
197*25c28e83SPiotr Jasiukajtis -7.385557567931505621170483708950557506819e-0003L,
198*25c28e83SPiotr Jasiukajtis 2.890459838416254326340844289785254883436e-0003L,
199*25c28e83SPiotr Jasiukajtis -1.193059036207136762877351596966718455737e-0003L,
200*25c28e83SPiotr Jasiukajtis 5.081914708100372836613371356529568937869e-0004L,
201*25c28e83SPiotr Jasiukajtis -2.289855016133600313131553005982542045338e-0004L,
202*25c28e83SPiotr Jasiukajtis 8.053454537980585879620331053833498511491e-0005L,
203*25c28e83SPiotr Jasiukajtis -9.574620532104845821243493405855672438998e-0005L,
204*25c28e83SPiotr Jasiukajtis -9.269085628207107155601445001196317715686e-0005L,
205*25c28e83SPiotr Jasiukajtis -2.183276779859490461716196344776208220180e-0004L,
206*25c28e83SPiotr Jasiukajtis -3.134834305597571096452454999737269668868e-0004L,
207*25c28e83SPiotr Jasiukajtis -3.973878894951937437018305986901392888619e-0004L,
208*25c28e83SPiotr Jasiukajtis -3.953352414899222799161275564386488057119e-0004L,
209*25c28e83SPiotr Jasiukajtis -3.136740932204038779362660900621212816511e-0004L,
210*25c28e83SPiotr Jasiukajtis -1.884502253819634073946130825196078627664e-0004L,
211*25c28e83SPiotr Jasiukajtis -8.192655799958926853585332542123631379301e-0005L,
212*25c28e83SPiotr Jasiukajtis -2.292183750010571062891605074281744854436e-0005L,
213*25c28e83SPiotr Jasiukajtis -3.223980628729716864927724265781406614294e-0006L,
214*25c28e83SPiotr Jasiukajtis };
215*25c28e83SPiotr Jasiukajtis
216*25c28e83SPiotr Jasiukajtis static const long double ap1[] = { /* 19 terms */
217*25c28e83SPiotr Jasiukajtis -0.0772156649015328606065120900824024296961L,
218*25c28e83SPiotr Jasiukajtis 3.224670334241132182362075833230047956465e-0001L,
219*25c28e83SPiotr Jasiukajtis -6.735230105319809513324605382963943777301e-0002L,
220*25c28e83SPiotr Jasiukajtis 2.058080842778454787900092126606252375465e-0002L,
221*25c28e83SPiotr Jasiukajtis -7.385551028673985266272518231365020063941e-0003L,
222*25c28e83SPiotr Jasiukajtis 2.890510330741523285681704570797770736423e-0003L,
223*25c28e83SPiotr Jasiukajtis -1.192753911703260971285304221165990244515e-0003L,
224*25c28e83SPiotr Jasiukajtis 5.096695247430420878696018188830886972245e-0004L,
225*25c28e83SPiotr Jasiukajtis -2.231547584535654004647639737841526025095e-0004L,
226*25c28e83SPiotr Jasiukajtis 9.945751278137201960636098805852315982919e-0005L,
227*25c28e83SPiotr Jasiukajtis -4.492623672777606053587919463929044226280e-0005L,
228*25c28e83SPiotr Jasiukajtis 2.050721258703289487603702670753053765201e-0005L,
229*25c28e83SPiotr Jasiukajtis -9.439485626565616989352750672499008021041e-0006L,
230*25c28e83SPiotr Jasiukajtis 4.374838162403994645138200419356844574219e-0006L,
231*25c28e83SPiotr Jasiukajtis -2.038979492862555348577006944451002161496e-0006L,
232*25c28e83SPiotr Jasiukajtis 9.536763152382263548086981191378885102802e-0007L,
233*25c28e83SPiotr Jasiukajtis -4.426111214332434049863595231916564014913e-0007L,
234*25c28e83SPiotr Jasiukajtis 1.911148847512947464234633846270287546882e-0007L,
235*25c28e83SPiotr Jasiukajtis -5.788673944861923038157839080272303519671e-0008L,
236*25c28e83SPiotr Jasiukajtis };
237*25c28e83SPiotr Jasiukajtis
238*25c28e83SPiotr Jasiukajtis static const long double ap2[] = { /* 19 terms */
239*25c28e83SPiotr Jasiukajtis -0.077215664901532860606428624449354836087L,
240*25c28e83SPiotr Jasiukajtis 3.224670334241132182271948744265855440139e-0001L,
241*25c28e83SPiotr Jasiukajtis -6.735230105319809467356126599005051676203e-0002L,
242*25c28e83SPiotr Jasiukajtis 2.058080842778453315716389815213496002588e-0002L,
243*25c28e83SPiotr Jasiukajtis -7.385551028673653323064118422580096222959e-0003L,
244*25c28e83SPiotr Jasiukajtis 2.890510330735923572088003424849289006039e-0003L,
245*25c28e83SPiotr Jasiukajtis -1.192753911629952368606185543945790688144e-0003L,
246*25c28e83SPiotr Jasiukajtis 5.096695239806718875364547587043220998766e-0004L,
247*25c28e83SPiotr Jasiukajtis -2.231547520600616108991867127392089144886e-0004L,
248*25c28e83SPiotr Jasiukajtis 9.945746913898151120612322833059416008973e-0005L,
249*25c28e83SPiotr Jasiukajtis -4.492599307461977003570224943054585729684e-0005L,
250*25c28e83SPiotr Jasiukajtis 2.050609891889165453592046505651759999090e-0005L,
251*25c28e83SPiotr Jasiukajtis -9.435329866734193796540515247917165988579e-0006L,
252*25c28e83SPiotr Jasiukajtis 4.362267138522223236241016136585565144581e-0006L,
253*25c28e83SPiotr Jasiukajtis -2.008556356653246579300491601497510230557e-0006L,
254*25c28e83SPiotr Jasiukajtis 8.961498103387207161105347118042844354395e-0007L,
255*25c28e83SPiotr Jasiukajtis -3.614187228330216282235692806488341157741e-0007L,
256*25c28e83SPiotr Jasiukajtis 1.136978988247816860500420915014777753153e-0007L,
257*25c28e83SPiotr Jasiukajtis -2.000532786387196664019286514899782691776e-0008L,
258*25c28e83SPiotr Jasiukajtis };
259*25c28e83SPiotr Jasiukajtis
260*25c28e83SPiotr Jasiukajtis static const long double ap3[] = { /* 19 terms */
261*25c28e83SPiotr Jasiukajtis -0.077215664901532859888521470795348856446L,
262*25c28e83SPiotr Jasiukajtis 3.224670334241131733364048614484228443077e-0001L,
263*25c28e83SPiotr Jasiukajtis -6.735230105319676541660495145259038151576e-0002L,
264*25c28e83SPiotr Jasiukajtis 2.058080842775975461837768839015444273830e-0002L,
265*25c28e83SPiotr Jasiukajtis -7.385551028347615729728618066663566606906e-0003L,
266*25c28e83SPiotr Jasiukajtis 2.890510327517954083379032008643080256676e-0003L,
267*25c28e83SPiotr Jasiukajtis -1.192753886919470728001821137439430882603e-0003L,
268*25c28e83SPiotr Jasiukajtis 5.096693728898932234814903769146577482912e-0004L,
269*25c28e83SPiotr Jasiukajtis -2.231540055048827662528594010961874258037e-0004L,
270*25c28e83SPiotr Jasiukajtis 9.945446210018649311491619999438833843723e-0005L,
271*25c28e83SPiotr Jasiukajtis -4.491608206598064519190236245753867697750e-0005L,
272*25c28e83SPiotr Jasiukajtis 2.047939071322271016498065052853746466669e-0005L,
273*25c28e83SPiotr Jasiukajtis -9.376824046522786006677541036631536790762e-0006L,
274*25c28e83SPiotr Jasiukajtis 4.259329829498149111582277209189150127347e-0006L,
275*25c28e83SPiotr Jasiukajtis -1.866064770421594266702176289764212873428e-0006L,
276*25c28e83SPiotr Jasiukajtis 7.462066721137579592928128104534957135669e-0007L,
277*25c28e83SPiotr Jasiukajtis -2.483546217529077735074007138457678727371e-0007L,
278*25c28e83SPiotr Jasiukajtis 5.915166576378161473299324673649144297574e-0008L,
279*25c28e83SPiotr Jasiukajtis -7.334139641706988966966252333759604701905e-0009L,
280*25c28e83SPiotr Jasiukajtis };
281*25c28e83SPiotr Jasiukajtis
282*25c28e83SPiotr Jasiukajtis static const long double ap4[] = { /* 19 terms */
283*25c28e83SPiotr Jasiukajtis -0.0772156649015326785569313252637238673675L,
284*25c28e83SPiotr Jasiukajtis 3.224670334241051435008842685722468344822e-0001L,
285*25c28e83SPiotr Jasiukajtis -6.735230105302832007479431772160948499254e-0002L,
286*25c28e83SPiotr Jasiukajtis 2.058080842553481183648529360967441889912e-0002L,
287*25c28e83SPiotr Jasiukajtis -7.385551007602909242024706804659879199244e-0003L,
288*25c28e83SPiotr Jasiukajtis 2.890510182473907253939821312248303471206e-0003L,
289*25c28e83SPiotr Jasiukajtis -1.192753098427856770847894497586825614450e-0003L,
290*25c28e83SPiotr Jasiukajtis 5.096659636418811568063339214203693550804e-0004L,
291*25c28e83SPiotr Jasiukajtis -2.231421144004355691166194259675004483639e-0004L,
292*25c28e83SPiotr Jasiukajtis 9.942073842343832132754332881883387625136e-0005L,
293*25c28e83SPiotr Jasiukajtis -4.483809261973204531263252655050701205397e-0005L,
294*25c28e83SPiotr Jasiukajtis 2.033260142610284888319116654931994447173e-0005L,
295*25c28e83SPiotr Jasiukajtis -9.153539544026646699870528191410440585796e-0006L,
296*25c28e83SPiotr Jasiukajtis 3.988460469925482725894144688699584997971e-0006L,
297*25c28e83SPiotr Jasiukajtis -1.609692980087029172567957221850825977621e-0006L,
298*25c28e83SPiotr Jasiukajtis 5.634916377249975825399706694496688803488e-0007L,
299*25c28e83SPiotr Jasiukajtis -1.560065465929518563549083208482591437696e-0007L,
300*25c28e83SPiotr Jasiukajtis 2.961350193868935325526962209019387821584e-0008L,
301*25c28e83SPiotr Jasiukajtis -2.834602215195368130104649234505033159842e-0009L,
302*25c28e83SPiotr Jasiukajtis };
303*25c28e83SPiotr Jasiukajtis
304*25c28e83SPiotr Jasiukajtis static long double
primary(long double s)305*25c28e83SPiotr Jasiukajtis primary(long double s) { /* assume |s|<=0.5 */
306*25c28e83SPiotr Jasiukajtis int i;
307*25c28e83SPiotr Jasiukajtis
308*25c28e83SPiotr Jasiukajtis i = (int) (8.0L * (s + 0.5L));
309*25c28e83SPiotr Jasiukajtis switch (i) {
310*25c28e83SPiotr Jasiukajtis case 0: return ch*s+s*poly(s, an4, 21);
311*25c28e83SPiotr Jasiukajtis case 1: return ch*s+s*poly(s, an3, 20);
312*25c28e83SPiotr Jasiukajtis case 2: return ch*s+s*poly(s, an2, 20);
313*25c28e83SPiotr Jasiukajtis case 3: return ch*s+s*poly(s, an1, 20);
314*25c28e83SPiotr Jasiukajtis case 4: return ch*s+s*poly(s, ap1, 19);
315*25c28e83SPiotr Jasiukajtis case 5: return ch*s+s*poly(s, ap2, 19);
316*25c28e83SPiotr Jasiukajtis case 6: return ch*s+s*poly(s, ap3, 19);
317*25c28e83SPiotr Jasiukajtis case 7: return ch*s+s*poly(s, ap4, 19);
318*25c28e83SPiotr Jasiukajtis }
319*25c28e83SPiotr Jasiukajtis /* NOTREACHED */
320*25c28e83SPiotr Jasiukajtis return (0.0L);
321*25c28e83SPiotr Jasiukajtis }
322*25c28e83SPiotr Jasiukajtis
323*25c28e83SPiotr Jasiukajtis static long double
poly(long double s,const long double * p,int n)324*25c28e83SPiotr Jasiukajtis poly(long double s, const long double *p, int n) {
325*25c28e83SPiotr Jasiukajtis long double y;
326*25c28e83SPiotr Jasiukajtis int i;
327*25c28e83SPiotr Jasiukajtis y = p[n-1];
328*25c28e83SPiotr Jasiukajtis for (i = n-2; i >= 0; i--) y = p[i]+s*y;
329*25c28e83SPiotr Jasiukajtis return (y);
330*25c28e83SPiotr Jasiukajtis }
331*25c28e83SPiotr Jasiukajtis
332*25c28e83SPiotr Jasiukajtis static const long double pt[] = {
333*25c28e83SPiotr Jasiukajtis 9.189385332046727417803297364056176804663e-0001L,
334*25c28e83SPiotr Jasiukajtis 8.333333333333333333333333333331286969123e-0002L,
335*25c28e83SPiotr Jasiukajtis -2.777777777777777777777777553194796036402e-0003L,
336*25c28e83SPiotr Jasiukajtis 7.936507936507936507927283071433584248176e-0004L,
337*25c28e83SPiotr Jasiukajtis -5.952380952380952362351042163192634108297e-0004L,
338*25c28e83SPiotr Jasiukajtis 8.417508417508395661774286645578379460131e-0004L,
339*25c28e83SPiotr Jasiukajtis -1.917526917525263651186066417934685675649e-0003L,
340*25c28e83SPiotr Jasiukajtis 6.410256409395203164659292973142293199083e-0003L,
341*25c28e83SPiotr Jasiukajtis -2.955065327248303301763594514012418438188e-0002L,
342*25c28e83SPiotr Jasiukajtis 1.796442830099067542945998615411893822886e-0001L,
343*25c28e83SPiotr Jasiukajtis -1.392413465829723742489974310411118662919e+0000L,
344*25c28e83SPiotr Jasiukajtis 1.339984238037267658352656597960492029261e+0001L,
345*25c28e83SPiotr Jasiukajtis -1.564707657605373662425785904278645727813e+0002L,
346*25c28e83SPiotr Jasiukajtis 2.156323807499211356127813962223067079300e+0003L,
347*25c28e83SPiotr Jasiukajtis -3.330486427626223184647299834137041307569e+0004L,
348*25c28e83SPiotr Jasiukajtis 5.235535072011889213611369254140123518699e+0005L,
349*25c28e83SPiotr Jasiukajtis -7.258160984602220710491988573430212593080e+0006L,
350*25c28e83SPiotr Jasiukajtis 7.316526934569686459641438882340322673357e+0007L,
351*25c28e83SPiotr Jasiukajtis -3.806450279064900548836571789284896711473e+0008L,
352*25c28e83SPiotr Jasiukajtis };
353*25c28e83SPiotr Jasiukajtis
354*25c28e83SPiotr Jasiukajtis static long double
polytail(long double s)355*25c28e83SPiotr Jasiukajtis polytail(long double s) {
356*25c28e83SPiotr Jasiukajtis long double t, z;
357*25c28e83SPiotr Jasiukajtis int i;
358*25c28e83SPiotr Jasiukajtis z = s*s;
359*25c28e83SPiotr Jasiukajtis t = pt[18];
360*25c28e83SPiotr Jasiukajtis for (i = 17; i >= 1; i--) t = pt[i]+z*t;
361*25c28e83SPiotr Jasiukajtis return (pt[0]+s*t);
362*25c28e83SPiotr Jasiukajtis }
363*25c28e83SPiotr Jasiukajtis
364*25c28e83SPiotr Jasiukajtis static long double
neg(long double z,int * signgamlp)365*25c28e83SPiotr Jasiukajtis neg(long double z, int *signgamlp) {
366*25c28e83SPiotr Jasiukajtis long double t, p;
367*25c28e83SPiotr Jasiukajtis
368*25c28e83SPiotr Jasiukajtis /*
369*25c28e83SPiotr Jasiukajtis * written by K.C. Ng, Feb 2, 1989.
370*25c28e83SPiotr Jasiukajtis *
371*25c28e83SPiotr Jasiukajtis * Since
372*25c28e83SPiotr Jasiukajtis * -z*G(-z)*G(z) = pi/sin(pi*z),
373*25c28e83SPiotr Jasiukajtis * we have
374*25c28e83SPiotr Jasiukajtis * G(-z) = -pi/(sin(pi*z)*G(z)*z)
375*25c28e83SPiotr Jasiukajtis * = pi/(sin(pi*(-z))*G(z)*z)
376*25c28e83SPiotr Jasiukajtis * Algorithm
377*25c28e83SPiotr Jasiukajtis * z = |z|
378*25c28e83SPiotr Jasiukajtis * t = sinpi(z); ...note that when z>2**112, z is an int
379*25c28e83SPiotr Jasiukajtis * and hence t=0.
380*25c28e83SPiotr Jasiukajtis *
381*25c28e83SPiotr Jasiukajtis * if (t == 0.0) return 1.0/0.0;
382*25c28e83SPiotr Jasiukajtis * if (t< 0.0) *signgamlp = -1; else t= -t;
383*25c28e83SPiotr Jasiukajtis * if (z<1.0e-40) ...tiny z
384*25c28e83SPiotr Jasiukajtis * return -log(z);
385*25c28e83SPiotr Jasiukajtis * else
386*25c28e83SPiotr Jasiukajtis * return log(pi/(t*z))-lgamma(z);
387*25c28e83SPiotr Jasiukajtis *
388*25c28e83SPiotr Jasiukajtis */
389*25c28e83SPiotr Jasiukajtis
390*25c28e83SPiotr Jasiukajtis t = sinpil(z); /* t := sin(pi*z) */
391*25c28e83SPiotr Jasiukajtis if (t == c0) /* return 1.0/0.0 = +INF */
392*25c28e83SPiotr Jasiukajtis return (c1/c0);
393*25c28e83SPiotr Jasiukajtis
394*25c28e83SPiotr Jasiukajtis z = -z;
395*25c28e83SPiotr Jasiukajtis if (z <= tiny)
396*25c28e83SPiotr Jasiukajtis p = -logl(z);
397*25c28e83SPiotr Jasiukajtis else
398*25c28e83SPiotr Jasiukajtis p = logl(pi/(fabsl(t)*z)) - __k_lgammal(z, signgamlp);
399*25c28e83SPiotr Jasiukajtis if (t < c0) *signgamlp = -1;
400*25c28e83SPiotr Jasiukajtis return (p);
401*25c28e83SPiotr Jasiukajtis }
402