master
1/**
2 * This file has no copyright assigned and is placed in the Public Domain.
3 * This file is part of the mingw-w64 runtime package.
4 * No warranty is given; refer to the file DISCLAIMER.PD within this package.
5 */
6/* log gamma(x+2), -.5 < x < .5 */
7static const float B[] = {
8 6.055172732649237E-004,
9 -1.311620815545743E-003,
10 2.863437556468661E-003,
11 -7.366775108654962E-003,
12 2.058355474821512E-002,
13 -6.735323259371034E-002,
14 3.224669577325661E-001,
15 4.227843421859038E-001
16};
17
18/* log gamma(x+1), -.25 < x < .25 */
19static const float C[] = {
20 1.369488127325832E-001,
21 -1.590086327657347E-001,
22 1.692415923504637E-001,
23 -2.067882815621965E-001,
24 2.705806208275915E-001,
25 -4.006931650563372E-001,
26 8.224670749082976E-001,
27 -5.772156501719101E-001
28};
29
30/* log( sqrt( 2*pi ) ) */
31static const float LS2PI = 0.91893853320467274178;
32#define MAXLGM 2.035093e36
33static const float PIINV = 0.318309886183790671538;
34
35#include "cephes_mconf.h"
36
37/* Reentrant version */
38/* Logarithm of gamma function */
39float __lgammaf_r(float x, int* sgngamf);
40
41float __lgammaf_r(float x, int* sgngamf)
42{
43 float p, q, w, z;
44 float nx, tx;
45 int i, direction;
46
47 *sgngamf = 1;
48#ifdef NANS
49 if (isnan(x))
50 return (x);
51#endif
52
53#ifdef INFINITIES
54 if (!isfinite(x))
55 return (INFINITY);
56#endif
57
58 if (x < 0.0)
59 {
60 q = -x;
61 w = __lgammaf_r(q, sgngamf); /* note this modifies sgngam! */
62 p = floorf(q);
63 if (p == q)
64 {
65lgsing:
66 _SET_ERRNO(EDOM);
67 mtherr("lgamf", SING);
68#ifdef INFINITIES
69 return (INFINITYF);
70#else
71 return( *sgngamf * MAXNUMF );
72#endif
73 }
74 i = p;
75 if ((i & 1) == 0)
76 *sgngamf = -1;
77 else
78 *sgngamf = 1;
79 z = q - p;
80 if (z > 0.5)
81 {
82 p += 1.0;
83 z = p - q;
84 }
85 z = q * sinf(PIF * z);
86 if (z == 0.0)
87 goto lgsing;
88 z = -logf(PIINV * z) - w;
89 return (z);
90 }
91
92 if (x < 6.5)
93 {
94 direction = 0;
95 z = 1.0;
96 tx = x;
97 nx = 0.0;
98 if (x >= 1.5)
99 {
100 while (tx > 2.5)
101 {
102 nx -= 1.0;
103 tx = x + nx;
104 z *=tx;
105 }
106 x += nx - 2.0;
107iv1r5:
108 p = x * polevlf( x, B, 7 );
109 goto cont;
110 }
111 if (x >= 1.25)
112 {
113 z *= x;
114 x -= 1.0; /* x + 1 - 2 */
115 direction = 1;
116 goto iv1r5;
117 }
118 if (x >= 0.75)
119 {
120 x -= 1.0;
121 p = x * polevlf( x, C, 7 );
122 q = 0.0;
123 goto contz;
124 }
125 while (tx < 1.5)
126 {
127 if (tx == 0.0)
128 goto lgsing;
129 z *=tx;
130 nx += 1.0;
131 tx = x + nx;
132 }
133 direction = 1;
134 x += nx - 2.0;
135 p = x * polevlf(x, B, 7);
136
137cont:
138 if (z < 0.0)
139 {
140 *sgngamf = -1;
141 z = -z;
142 }
143 else
144 {
145 *sgngamf = 1;
146 }
147 q = logf(z);
148 if (direction)
149 q = -q;
150contz:
151 return( p + q );
152 }
153
154 if (x > MAXLGM)
155 {
156 _SET_ERRNO(ERANGE);
157 mtherr("lgamf", OVERFLOW);
158#ifdef INFINITIES
159 return (*sgngamf * INFINITYF);
160#else
161 return (*sgngamf * MAXNUMF);
162#endif
163
164 }
165
166/* Note, though an asymptotic formula could be used for x >= 3,
167 * there is cancellation error in the following if x < 6.5. */
168 q = LS2PI - x;
169 q += (x - 0.5) * logf(x);
170
171 if (x <= 1.0e4)
172 {
173 z = 1.0/x;
174 p = z * z;
175 q += (( 6.789774945028216E-004 * p
176 - 2.769887652139868E-003 ) * p
177 + 8.333316229807355E-002 ) * z;
178 }
179 return (q);
180}
181
182/* This is the C99 version */
183float lgammaf(float x)
184{
185 return (__lgammaf_r(x, &signgam));
186}
187