← Back to branch summary

~mingw-w64/mingw-w64/experimental

« back to all changes in this revision

Viewing changes to ros-privexp/mingw-w64-crt/math/tgamma.c

  • Committer: NightStrike
  • Date: 2010-08-11 22:20:57 UTC
  • Revision ID: svn-v4:4407c894-4637-0410-b4f5-ada5f102cad1:experimental:3266
Branch for adding option for supporting ros

Show diffs side-by-side

added added

removed removed

Lines of Context:
ros-privexp/mingw-w64-crt/math/tgamma.c
 
1
/**
 
2
 * This file has no copyright assigned and is placed in the Public Domain.
 
3
 * This file is part of the w64 mingw-runtime package.
 
4
 * No warranty is given; refer to the file DISCLAIMER.PD within this package.
 
5
 */
 
6
#include "cephes_mconf.h"
 
7
 
 
8
#ifdef UNK
 
9
static const uD P[7] = {
 
10
  { { 1.60119522476751861407E-4 } },
 
11
  { { 1.19135147006586384913E-3 } },
 
12
  { { 1.04213797561761569935E-2 } },
 
13
  { { 4.76367800457137231464E-2 } },
 
14
  { { 2.07448227648435975150E-1 } },
 
15
  { { 4.94214826801497100753E-1 } },
 
16
  { { 9.99999999999999996796E-1 } }
 
17
};
 
18
static const uD Q[8] = {
 
19
  { { -2.31581873324120129819E-5 } },
 
20
  { { 5.39605580493303397842E-4 } },
 
21
  { { -4.45641913851797240494E-3 } },
 
22
  { { 1.18139785222060435552E-2 } },
 
23
  { { 3.58236398605498653373E-2 } },
 
24
  { { -2.34591795718243348568E-1 } },
 
25
  { { 7.14304917030273074085E-2 } },
 
26
  { { 1.00000000000000000320E0 } }
 
27
};
 
28
#define MAXGAM 171.624376956302725
 
29
static const double LOGPI = 1.14472988584940017414;
 
30
#endif
 
31
 
 
32
#ifdef IBMPC
 
33
static const uD P[7] = {
 
34
  { { 0x2153,0x3998,0xfcb8,0x3f24 } },
 
35
  { { 0xbfab,0xe686,0x84e3,0x3f53 } },
 
36
  { { 0x14b0,0xe9db,0x57cd,0x3f85 } },
 
37
  { { 0x23d3,0x18c4,0x63d9,0x3fa8 } },
 
38
  { { 0x7d31,0xdcae,0x8da9,0x3fca } },
 
39
  { { 0xe312,0x3993,0xa137,0x3fdf } },
 
40
  { { 0x0000,0x0000,0x0000,0x3ff0 } }
 
41
};
 
42
static const uD Q[8] = {
 
43
  { { 0xd3af,0x8400,0x487a,0xbef8 } },
 
44
  { { 0x2573,0x2915,0xae8a,0x3f41 } },
 
45
  { { 0xb44a,0xe750,0x40e4,0xbf72 } },
 
46
  { { 0xb117,0x5b1b,0x31ed,0x3f88 } },
 
47
  { { 0xde67,0xe33f,0x5779,0x3fa2 } },
 
48
  { { 0x87c2,0x9d42,0x071a,0xbfce } },
 
49
  { { 0x3c51,0xc9cd,0x4944,0x3fb2 } },
 
50
  { { 0x0000,0x0000,0x0000,0x3ff0 } }
 
51
};
 
52
#define MAXGAM 171.624376956302725
 
53
#endif 
 
54
 
 
55
#ifdef MIEEE
 
