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