1 /******************************************************************************
2 Copyright (c) 2007-2011, Intel Corp.
3 All rights reserved.
4
5 Redistribution and use in source and binary forms, with or without
6 modification, are permitted provided that the following conditions are met:
7
8 * Redistributions of source code must retain the above copyright notice,
9 this list of conditions and the following disclaimer.
10 * Redistributions in binary form must reproduce the above copyright
11 notice, this list of conditions and the following disclaimer in the
12 documentation and/or other materials provided with the distribution.
13 * Neither the name of Intel Corporation nor the names of its contributors
14 may be used to endorse or promote products derived from this software
15 without specific prior written permission.
16
17 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
18 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
21 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23 SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24 INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25 CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26 ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
27 THE POSSIBILITY OF SUCH DAMAGE.
28 ******************************************************************************/
29
30 #ifndef BASE_NAME
31 # define BASE_NAME LGAMMA_BASE_NAME
32 #endif
33
34 #include "dpml_private.h"
35
36 #if (F_NAME_SUFFIX == DPML_NULL_MACRO_TOKEN)
37 # define SIGNGAM_NAME signgam
38 #else
39 # if (OP_SYSTEM == vms)
40 # define SIGNGAM_NAME PASTE_3(F_NAME_PREFIX,signgam, F_NAME_SUFFIX)
41 # else
42 # define SIGNGAM_NAME PASTE_2(signgam, F_NAME_SUFFIX)
43 # endif
44 #endif
45
46 #if HACK_GAMMAS_INLINE
47 # define SIGNGAM SIGNGAM_NAME
48 #else
49 # define SIGNGAM *SIGNGAM_NAME
50 #endif
51
52
53 /*
54 * Lgamma(x) is defined as the log(|gamma(x)|), where gamma(x) is defined
55 * for positive x as
56 *
57 * gamma(x) = integral{ 0 to infinity | t^(x-1)e^t dt }
58 *
59 * From the definition of gamma(x) it follows that
60 *
61 * x*gamma(x) = gamma(x+1) (1)
62 *
63 * and the limit as x --> +0 of gamma(x) = +infinity. Equation (1) can be
64 * used to extend gamma(x) to negitive numbers by recursively applying:
65 *
66 * gamma(-x) = gamma(1 - x)/(-x) (2)
67 *
68 * Since gamma(0) = + infinity, it follows that gamma(n) is undefined for
69 * any non-positive integer. An alternative extension of gamma to negative
70 * arguments is the reflection fomula
71 *
72 * gamma(-x) = -pi/(sin(pi*x)*gamma(1 + x)) (3)
73 *
74 * Evalutation of lgamma(x) suffers potential loss of significance at
75 * its zeros or alternatively, when |gamma(x)| = 1. From the definition
76 * of gamma and (1) we see that |gamma(x)| = 1 for positive x only at
77 * x = 1 and 2. From equation (2), we see that |gamma(x)| = 1 when x
78 * is a negative integer +/- epsilon, where epsilon is on the order
79 * of 1/n!.
80 *
81 * Computation of lgamma(x) is based on two identities:
82 *
83 * zeta(2)-1 zeta(3)-1
84 * lgamma(1+x) = (1-G)x - ln(1+x) + ---------x^2 - ---------x^3 + ...
85 * 2 3
86 *
87 * zeta(n)-1
88 * ---------(-x)^n ... (4)
89 * n
90 *
91 * = -ln(1+x) + x*Q(x)
92 *
93 * where G is Euler's constant and zeta(n) is the Reimann zeta function:
94 *
95 * 1 1 1
96 * zeta(n) = 1 + --- + --- + --- + ...
97 * 2^n 3^n 4^n
98 *
99 * and Stirlings asymtotic approximation to gamma(x):
100 *
101 * 1 1 1 1
102 * lgamma(x) ~ ---ln(2*pi) - x + (x - ---)*ln(x) + ---*phi(---) (5)
103 * 2 2 x x^2
104 *
105 * 1 B(2) B(4) B(6) B(8)
106 * phi(---) = ----- - ------- + ------- - ------- ..... (6)
107 * x^2 2*1 4*3*x^2 6*5*x^4 8*7*x^6
108 *
109 * where B(n) is the n-th bernoulli number.
110 */
111
112 # define TMP_FILE ADD_EXTENSION(BUILD_FILE_NAME,tmp)
113 # define RND_TO_FMT(x) bround(x, F_PRECISION)
114 # define PRINT_TABLE_ENTRY(a) PRINT_1_F_TYPE_ENTRY(a, offset)
115 # define PRINT_TABLE_VALUE_F_DEFINE(n) \
116 PRINT_TABLE_VALUE_DEFINE(n, TABLE_NAME, offset, F_TYPE)
117 # define F_PRINT_A_DEFINE(name) PRINT_TABLE_ADDRESS_DEFINE(name, \
118 TABLE_NAME, offset, F_TYPE)
119
120 #ifndef MAKE_INCLUDE
121
122 # include STR(BUILD_FILE_NAME)
123
124 #else
125
126 @divert divertText
127
128 precision = ceil(2*F_PRECISION/MP_RADIX_BITS) + 2;
129
130 /*
131 * The following macro defines an mphoc routine that finds a root of
132 * the function f between x_0 and x_1 to precsion p and returns the
133 * result in y.
134 */
135
136 # define FIND_ROOT(x_0, x_1, f, p, y) \
137 y_0 = f(x_0); y_1 = f(x_1); \
138 if (y_0 * y_1 > 0) \
139 { \
140 printf("Invalid input to FIND_ROOT\n"); \
141 exit; \
142 } \
143 while (1) \
144 { \
145 delta = x_1 - x_0; \
146 if (bexp(x_1) - bexp(delta) > p) \
147 break; \
148 x = x_1 - y_1*(delta/(y_1 - y_0)); \
149 x_0 = x_1; x_1 = x; \
150 y_0 = y_1; y_1 = f(x); \
151 } \
152 y = x_1
153
154
155 START_STATIC_TABLE(TABLE_NAME, offset);
156 TABLE_COMMENT("Miscelaneous constants");
157
158 /*
159 * lgamma(x) will overflow for large positive values of x. Note that
160 * For large negative values of x, x is a negative integer, and hence
161 * the function is not defined. To compute the overflow threshold
162 * we need to solve the equatation lgamma(x) = MP_MAX_FLOAT + 1/2 lsb
163 * and rounding down the result to working precision. We do this using
164 * the macro FIND_ROOT defined above with f = lgamma(x) - MP_MAX_FLOAT +
165 * 1/2 lsb. To find the starting values we note that lgamma is ~ x*log(x)
166 * and assume that x = 2^k/k*log(2). Then x*log(x) ~ MP_MAX_FLOAT when
167 * k = F_MAX_BIN_EXP
168 */
169
170 k = F_MAX_BIN_EXP;
171 x0 = 2^k/(k*log(2));
172 x1 = 3*x0;
173
174 c = MP_MAX_FLOAT + 2^(bexp(MP_MAX_FLOAT) - F_PRECISION);
f()175 function f() { return lgamma($1) - MP_MAX_FLOAT; }
176
177 FIND_ROOT(x0, x1, f, F_PRECISION + 1, y);
178 y = bchop(y, F_PRECISION);
179
180 PRINT_TABLE_VALUE_F_DEFINE(OVERFLOW_THRESHOLD);
181 PRINT_TABLE_ENTRY(y);
182
183 /*
184 * For large values of x, it is most efficient to use equation (5).
185 * When x is very large, 1/x^2 will underflow. However, long before
186 * the underflow threshold is reached, (1/x)*phi(1/x^2) will become
187 * insignificant when compared with the other terms in (5).
188 * Consequently, we should stop computing z(x) = (1/x)*phi(1/x^2) when
189 * x is big enough. This is more efficient and avoids the underflow.
190 *
191 * z(x) will be insignificant when z(x)/lgamma(x) < 1/2^(F_PRECISION + 1),
192 * or when
193 *
194 * (1 - 2^-(F_PRECISION+1))*lgamma(x)-.5*ln(2*pi)+x-(x-.5)*ln(x) < 0
195 *
196 * Using the macro, FIND_ROOT, we determine an x that satisfies the above.
197 */
198
199 a = 1 - 1/2^(F_PRECISION + 1);
200 b = .5*log(2*pi);
g()201 function g()
202 {
203 s = a*lgamma($1);
204 t = ($1 - .5)*log($1) - $1 + b;
205 return s - t;
206 }
207
208 k = .5*(F_PRECISION + 1 - log2(12.));
209 x0 = 2^k/sqrt(k*log(2));
210 x1 = x0 + x0;
211
212 FIND_ROOT(x0, x1, g, F_PRECISION + 1, real_big);
213 PRINT_TABLE_VALUE_F_DEFINE(REAL_BIG);
214 PRINT_TABLE_ENTRY(real_big);
215
216 /*
217 * Using equation (5) requires the constant .5*ln(2*pi)
218 */
219 y = .5*log(2*pi);
220 PRINT_TABLE_VALUE_F_DEFINE(HALF_LN_2_PI);
221 PRINT_TABLE_ENTRY(y);
222
223 /*
224 * For suitably large negative x, we would like to a computation
225 * based on equation (3).
226 *
227 * lgamma(-x) = ln|gamma(-x)|
228 * = ln|-pi/(sin(pi*x)*x*gamma(x))|
229 * = ln(pi) - ln|sin(pi*x)| - ln(x) - ln(gamma(x))
230 * = ln(pi) - ln|sin(pi*x)| - ln(x) - lgamma(x)
231 * = ln(pi) - ln|sin(pi*x)| - ln(x) - lgamma(x)
232 *
233 * combined with (5) this gives:
234 *
235 * lgamma(-x) ~ ln(pi) - ln|sin(pi*x)| - ln(x) -
236 * [.5*ln(2*pi) - x + (x - .5)*ln(x) + phi(x)/x]
237 * ~ .5*ln(pi/2) - ln|sin(pi*x)| + x - (x + .5)*ln(x) - phi(x)/x
238 *
239 * Consequently, we also need the constants .5*ln(pi/2) and pi
240 */
241
242 y = .5*log(pi/2);
243 PRINT_TABLE_VALUE_F_DEFINE(HALF_LN_PI_OVER_2);
244 PRINT_TABLE_ENTRY(y);
245
246 y = pi;
247 PRINT_TABLE_VALUE_F_DEFINE(PI);
248 PRINT_TABLE_ENTRY(y);
249
250 /*
251 * When x is not large, the computation of lgamma is based on equations
252 * (1) and (2). Specifically, let
253 *
254 * lgamma(n+x) = log(F(n,x)) + x*Q(x)
255 *
256 * where Q(x) is defined by equation (4). From equation (1) it follows
257 * that
258 *
259 * lgamma(n+1+x) = log((n+x)*gamma(n+x)
260 * = log(n+x) + lgamma(n+x)
261 * = log(n+x) + log(F(n,x)) + x*Q(x)
262 * = log[(n+x)*F(n,x)] + x*Q(x)
263 *
264 * From the above and equation (4) it follows that F(1,x) = 1+x and
265 * F(n+1, x) = (n+x)*F(n,x). Note the F(n,x) is define for both
266 * negative and positive integers.
267 *
268 * Since we know the range of our x value for this evaluation we can
269 * increase the accuracy of the computation of x*Q(x) by performing
270 * the following transformation:
271 *
272 * Given Q(x) = p(x)/q(x), define R(X) as
273 *
274 * Q(x) = 1/2 - R(x)
275 *
276 * This yields
277 *
278 * R(x) = p(x) - q(x)/2
279 * -------------
280 * q(x)
281 *
282 * Now x*Q(x) can be computed as x*(1/2 - R(x)), or rather x*(1/2) - x*R(x)
283 * which forces, x*(1/2), the most significant term, to be exact.
284 *
285 *
286 * NOTE: We need coefficients for Q and phi. From (4) we obtain
287 * Q by approximating
288 *
289 * (lgamma(1+x) + ln(1 + x))/x
290 *
291 * on the interval [-.5, .5]. A rational approximation for Q
292 * has competative performance on ALPHA with a polynomial
293 * approximation.
294 *
295 * From (5) we obtain phi by approximating
296 *
297 * x * [lgamma(x) - .5*ln(2*pi) + x - (x - .5)*ln(x)]
298 *
299 * on the range [8,max_val], where max_val is the largest
300 * value of x which will be evaluated by phi (i.e. for X>x, phi(x)
301 * is insignificant to the other terms of the sum in (5).
302 */
303
304 old_precision = precision;
305 precision = ceil(2*F_PRECISION/MP_RADIX_BITS) + 4;
306
lgamma_approx()307 function lgamma_approx()
308 {
309 if ($1 == 0)
310 return (1 - euler_gamma);
311 else {
312 /* logx1(x) is more accurate than ln(x) for |x| < 1/MP_RADIX */
313 if ( abs($1) < (1 / MP_RADIX))
314 return (lgamma(1+$1) + logx1($1))/$1;
315 else
316 return (lgamma(1+$1) + ln(1 + $1))/$1;
317 }
318 }
319
320 /* To shorten our search time we'll make some initial estimates based
321 on experience. (These estimates are on the low side to assure we
322 don't over step the optimal degree)
323 */
324 #if (F_PRECISION == 24)
325 degree = 3;
326 #elif (F_PRECISION == 53)
327 degree = 5;
328 #elif (F_PRECISION == 113)
329 degree = 11;
330 #else
331 degree = 0;
332 #endif
333
334 tol = 0;
335 while (tol < (F_PRECISION + 1 + 3)) {
336 den_degree = num_degree = ++degree;
337 tol = remes( REMES_STATIC + REMES_LINEAR_ARG + REMES_RELATIVE_WEIGHT,
338 -0.5, 0.5, lgamma_approx,
339 num_degree, den_degree, &rational_coefs);
340 }
341
342 precision = old_precision;
343
344
345 /* Extract denominator coefficients */
346 first_den_coef = num_degree + 1;
347 for (i = 0; i <= den_degree; i++)
348 q[i] = rational_coefs[i + first_den_coef];
349
350 /* Extract numerator coefficients */
351 for (i = 0; i <= num_degree; i++)
352 p[i] = rational_coefs[i] - q[i]/2;
353
354
355
356
357 /* Generate constants for Phi */
358
359 old_precision = precision;
360 precision = ceil(2*F_PRECISION/MP_RADIX_BITS) + 4;
361
362 half_ln_of_2pi = .5*ln(2*pi);
lgamma_asym_approx()363 function lgamma_asym_approx()
364 {
365 x = $1;
366 if (x == 0)
367 return (1/12); /* B2(0)/2 where B2(x) = x^2 - x + 1/6 */
368 else
369 return x*(lgamma(x) - half_ln_of_2pi + x - (x - .5)*ln(x));
370 }
371
372 max_arg = real_big;
373 #if QUAD_PRECISION
374 max_arg = 100000; /* This is temporary until mp_remes is corrected */
375 #endif
376 remes(REMES_FIND_POLYNOMIAL+ REMES_RELATIVE_WEIGHT+ REMES_RECIP_SQUARE_ARG,
377 8.0, max_arg, lgamma_asym_approx, (F_PRECISION + 1),
378 &poly_degree, &r);
379 precision = old_precision;
380
381
382
383 #define PRINT_COEFS(n,p) for (i = 0; i <= n; i++) \
384 { PRINT_TABLE_ENTRY(p[i]); }
385
386 TABLE_COMMENT("Rational Coefficents for Q(x)");
387 F_PRINT_A_DEFINE(P_COEFS);
388 PRINT_COEFS(num_degree, p);
389 printf("\n");
390
391 F_PRINT_A_DEFINE(Q_COEFS);
392 PRINT_COEFS(den_degree, q);
393
394 TABLE_COMMENT("Polynomial Coefficents phi(x)");
395 F_PRINT_A_DEFINE(PHI_COEFS);
396 PRINT_COEFS(poly_degree, r);
397
398 END_TABLE;
399
400 /*
401 * Print out defines for polynomial and rational approximations
402 */
403
404 printf("#define PHI(a,u) u = a*a; u = a*POLY%i(PHI_COEFS, u)\n",
405 poly_degree);
406 printf("#define Q(x) (POLY%i(P_COEFS, x)/POLY%i(Q_COEFS, x))\n",
407 num_degree, den_degree);
408
409 @end_divert
410 @eval my $outText = MphocEval( GetStream( "divertText" ) ); \
411 my $defineText = Egrep( "#define", $outText, \$tableText ); \
412 my $headerText = GetHeaderText( STR(BUILD_FILE_NAME), \
413 "Definitions and constants for " . \
414 STR(F_ENTRY_NAME), __FILE__); \
415 print "$headerText\n\n$tableText\n\n$defineText";
416 #endif
417
418
419 #if IEEE_FLOATING
420 # define SCREEN_SPECIAL_ARGS(x,i) GET_EXP_WORD(x, i); \
421 if (F_EXP_WORD_IS_ABNORMAL(i)) \
422 goto special_args
423 #else
424 # define SCREEN_SPECIAL_ARGS(x,i)
425 #endif
426
427 #define NEG_2_POW_F_PRECISION ALIGN_W_EXP_FIELD(F_PRECISION + F_EXP_BIAS - F_NORM) + \
428 F_SIGN_BIT_MASK
429
430 #ifdef F_COPY_SIGN_FAST
431 # define F_SET_SIGN(val, sign, res) F_COPY_SIGN(val, sign, result)
432 #else
433 # define F_SET_SIGN(val, sign, res) res = val; if ((sign) < 0) F_NEGATE(res)
434 #endif
435
436
437 #if DO_LGAMMA
438
439 int SIGNGAM_NAME = 0;
440 # define _F_ENTRY_NAME F_LGAMMA_NAME
441 # define OPT_PTR_ARG
442 # define USE_CALL !HACK_GAMMAS_INLINE
443
444 #elif DO_GAMMA
445
446 extern int SIGNGAM_NAME;
447 # define _F_ENTRY_NAME F_GAMMA_NAME
448 # define OPT_PTR_ARG
449 # define USE_CALL !HACK_GAMMAS_INLINE
450
451 #else
452
453 # define _F_ENTRY_NAME F_RT_LGAMMA_NAME
454 # define OPT_PTR_ARG , int *SIGNGAM_NAME
455 # undef HACK_GAMMAS_INLINE
456 # define USE_CALL 0
457
458 #endif
459
460 #if !defined F_ENTRY_NAME
461 # define F_ENTRY_NAME _F_ENTRY_NAME
462 #endif
463
464 F_F_PROTO( F_LN_NAME ) ;
465 F_F_PROTO( F_SIN_NAME ) ;
466
467 F_TYPE
F_ENTRY_NAME(F_TYPE x OPT_PTR_ARG)468 F_ENTRY_NAME(F_TYPE x OPT_PTR_ARG)
469 {
470 F_TYPE y;
471 WORD i;
472
473 #if USE_CALL
474
475 F_FpI_PROTO( F_RT_LGAMMA_NAME ) ;
476
477 y = F_RT_LGAMMA_NAME(x, &i);
478 SIGNGAM_NAME = i;
479 return y;
480
481 #else
482
483 EXCEPTION_RECORD_DECLARATION
484 F_TYPE s, t;
485
486 SIGNGAM = 1;
487
488 /* screen for NaNs, infinities, zeros & denorms */
489 SCREEN_SPECIAL_ARGS(x, i);
490
491 /*
492 * Initialize the SIGNGAM to 1 and send large arguments to asymtotic
493 * region. Note the choice of asymtotic region being |x| >= 8 is
494 * fairly arbitrary and need not be symetric. As the lower bound of
495 * the asymtotic region increases, the more multiplies are performed
496 * in computing F(n,x). Eventually, it is faster to use the asymtotic
497 * approximations. Experimentally, it appears that the asymtotic
498 * regions are not as accurate. However, that might be caused by a
499 * sloppy implementation in that region.
500 *
501 * For the non-asymtotic region, we need to compute rint(x). Get 1/2
502 * with the correct sign now.
503 */
504
505 F_SET_SIGN((F_TYPE) .5, x, y);
506
507 if (x >= (F_TYPE) 8.)
508 goto pos_asymtotic;
509 if (x <= (F_TYPE) -8.)
510 goto neg_asymtotic;
511
512 /* For small x, get i = rint(x), y = x - i */
513
514 i = (WORD)(x + y);
515 y = x - (F_TYPE) i;
516 t = (F_TYPE) 1.;
517
518 /*
519 * Compute F(n,x) and take its log. In most cases this switch statement
520 * is faster than a loop.
521 */
522
523 switch (i)
524 {
525 case -8:
526 t *= (y - 8);
527 /* Fall through */
528
529 case -7:
530 t *= (y - 7);
531 /* Fall through */
532
533 case -6:
534 t *= (y - 6);
535 /* Fall through */
536
537 case -5:
538 t *= (y - 5);
539 /* Fall through */
540
541 case -4:
542 t *= (y - 4);
543 /* Fall through */
544
545 case -3:
546 t *= (y - 3);
547 /* Fall through */
548
549 case -2:
550 t *= (y - 2);
551 /* Fall through */
552
553 case -1:
554 t *= (y - 1);
555 /* Fall through */
556
557 case 0:
558 /*
559 * Since all of the negative cases come through here, we need
560 * to check for integer values and set signgam correctly;
561 */
562 t *= (y*(y+1));
563 if (y == 0) goto non_pos_int;
564 if (t < 0)
565 SIGNGAM = -1;
566 F_ABS(t, t);
567 t = - F_LN_NAME(t);
568 goto pos_eval;
569
570 case 1:
571 t = - F_LN_NAME(x);
572 goto pos_eval;
573
574 case 2:
575 t = 0;
576 goto pos_eval;
577
578 case 8:
579 t *= (x-6);
580 /* Fall through */
581
582 case 7:
583 t *= (x-5);
584 /* Fall through */
585
586 case 6:
587 t *= (x-4);
588 /* Fall through */
589
590 case 5:
591 t *= (x-3);
592 /* Fall through */
593
594 case 4:
595 t *= (x - 2);
596 /* Fall through */
597
598 case 3:
599 t *= (x - 1);
600 t = F_LN_NAME(t);
601 goto pos_eval;
602 }
603
604 pos_eval:
605 /*
606 * OK - just need to compute rational approximation and we're done.
607 */
608 t = t + (y*0.5 + y*Q(y));
609 return t;
610
611 pos_asymtotic:
612
613 /*
614 * In this region we compute lgamma using an asymtotic expansion.
615 * If x is really big, we don't need phi(x), so we can skip it.
616 */
617 t = HALF_LN_2_PI;
618 if (x > REAL_BIG) goto skip_poly;
619 y = 1/x;
620 PHI(y, s);
621 t += s;
622
623 add_in_log:
624 y = F_LN_NAME(x);
625 s = x * (y - 1);
626 s -= (F_TYPE) .5 * y;
627 t += s;
628
629 return t;
630
631 skip_poly:
632 /* If x is reaally, really big, result will overflow */
633 if (x <= OVERFLOW_THRESHOLD)
634 goto add_in_log;
635
636 GET_EXCEPTION_RESULT_1(LGAMMA_OVERFLOW, x, t);
637 return t;
638
639 neg_asymtotic:
640
641 /*
642 * Here we are dealing with large negative arguments we need to
643 * determine an integer n, such that n <= x < n+1. The parity
644 * of n determines whether SIGNGAM is + or - 1. Also, we are
645 * going to compute log(|sin(pi*x)|). If we can find and integer
646 * k such that k = rint(x) and define y = x - k, then log(|sin(pi*x)|)
647 * = log(sin(|y|*pi)). We begin by using "+ big - big" to determine
648 * k and y.
649 */
650
651 x = -x;
652 s = F_POW_2(F_PRECISION - 1);
653 if (x >= s)
654 /* x is so big that it must be an integer */
655 goto non_pos_int;
656 y = x + s;
657
658 /*
659 * get the low fraction bits of y. These are the same as the low
660 * bits of k
661 */
662 GET_LO_FRAC_WORD(y,i);
663 i = PDP_SHUFFLE(i);
664 t = y - s;
665 y = x - t;
666
667 /* Figure out n so we can set signgam correctly, and get |y| */
668 if (y < 0)
669 {
670 i--;
671 t--;
672 y = -y;
673 }
674 if (x == t) goto non_pos_int;
675 SIGNGAM = ((i + i) & 2) - 1;
676
677 /* OK compute aymtotic polynomial approximation for lgamma(|x|) */
678 s = ((F_TYPE) 1.)/x;
679 PHI(s, t);
680 t = HALF_LN_PI_OVER_2 - t;
681
682 /* Get log(|sin(pi*x)|) and remainder of asymtotic approximation */
683 s = F_LN_NAME(F_SIN_NAME(y*PI));
684 t = x + (t - s);
685 s = (x + (F_TYPE) .5)*F_LN_NAME(x);
686 t = t - s;
687 return t;
688
689
690
691 special_args:
692
693 #if IEEE_FLOATING
694
695 /* Note: The code below assumes that SIGNGAM has already been set to 1.
696 Thus, we only bother to set it here when gamma(x) is known to be
697 negative.
698 */
699
700 F_CLASSIFY(x, i);
701 switch (i)
702 {
703 case F_C_POS_INF:
704 GET_EXCEPTION_RESULT_1(LGAMMA_POS_INF, x, t);
705 break;
706
707 case F_C_NEG_INF:
708 GET_EXCEPTION_RESULT_1(LGAMMA_NEG_INF, x, t);
709 break;
710
711 case F_C_QUIET_NAN:
712 case F_C_SIG_NAN:
713 t = x;
714 break;
715
716 case F_C_NEG_ZERO:
717 SIGNGAM = -1;
718 /* fall through */
719
720 case F_C_POS_ZERO:
721 GET_EXCEPTION_RESULT_1(LGAMMA_OF_ZERO, x, t);
722 break;
723
724 default: /* +-denorm */
725 if (F_C_IS_NEG_CLASS(i)) {
726 SIGNGAM = -1;
727 F_ABS(x, x);
728 }
729 t = -F_LN_NAME(x);
730
731 }
732 return t;
733
734 #endif
735
736 non_pos_int:
737 GET_EXCEPTION_RESULT_1(LGAMMA_NON_POS_INT, -x, t);
738 return t;
739
740 #endif
741 }
742
743