home *** CD-ROM | disk | FTP | other *** search
/ Il CD di internet / CD.iso / SOURCE / D / LIBC / LIBC-4.6 / LIBC-4 / libc-linux / sysdeps / linux / i386 / math / lgamma.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-11-25  |  3.4 KB  |  156 lines
  1. /*
  2.  * Copyright (c) 1985 Regents of the University of California.
  3.  * All rights reserved.  The Berkeley software License Agreement
  4.  * specifies the terms and conditions for redistribution.
  5.  */
  6.  
  7. #if 0
  8. #ifndef lint
  9. static char sccsid[] = "@(#)lgamma.c    5.3 (Berkeley) 9/22/88";
  10. #endif /* not lint */
  11. #endif
  12.  
  13. /*
  14.     C program for floating point log Gamma function
  15.  
  16.     lgamma(x) computes the log of the absolute
  17.     value of the Gamma function.
  18.     The sign of the Gamma function is returned in the
  19.     external quantity signgam.
  20.  
  21.     The coefficients for expansion around zero
  22.     are #5243 from Hart & Cheney; for expansion
  23.     around infinity they are #5404.
  24.  
  25.     Calls log, floor and sin.
  26. */
  27.  
  28. #if 0
  29. #include "mathimpl.h"
  30. #else
  31. #include <math.h>
  32. #endif
  33. #if defined(vax)||defined(tahoe)
  34. #include <errno.h>
  35. #endif    /* defined(vax)||defined(tahoe) */
  36.  
  37. int signgam;
  38.  
  39. static const double goobie = 0.9189385332046727417803297;  /* log(2*pi)/2 */
  40. static const double pi       = 3.1415926535897932384626434;
  41.  
  42. #define M 6
  43. #define N 8
  44. static const double p1[] = {
  45.     0.83333333333333101837e-1,
  46.     -.277777777735865004e-2,
  47.     0.793650576493454e-3,
  48.     -.5951896861197e-3,
  49.     0.83645878922e-3,
  50.     -.1633436431e-2,
  51. };
  52. static const double p2[] = {
  53.     -.42353689509744089647e5,
  54.     -.20886861789269887364e5,
  55.     -.87627102978521489560e4,
  56.     -.20085274013072791214e4,
  57.     -.43933044406002567613e3,
  58.     -.50108693752970953015e2,
  59.     -.67449507245925289918e1,
  60.     0.0,
  61. };
  62. static const double q2[] = {
  63.     -.42353689509744090010e5,
  64.     -.29803853309256649932e4,
  65.     0.99403074150827709015e4,
  66.     -.15286072737795220248e4,
  67.     -.49902852662143904834e3,
  68.     0.18949823415702801641e3,
  69.     -.23081551524580124562e2,
  70.     0.10000000000000000000e1,
  71. };
  72.  
  73. static double pos(double), neg(double), asym(double);
  74.  
  75. double
  76. lgamma(double arg)
  77. {
  78.  
  79.     signgam = 1;
  80.     if(arg <= 0.) return(neg(arg));
  81.     if(arg > 8.) return(asym(arg));
  82.     return(log(pos(arg)));
  83. }
  84.  
  85. static double
  86. asym(double arg)
  87. {
  88.     double n, argsq;
  89.     int i;
  90.  
  91.     /* Avoid overflow of arg*arg.  */
  92.     n = 0.0;
  93.     if(arg > 1.0e8)
  94.       goto noasym;
  95.  
  96.     argsq = 1./(arg*arg);
  97.     for(i=M-1; i>=0; i--){
  98.         n = n*argsq + p1[i];
  99.     }
  100. noasym:
  101.     return((arg-.5)*log(arg) - arg + goobie + n/arg);
  102. }
  103.  
  104. static double
  105. neg(double arg)
  106. {
  107.     double t;
  108.  
  109.     arg = -arg;
  110.      /*
  111.       * to see if arg were a true integer, the old code used the
  112.       * mathematically correct observation:
  113.       * sin(n*pi) = 0 <=> n is an integer.
  114.       * but in finite precision arithmetic, sin(n*PI) will NEVER
  115.       * be zero simply because n*PI is a rational number.  hence
  116.       *    it failed to work with our newer, more accurate sin()
  117.       * which uses true pi to do the argument reduction...
  118.       *    temp = sin(pi*arg);
  119.       */
  120.     t = floor(arg);
  121.     if (arg - t  > 0.5e0)
  122.         t += 1.e0;                /* t := integer nearest arg */
  123. #if defined(vax)||defined(tahoe)
  124.     if (arg == t) {
  125.         return(infnan(ERANGE));        /* +INF */
  126.     }
  127. #endif    /* defined(vax)||defined(tahoe) */
  128.     signgam = (int) (t - 2*floor(t/2));    /* signgam =  1 if t was odd, */
  129.                         /*            0 if t was even */
  130.     signgam = signgam - 1 + signgam;    /* signgam =  1 if t was odd, */
  131.                         /*           -1 if t was even */
  132.     t = arg - t;                /*  -0.5 <= t <= 0.5 */
  133.     if (t < 0.e0) {
  134.         t = -t;
  135.         signgam = -signgam;
  136.     }
  137.     return(-log(arg*pos(arg)*sin(pi*t)/pi));
  138. }
  139.  
  140. static double
  141. pos(double arg)
  142. {
  143.     double n, d, s;
  144.     register i;
  145.  
  146.     if(arg < 2.) return(pos(arg+1.)/arg);
  147.     if(arg > 3.) return((arg-1.)*pos(arg-1.));
  148.  
  149.     s = arg - 2.;
  150.     for(n=0,d=0,i=N-1; i>=0; i--){
  151.         n = n*s + p2[i];
  152.         d = d*s + q2[i];
  153.     }
  154.     return(n/d);
  155. }
  156.