home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 2.ddi / SPECFUN.DI$ / GAMMALN.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  4.8 KB  |  154 lines
  1. function res = gammaln(x)
  2. %GAMMALN The loggamma function.
  3. %    Y = GAMMALN(X) is LOG(GAMMA(X)), avoiding underflow and overflow.
  4. %    See also GAMMA.
  5.  
  6. %    C. Moler, 2-1-91.
  7. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  8.  
  9. %   This is based on a FORTRAN program by W. J. Cody,
  10. %   Argonne National Laboratory, NETLIB/SPECFUN, June 16, 1988.
  11. %
  12. % References:
  13. %
  14. %  1) W. J. Cody and K. E. Hillstrom, 'Chebyshev Approximations for
  15. %     the Natural Logarithm of the Gamma Function,' Math. Comp. 21,
  16. %     1967, pp. 198-203.
  17. %
  18. %  2) K. E. Hillstrom, ANL/AMD Program ANLC366S, DGAMMA/DLGAMA, May,
  19. %     1969.
  20. % 
  21. %  3) Hart, Et. Al., Computer Approximations, Wiley and sons, New
  22. %     York, 1968.
  23. %
  24.  
  25.      d1 = -5.772156649015328605195174e-1;
  26.      p1 = [4.945235359296727046734888e0; 2.018112620856775083915565e2; 
  27.            2.290838373831346393026739e3; 1.131967205903380828685045e4; 
  28.            2.855724635671635335736389e4; 3.848496228443793359990269e4; 
  29.            2.637748787624195437963534e4; 7.225813979700288197698961e3];
  30.      q1 = [6.748212550303777196073036e1; 1.113332393857199323513008e3; 
  31.            7.738757056935398733233834e3; 2.763987074403340708898585e4; 
  32.            5.499310206226157329794414e4; 6.161122180066002127833352e4; 
  33.            3.635127591501940507276287e4; 8.785536302431013170870835e3];
  34.      d2 = 4.227843350984671393993777e-1;
  35.      p2 = [4.974607845568932035012064e0; 5.424138599891070494101986e2; 
  36.            1.550693864978364947665077e4; 1.847932904445632425417223e5; 
  37.            1.088204769468828767498470e6; 3.338152967987029735917223e6; 
  38.            5.106661678927352456275255e6; 3.074109054850539556250927e6];
  39.      q2 = [1.830328399370592604055942e2; 7.765049321445005871323047e3; 
  40.            1.331903827966074194402448e5; 1.136705821321969608938755e6; 
  41.            5.267964117437946917577538e6; 1.346701454311101692290052e7; 
  42.            1.782736530353274213975932e7; 9.533095591844353613395747e6];
  43.      d4 = 1.791759469228055000094023e0;
  44.      p4 = [1.474502166059939948905062e4; 2.426813369486704502836312e6; 
  45.            1.214755574045093227939592e8; 2.663432449630976949898078e9; 
  46.            2.940378956634553899906876e10; 1.702665737765398868392998e11; 
  47.            4.926125793377430887588120e11; 5.606251856223951465078242e11];
  48.      q4 = [2.690530175870899333379843e3; 6.393885654300092398984238e5; 
  49.            4.135599930241388052042842e7; 1.120872109616147941376570e9; 
  50.            1.488613728678813811542398e10; 1.016803586272438228077304e11; 
  51.            3.417476345507377132798597e11; 4.463158187419713286462081e11];
  52.      c = [-1.910444077728e-03; 8.4171387781295e-04; 
  53.           -5.952379913043012e-04; 7.93650793500350248e-04; 
  54.           -2.777777777777681622553e-03; 8.333333333333333331554247e-02; 
  55.            5.7083835261e-03];
  56.  
  57.    if any((imag(x) ~= 0) | (x < 0))
  58.       error('Input arguments must be real and nonnegative.')
  59.    end
  60.    res = x;
  61. %
  62. %  0 <= x <= eps
  63. %
  64.    k = find(x <= eps);
  65.    if any(k)
  66.       res(k) = -log(x(k));
  67.    end
  68. %
  69. %  eps < x <= 0.5
  70. %
  71.    k = find((x > eps) & (x <= 0.5));
  72.    if any(k)
  73.       y = x(k);
  74.       xden = ones(size(k));
  75.       xnum = 0;
  76.       for i = 1:8
  77.          xnum = xnum .* y + p1(i);
  78.          xden = xden .* y + q1(i);
  79.       end
  80.       res(k) = -log(y) + (y .* (d1 + y .* (xnum ./ xden)));
  81.    end
  82. %
  83. %  0.5 < x <= 0.6796875
  84. %
  85.    k = find((x > 0.5) & (x <= 0.6796875));
  86.    if any(k)
  87.       xm1 = (x(k) - 0.5) - 0.5;
  88.       xden = ones(size(k));
  89.       xnum = 0;
  90.       for i = 1:8
  91.          xnum = xnum .* xm1 + p2(i);
  92.          xden = xden .* xm1 + q2(i);
  93.       end
  94.       res(k) = -log(x(k)) + xm1 .* (d2 + xm1 .* (xnum ./ xden));
  95.    end
  96. %
  97. %  0.6796875 < x <= 1.5
  98. %
  99.    k = find((x > 0.6796875) & (x <= 1.5));
  100.    if any(k)
  101.       xm1 = (x(k) - 0.5) - 0.5;
  102.       xden = ones(size(k));
  103.       xnum = 0;
  104.       for i = 1:8
  105.          xnum = xnum .* xm1 + p1(i);
  106.          xden = xden .* xm1 + q1(i);
  107.       end
  108.       res(k) = xm1 .* (d1 + xm1 .* (xnum ./ xden));
  109.    end
  110. %
  111. %  1.5 < x <= 4
  112. %
  113.    k = find((x > 1.5) & (x <= 4));
  114.    if any(k)
  115.       xm2 = x(k) - 2;
  116.       xden = ones(size(k));
  117.       xnum = 0;
  118.       for i = 1:8
  119.          xnum = xnum .* xm2 + p2(i);
  120.          xden = xden .* xm2 + q2(i);
  121.       end
  122.       res(k) = xm2 .* (d2 + xm2 .* (xnum ./ xden));
  123.    end
  124. %
  125. %  4 < x <= 12
  126. %
  127.    k = find((x > 4) & (x <= 12));
  128.    if any(k)
  129.       xm4 = x(k) - 4;
  130.       xden = -ones(size(k));
  131.       xnum = 0;
  132.       for i = 1:8
  133.          xnum = xnum .* xm4 + p4(i);
  134.          xden = xden .* xm4 + q4(i);
  135.       end
  136.       res(k) = d4 + xm4 .* (xnum ./ xden);
  137.    end
  138. %
  139. %  x > 12
  140. %
  141.    k = find(x > 12);
  142.    if any(k)
  143.       y = x(k);
  144.       r = c(7)*ones(size(k));
  145.       ysq = y .* y;
  146.       for i = 1:6
  147.          r = r ./ ysq + c(i);
  148.       end
  149.       r = r ./ y;
  150.       corr = log(y);
  151.       spi = 0.9189385332046727417803297;
  152.       res(k) = r + spi - 0.5*corr + y .* (corr-1);
  153.    end
  154.