home *** CD-ROM | disk | FTP | other *** search
/ InfoMagic Source Code 1993 July / THE_SOURCE_CODE_CD_ROM.iso / bsd_srcs / lib / libm / common_source / lgamma.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-12-16  |  8.3 KB  |  309 lines
  1. /*-
  2.  * Copyright (c) 1992 The Regents of the University of California.
  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
  7.  * are met:
  8.  * 1. Redistributions of source code must retain the above copyright
  9.  *    notice, this list of conditions and the following disclaimer.
  10.  * 2. 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.  * 3. All advertising materials mentioning features or use of this software
  14.  *    must display the following acknowledgement:
  15.  *    This product includes software developed by the University of
  16.  *    California, Berkeley and its contributors.
  17.  * 4. Neither the name of the University nor the names of its contributors
  18.  *    may be used to endorse or promote products derived from this software
  19.  *    without specific prior written permission.
  20.  *
  21.  * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
  22.  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  23.  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  24.  * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
  25.  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  26.  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  27.  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  28.  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  29.  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  30.  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  31.  * SUCH DAMAGE.
  32.  */
  33.  
  34. #ifndef lint
  35. static char sccsid[] = "@(#)lgamma.c    5.11 (Berkeley) 12/16/92";
  36. #endif /* not lint */
  37.  
  38. /*
  39.  * Coded by Peter McIlroy, Nov 1992;
  40.  *
  41.  * The financial support of UUNET Communications Services is greatfully
  42.  * acknowledged.
  43.  */
  44.  
  45. #include <math.h>
  46. #include <errno.h>
  47.  
  48. #include "mathimpl.h"
  49.  
  50. /* Log gamma function.
  51.  * Error:  x > 0 error < 1.3ulp.
  52.  *       x > 4, error < 1ulp.
  53.  *       x > 9, error < .6ulp.
  54.  *        x < 0, all bets are off. (When G(x) ~ 1, log(G(x)) ~ 0)
  55.  * Method:
  56.  *    x > 6:
  57.  *        Use the asymptotic expansion (Stirling's Formula)
  58.  *    0 < x < 6:
  59.  *        Use gamma(x+1) = x*gamma(x) for argument reduction.
  60.  *        Use rational approximation in
  61.  *        the range 1.2, 2.5
  62.  *        Two approximations are used, one centered at the
  63.  *        minimum to ensure monotonicity; one centered at 2
  64.  *        to maintain small relative error.
  65.  *    x < 0:
  66.  *        Use the reflection formula,
  67.  *        G(1-x)G(x) = PI/sin(PI*x)
  68.  * Special values:
  69.  *    non-positive integer    returns +Inf.
  70.  *    NaN            returns NaN
  71. */
  72. static int endian;
  73. #if defined(vax) || defined(tahoe)
  74. #define _IEEE        0
  75. /* double and float have same size exponent field */
  76. #define TRUNC(x)    x = (double) (float) (x)
  77. #else
  78. #define _IEEE        1
  79. #define TRUNC(x)    *(((int *) &x) + endian) &= 0xf8000000
  80. #define infnan(x)    0.0
  81. #endif
  82.  
  83. extern double log1p(double);
  84. static double small_lgam(double);
  85. static double large_lgam(double);
  86. static double neg_lgam(double);
  87. static double zero = 0.0, one = 1.0;
  88. int signgam;
  89.  
  90. #define UNDERFL (1e-1020 * 1e-1020)
  91.  
  92. #define LEFT    (1.0 - (x0 + .25))
  93. #define RIGHT    (x0 - .218)
  94. /*
  95. /* Constants for approximation in [1.244,1.712]
  96. */
  97. #define x0    0.461632144968362356785
  98. #define x0_lo    -.000000000000000015522348162858676890521
  99. #define a0_hi    -0.12148629128932952880859
  100. #define a0_lo    .0000000007534799204229502
  101. #define r0    -2.771227512955130520e-002
  102. #define r1    -2.980729795228150847e-001
  103. #define r2    -3.257411333183093394e-001
  104. #define r3    -1.126814387531706041e-001
  105. #define r4    -1.129130057170225562e-002
  106. #define r5    -2.259650588213369095e-005
  107. #define s0     1.714457160001714442e+000
  108. #define s1     2.786469504618194648e+000
  109. #define s2     1.564546365519179805e+000
  110. #define s3     3.485846389981109850e-001
  111. #define s4     2.467759345363656348e-002
  112. /*
  113.  * Constants for approximation in [1.71, 2.5]
  114. */
  115. #define a1_hi    4.227843350984671344505727574870e-01
  116. #define a1_lo    4.670126436531227189e-18
  117. #define p0    3.224670334241133695662995251041e-01
  118. #define p1    3.569659696950364669021382724168e-01
  119. #define p2    1.342918716072560025853732668111e-01
  120. #define p3    1.950702176409779831089963408886e-02
  121. #define p4    8.546740251667538090796227834289e-04
  122. #define q0    1.000000000000000444089209850062e+00
  123. #define q1    1.315850076960161985084596381057e+00
  124. #define q2    6.274644311862156431658377186977e-01
  125. #define q3    1.304706631926259297049597307705e-01
  126. #define q4    1.102815279606722369265536798366e-02
  127. #define q5    2.512690594856678929537585620579e-04
  128. #define q6    -1.003597548112371003358107325598e-06
  129. /*
  130.  * Stirling's Formula, adjusted for equal-ripple. x in [6,Inf].
  131. */
  132. #define lns2pi    .418938533204672741780329736405
  133. #define pb0     8.33333333333333148296162562474e-02
  134. #define pb1    -2.77777777774548123579378966497e-03
  135. #define pb2     7.93650778754435631476282786423e-04
  136. #define pb3    -5.95235082566672847950717262222e-04
  137. #define pb4     8.41428560346653702135821806252e-04
  138. #define pb5    -1.89773526463879200348872089421e-03
  139. #define pb6     5.69394463439411649408050664078e-03
  140. #define pb7    -1.44705562421428915453880392761e-02
  141.  
  142. double
  143. lgamma(double x)
  144. {
  145.     double r;
  146.  
  147.     signgam = 1;
  148.     endian = ((*(int *) &one)) ? 1 : 0;
  149.  
  150.     if (!finite(x))
  151.         if (_IEEE)
  152.             return (x+x);
  153.         else return (infnan(EDOM));
  154.  
  155.     if (x > 6 + RIGHT) {
  156.         r = large_lgam(x);
  157.         return (r);
  158.     } else if (x > 1e-16)
  159.         return (small_lgam(x));
  160.     else if (x > -1e-16) {
  161.         if (x < 0)
  162.             signgam = -1, x = -x;
  163.         return (-log(x));
  164.     } else
  165.         return (neg_lgam(x));
  166. }
  167.  
  168. static double
  169. large_lgam(double x)
  170. {
  171.     double z, p, x1;
  172.     int i;
  173.     struct Double t, u, v;
  174.     u = log__D(x);
  175.     u.a -= 1.0;
  176.     if (x > 1e15) {
  177.         v.a = x - 0.5;
  178.         TRUNC(v.a);
  179.         v.b = (x - v.a) - 0.5;
  180.         t.a = u.a*v.a;
  181.         t.b = x*u.b + v.b*u.a;
  182.         if (_IEEE == 0 && !finite(t.a))
  183.             return(infnan(ERANGE));
  184.         return(t.a + t.b);
  185.     }
  186.     x1 = 1./x;
  187.     z = x1*x1;
  188.     p = pb0+z*(pb1+z*(pb2+z*(pb3+z*(pb4+z*(pb5+z*(pb6+z*pb7))))));
  189.                     /* error in approximation = 2.8e-19 */
  190.  
  191.     p = p*x1;            /* error < 2.3e-18 absolute */
  192.                     /* 0 < p < 1/64 (at x = 5.5) */
  193.     v.a = x = x - 0.5;
  194.     TRUNC(v.a);            /* truncate v.a to 26 bits. */
  195.     v.b = x - v.a;
  196.     t.a = v.a*u.a;            /* t = (x-.5)*(log(x)-1) */
  197.     t.b = v.b*u.a + x*u.b;
  198.     t.b += p; t.b += lns2pi;    /* return t + lns2pi + p */
  199.     return (t.a + t.b);
  200. }
  201.  
  202. static double
  203. small_lgam(double x)
  204. {
  205.     int x_int;
  206.     double y, z, t, r = 0, p, q, hi, lo;
  207.     struct Double rr;
  208.     x_int = (x + .5);
  209.     y = x - x_int;
  210.     if (x_int <= 2 && y > RIGHT) {
  211.         t = y - x0;
  212.         y--; x_int++;
  213.         goto CONTINUE;
  214.     } else if (y < -LEFT) {
  215.         t = y +(1.0-x0);
  216. CONTINUE:
  217.         z = t - x0_lo;
  218.         p = r0+z*(r1+z*(r2+z*(r3+z*(r4+z*r5))));
  219.         q = s0+z*(s1+z*(s2+z*(s3+z*s4)));
  220.         r = t*(z*(p/q) - x0_lo);
  221.         t = .5*t*t;
  222.         z = 1.0;
  223.         switch (x_int) {
  224.         case 6:    z  = (y + 5);
  225.         case 5:    z *= (y + 4);
  226.         case 4:    z *= (y + 3);
  227.         case 3:    z *= (y + 2);
  228.             rr = log__D(z);
  229.             rr.b += a0_lo; rr.a += a0_hi;
  230.             return(((r+rr.b)+t+rr.a));
  231.         case 2: return(((r+a0_lo)+t)+a0_hi);
  232.         case 0: r -= log1p(x);
  233.         default: rr = log__D(x);
  234.             rr.a -= a0_hi; rr.b -= a0_lo;
  235.             return(((r - rr.b) + t) - rr.a);
  236.         }
  237.     } else {
  238.         p = p0+y*(p1+y*(p2+y*(p3+y*p4)));
  239.         q = q0+y*(q1+y*(q2+y*(q3+y*(q4+y*(q5+y*q6)))));
  240.         p = p*(y/q);
  241.         t = (double)(float) y;
  242.         z = y-t;
  243.         hi = (double)(float) (p+a1_hi);
  244.         lo = a1_hi - hi; lo += p; lo += a1_lo;
  245.         r = lo*y + z*hi;    /* q + r = y*(a0+p/q) */
  246.         q = hi*t;
  247.         z = 1.0;
  248.         switch (x_int) {
  249.         case 6:    z  = (y + 5);
  250.         case 5:    z *= (y + 4);
  251.         case 4:    z *= (y + 3);
  252.         case 3:    z *= (y + 2);
  253.             rr = log__D(z);
  254.             r += rr.b; r += q;
  255.             return(rr.a + r);
  256.         case 2:    return (q+ r);
  257.         case 0: rr = log__D(x);
  258.             r -= rr.b; r -= log1p(x);
  259.             r += q; r-= rr.a;
  260.             return(r);
  261.         default: rr = log__D(x);
  262.             r -= rr.b;
  263.             q -= rr.a;
  264.             return (r+q);
  265.         }
  266.     }
  267. }
  268.  
  269. static double
  270. neg_lgam(double x)
  271. {
  272.     int xi;
  273.     double y, z, one = 1.0, zero = 0.0;
  274.     extern double gamma();
  275.  
  276.     /* avoid destructive cancellation as much as possible */
  277.     if (x > -170) {
  278.         xi = x;
  279.         if (xi == x)
  280.             if (_IEEE)
  281.                 return(one/zero);
  282.             else
  283.                 return(infnan(ERANGE));
  284.         y = gamma(x);
  285.         if (y < 0)
  286.             y = -y, signgam = -1;
  287.         return (log(y));
  288.     }
  289.     z = floor(x + .5);
  290.     if (z == x) {        /* convention: G(-(integer)) -> +Inf */
  291.         if (_IEEE)
  292.             return (one/zero);
  293.         else
  294.             return (infnan(ERANGE));
  295.     }
  296.     y = .5*ceil(x);
  297.     if (y == ceil(y))
  298.         signgam = -1;
  299.     x = -x;
  300.     z = fabs(x + z);    /* 0 < z <= .5 */
  301.     if (z < .25)
  302.         z = sin(M_PI*z);
  303.     else
  304.         z = cos(M_PI*(0.5-z));
  305.     z = log(M_PI/(z*x));
  306.     y = large_lgam(x);
  307.     return (z - y);
  308. }
  309.