56
static const uD P[7] = {
 
57
  { { 0x3f24,0xfcb8,0x3998,0x2153 } },
 
58
  { { 0x3f53,0x84e3,0xe686,0xbfab } },
 
59
  { { 0x3f85,0x57cd,0xe9db,0x14b0 } },
 
60
  { { 0x3fa8,0x63d9,0x18c4,0x23d3 } },
 
61
  { { 0x3fca,0x8da9,0xdcae,0x7d31 } },
 
62
  { { 0x3fdf,0xa137,0x3993,0xe312 } },
 
63
  { { 0x3ff0,0x0000,0x0000,0x0000 } }
 
64
};
 
65
static const unsigned short Q[8] = {
 
66
  { { 0xbef8,0x487a,0x8400,0xd3af } },
 
67
  { { 0x3f41,0xae8a,0x2915,0x2573 } },
 
68
  { { 0xbf72,0x40e4,0xe750,0xb44a } },
 
69
  { { 0x3f88,0x31ed,0x5b1b,0xb117 } },
 
70
  { { 0x3fa2,0x5779,0xe33f,0xde67 } },
 
71
  { { 0xbfce,0x071a,0x9d42,0x87c2 } },
 
72
  { { 0x3fb2,0x4944,0xc9cd,0x3c51 } },
 
73
  { { 0x3ff0,0x0000,0x0000,0x0000 } }
 
74
};
 
75
#define MAXGAM 171.624376956302725
 
76
#endif 
 
77
 
 
78
/* Stirling's formula for the gamma function */
 
79
#if UNK
 
80
static const uD STIR[5] = {
 
81
  { { 7.87311395793093628397E-4 } },
 
82
  { { -2.29549961613378126380E-4 } },
 
83
  { { -2.68132617805781232825E-3 } },
 
84
  { { 3.47222221605458667310E-3 } },
 
85
  { { 8.33333333333482257126E-2 } }
 
86
};
 
87
#define MAXSTIR 143.01608
 
88
static const double SQTPI = 2.50662827463100050242E0;
 
89
#endif
 
90
#if IBMPC
 
91
static const uD STIR[5] = {
 
92
  { { 0x7293,0x592d,0xcc72,0x3f49 } },
 
93
  { { 0x1d7c,0x27e6,0x166b,0xbf2e } },
 
94
  { { 0x4fd7,0x07d4,0xf726,0xbf65 } },
 
95
  { { 0xc5fd,0x1b98,0x71c7,0x3f6c } },
 
96
  { { 0x5986,0x5555,0x5555,0x3fb5 } }
 
97
};
 
98
#define MAXSTIR 143.01608
 
99
 
 
100
static const union
 
101
{
 
102
  unsigned short s[4];
 
103
  double d;
 
104
} sqt = {{0x2706,0x1ff6,0x0d93,0x4004}};
 
105
#define SQTPI (sqt.d)
 
106
#endif
 
107
#if MIEEE
 
108
static const uD STIR[5] = {
 
109
  { { 0x3f49,0xcc72,0x592d,0x7293 } },
 
110
  { { 0xbf2e,0x166b,0x27e6,0x1d7c } },
 
111
  { { 0xbf65,0xf726,0x07d4,0x4fd7 } },
 
112
  { { 0x3f6c,0x71c7,0x1b98,0xc5fd } },
 
113
  { { 0x3fb5,0x5555,0x5555,0x5986 } }
 
114
};
 
115
#define MAXSTIR 143.01608
 
116
static const uD SQT = {
 
117
  { { 0x4004,0x0d93,0x1ff6,0x2706 } }
 
118
};
 
119
#define SQTPI SQT.d
 
120
#endif
 
121
 
 
122
static double stirf (double);
 
123
 
 
124
/* Gamma function computed by Stirling's formula.
 
125
 * The polynomial STIR is valid for 33 <= x <= 172.
 
126
 */
 
127
static double stirf(double x)
 
128
{
 
129
        double y, w, v;
 
130
 
 
131
        w = 1.0/x;
 
132
        w = 1.0 + w * polevl(w, STIR, 4);
 
133
        y = exp(x);
 
134
        if (x > MAXSTIR)
 
135
        { /* Avoid overflow in pow() */
 
136
                v = pow(x, 0.5 * x - 0.25);
 
137
                y = v * (v / y);
 
138
        }
 
139
        else
 
140
        {
 
141
                y = pow(x, x - 0.5) / y;
 
142
        }
 
143
        y = SQTPI * y * w;
 
144
        return (y);
 
145
}
 
