Mercurial > hg > octave-lyh

comparison scripts/specfun/gammai.m @ 1586:54350c98055c

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1995-10-31 05:22:53 by jwe]
author jwe
date Tue, 31 Oct 1995 05:24:40 +0000
parents 100413a7e8a2
children 5d29638dd524
comparison
equal deleted inserted replaced
1585:100413a7e8a2 1586:54350c98055c
78 endif 78 endif
79 if (any (x < 0)) 79 if (any (x < 0))
80 error ("gammai: all entries of x must be nonnegative"); 80 error ("gammai: all entries of x must be nonnegative");
81 endif 81 endif
82 82
83 y = zeros(1, n); 83 y = zeros (1, n);
84 84
85 # For x < a + 1, use summation. The below choice of k should ensure 85 # For x < a + 1, use summation. The below choice of k should ensure
86 # that the overall error is less than eps ... 86 # that the overall error is less than eps ...
87 87
88 S = find ((x > 0) & (x < a + 1)); 88 S = find ((x > 0) & (x < a + 1));
89 s = length (S); 89 s = length (S);
90 if (s > 0) 90 if (s > 0)
91 k = ceil (- max ([a(S), x(S)]) * log (eps)); 91 k = ceil (- max ([a(S), x(S)]) * log (eps));
92 K = (1:k)'; 92 K = (1:k)';
93 M = ones(k, 1); 93 M = ones (k, 1);
94 A = cumprod((M * x(S)) ./ (M * a(S) + K * ones(1, s))); 94 A = cumprod ((M * x(S)) ./ (M * a(S) + K * ones(1, s)));
95 y(S) = exp (-x(S) + a(S) .* log (x(S))) .* (1 + sum (A)) ./ gamma (a(S)+1); 95 y(S) = exp (-x(S) + a(S) .* log (x(S))) .* (1 + sum (A)) ./ gamma (a(S)+1);
96 endif 96 endif
97 97
98 # For x >= a + 1, use the continued fraction. 98 # For x >= a + 1, use the continued fraction.
99 # Note, however, that this converges MUCH slower than the series 99 # Note, however, that this converges MUCH slower than the series
100 # expansion for small a and x not too large! 100 # expansion for small a and x not too large!
101 101
102 S = find((x >= a + 1) & (x < Inf)); 102 S = find ((x >= a + 1) & (x < Inf));
103 s = length(S); 103 s = length (S);
104 if (s > 0) 104 if (s > 0)
105 u = [zeros(1, s); ones(1, s)]; 105 t1 = zeros (1, s);
106 v = [ones(1, s); x(S)]; 106 t2 = ones (1, s);
107 u = [t1; t2];
108 v = [t2; x(S)];
107 c_old = 0; 109 c_old = 0;
108 c_new = v(1,:) ./ v(2,:); 110 c_new = v(1,:) ./ v(2,:);
109 n = 1; 111 n = 1;
110 while (max (abs (c_old ./ c_new - 1)) > 10 * eps) 112 while (max (abs (c_old ./ c_new - 1)) > 10 * eps)
111 c_old = c_new; 113 c_old = c_new;