overby.me
this repo has no description
1/*
2 * Copyright (c) 2002 Apple Computer, Inc. All rights reserved.
3 *
4 * @APPLE_LICENSE_HEADER_START@
5 *
6 * The contents of this file constitute Original Code as defined in and
7 * are subject to the Apple Public Source License Version 1.1 (the
8 * "License"). You may not use this file except in compliance with the
9 * License. Please obtain a copy of the License at
10 * http://www.apple.com/publicsource and read it before using this file.
11 *
12 * This Original Code and all software distributed under the License are
13 * distributed on an "AS IS" basis, WITHOUT WARRANTY OF ANY KIND, EITHER
14 * EXPRESS OR IMPLIED, AND APPLE HEREBY DISCLAIMS ALL SUCH WARRANTIES,
15 * INCLUDING WITHOUT LIMITATION, ANY WARRANTIES OF MERCHANTABILITY,
16 * FITNESS FOR A PARTICULAR PURPOSE OR NON-INFRINGEMENT. Please see the
17 * License for the specific language governing rights and limitations
18 * under the License.
19 *
20 * @APPLE_LICENSE_HEADER_END@
21 */
22// GammaDD.c
23//
24// Double-double Function Library
25// Copyright: � 1996 by Apple Computer, Inc., all rights reserved
26//
27// long double gammal( long double x );
28// long double lgammal( long double x );
29//
30
31#include "math.h"
32#include "fp_private.h"
33#if (__WANT_LONG_DOUBLE_FORMAT__ - 0L == 128L)
34#include "DD.h"
35
36// Requires the following functions:
37
38void _d3div(double a, double b, double c, double d, double e, double f,
39 double *high, double *mid, double *low);
40void _d3mul(double a, double b, double c, double d, double e, double f,
41 double *high, double *mid, double *low);
42void _d3add(double a, double b, double c, double d, double e, double f,
43 double *high, double *mid, double *low);
44
45long double _LogInnerLD(double, double, double, double *, int);
46long double _ExpInnerLD(double alpha, double beta, double gamma, double *extra,
47 int exptype);
48long double rintl(long double x);
49long double sinl(long double x);
50
51typedef double DD[2];
52
53// Coefficients for asymptotic formula:
54// real*8 bern(1:2,1:11)
55// FORTRAN: bern(i,j) C: bern[j][i-1]
56// DD *bern = (DD *)(bernData - 4);
57
58static const uint32_t bernData[] __attribute__ ((aligned(8))) = {
59 0x3FB55555, 0x55555555, 0x3C555555, 0x55555555,
60 0xBF66C16C, 0x16C16C17, 0x3BFF49F4, 0x9F49F49F,
61 0x3F4A01A0, 0x1A01A01A, 0x3B8A01A0, 0x1A01A01A,
62 0xBF438138, 0x13813814, 0x3BEFB1FB, 0x1FB1FB20,
63 0x3F4B951E, 0x2B18FF23, 0x3BE5C3A9, 0xCE01B952,
64 0xBF5F6AB0, 0xD9993C7D, 0x3BFF8255, 0x3C999B0E,
65 0x3F7A41A4, 0x1A41A41A, 0x3C106906, 0x90690690,
66 0xBF9E4286, 0xCB0F5398, 0x3C21EFCD, 0xAB896745,
67 0x3FC6FE96, 0x381E0680, 0xBC279E24, 0x05A71F88,
68 0xBFF64767, 0x01181F3A, 0x3C724246, 0x319DA678,
69 0x402ACE44, 0x322CE006, 0xBCC62C2B, 0x1BBCDD32
70};
71
72// Coefficients of approximating polynomial for LnGamma(2+x)
73// real*8 zeta(1:2,2:53)
74// FORTRAN: zeta(i,j) C: zeta[j][i-1]
75// DD *zeta = (DD *)(zetaData - 8);
76
77static const uint32_t zetaData[] __attribute__ ((aligned(8))) = {
78 0x3FD4A34C, 0xC4A60FA6, 0x3C71873D, 0x8912200C,
79 0xBFB13E00, 0x1A557607, 0x3C5FB68B, 0xE2F8821F,
80 0x3F951322, 0xAC7D8483, 0x3C3AFC89, 0x088CB729,
81 0xBF7E404F, 0xC218F5F2, 0x3C1E4A62, 0x7CF1EB34,
82 0x3F67ADD6, 0xEADB6C30, 0xBBF5B782, 0x8C7FD7F4,
83 0xBF538AC5, 0xC2BF8E08, 0x3BE8A4C1, 0xCFD9CEC8,
84 0x3F40B36A, 0xF86396E9, 0xBBE0698D, 0x6C892967,
85 0xBF2D3FD4, 0xC76D2FC8, 0x3BBC7C55, 0xCFCCBB83,
86 0x3F1A127B, 0x0F17D65A, 0x3BA9D309, 0xAA700268,
87 0xBF078DE5, 0xBD7C81EF, 0x3B7A2054, 0x1CDE47A6,
88 0x3EF580DC, 0xEE66EB02, 0x3B826057, 0x4B258F72,
89 0xBEE3CBC9, 0x63CE2243, 0x3B8EA56E, 0x6C7D5329,
90 0x3ED2597A, 0x39F34AAC, 0xBB7BF911, 0x462A7D81,
91 0xBEC11B2E, 0xB7679541, 0xBB4C76B0, 0xE65AC63A,
92 0x3EB0064C, 0xDEB22F0F, 0x3B4D0156, 0xAFFDBC11,
93 0xBE9E2600, 0xD93CFD2F, 0x3B3130AC, 0x39E5C106,
94 0x3E8C76BB, 0xB3F07A4D, 0x3B2D9A2B, 0x77769B52,
95 0xBE7AF5A6, 0xCBBF8A97, 0xBB195F22, 0x7E96D83E,
96 0x3E699B93, 0xC2070B0F, 0x3B003271, 0x64736428,
97 0xBE5862C7, 0x34DF3EAC, 0xBAFB3280, 0x2BEC0DA0,
98 0x3E47469D, 0xACCFADCD, 0xBAE369D3, 0x88CEBAA9,
99 0xBE36434A, 0x8447AEAD, 0xBA8AF72E, 0xDF876FCD,
100 0x3E2555A8, 0x77FFD2C3, 0xBAC87506, 0x5F26A43B,
101 0xBE147B16, 0x79258D0E, 0xBAB04F36, 0xE0E854E4,
102 0x3E03B15D, 0x2B2FC10C, 0xBA9D79F6, 0xFEEEB28B,
103 0xBDF2F69A, 0x9FABE3E0, 0x3A9A162A, 0xB374C789,
104 0x3DE24932, 0xA337434C, 0x3A806082, 0x9C24508F,
105 0xBDD1A7C2, 0x6EC2523C, 0xBA74F4EB, 0xDB4A04B5,
106 0x3DC11116, 0xE693ED98, 0xBA6C7034, 0xD49E7FC7,
107 0xBDB08424, 0xCBC543D8, 0xBA440EF8, 0x20DBC9EA,
108 0x3DA00002, 0x6E3F644F, 0x3A43546A, 0x6054C889,
109 0xBD8F07C5, 0x14FC9F0A, 0xB9F75B6B, 0xE545AC09,
110 0x3D7E1E20, 0x26AAFCD8, 0xBA162A85, 0x86538620,
111 0xBD6D41D5, 0x6E5EE2E2, 0x39F43894, 0xD27CED5E,
112 0x3D5C71C7, 0xF6F10E37, 0xB9F01074, 0x764D33F2,
113 0xBD4BACF9, 0xA27BC89B, 0x39D4A5A2, 0x15E0508E,
114 0x3D3AF287, 0x18A10D6E, 0x39C40D7F, 0x1B842CC3,
115 0xBD2A41A4, 0x5603E5B6, 0x39C62BE9, 0xCF212D90,
116 0x3D199999, 0xC0716EE9, 0xB9B39E10, 0xF90435CB,
117 0xBD08F9C1, 0xA8DF9D78, 0x3989DA56, 0xD4471922,
118 0x3CF86186, 0x28D28905, 0xB989D7D4, 0xEE5A8893,
119 0xBCE7D05F, 0x4C31C560, 0xB9871BBA, 0x0B7CC339,
120 0x3CD745D1, 0x7B56BA4A, 0x3979D38B, 0xC00D6F02,
121 0xBCC6C16C, 0x1B4D6456, 0xB96AED17, 0x2E5C90F2,
122 0x3CB642C8, 0x5C023D9D, 0xB95DE052, 0x190D755D,
123 0xBCA5C988, 0x2D825E9D, 0x3949723F, 0x1BF240B7,
124 0x3C955555, 0x5698A866, 0x392CF5C8, 0x64970970,
125 0xBC84E5E0, 0xA8022BC9, 0x39128B9D, 0xC88F5BE2,
126 0x3C747AE1, 0x4838081F, 0xB91DF461, 0x306465C0,
127 0xBC641414, 0x146E3E31, 0xB8FE4773, 0xEA1309D6,
128 0x3C53B13B, 0x13EC2F3A, 0x38C41C5B, 0x07A32CA9,
129 0xBC43521C, 0xFB520859, 0xB8E225B1, 0x0AA3CE51
130};
131
132static const double zero = 0.0;
133static const Double pi = {{0x400921FB, 0x54442D18}};
134static const Double pib = {{0x3CA1A626, 0x33145C07}};
135static const Double pic = {{0xB92F1976, 0xB7ED8FBC}};
136static const Double infinity = {{0x7ff00000, 0x00000000}};
137static const Double scaleup = {{0x4ff00000, 0x00000000}};
138static const Double tiny = {{0x0E000000, 0x00000000}};
139static const Double sclog = {{0x40662E42, 0xFEFA39EF}};
140static const Double sclog2 = {{0x3CFABC9E, 0x3B39803F}};
141static const Double sclog3 = {{0x3987B57A, 0x079A1934}};
142
143static const Double bump = {{0x3FED67F1, 0xC864BEB5}}; // 0.5 ln (2 Pi)
144static const Double bump2 = {{0xBC865B5A, 0x1B7FF5DF}};
145static const Double bump3 = {{0xB91B7F70, 0xC13DC168}};
146
147static const Double ec = {{0x3FDB0EE6, 0x072093CE}}; // 1 - Euler's constant
148static const Double ec2 = {{0x3C56CB90, 0x701FBFAB}};
149static const Double ec3 = {{0x38F34A95, 0xE3133C51}};
150
151static const Double piln = {{0x3FF250D0, 0x48E7A1BD}}; // log(pi)
152static const Double pilnb = {{0x3C67ABF2, 0xAD8D5088}};
153static const Double pilnc = {{0xB8E6CCF4, 0x32447B52}};
154
155static const Double zeta32 = {{0xB914C685, 0x28DDC956}}; // Extra word of precision
156 // for (zeta(2,1),zeta(2,2))
157
158static long double GammaCommon( double dhead, double dtail, int igamma )
159{
160 int signgam, ireflect, increase, i;
161
162 double int1, int2, arg, arg2, arg3, argt, arg2t, anew;
163 double factor, factor1, factor2, factor3, f1, f2, f3;
164 double prod1, prod2, pextra, pextra1, pextra2, sum1, sum2, sextra;
165 double sum21, sum22, sextra2, sum31, sum32, sextra3;
166 double diff1, diff2, dextra, sarg1, sarg2, saextra;
167 double denom1, denom2, sum, suml, sumt, recip1, recip2, recipsq1, recipsq2;
168 double extra, extra2, extra3, extrax, extray, extralg, gmextra;
169 double temp, eexp, diff, z, z2, y, y2, y3, asize, c;
170
171 DblDbl gamdd, tempdd, intdd, xhintdd, sinargdd;
172
173 DD *bern = (DD *)(bernData - 4);
174 DD *zeta = (DD *)(zetaData - 8);
175
176 signgam = 1;
177
178 if (dhead <= 0) { // special for neg arguments
179 tempdd.d[0] = dhead;
180 tempdd.d[1] = dtail;
181 intdd.f = rintl(tempdd.f); // get integer part of argument
182 int1 = intdd.d[0];
183 int2 = intdd.d[1];
184
185 if (isinf(dhead) || (dhead - int1 + dtail - int2) == 0) {
186
187 //*********************************************************************
188 // Argument is a non-positive integer. Gamma is not defined. *
189 // Return a NaN. *
190 //*********************************************************************
191 gamdd.d[0] = INFINITY; // per IEC, [was: zero/zero;]
192 gamdd.d[1] = zero;
193 return gamdd.f;
194 }
195 if (dhead > -20.0) {
196 ireflect = 0; // reflection formula not used for
197 // modest sized neg numbers
198 }
199 else {
200 ireflect = 1; // signal use of reflection formula
201 dhead = -dhead; // change sign of argument
202 dtail = -dtail;
203 }
204 }
205 else // else, positive arguments --
206 ireflect = 0; // signal reflection formula not used
207
208 if (isinf(dhead) || (dhead + dtail) != (dhead + dtail)) {
209 gamdd.d[0] = dhead + dtail;
210 gamdd.d[1] = zero;
211 return gamdd.f;
212 }
213
214 arg = dhead; // working copy of argument
215 arg2 = dtail;
216 arg3 = 0.0;
217 factor1 = 1.0;
218 factor2 = 0.0;
219 factor3 = 0.0;
220
221 if (arg > 18.0) { // use asymptotic formula
222 while (arg < 31.0) {
223 _d3mul( factor1, factor2, factor3, arg, arg2, arg3,
224 &factor1, &factor2, &factor3 );
225 arg = arg + 1.0;
226 } // argument in range for asympt. formula
227 while (arg < 35.0) {
228 _d3mul( factor1, factor2, factor3, arg, arg2, arg3,
229 &factor1, &factor2, &factor3 );
230 argt = __FADD( arg, 1.0 );
231 arg2t = __FADD( (arg - argt + 1.0), arg2 );
232 arg3 = ((arg - argt + 1.0) - arg2t + arg2) + arg3;
233 arg = argt;
234 arg2 = arg2t;
235 }
236 tempdd.f = _LogInnerLD(arg, arg2, 0.0, &f3, 1); // ln x
237 f3 = f3 + arg3/arg;
238 f1 = tempdd.d[0];
239 f2 = tempdd.d[1];
240
241 _d3mul(f1, f2, f3, (arg - 0.5), arg2, arg3,
242 &prod1, &prod2, &c); // (x-.5)ln x
243 _d3add(prod1, prod2, c, -arg, -arg2, -arg3,
244 &sum1, &sum2, &sextra); // (x-.5)ln x-x
245 _d3add(sum1, sum2, sextra, bump.f, bump2.f, bump3.f,
246 &sum21, &sum22, &sextra2);
247
248 //**************************************************************
249 // sum21, sum22, sextra2 represent (x-.5)ln x-x+.5 ln(2 Pi) *
250 //**************************************************************
251
252 _d3div(1.0, 0.0, 0.0, arg, arg2, arg3, &recip1, &recip2, &extra);
253
254 // argument for asymptotic formula
255 _d3mul(recip1, recip2, extra, recip1, recip2, extra,
256 &recipsq1, &recipsq2, &extra2);
257
258 sum = bern[11][0];
259 suml = bern[11][1];
260 for (i = 10; i >= 1; --i) {
261 temp = __FMADD( sum, recipsq1, bern[i][0] ); // bern[i][0] + sum*recipsq1;
262 /* suml = bern[i][0] - temp + sum*recipsq1 +
263 (bern[i][1] + sum*recipsq2 + suml*recipsq1); */
264 suml = __FMADD( sum, recipsq1, bern[i][0] - temp ) +
265 __FMADD( suml, recipsq1, __FMADD( sum, recipsq2, bern[i][1] ) );
266 sum = temp;
267 } // finish summation of asymptotic series
268
269 _d3mul(sum, suml, 0.0, recip1, recip2, extra, &prod1, &prod2, &extra3);
270 _d3add(sum21, sum22, sextra2, prod1, prod2, extra3, &sum31, &sum32, &sextra3);
271
272 //******************************************************************************
273 // At this point, log(gamma*factor) is computed as (sum31, sum32, sextra3). *
274 //******************************************************************************
275
276 if ((igamma == 1) && (ireflect != 1)) {
277 // Gamma entry point (without use of reflection formula)?
278 tempdd.f = _ExpInnerLD(sum31, sum32, sextra3, &eexp, 4);
279 if (factor1 == 1.0) {
280 gamdd.f = tempdd.f;
281 gmextra = eexp;
282 }
283 else {
284 _d3div(tempdd.d[0], tempdd.d[1], eexp, factor1, factor2, factor3,
285 &(gamdd.d[0]), &(gamdd.d[1]), &gmextra);
286 }
287 }
288 else { // Computing log(gamma(x))
289 factor = factor1;
290 if (factor == 1.0) {
291 gamdd.d[0] = sum31;
292 gamdd.d[1] = sum32;
293 gmextra = sextra3;
294 }
295 else {
296 tempdd.f = _LogInnerLD(factor1, factor2, 0.0, &extrax, 1);
297 // computing log gamma.
298 extray = -(extrax + factor3/factor);
299 _d3add(sum31, sum32, sextra3, -tempdd.d[0], -tempdd.d[1], extray,
300 &(gamdd.d[0]), &(gamdd.d[1]), &gmextra);
301 }
302 }
303 }
304 else { // use formula for interval [0.5,1.5]
305 arg = dhead;
306 arg2 = dtail; // working copy of argument
307 arg3 = 0.0;
308 increase = 0; // signal-> no scaling for result
309 if (arg < 1.5) { // scale up argument by recursion formula
310 increase = -1; // signal-> result to be divided by factor
311 factor1 = arg;
312 factor2 = arg2;
313 factor3 = 0.0; // factor is the old argument
314 while (arg < 1.5) {
315 _d3add(arg, arg2, arg3, 1.0, 0.0, 0.0, &arg, &arg2, &arg3);
316 if (arg < 1.5) {
317 _d3mul(factor1, factor2, factor3, arg, arg2, arg3,
318 &factor1, &factor2, &factor3);
319 }
320 }
321 if (factor1 < 0.0) {
322 signgam = -1; // special case of negative arguments
323 }
324 }
325 else if (arg > 2.5) {
326 increase = +1; // signal-> result must be mult by factor
327 factor1 = 1.0;
328 factor2 = 0.0;
329 factor3 = 0.0;
330 while (arg > 2.5) { // reduce argument by 1
331 arg = arg - 1.0; // there may be room for bits to
332 anew = __FADD( arg, arg2 ); // shift from low order word to
333 arg2 = arg - anew + arg2; // high order word
334 arg = anew;
335 _d3mul(factor1, factor2, factor3, arg, arg2, 0.0,
336 &factor1, &factor2, &factor3);
337 }
338 arg3 = 0.0;
339 }
340 diff = __FSUB( arg, 2.0 ); // series in terms of
341 z = __FADD( (arg - (diff + 2.0)), arg2 ); // (diff,z,x2)=arg-2
342 z2 = (arg - (diff + 2.0)) - z + arg2 + arg3;
343 y = __FADD( diff, z );
344 y2 = (diff - y + z) + z2;
345 y3 = (diff - y + z) - y2 + z2;
346 sum = zeta[53][0];
347 suml = zeta[53][1];
348 for (i = 52; i >= 40; --i) {
349 sum = __FMADD( sum, y, zeta[i][0] ); // zeta[i][0] + sum*y;
350 }
351 sumt = __FMADD( sum, y, zeta[39][0] ); // zeta[39][0] + sum*y;
352 suml = __FMADD( sum, y, zeta[39][0] - sumt) + __FMADD( sum, y2, zeta[39][1] ); // zeta[39][0] - sumt + sum*y + (zeta[39][1] + sum*y2);
353 sum = sumt;
354 for (i = 38; i >= 3; --i) {
355 temp = __FMADD( sum, y, zeta[i][0] ); // zeta[i][0] + sum*y;
356 // suml = (zeta[i][0] - temp + sum*y) + zeta[i][1] + (sum*y2 + suml*y);
357 suml = __FMADD( sum, y, zeta[i][0] - temp ) + zeta[i][1] + __FMADD( sum, y2, suml*y );
358 sum = temp;
359 }
360 _d3mul(sum, suml, 0.0, y, y2, y3, &prod1, &prod2, &pextra2);
361 _d3add(prod1, prod2, pextra2, zeta[2][0], zeta[2][1], zeta32.f,
362 &sum1, &sum2, &sextra);
363 _d3mul(sum1, sum2, sextra, y, y2, y3,
364 &prod1, &prod2, &pextra1); // multiply sum by z
365 _d3add(prod1, prod2, pextra1, ec.f, ec2.f, ec3.f,
366 &sum1, &sum2, &sextra); // add linear term
367 _d3mul(sum1, sum2, sextra, y, y2, y3,
368 &prod1, &prod2, &pextra); // final mult. by z.
369
370 if (igamma == 1) { // a Gamma entry
371 tempdd.f = _ExpInnerLD(prod1, prod2, pextra, &eexp, 4);
372 if (increase == 1) {
373 _d3mul(tempdd.d[0], tempdd.d[1], eexp, factor1, factor2, factor3,
374 &(gamdd.d[0]), &(gamdd.d[1]), &gmextra);
375 }
376 else if (increase == -1) {
377 _d3div(tempdd.d[0], tempdd.d[1], eexp, factor1, factor2, factor3,
378 &(gamdd.d[0]), &(gamdd.d[1]), &gmextra);
379 }
380 else {
381 gamdd.f = tempdd.f;
382 gmextra = eexp;
383 }
384 }
385 else { // entry was for log gamma
386 if (increase == 0) {
387 gamdd.d[0] = prod1;
388 gamdd.d[1] = prod2;
389 gmextra = pextra;
390 }
391 else {
392 if (signgam < 0) {
393 factor1 = -factor1;
394 factor2 = -factor2;
395 factor3 = -factor3;
396 }
397 factor = factor1;
398 tempdd.f = _LogInnerLD(factor1, factor2, 0.0, &extra, 1);
399 extra = extra + factor3/factor;
400 if (increase == -1) { // change sign of log
401 tempdd.d[0] = -tempdd.d[0];
402 tempdd.d[1] = -tempdd.d[1];
403 extra = -extra;
404 }
405 _d3add(prod1, prod2, pextra, tempdd.d[0], tempdd.d[1], extra,
406 &(gamdd.d[0]), &(gamdd.d[1]), &gmextra);
407 }
408 }
409 }
410
411 if (gamdd.d[0] != gamdd.d[0]) {
412 gamdd.d[0] = infinity.f;
413 gamdd.d[1] = zero;
414 gmextra = 0.0;
415 }
416
417 if (ireflect == 1) { // original argument less than 0.0
418 arg = -dhead; // recover original argument
419 arg2 = -dtail;
420 // reduce argument for computation
421 _d3add(arg, arg2, 0.0, -int1, -int2, 0.0, &diff1, &diff2, &dextra);
422 asize = fabs(diff1);
423
424 if (asize < tiny.f) { // For arguments very close
425 diff1 *= scaleup.f; // to an integer, rescale,
426 diff2 *= scaleup.f; // so that sin can be computed
427 dextra *= scaleup.f; // to a full 107+ bits
428 } // of sin (Pi x)
429
430 _d3mul(diff1, diff2, dextra, pi.f, pib.f, pic.f, &sarg1, &sarg2, &saextra);
431 xhintdd.f = 2.0*rintl(0.5*intdd.f);
432
433 tempdd.d[0] = sarg1;
434 tempdd.d[1] = sarg2;
435 sinargdd.f = sinl(tempdd.f); // sin of argument
436
437 if ((xhintdd.d[0] - intdd.d[0] + xhintdd.d[1] - intdd.d[1]) != 0.0) {
438 sinargdd.f = -sinargdd.f;
439 }
440
441 if (sinargdd.d[0] < 0.0) {
442 sinargdd.f = -sinargdd.f; // force sin(pi x) to have plus sign
443 signgam = -signgam; // show gamma(x) has negative sign
444 }
445
446 if (fabs(gamdd.d[0]) == infinity.f){ // result will underflow
447 if(igamma == 1) {
448 gamdd.d[0] = zero;
449 gamdd.d[1] = zero;
450 }
451 else {
452 gamdd.d[0] = -infinity.f; // gamma underflows, so
453 gamdd.d[1] = zero; // lgamma overflows to -INF
454 }
455 return gamdd.f;
456 }
457 _d3mul(dhead, dtail, 0.0, sinargdd.d[0], sinargdd.d[1], 0.0,
458 &prod1, &prod2, &pextra); // x sin(pi x)
459
460 tempdd.f = _LogInnerLD(prod1, prod2, 0.0, &extralg, 1); // log (x sin(pi x))
461 extralg = extralg + pextra/prod1;
462 _d3add(tempdd.d[0], tempdd.d[1], extralg, gamdd.d[0], gamdd.d[1], gmextra,
463 &denom1, &denom2, &dextra);
464 _d3add(piln.f, pilnb.f, pilnc.f, -denom1, -denom2, -dextra,
465 &(gamdd.d[0]), &(gamdd.d[1]), &gmextra);
466
467 if (asize < tiny.f) {
468 // Compensate for scaling of argument to sin(pi x)
469 _d3add(gamdd.d[0], gamdd.d[1], gmextra, sclog.f, sclog2.f, sclog3.f,
470 &(gamdd.d[0]), &(gamdd.d[1]), &gmextra);
471 }
472 if (igamma == 1) { // we really want gamma itself ?
473 tempdd.f = _ExpInnerLD(gamdd.d[0], gamdd.d[1], gmextra, &eexp, 4);
474
475 if (signgam == 1) {
476 gamdd.f = tempdd.f;
477 }
478 else {
479 gamdd.f = -tempdd.f;
480 }
481 }
482 }
483 return gamdd.f;
484}
485
486long double tgammal( long double x )
487{
488 double head, fenv;
489 DblDbl t;
490
491 if (x == 0.0L)
492 return copysignl( INFINITY, x );
493
494 FEGETENVD(fenv);
495 FESETENVD(0.0);
496
497 t.f = x;
498
499 t.f = GammaCommon(t.d[0], t.d[1], 1);
500
501 head = __FADD( t.d[0], t.d[1] );
502 t.d[1] = (t.d[0] - head) + t.d[1]; // cannonize tail
503 t.d[0] = head; // cannonize head
504
505 FESETENVD(fenv);
506 return t.f;
507}
508
509long double lgammal( long double x )
510{
511 double head, fenv;
512 DblDbl t;
513
514 FEGETENVD(fenv);
515 FESETENVD(0.0);
516
517 t.f = x;
518
519 t.f = GammaCommon(t.d[0], t.d[1], 0);
520
521 head = __FADD( t.d[0], t.d[1] );
522 t.d[1] = (t.d[0] - head) + t.d[1]; // cannonize tail
523 t.d[0] = head; // cannonize head
524
525 FESETENVD(fenv);
526 return t.f;
527}
528#endif