overby.me
this repo has no description
1/*
2 * tgamma.c
3 *
4 * by Ian Ollmann, drawing heavily from work by Steve Peters, Jeff Kidder, and Ali Sazegari.
5 *
6 * Copyright (c) 2007 Apple Inc. All Rights Resterved.
7 */
8
9#include <math.h>
10#include <stdint.h>
11
12#ifdef ARMLIBM_SET_FLAGS
13#include "required_arithmetic.h"
14#endif
15
16static const double LogSqrt2pi = 0.9189385332046727417803297e+0;
17static const double pi = 3.1415926535897932384626434e+0;
18static const double eps = 0x1.0p-52;
19static const double xbig = 171.624e+0;
20
21
22static const double p[8] = { -1.71618513886549492533811e+0,
23 2.47656508055759199108314e+1,
24 -3.79804256470945635097577e+2,
25 6.29331155312818442661052e+2,
26 8.66966202790413211295064e+2,
27 -3.14512729688483675254357e+4,
28 -3.61444134186911729807069e+4,
29 6.64561438202405440627855e+4 };
30
31static const double q[8] = { -3.08402300119738975254353e+1,
32 3.15350626979604161529144e+2,
33 -1.01515636749021914166146e+3,
34 -3.10777167157231109440444e+3,
35 2.25381184209801510330112e+4,
36 4.75584627752788110767815e+3,
37 -1.34659959864969306392456e+5,
38 -1.15132259675553483497211e+5 };
39
40static const double c[7] = { 0xa.aaaaaaaaaaaaaabp-7,
41 -0xb.60b60b60b60b60bp-12,
42 0xd.00d00d00d00d00dp-14,
43 -0x9.c09c09c09c09c0ap-14,
44 0xd.ca8f158c7f91ab8p-14,
45 -0xf.b5586ccc9e3e41p-13,
46 0xd.20d20d20d20d20dp-11 };
47
48double tgamma( double x )
49{
50 union{ double d; uint64_t u; }u = {x};
51
52 // deal with 0, NaN, Inf
53 if( (u.u & 0x7fffffffffffffffULL) -1ULL >= 0x7ff0000000000000ULL -1ULL )
54 {
55 // 0
56 if( x == 0.0 ) {
57#ifdef ARMLIBM_SET_FLAGS
58 return required_divide_double( 1.0, x );
59#else
60 return __builtin_copysign(__builtin_inf(), x);
61#endif
62 }
63
64 //NaN
65 if( x != x )
66 return x + x;
67
68 //Inf
69 if( x < 0.0 ) {
70#ifdef ARMLIBM_SET_FLAGS
71 return required_multiply_double( x, 0.0 ); // -Inf: set invalid, return NaN
72#else
73 return __builtin_nan("");
74#endif
75 }
76
77 return x; // Inf
78 }
79
80 int n = 0;
81 int parity = 0;
82 double y = x;
83 double fact = 1.0;
84
85 // negative x
86 if( x <= 0.0 )
87 {
88 y = -x;
89 if( y < 0x1.0p-1022 )
90 return 1.0 / x;
91
92 // y1 = trunc ( y ); //fast trunc, no edge cases
93 double iy = y;
94 int64_t ll = 0;
95 if( y < 0x1.0p54 )
96 {
97 ll = y;
98 iy = ll;
99 }
100
101 double fraction = y - iy;
102
103 if( fraction != 0.0 )
104 {
105 parity = (int) (ll & 1);
106 fact = -pi / sin( pi * fraction );
107 y += 1.0;
108 }
109 else {
110#ifdef ARMLIBM_SET_FLAGS
111 return required_add_double( -__builtin_inf(), __builtin_inf() ); // set invalid, return NaN
112#else
113 return __builtin_nan("");
114#endif
115 }
116 }
117
118 double result;
119 int i;
120
121 if( y < eps )
122 {
123 result = 1.0 / y;
124 }
125 else if( y < 12.0 )
126 {
127 double y1 = y;
128 double z;
129 if( y < 1.0)
130 {
131 z = y;
132 y += 1.0;
133 }
134 else
135 {
136 n = (int)y - 1;
137 y -= (double) n;
138 z = y - 1.0;
139 }
140
141 double numerator = 0.0;
142 double denomenator = 1.0;
143 for( i = 0; i < 8; i++ )
144 {
145 numerator = (numerator + p[i]) * z;
146 denomenator = denomenator * z + q[i];
147 }
148 result = numerator / denomenator + 1.0;
149 if( y1 < y )
150 result /= y1;
151 else if( y1 > y )
152 {
153 for( i = 0; i < n; i++ )
154 {
155 result *= y;
156 y += 1.0;
157 }
158 }
159 }
160 else
161 { // y large
162 if( x < xbig )
163 {
164 double yy = y * y;
165 double sum = c[6];
166 for( i = 5; i >= 0; i-- )
167 sum = sum / yy + c[i];
168 sum = sum / y - y + LogSqrt2pi;
169 sum += (y - 0.5) * log( y );
170 result = exp( sum );
171 }
172 else
173 return x * 0x1.0p1023; //set overflow, return inf
174 }
175
176 if( parity )
177 result = -result;
178
179 if( fact != 1.0 )
180 result = fact / result;
181
182 return result;
183}