5 * Common Development and Distribution License (the "License").
6 * You may not use this file except in compliance with the License.
7 *
8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9 * or http://www.opensolaris.org/os/licensing.
10 * See the License for the specific language governing permissions
11 * and limitations under the License.
12 *
13 * When distributing Covered Code, include this CDDL HEADER in each
14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
15 * If applicable, add the following below this CDDL HEADER, with the
16 * fields enclosed by brackets "[]" replaced with your own identifying
17 * information: Portions Copyright [yyyy] [name of copyright owner]
18 *
19 * CDDL HEADER END
20 */
21
22 /*
23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
24 */
25 /*
26 * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
27 * Use is subject to license terms.
28 */
29
30 /*
31 * long double __k_lgammal(long double x, int *signgamlp);
32 * K.C. Ng, August, 1989.
33 *
34 * We choose [1.5,2.5] to be the primary interval. Our algorithms
35 * are mainly derived from
36 *
37 *
38 * zeta(2)-1 2 zeta(3)-1 3
39 * lgamma(2+s) = s*(1-euler) + --------- * s - --------- * s + ...
40 * 2 3
41 *
42 *
43 * 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).
46 * When s is really tiny (like roundoff), lgamma(1+s) ~ s(1-enler)
47 * Hence lgamma(s) ~ -log(s) for tiny s
48 *
49 */
50
51 #include "libm.h"
52 #include "longdouble.h"
53
54 static long double neg(long double, int *);
55 static long double poly(long double, const long double *, int);
56 static long double polytail(long double);
57 static long double primary(long double);
58
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;
70
71 long double
72 __k_lgammal(long double x, int *signgamlp) {
73 long double t, y;
74 int i;
75
76 /* purge off +-inf, NaN and negative arguments */
77 if (!finitel(x))
78 return (x*x);
79 *signgamlp = 1;
80 if (signbitl(x))
81 return (neg(x, signgamlp));
82
83 /* for x < 8.0 */
84 if (x < 8.0L) {
85 y = anintl(x);
86 i = (int) y;
87 switch (i) {
88 case 0:
89 if (x < 1.0e-40L)
90 return (-logl(x));
91 else
92 return (primary(x)-log1pl(x))-logl(x);
93 case 1:
94 return (primary(x-y)-logl(x));
95 case 2:
96 return (primary(x-y));
97 case 3:
98 return (primary(x-y)+logl(x-c1));
99 case 4:
100 return (primary(x-y)+logl((x-c1)*(x-c2)));
101 case 5:
102 return (primary(x-y)+logl((x-c1)*(x-c2)*(x-c3)));
103 case 6:
104 return (primary(x-y)+logl((x-c1)*(x-c2)*(x-c3)*(x-c4)));
105 case 7:
106 return (primary(x-y)+logl((x-c1)*(x-c2)*(x-c3)*(x-c4)*(x-c5)));
107 case 8:
108 return primary(x-y)+
109 logl((x-c1)*(x-c2)*(x-c3)*(x-c4)*(x-c5)*(x-c6));
110 }
111 }
112
113 /* 8.0 <= x < 1.0e40 */
114 if (x < 1.0e40L) {
115 t = logl(x);
116 return (x*(t-c1)-(ch*t-polytail(c1/x)));
117 }
118
119 /* 1.0e40 <= x <= inf */
120 return (x*(logl(x)-c1));
121 }
122
123 static const long double an1[] = { /* 20 terms */
124 -0.0772156649015328606065120900824024309741L,
125 3.224670334241132182362075833230130289059e-0001L,
126 -6.735230105319809513324605383668929964120e-0002L,
127 2.058080842778454787900092432928910226297e-0002L,
128 -7.385551028673985266273054086081102125704e-0003L,
129 2.890510330741523285758867304409628648727e-0003L,
130 -1.192753911703260976581414338096267498555e-0003L,
131 5.096695247430424562831956662855697824035e-0004L,
132 -2.231547584535777978926798502084300123638e-0004L,
133 9.945751278186384670278268034322157947635e-0005L,
134 -4.492623673665547726647838474125147631082e-0005L,
135 2.050721280617796810096993154281561168706e-0005L,
136 -9.439487785617396552092393234044767313568e-0006L,
137 4.374872903516051510689234173139793159340e-0006L,
138 -2.039156676413643091040459825776029327487e-0006L,
139 9.555777181318621470466563543806211523634e-0007L,
140 -4.468344919709630637558538313482398989638e-0007L,
285 -6.735230105302832007479431772160948499254e-0002L,
286 2.058080842553481183648529360967441889912e-0002L,
287 -7.385551007602909242024706804659879199244e-0003L,
288 2.890510182473907253939821312248303471206e-0003L,
289 -1.192753098427856770847894497586825614450e-0003L,
290 5.096659636418811568063339214203693550804e-0004L,
291 -2.231421144004355691166194259675004483639e-0004L,
292 9.942073842343832132754332881883387625136e-0005L,
293 -4.483809261973204531263252655050701205397e-0005L,
294 2.033260142610284888319116654931994447173e-0005L,
295 -9.153539544026646699870528191410440585796e-0006L,
296 3.988460469925482725894144688699584997971e-0006L,
297 -1.609692980087029172567957221850825977621e-0006L,
298 5.634916377249975825399706694496688803488e-0007L,
299 -1.560065465929518563549083208482591437696e-0007L,
300 2.961350193868935325526962209019387821584e-0008L,
301 -2.834602215195368130104649234505033159842e-0009L,
302 };
303
304 static long double
305 primary(long double s) { /* assume |s|<=0.5 */
306 int i;
307
308 i = (int) (8.0L * (s + 0.5L));
309 switch (i) {
310 case 0: return ch*s+s*poly(s, an4, 21);
311 case 1: return ch*s+s*poly(s, an3, 20);
312 case 2: return ch*s+s*poly(s, an2, 20);
313 case 3: return ch*s+s*poly(s, an1, 20);
314 case 4: return ch*s+s*poly(s, ap1, 19);
315 case 5: return ch*s+s*poly(s, ap2, 19);
316 case 6: return ch*s+s*poly(s, ap3, 19);
317 case 7: return ch*s+s*poly(s, ap4, 19);
318 }
319 /* NOTREACHED */
320 return (0.0L);
321 }
322
323 static long double
324 poly(long double s, const long double *p, int n) {
325 long double y;
326 int i;
327 y = p[n-1];
328 for (i = n-2; i >= 0; i--) y = p[i]+s*y;
329 return (y);
330 }
331
332 static const long double pt[] = {
333 9.189385332046727417803297364056176804663e-0001L,
334 8.333333333333333333333333333331286969123e-0002L,
335 -2.777777777777777777777777553194796036402e-0003L,
336 7.936507936507936507927283071433584248176e-0004L,
337 -5.952380952380952362351042163192634108297e-0004L,
338 8.417508417508395661774286645578379460131e-0004L,
339 -1.917526917525263651186066417934685675649e-0003L,
340 6.410256409395203164659292973142293199083e-0003L,
341 -2.955065327248303301763594514012418438188e-0002L,
342 1.796442830099067542945998615411893822886e-0001L,
343 -1.392413465829723742489974310411118662919e+0000L,
344 1.339984238037267658352656597960492029261e+0001L,
345 -1.564707657605373662425785904278645727813e+0002L,
346 2.156323807499211356127813962223067079300e+0003L,
347 -3.330486427626223184647299834137041307569e+0004L,
348 5.235535072011889213611369254140123518699e+0005L,
349 -7.258160984602220710491988573430212593080e+0006L,
350 7.316526934569686459641438882340322673357e+0007L,
351 -3.806450279064900548836571789284896711473e+0008L,
352 };
353
354 static long double
355 polytail(long double s) {
356 long double t, z;
357 int i;
358 z = s*s;
359 t = pt[18];
360 for (i = 17; i >= 1; i--) t = pt[i]+z*t;
361 return (pt[0]+s*t);
362 }
363
364 static long double
365 neg(long double z, int *signgamlp) {
366 long double t, p;
367
368 /*
369 * written by K.C. Ng, Feb 2, 1989.
370 *
371 * Since
372 * -z*G(-z)*G(z) = pi/sin(pi*z),
373 * we have
374 * G(-z) = -pi/(sin(pi*z)*G(z)*z)
375 * = pi/(sin(pi*(-z))*G(z)*z)
376 * Algorithm
377 * z = |z|
378 * t = sinpi(z); ...note that when z>2**112, z is an int
379 * and hence t=0.
380 *
381 * if (t == 0.0) return 1.0/0.0;
382 * if (t< 0.0) *signgamlp = -1; else t= -t;
383 * if (z<1.0e-40) ...tiny z
384 * return -log(z);
385 * else
386 * return log(pi/(t*z))-lgamma(z);
387 *
388 */
389
390 t = sinpil(z); /* t := sin(pi*z) */
391 if (t == c0) /* return 1.0/0.0 = +INF */
392 return (c1/c0);
393
394 z = -z;
395 if (z <= tiny)
396 p = -logl(z);
397 else
398 p = logl(pi/(fabsl(t)*z)) - __k_lgammal(z, signgamlp);
399 if (t < c0) *signgamlp = -1;
400 return (p);
401 }
|
5 * Common Development and Distribution License (the "License").
6 * You may not use this file except in compliance with the License.
7 *
8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9 * or http://www.opensolaris.org/os/licensing.
10 * See the License for the specific language governing permissions
11 * and limitations under the License.
12 *
13 * When distributing Covered Code, include this CDDL HEADER in each
14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
15 * If applicable, add the following below this CDDL HEADER, with the
16 * fields enclosed by brackets "[]" replaced with your own identifying
17 * information: Portions Copyright [yyyy] [name of copyright owner]
18 *
19 * CDDL HEADER END
20 */
21
22 /*
23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
24 */
25
26 /*
27 * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
28 * Use is subject to license terms.
29 */
30
31 /*
32 * long double __k_lgammal(long double x, int *signgamlp);
33 * K.C. Ng, August, 1989.
34 *
35 * We choose [1.5,2.5] to be the primary interval. Our algorithms
36 * are mainly derived from
37 *
38 *
39 * zeta(2)-1 2 zeta(3)-1 3
40 * lgamma(2+s) = s*(1-euler) + --------- * s - --------- * s + ...
41 * 2 3
42 *
43 *
44 * Note 1. Since gamma(1+s)=s*gamma(s), hence
45 * lgamma(1+s) = log(s) + lgamma(s), or
46 * lgamma(s) = lgamma(1+s) - log(s).
47 * When s is really tiny (like roundoff), lgamma(1+s) ~ s(1-enler)
48 * Hence lgamma(s) ~ -log(s) for tiny s
49 *
50 */
51
52 #include "libm.h"
53 #include "longdouble.h"
54
55 static long double neg(long double, int *);
56 static long double poly(long double, const long double *, int);
57 static long double polytail(long double);
58 static long double primary(long double);
59
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
71 long double
72 __k_lgammal(long double x, int *signgamlp)
73 {
74 long double t, y;
75 int i;
76
77 /* purge off +-inf, NaN and negative arguments */
78 if (!finitel(x))
79 return (x * x);
80
81 *signgamlp = 1;
82
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 }
120 }
121
122 /* 8.0 <= x < 1.0e40 */
123 if (x < 1.0e40L) {
124 t = logl(x);
125 return (x * (t - c1) - (ch * t - polytail(c1 / x)));
126 }
127
128 /* 1.0e40 <= x <= inf */
129 return (x * (logl(x) - c1));
130 }
131
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,
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,
311 };
312
313 static long double
314 primary(long double s) /* assume |s|<=0.5 */
315 {
316 int i;
317
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));
337 }
338
339 /* NOTREACHED */
340 return (0.0L);
341 }
342
343 static long double
344 poly(long double s, const long double *p, int n)
345 {
346 long double y;
347 int i;
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);
355 }
356
357 static const long double pt[] = {
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,
377 };
378
379 static long double
380 polytail(long double s)
381 {
382 long double t, z;
383 int i;
384
385 z = s * s;
386 t = pt[18];
387
388 for (i = 17; i >= 1; i--)
389 t = pt[i] + z * t;
390
391 return (pt[0] + s * t);
392 }
393
394 static long double
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 */
422
423 t = sinpil(z); /* t := sin(pi*z) */
424
425 if (t == c0) /* return 1.0/0.0 = +INF */
426 return (c1 / c0);
427
428 z = -z;
429
430 if (z <= tiny)
431 p = -logl(z);
432 else
433 p = logl(pi / (fabsl(t) * z)) - __k_lgammal(z, signgamlp);
434
435 if (t < c0)
436 *signgamlp = -1;
437
438 return (p);
439 }
|