1 /* -------------------------------------------------------------- */
2 /* (C)Copyright 2007,2008, */
3 /* International Business Machines Corporation */
4 /* All Rights Reserved. */
5 /* */
6 /* Redistribution and use in source and binary forms, with or */
7 /* without modification, are permitted provided that the */
8 /* following conditions are met: */
9 /* */
10 /* - Redistributions of source code must retain the above copyright*/
11 /* notice, this list of conditions and the following disclaimer. */
12 /* */
13 /* - Redistributions in binary form must reproduce the above */
14 /* copyright notice, this list of conditions and the following */
15 /* disclaimer in the documentation and/or other materials */
16 /* provided with the distribution. */
17 /* */
18 /* - Neither the name of IBM Corporation nor the names of its */
19 /* contributors may be used to endorse or promote products */
20 /* derived from this software without specific prior written */
21 /* permission. */
22 /* */
23 /* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND */
24 /* CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, */
25 /* INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF */
26 /* MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE */
27 /* DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR */
28 /* CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, */
29 /* SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT */
30 /* NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; */
31 /* LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) */
32 /* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN */
33 /* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR */
34 /* OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, */
35 /* EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. */
36 /* -------------------------------------------------------------- */
37 /* PROLOG END TAG zYx */
38 #ifdef __SPU__
39
40 #ifndef _TGAMMAD2_H_
41 #define _TGAMMAD2_H_ 1
42
43 #include <spu_intrinsics.h>
44 #include "simdmath.h"
45
46 #include "recipd2.h"
47 #include "truncd2.h"
48 #include "expd2.h"
49 #include "logd2.h"
50 #include "divd2.h"
51 #include "sind2.h"
52 #include "powd2.h"
53
54
55 /*
56 * FUNCTION
57 * vector double _tgammad2(vector double x)
58 *
59 * DESCRIPTION
60 * _tgammad2
61 *
62 * This is an interesting function to approximate fast
63 * and accurately. We take a fairly standard approach - break
64 * the domain into 5 separate regions:
65 *
66 * 1. [-infinity, 0) - use
67 * 2. [0, 1) - push x into [1,2), then adjust the
68 * result.
69 * 3. [1, 2) - use a rational approximation.
70 * 4. [2, 10) - pull back into [1, 2), then adjust
71 * the result.
72 * 5. [10, +infinity] - use Stirling's Approximation.
73 *
74 *
75 * Special Cases:
76 * - tgamma(+/- 0) returns +/- infinity
77 * - tgamma(negative integer) returns NaN
78 * - tgamma(-infinity) returns NaN
79 * - tgamma(infinity) returns infinity
80 *
81 */
82
83
84 /*
85 * Coefficients for Stirling's Series for Gamma()
86 */
87 /* 1/ 1 */
88 #define STIRLING_00 1.000000000000000000000000000000000000E0
89 /* 1/ 12 */
90 #define STIRLING_01 8.333333333333333333333333333333333333E-2
91 /* 1/ 288 */
92 #define STIRLING_02 3.472222222222222222222222222222222222E-3
93 /* -139/ 51840 */
94 #define STIRLING_03 -2.681327160493827160493827160493827160E-3
95 /* -571/ 2488320 */
96 #define STIRLING_04 -2.294720936213991769547325102880658436E-4
97 /* 163879/ 209018880 */
98 #define STIRLING_05 7.840392217200666274740348814422888497E-4
99 /* 5246819/ 75246796800 */
100 #define STIRLING_06 6.972813758365857774293988285757833083E-5
101 /* -534703531/ 902961561600 */
102 #define STIRLING_07 -5.921664373536938828648362256044011874E-4
103 /* -4483131259/ 86684309913600 */
104 #define STIRLING_08 -5.171790908260592193370578430020588228E-5
105 /* 432261921612371/ 514904800886784000 */
106 #define STIRLING_09 8.394987206720872799933575167649834452E-4
107 /* 6232523202521089/ 86504006548979712000 */
108 #define STIRLING_10 7.204895416020010559085719302250150521E-5
109 /* -25834629665134204969/ 13494625021640835072000 */
110 #define STIRLING_11 -1.914438498565477526500898858328522545E-3
111 /* -1579029138854919086429/ 9716130015581401251840000 */
112 #define STIRLING_12 -1.625162627839158168986351239802709981E-4
113 /* 746590869962651602203151/ 116593560186976815022080000 */
114 #define STIRLING_13 6.403362833808069794823638090265795830E-3
115 /* 1511513601028097903631961/ 2798245444487443560529920000 */
116 #define STIRLING_14 5.401647678926045151804675085702417355E-4
117 /* -8849272268392873147705987190261/ 299692087104605205332754432000000 */
118 #define STIRLING_15 -2.952788094569912050544065105469382445E-2
119 /* -142801712490607530608130701097701/ 57540880724084199423888850944000000 */
120 #define STIRLING_16 -2.481743600264997730915658368743464324E-3
121
122
123 /*
124 * Rational Approximation Coefficients for the
125 * domain [1, 2).
126 */
127 #define TGD2_P00 -1.8211798563156931777484715e+05
128 #define TGD2_P01 -8.7136501560410004458390176e+04
129 #define TGD2_P02 -3.9304030489789496641606092e+04
130 #define TGD2_P03 -1.2078833505605729442322627e+04
131 #define TGD2_P04 -2.2149136023607729839568492e+03
132 #define TGD2_P05 -7.2672456596961114883015398e+02
133 #define TGD2_P06 -2.2126466212611862971471055e+01
134 #define TGD2_P07 -2.0162424149396112937893122e+01
135
136 #define TGD2_Q00 1.0000000000000000000000000
137 #define TGD2_Q01 -1.8212849094205905566923320e+05
138 #define TGD2_Q02 -1.9220660507239613798446953e+05
139 #define TGD2_Q03 2.9692670736656051303725690e+04
140 #define TGD2_Q04 3.0352658363629092491464689e+04
141 #define TGD2_Q05 -1.0555895821041505769244395e+04
142 #define TGD2_Q06 1.2786642579487202056043316e+03
143 #define TGD2_Q07 -5.5279768804094054246434098e+01
144
_tgammad2(vector double x)145 static __inline vector double _tgammad2(vector double x)
146 {
147 vector double signbit = spu_splats(-0.0);
148 vector double zerod = spu_splats(0.0);
149 vector double halfd = spu_splats(0.5);
150 vector double oned = spu_splats(1.0);
151 vector double ninep9d = (vec_double2)spu_splats(0x4023FFFFFFFFFFFFull);
152 vector double twohd = spu_splats(200.0);
153 vector double pi = spu_splats(SM_PI);
154 vector double sqrt2pi = spu_splats(2.50662827463100050241576528481);
155 vector double inf = (vector double)spu_splats(0x7FF0000000000000ull);
156 vector double nan = (vector double)spu_splats(0x7FF8000000000000ull);
157
158
159 vector double xabs;
160 vector double xscaled;
161 vector double xtrunc;
162 vector double xinv;
163 vector double nresult;
164 vector double rresult; /* Rational Approx result */
165 vector double sresult; /* Stirling's result */
166 vector double result;
167 vector double pr,qr;
168
169 vector unsigned long long gt0 = spu_cmpgt(x, zerod);
170 vector unsigned long long gt1 = spu_cmpgt(x, oned);
171 vector unsigned long long gt9p9 = spu_cmpgt(x, ninep9d);
172 vector unsigned long long gt200 = spu_cmpgt(x, twohd);
173
174
175 xabs = spu_andc(x, signbit);
176
177 /*
178 * For x in [0, 1], add 1 to x, use rational
179 * approximation, then use:
180 *
181 * gamma(x) = gamma(x+1)/x
182 *
183 */
184 xabs = spu_sel(spu_add(xabs, oned), xabs, gt1);
185 xtrunc = _truncd2(xabs);
186
187
188 /*
189 * For x in [2, 10):
190 */
191 xscaled = spu_add(oned, spu_sub(xabs, xtrunc));
192
193 /*
194 * For x in [1,2), use a rational approximation.
195 */
196 pr = spu_madd(xscaled, spu_splats(TGD2_P07), spu_splats(TGD2_P06));
197 pr = spu_madd(pr, xscaled, spu_splats(TGD2_P05));
198 pr = spu_madd(pr, xscaled, spu_splats(TGD2_P04));
199 pr = spu_madd(pr, xscaled, spu_splats(TGD2_P03));
200 pr = spu_madd(pr, xscaled, spu_splats(TGD2_P02));
201 pr = spu_madd(pr, xscaled, spu_splats(TGD2_P01));
202 pr = spu_madd(pr, xscaled, spu_splats(TGD2_P00));
203
204 qr = spu_madd(xscaled, spu_splats(TGD2_Q07), spu_splats(TGD2_Q06));
205 qr = spu_madd(qr, xscaled, spu_splats(TGD2_Q05));
206 qr = spu_madd(qr, xscaled, spu_splats(TGD2_Q04));
207 qr = spu_madd(qr, xscaled, spu_splats(TGD2_Q03));
208 qr = spu_madd(qr, xscaled, spu_splats(TGD2_Q02));
209 qr = spu_madd(qr, xscaled, spu_splats(TGD2_Q01));
210 qr = spu_madd(qr, xscaled, spu_splats(TGD2_Q00));
211
212 rresult = _divd2(pr, qr);
213 rresult = spu_sel(_divd2(rresult, x), rresult, gt1);
214
215 /*
216 * If x was in [2,10) and we pulled it into [1,2), we need to push
217 * it back out again.
218 */
219 rresult = spu_sel(rresult, spu_mul(rresult, xscaled), spu_cmpgt(x, xscaled)); /* [2,3) */
220 xscaled = spu_add(xscaled, oned);
221 rresult = spu_sel(rresult, spu_mul(rresult, xscaled), spu_cmpgt(x, xscaled)); /* [3,4) */
222 xscaled = spu_add(xscaled, oned);
223 rresult = spu_sel(rresult, spu_mul(rresult, xscaled), spu_cmpgt(x, xscaled)); /* [4,5) */
224 xscaled = spu_add(xscaled, oned);
225 rresult = spu_sel(rresult, spu_mul(rresult, xscaled), spu_cmpgt(x, xscaled)); /* [5,6) */
226 xscaled = spu_add(xscaled, oned);
227 rresult = spu_sel(rresult, spu_mul(rresult, xscaled), spu_cmpgt(x, xscaled)); /* [6,7) */
228 xscaled = spu_add(xscaled, oned);
229 rresult = spu_sel(rresult, spu_mul(rresult, xscaled), spu_cmpgt(x, xscaled)); /* [7,8) */
230 xscaled = spu_add(xscaled, oned);
231 rresult = spu_sel(rresult, spu_mul(rresult, xscaled), spu_cmpgt(x, xscaled)); /* [8,9) */
232 xscaled = spu_add(xscaled, oned);
233 rresult = spu_sel(rresult, spu_mul(rresult, xscaled), spu_cmpgt(x, xscaled)); /* [9,10) */
234
235
236 /*
237 * For x >= 10, we use Stirling's Approximation
238 */
239 vector double sum;
240 xinv = _recipd2(xabs);
241 sum = spu_madd(xinv, spu_splats(STIRLING_16), spu_splats(STIRLING_15));
242 sum = spu_madd(sum, xinv, spu_splats(STIRLING_14));
243 sum = spu_madd(sum, xinv, spu_splats(STIRLING_13));
244 sum = spu_madd(sum, xinv, spu_splats(STIRLING_12));
245 sum = spu_madd(sum, xinv, spu_splats(STIRLING_11));
246 sum = spu_madd(sum, xinv, spu_splats(STIRLING_10));
247 sum = spu_madd(sum, xinv, spu_splats(STIRLING_09));
248 sum = spu_madd(sum, xinv, spu_splats(STIRLING_08));
249 sum = spu_madd(sum, xinv, spu_splats(STIRLING_07));
250 sum = spu_madd(sum, xinv, spu_splats(STIRLING_06));
251 sum = spu_madd(sum, xinv, spu_splats(STIRLING_05));
252 sum = spu_madd(sum, xinv, spu_splats(STIRLING_04));
253 sum = spu_madd(sum, xinv, spu_splats(STIRLING_03));
254 sum = spu_madd(sum, xinv, spu_splats(STIRLING_02));
255 sum = spu_madd(sum, xinv, spu_splats(STIRLING_01));
256 sum = spu_madd(sum, xinv, spu_splats(STIRLING_00));
257
258 sum = spu_mul(sum, sqrt2pi);
259 sum = spu_mul(sum, _powd2(x, spu_sub(x, halfd)));
260 sresult = spu_mul(sum, _expd2(spu_or(x, signbit)));
261
262 /*
263 * Choose rational approximation or Stirling's result.
264 */
265 result = spu_sel(rresult, sresult, gt9p9);
266
267
268 result = spu_sel(result, inf, gt200);
269
270 /* For x < 0, use:
271 *
272 * gamma(x) = pi/(x*gamma(-x)*sin(x*pi))
273 * or
274 * gamma(x) = pi/(gamma(1 - x)*sin(x*pi))
275 */
276 nresult = _divd2(pi, spu_mul(x, spu_mul(result, _sind2(spu_mul(x, pi)))));
277 result = spu_sel(nresult, result, gt0);
278
279 /*
280 * x = non-positive integer, return NaN.
281 */
282 result = spu_sel(result, nan, spu_andc(spu_cmpeq(x, xtrunc), gt0));
283
284
285 return result;
286 }
287
288 #endif /* _TGAMMAD2_H_ */
289 #endif /* __SPU__ */
290