home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 2.ddi / SPECFUN.DI$ / GAMMA.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  3.5 KB  |  125 lines
  1. function res = gamma(x,a)
  2. %GAMMA    The gamma function.
  3. %    Y = GAMMA(X) evaluates the gamma function at all the elements of X.
  4. %    X must be real.
  5. %    gamma(x) = integral from 0 to inf of t^(x-1) exp(-t) dt.
  6. %    gamma(n+1) = n! = n factorial = prod(1:n).
  7. %
  8. %    See also GAMMALN, GAMMAINC.
  9.  
  10. %    C. B. Moler, 5-7-91, 11-4-92.
  11. %    Ref: Abramowitz & Stegun, Handbook of Mathemtical Functions, sec. 6.1.
  12. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  13.  
  14. %   This is based on a FORTRAN program by W. J. Cody,
  15. %   Argonne National Laboratory, NETLIB/SPECFUN, October 12, 1989.
  16. %
  17. % References: "An Overview of Software Development for Special
  18. %              Functions", W. J. Cody, Lecture Notes in Mathematics,
  19. %              506, Numerical Analysis Dundee, 1975, G. A. Watson
  20. %              (ed.), Springer Verlag, Berlin, 1976.
  21. %
  22. %              Computer Approximations, Hart, Et. Al., Wiley and
  23. %              sons, New York, 1968.
  24. %
  25. if nargin == 2
  26.     % For backward compatibility with MATLAB 3.5
  27.     res = gammainc(a,x);
  28.     return
  29. end
  30.  
  31.       p = [-1.71618513886549492533811e+0; 2.47656508055759199108314e+1,
  32.            -3.79804256470945635097577e+2; 6.29331155312818442661052e+2,
  33.             8.66966202790413211295064e+2; -3.14512729688483675254357e+4,
  34.            -3.61444134186911729807069e+4; 6.64561438202405440627855e+4];
  35.       q = [-3.08402300119738975254353e+1; 3.15350626979604161529144e+2,
  36.            -1.01515636749021914166146e+3; -3.10777167157231109440444e+3,
  37.             2.25381184209801510330112e+4; 4.75584627752788110767815e+3,
  38.            -1.34659959864969306392456e+5; -1.15132259675553483497211e+5];
  39.       c = [-1.910444077728e-03; 8.4171387781295e-04,
  40.            -5.952379913043012e-04; 7.93650793500350248e-04,
  41.            -2.777777777777681622553e-03; 8.333333333333333331554247e-02,
  42.             5.7083835261e-03];
  43.  
  44.    if any(imag(x) ~= 0)
  45.       error('Input argument must be real.')
  46.    end
  47.    res = zeros(size(x));
  48.    xn = zeros(size(x));
  49. %
  50. %  Catch negative x.
  51. %
  52.    kneg = find(x <= 0);
  53.    if any(kneg)
  54.       y = -x(kneg);
  55.       y1 = fix(y);
  56.       res(kneg) = y - y1;
  57.       fact = -pi ./ sin(pi*res(kneg)) .* (1 - 2*rem(y1,2));
  58.       x(kneg) = y + 1;
  59.    end
  60. %
  61. %  x is now positive.
  62. %  Map x in interval [0,1] to [1,2]
  63. %
  64.    k1 = find(x < 1);
  65.    x1 = x(k1);
  66.    x(k1) = x1 + 1;
  67. %
  68. %  Map x in interval [1,12] to [1,2]
  69. %
  70.    k = find(x < 12);
  71.    xn(k) = fix(x(k)) - 1;
  72.    x(k) = x(k) - xn(k);
  73. %
  74. %  Evaluate approximation for 1 < x < 2
  75. %
  76.    if any(k)
  77.       z = x(k) - 1;
  78.       xnum = 0*z;
  79.       xden = xnum + 1;
  80.       for i = 1:8
  81.          xnum = (xnum + p(i)) .* z;
  82.          xden = xden .* z + q(i);
  83.       end
  84.       res(k) = xnum ./ xden + 1;
  85.    end
  86. %
  87. %  Adjust result for case  0.0 < x < 1.0
  88. %
  89.    res(k1) = res(k1) ./ x1;
  90. %
  91. %  Adjust result for case  2.0 < x < 12.0
  92. %
  93.    for j = 1:max(max(xn))
  94.       k = find(xn);
  95.       res(k) = res(k) .* x(k);
  96.       x(k) = x(k) + 1;
  97.       xn(k) = xn(k) - 1;
  98.    end
  99. %
  100. %  Evaluate approximation for x >= 12
  101. %
  102.    k = find(x >= 12);
  103.    if any(k)
  104.       y = x(k);
  105.       ysq = y .* y;
  106.       sum = c(7);
  107.       for i = 1:6
  108.          sum = sum ./ ysq + c(i);
  109.       end
  110.       spi = 0.9189385332046727417803297;
  111.       sum = sum ./ y - y + spi;
  112.       sum = sum + (y-0.5).*log(y);
  113.       res(k) = exp(sum);
  114.    end
  115. %
  116. %  Final adjustments.
  117. %
  118.    if any(~finite(x))
  119.       k = find(isinf(x)); res(k) = Inf*ones(size(k));
  120.       k = find(isnan(x)); res(k) = NaN*ones(size(k));
  121.    end
  122.    if any(kneg)
  123.       res(kneg) = fact ./ res(kneg);
  124.    end
  125.