# v_hypergeom1f1

## PURPOSE

V_HYPERGEOM1F1 Confluent hypergeometric function, 1F1 a.k.a Kummer's M function [h,l]=(a,b,z,tol,maxj)

## SYNOPSIS

function [h,l,alg]=v_hypergeom1f1(a,b,z,tol,maxj,th)

## DESCRIPTION

    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.

## CROSS-REFERENCE INFORMATION

This function calls:

- v_gammalns V_GAMMALNS Log of Gamma(x) for positive or negative real x [y,s]=(x)

This function is called by:

## SOURCE CODE

    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 Lesser General Public License as published by
    0053 %   the Free Software Foundation; either version 3 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 Lesser General Public License for more details.
    0060 %
    0061 %   You can obtain a copy of the GNU Lesser General Public License from
    0062 %   https://www.gnu.org/licenses/ .
    0063 %
    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