146
 
 
147
 
 
148
double __tgamma_r(double x, int *sgngam);
 
149
 
 
150
double __tgamma_r(double x, int *sgngam)
 
151
{
 
152
        double p, q, z;
 
153
        int i;
 
154
 
 
155
        *sgngam = 1;
 
156
#ifdef NANS
 
157
        if (isnan(x))
 
158
                return (x);
 
159
#endif
 
160
#ifdef INFINITIES
 
161
#ifdef NANS
 
162
        if (x == INFINITY)
 
163
                return (x);
 
164
        if (x == -INFINITY)
 
165
                return (NAN);
 
166
#else
 
167
        if (!isfinite(x))
 
168
                return (x);
 
169
#endif
 
170
#endif
 
171
        q = fabs(x);
 
172
 
 
173
        if (q > 33.0)
 
174
        {
 
175
                if (x < 0.0)
 
176
                {
 
177
                        p = floor(q);
 
178
                        if (p == q)
 
179
                        {
 
180
gsing:
 
181
                                _SET_ERRNO(EDOM);
 
182
                                mtherr("tgamma", SING);
 
183
#ifdef INFINITIES
 
184
                                return (INFINITY);
 
185
#else
 
186
                                return (MAXNUM);
 
187
#endif
 
188
                        }
 
189
                        i = p;
 
190
                        if ((i & 1) == 0)
 
191
                                *sgngam = -1;
 
192
                        z = q - p;
 
193
                        if (z > 0.5)
 
194
                        {
 
195
                                p += 1.0;
 
196
                                z = q - p;
 
197
                        }
 
198
                        z = q * sin(PI * z);
 
199
                        if (z == 0.0)
 
200
                        {
 
201
                                _SET_ERRNO(ERANGE);
 
202
                                mtherr("tgamma", OVERFLOW);
 
203
#ifdef INFINITIES
 
204
                                return (*sgngam * INFINITY);
 
205
#else
 
206
                                return (*sgngam * MAXNUM);
 
207
#endif
 
208
                        }
 
209
                        z = fabs(z);
 
210
                        z = PI/(z * stirf(q));
 
211
                }
 
212
                else
 
213
                {
 
214
                        z = stirf(x);
 
215
                }
 
216
                return (*sgngam * z);
 
217
        }
 
218
 
 
219
        z = 1.0;
 
220
        while (x >= 3.0)
 
221
        {
 
222
                x -= 1.0;
 
223
                z *= x;
 
224
        }
 
225
 
 
226
        while (x < 0.0)
 
227
        {
 
228
                if (x > -1.E-9)
 
229
                        goto Small;
 
230
                z /= x;
 
231
                x += 1.0;
 
232
        }
 
233
 
 
234
        while (x < 2.0)
 
235
        {
 
236
                if (x < 1.e-9)
 
237
                        goto Small;
 
238
                z /= x;
 
239
                x += 1.0;
 
240
        }
 
241
 
 
242
        if (x == 2.0)
 
243
                return (z);
 
244
 
 
245
        x -= 2.0;
 
246
        p = polevl( x, P, 6 );
 
247
        q = polevl( x, Q, 7 );
 
248
        return (z * p / q);
 
249
 
 
250
Small:
 
251
        if (x == 0.0)
 
252
        {
 
253
                goto gsing;
 
254
        }
 
255
        else
 
256
                return (z/((1.0 + 0.5772156649015329 * x) * x));
 
257
}
 
258
 
 
259
/* This is the C99 version */
 
260
double tgamma(double x)
 
261
{
 
262
        int local_sgngam = 0;
 
263
        return (__tgamma_r(x, &local_sgngam));
 
264
}
 
265