V_HYPERGEOM1F1 Confluent hypergeometric function, 1F1 a.k.a Kummer's M function [h,l]=(a,b,z,tol,maxj) Inputs: a,b,z input arguments a and b must be real scalars, z can be a real matrix The function is undefined if b is a non-positive integer tol Optional tolerance [default 1e-10] maxj Optional iteration limit [default 500] th Threshold for choice of algorithm [default 30] Outputs: h output result (size of z): 1F1(a;b;z) = M(a;b;z) l actual number of iterations (size of z) alg algorithm used This routine is functionally equivalent to the symbolic toolbox function hypergeom(a,b,z) but much faster. The function calculated is the solution to z M'' + (b-z) M' - a M = 0 where M' and M'' are 1st and 2nd derivatives of M with respect to z. This routine is closely based on taylorb1f1.m from [4] which is explained in Sec 3.2 of [3]. Special cases and relations [2]: M(a;b;z) = Inf for integer b<=0 M(a;b;0) = 1 M(0;b;z) = 1 M(a;a;z) = exp(z) M(1;2;2z) = exp(z)*sinh(z)/z M(a;b;z) = exp(z)*M(b-a;b;-z) (13.2.39 from [2]) M(a;a+1;-z) = exp(z)*M(1;a+1;z) = a*gamma(z,a)/z^a M(a;b;z) = M(a-1;b;z) + M(a;b+1;z)*z/b (b-a)M(a-1,b,z)+(2a-b+z)M(a,b,z)-aM(a+1,b,z)=0 [1] 13.4.1 b(b-1)M(a,b-1,z)+b(1-b-z)M(a,b,z)+z(b-a)(a,b+1,z)=0 [1] 13.4.2 References: [1] M. Abramowitz and I. A. Stegun, editors. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. Dover, New York, 9 edition, 1964. ISBN 0-486-61272-4. [2] F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark, editors. NIST Handbook of Mathematical Functions. CUP, 2010. ISBN 978-0-521-14063-8. [3] J. Pearson. Computation of hypergeometric functions. Master's thesis, Oxford University, 2009. URL http://people.maths.ox.ac.uk/porterm/research/-pearson_final.pdf. [4] J. Pearson. Computation of hypergeometric functions. Matlab code, Oxford University, 2009. URL http://people.maths.ox.ac.uk/porterm/research/-hypergeometricpackage.zip.
0001 function [h,l,alg]=v_hypergeom1f1(a,b,z,tol,maxj,th) 0002 %V_HYPERGEOM1F1 Confluent hypergeometric function, 1F1 a.k.a Kummer's M function [h,l]=(a,b,z,tol,maxj) 0003 % 0004 % Inputs: a,b,z input arguments 0005 % a and b must be real scalars, z can be a real matrix 0006 % The function is undefined if b is a non-positive integer 0007 % tol Optional tolerance [default 1e-10] 0008 % maxj Optional iteration limit [default 500] 0009 % th Threshold for choice of algorithm [default 30] 0010 % 0011 % Outputs: h output result (size of z): 1F1(a;b;z) = M(a;b;z) 0012 % l actual number of iterations (size of z) 0013 % alg algorithm used 0014 % 0015 % This routine is functionally equivalent to the symbolic toolbox function hypergeom(a,b,z) but much faster. 0016 % The function calculated is the solution to z M'' + (b-z) M' - a M = 0 where 0017 % M' and M'' are 1st and 2nd derivatives of M with respect to z. 0018 % This routine is closely based on taylorb1f1.m from [4] which is explained in Sec 3.2 of [3]. 0019 % 0020 % Special cases and relations [2]: 0021 % M(a;b;z) = Inf for integer b<=0 0022 % M(a;b;0) = 1 0023 % M(0;b;z) = 1 0024 % M(a;a;z) = exp(z) 0025 % M(1;2;2z) = exp(z)*sinh(z)/z 0026 % M(a;b;z) = exp(z)*M(b-a;b;-z) (13.2.39 from [2]) 0027 % M(a;a+1;-z) = exp(z)*M(1;a+1;z) = a*gamma(z,a)/z^a 0028 % M(a;b;z) = M(a-1;b;z) + M(a;b+1;z)*z/b 0029 % (b-a)M(a-1,b,z)+(2a-b+z)M(a,b,z)-aM(a+1,b,z)=0 [1] 13.4.1 0030 % b(b-1)M(a,b-1,z)+b(1-b-z)M(a,b,z)+z(b-a)(a,b+1,z)=0 [1] 13.4.2 0031 % 0032 % References: 0033 % [1] M. Abramowitz and I. A. Stegun, editors. 0034 % Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. 0035 % Dover, New York, 9 edition, 1964. ISBN 0-486-61272-4. 0036 % [2] F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark, editors. 0037 % NIST Handbook of Mathematical Functions. CUP, 2010. ISBN 978-0-521-14063-8. 0038 % [3] J. Pearson. Computation of hypergeometric functions. Master's thesis, Oxford University, 2009. 0039 % URL http://people.maths.ox.ac.uk/porterm/research/-pearson_final.pdf. 0040 % [4] J. Pearson. Computation of hypergeometric functions. Matlab code, Oxford University, 2009. 0041 % URL http://people.maths.ox.ac.uk/porterm/research/-hypergeometricpackage.zip. 0042 % 0043 0044 % Copyright (C) Mike Brookes 2016 0045 % Version: $Id: v_hypergeom1f1.m 10865 2018-09-21 17:22:45Z dmb $ 0046 % 0047 % VOICEBOX is a MATLAB toolbox for speech processing. 0048 % Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html 0049 % 0050 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0051 % This program is free software; you can redistribute it and/or modify 0052 % it under the terms of the GNU General Public License as published by 0053 % the Free Software Foundation; either version 2 of the License, or 0054 % (at your option) any later version. 0055 % 0056 % This program is distributed in the hope that it will be useful, 0057 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0058 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0059 % GNU General Public License for more details. 0060 % 0061 % You can obtain a copy of the GNU General Public License from 0062 % http://www.gnu.org/copyleft/gpl.html or by writing to 0063 % Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. 0064 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0065 if nargin<4 || isempty(tol) 0066 tol=1e-10; 0067 end 0068 if nargin<5 || isempty(maxj) 0069 maxj=500; 0070 end 0071 if nargin<6 || isempty(th) 0072 th=30; 0073 end 0074 a1=a-1; % define some useful constants 0075 b1=b-1; 0076 ba=b-a; 0077 h=zeros(size(z)); 0078 l=zeros(size(z)); 0079 nz=numel(z); 0080 for i=1:nz 0081 y=z(i); 0082 q=0; % break criterion initially false 0083 if abs(y)<th % small |y| 0084 % compute using the series from 13.2.2 in [2] 0085 % sum(j=0:inf,prod(a+(0:j-1)) y^j /(prod(b+(0:j-1)) j!)) 0086 alg=1; 0087 d=1; 0088 g=1; 0089 jlim=0; % find j beyond which terms definitely decrease 0090 k=(b1+y)^2-4*(a1)*y; 0091 if k>=0 0092 jlim=max(jlim,0.5*(sqrt(k)-(b1+y))); 0093 end 0094 k=(b1-y)^2+4*(a1)*y; 0095 if k>=0 0096 jlim=max(jlim,0.5*(sqrt(k)-(b1-y))); 0097 end 0098 jlim=min(maxj,jlim); 0099 for j=1:maxj 0100 d=d*y*(a1+j)/(j*(b1+j)); % d = A(j)-A(j-1) = (A(j-1)-A(j-2))*r(j)*z from [4] 0101 g=g+d; % g = A(j) from [4] 0102 p=abs(d)<tol*abs(g); 0103 if q && p && j>=jlim 0104 break 0105 end 0106 q=p; 0107 end 0108 elseif y>0 % big positive y 0109 % see 13.7.1 combined with 13.2.3 [2]. Valid for large y unless a is a non-positive integer 0110 alg=2; 0111 d=1; 0112 g=1; 0113 jlim=1; % find j beyond which terms definitely increase 0114 k=(ba-1-a+y)^2+4*a*(ba-1); 0115 if k>=0 0116 jlim=max(jlim,0.5*(sqrt(k)-(ba-a-1+y))); 0117 end 0118 k=(ba-a-1-y)^2+4*a*(ba-1); 0119 if k>=0 0120 jlim=max(jlim,0.5*(sqrt(k)-(ba-a-1-y))); 0121 end 0122 jlim=min(maxj,jlim); 0123 for j=1:jlim 0124 d=d*(ba-1+j)*(j-a)/(j*y); 0125 g=g+d; 0126 p=abs(d)<tol*abs(g); 0127 if q && p 0128 break 0129 end 0130 q=p; 0131 end 0132 [gl,gs]=v_gammalns([a b]); 0133 g=gs(1)*gs(2)*exp(y+gl(2)-gl(1)+(a-b)*log(y))*g; 0134 else % big negative y 0135 % Combine M(a;b;z) = exp(z)*M(b-a;b;-z) (13.2.39 from [2]) with algorithm 2 0136 alg=3; 0137 d=1; 0138 g=1; 0139 jlim=1; % find j beyond which terms definitely increase 0140 k=(a1-ba+y)^2+4*a1*ba; 0141 if k>=0 0142 jlim=max(jlim,0.5*(sqrt(k)-(a1-ba+y))); 0143 end 0144 k=(a1-ba-y)^2+4*a1*ba; 0145 if k>=0 0146 jlim=max(jlim,0.5*(sqrt(k)-(a1-ba-y))); 0147 end 0148 jlim=min(maxj,jlim); 0149 for j=1:jlim 0150 d=d*(a1+j)*(ba-j)/(j*y); 0151 g=g+d; 0152 p=abs(d)<tol*abs(g); 0153 if q && p 0154 break 0155 end 0156 q=p; 0157 end 0158 [gl,gs]=v_gammalns([ba b]); 0159 g=gs(1)*gs(2)*exp(gl(2)-gl(1)-a*log(-y))*g; 0160 end 0161 h(i)=g; 0162 l(i)=j; 0163 end