Bessels and Gammas (new)

[Eyal Doron]

Hi all,

Here is an update of the Bessel function/Gamma function package. It fixes
the following:

*) jybess, cgamma and cgammaln failed on vector inputs.
*) cgamma and cgammaln failed on complex arguments with a zero imaginary
   part.
*) Removed dependence on "implicit_str_to_num_ok" and "do_fortran_indexing"
   settings in Octave.
*) jybess returned wrong results on complex orders with negative real parts.
*) Worked around an Octave indexing bug (I hope).

Files in this release are:

1) jybess.m   - Bessel and Newmann functions for real or complex orders
                and arguments.
2) cgamma.m   - Gamma function for complex arguments. Required for jybess 
                if complex orders are evaluated.
3) cgammaln.m - log-gamma function for complex arguments.
4) cdigamma.m - Digamma (psi) function for complex arguments.

Again, please regard these files as beta, and let me know any comments.

Good luck,
      Eyal Doron


-----------------------------jybess.m---------------------------
function [Jn,Yn] = jybess(n,x,ToDo);
% FUNCTION [Jn,Yn] = jybess(n,x,ToDo):
% Returns the Bessel and Newmann functions for complex arguments and order.
%
% n    - Maximal order, real or complex scalar.
% x    - Real or complex argument vector.
% ToDo - Optional string of command characters (case insensitive):
%        'S' - return both negative and positive orders (see below).
%        'Y' - if nargout<2, return Y_n(x) instead of J_n(x).
%        'L' - return only the highest computed order.
% -----------------------------------------
% Jn,Yn - Bessel functions of the first and (optionally) second kind.
%
% Jn, Yn are length(x) X (# orders) matrices. The orders returned are:
%   -- i*ni + [frac(nr):sign(nr):nr], or
%   -- i*ni + [-|nr|:-frac(|nr|), frac(|nr|):|nr|] if 'S'.
% where nr=real(n) and ni=imag(n).
%
% Note: Support for complex n is only partial; The algorithm doesn't work very
% well for large imag(n) and |x/n|>~1 .

% Algorithm:
% ----------
% 1) Calculate J_n by backwards recursion from a calculated starting point.
%    Minimal order is \nu=|n|-floor(|n|).
% 2) Normalization: Normalize J_\nu(x) by:
%    a) for integer n and |\Im(x)|<log(10), use
%                         1 + 2 sum_{k=1}^\infty J_{2k}(x) = 1 
%    b) for integer n and |\Im(x)|>log(10), use
%                         1 + 2 sum_{k=1}^\infty (-)^k J_{2k}(x) = cos(x)
%    c) for half-integer n, use explicit formula for the spherical Bessels.
%    d) for complex n and |\Im(x)|<5, use
%       sum_{k=0}^\infty (\nu+2k)\Gamma(\nu+k)/k! J_{\nu+2k}(x) = (x/2)^\nu
%    e) for complex n and |\Im(x)|>5, use
%       e^{i\phi\nu} \sum_{k=0}^\infty [r/2*(1-e^{2i\phi})]^k/k! J_{\nu+k}(r)
%                   = J_\nu(r e^i\phi)
%       This requires a recursive call to obtain the J_{\nu+k}(r).
% 3) If Y_n, or J_n for negative non-integer n, are required, calculate
%    the Y_n(x) as follows:
% 4) Evaluate Y_\nu:
%    a) for nu=0, use
%      2/\pi(\log(x/2)+\gamma)J_0(x) - 4/\pi\sum_{k=1}^\infty (-)^k J_{2k}(x)/k
%                  = Y_0(x)
%    b) for nu=0.5, use explicit formula for the spherical Bessels.
%    c) otherwise, do the following:
%        i) Calculate J_{1-\nu}(x) and J_{2-\nu}(x) using a recursive call.
%       ii) Calculate J_{-\nu}(x) using one backwards recursion step.
%      iii) Use [J_\nu(x)\cos(\nu\pi)-J_{-\nu}(x)]/\sin(\nu\pi) = Y_\nu(x)
% 5) For real arguments, obtain the rest of the Y_{\nu+k}(x) using forward
%    recurrence, Y_{\nu+1}(x)={2\nu\over x}Y_\nu(x) - Y_{\nu-1}(x) .
% 6) Obtain the rest of the Y_{\nu+k}(x) using the Wronskian relation
%    J_{\nu+1}(x)Y_\nu(x) - J_\nu(x)Y_{\nu+1}(x) = 2/(\pi x)
%    This turns out to be much more stable for complex arguments than forward
%    recurrence technique, but fails at zeros of the J_\nu(x) .
% 7) If J for negative orders is required:
%    a) for nu=0, use J_n(x) = (-)^n J_{-n}(x)
%    b) otherwise, use 4.c.iii
% 8) If Y for negative orders is required:
%    a) for nu=0, use Y_n(x) = (-)^n Y_{-n}(x)
%    b) otherwise, use 4.c.iii
%
% For complex orders, jybess also requires the complex Gamma function cgamma.

if nargin==0
  help jybess
  return
end
if nargin<2
  error('Not enough input parameters!');
end
if max(max(size(n)))>1
  error('Max order must be a scalar!')
end
sym=0; Return_Y=(nargout==2); Return_Last=0;
if nargin>2
  if exist('OCTAVE_VERSION')
    ToDo=toascii(ToDo);
    sym=any(ToDo==1) | any(ToDo==toascii('s')) | any(ToDo==toascii('S'));
    Return_Y=Return_Y | any(ToDo==toascii('y') | ToDo==toascii('Y'));
    Return_Last=any(ToDo==toascii('l') | ToDo==toascii('L'));
  else
    sym=any(ToDo==1) | any(ToDo=='s') | any(ToDo=='S');
    Return_Y=Return_Y | any(ToDo=='y' | ToDo=='Y');
    Return_Last=any(ToDo=='l' | ToDo=='L');
  end
end
sym=(sym~=0);
if sym & real(n)<0, n=-real(n)+i*imag(n); end
orig_n=n;
n_is_negative=(real(n)<0);
if n_is_negative, n=-real(n)+i*imag(n); end
nu=n-floor(real(n)); n=floor(real(n));
calc_Y=(Return_Y | sym | n_is_negative);
Acc_Jnorm=(nu~=0.5);
Acc_Ynorm=calc_Y & nu==0;

xlen=prod(size(x)); x=reshape(x,1,xlen);
z = x==0; x = x + z;            % Temporarily replace x=0 with x=1

tiny = 16^(-250);

%                             Starting index for backwards recurrence.
c = [ 0.9507,    1.4208,   14.1850;
      0.9507,    1.4208,   14.1850;
      0.9507,    1.4208,   14.1850;
      0.7629,    1.4222,   13.9554;
      0.7369,    1.4289,   13.1756;
      0.7674,    1.4311,   12.4523;
      0.8216,    1.4354,   11.2121;
      0.8624,    1.4397,    9.9718;
      0.8798,    1.4410,    8.9217;
      0.9129,    1.4360,    7.8214;
      0.9438,    1.5387,    6.5014;
      0.9609,    1.5216,    5.7256;
      0.9693,    1.5377,    5.3565;
      0.9823,    1.5220,    4.5659;
      0.9934,    1.5049,    3.7902;
      0.9985,    1.4831,    3.2100;
      1.0006,    1.4474,    3.0239;
      0.9989,    1.4137,    2.8604;
      0.9959,    1.3777,    2.7760;
      1.0005,    1.3500,    2.3099]';
j = 1+min(abs(n),19);
m = c(1,j).*max(3,j) + c(2,j).*(max(1,abs(x))-1) + ...
    c(3,j)./(1-log(min(1,abs(x))));
m=max([m; abs(n)+10+0*m]);
m = 4*ceil(m/4);                      % Make sure rem(m,4)=0.

%                                       Prevent underflow
logtiny=log(tiny); logx=log(abs(x)/2)+1;
II=1:xlen;
while ~isempty(II)
  est=-log(2*pi*m(II))/2+m(II).*(logx(II)-log(abs(m(II)+nu)))<logtiny;
  II=II(find(est)); m(II)=m(II)-4;
end;
mm = max(max(m));

% Normalization summation coefficients

if nu==0
  Nrm=[1; 2*ones(ceil(mm/2),1)]; 
  UseCos=find(abs(imag(x))>log(10));
  if ~isempty(UseCos)                  % Use norm. to cos(z) instead of 1.
    Nrm=Nrm*ones(1,xlen);
    [r,c]=size(Nrm);
    Nrm(2:2:r,UseCos)=-Nrm(2:2:r,UseCos);
  end
elseif nu~=0.5               % nu=0.5 doesn't use normalizations
  k=(1:ceil(mm/2)).';
  if imag(nu)==0
    Nrm=cumprod([nu*gamma(nu); (nu+2*k).*(nu+k-1)./k./(nu+2*k-2)]);
  else                                    %     Needs complex gamma function
    Nrm=cumprod([nu*cgamma(nu); (nu+2*k).*(nu+k-1)./k./(nu+2*k-2)]);
  end
  Nrm=Nrm*ones(1,xlen);
end

%                           Backwards recursion for the Jn
Jn=zeros(n+1,xlen);
k = mm;
bkp1 = 0*x;
bk = tiny*(m==k);
if Acc_Jnorm, t = Nrm(k/2+1,:).*bk; end
if Acc_Ynorm, y = 2*bk/k; end
sig = 1;
for k = mm-1:-1:0
  bkp2 = bkp1;
  bkp1 = bk;
  bk = 2*(k+1+nu)*bkp1./x - bkp2 + tiny*(m==k);
  if k<=n
    Jn(k+1,:)=bk;
  end
  if (floor(k/2)*2-k==0)
     if Acc_Jnorm, t = t + Nrm(k/2+1,:).*bk; end
     if Acc_Ynorm & k>0
       sig = -sig;
       y = y + sig*2*bk/k;
     end
  end
end
clear Nrm bkp2 sig

if nu==0.5               % Use explicit formulae for the spherical Bessels
  si=sin(x);
  II=(abs(si)>0.1);
  F1=find(II); F2=find(1-II);
  nrm=ones(1,xlen);
  SpFact=sqrt(2*x/pi);
  if ~isempty(F1)
    nrm(F1)=SpFact(F1).*(si(F1)./x(F1))./Jn(1,F1);
  end
  if ~isempty(F2)
    nrm(F2)=SpFact(F2).*(si(F2)./x(F2).^2-cos(x(F2))./x(F2))./bkp1(F2);
  end
elseif nu==0
  nrm=ones(size(x));
  if ~isempty(UseCos)
    nrm(UseCos)=cos(x(UseCos));
  end
  abst=abs(t); nrm=nrm./abst./(t./abst);
else
  II=(abs(imag(x))<5);
  F1=find(II); F2=find(1-II);
  nrm=ones(1,xlen);
  if ~isempty(F1)
    nrm(F1)=(x(F1)/2).^nu ./t(F1);
  end
  if ~isempty(F2)               % Use mult. Theorem to calc J_0
    r=abs(x(F2)); ph=angle(x(F2)); rn=max(ceil(r));
    MaxOrder=ceil(18+1.25*rn);
    Jtmp=jybess(nu+MaxOrder,r);
    k=(0:MaxOrder).';
    Kg=cumsum([0;log(1:MaxOrder).']);          % log(k!)
    v=i*ph-i*pi/2+log(imag(x(F2)));
    
nrm(F2)=exp(i*nu*ph).*sum(exp(k*v-(Kg*ones(1,length(F2)))).*Jtmp.')./Jn(1,F2);
  end
end
Jn=Jn.*(ones(n+1,1)*nrm);           % Normalizing condition.

%                Restore results for x = 0; j0(0) = 1, jn(0) = 0 for n > 0.
if any(z)
  Ii=find(z); LI=max(size(Ii));
  Jn(:,Ii)=[ones(1,LI)*(nu==0);zeros(n,LI)];
end
Jn=Jn.';

if calc_Y                                    % Also Yn
%                         First, get Y_0(x)
  if n==0
    J1 = bkp1.*nrm;
  else
    J1=Jn(:,2).';
  end
  if nu==0                          % Use dependence on sums of J
    gam = 0.57721566490153286;     % Euler number
    lx=log(x/2);
    II=find(imag(x)==0 & real(x)<0);
    if ~isempty(II)
      lx(II)=log(abs(x(II))/2);
    end
    y = 2/pi*(lx + gam).*bk - 4/pi*y;
    Y0 = y.*nrm;
  elseif nu==0.5                    % explicit formula
    Y0=-cos(x)./x.*SpFact;
  else                              % Use linear relations between J and Y
    if real(nu)==0
      Jminus=jybess(-nu,x);
    else                            % One step of the recursion for the J
      Jtmp=jybess(2-nu,x);
      Jminus=2*(1-nu)./x.'.*Jtmp(:,1)-Jtmp(:,2);
    end
    Y0=(Jn(:,1).'*cos(pi*nu)-Jminus.')/sin(pi*nu);
  end
  
  Yn=[Y0.', zeros(xlen,n)];
  if n>0
    Yn(:,2)=(Jn(:,2).*Yn(:,1)-2/pi./x.')./Jn(:,1);
  end
  if n>1
    II=(imag(x)==0) * (imag(nu)==0);
    if any(II)                        % Iterate using forward recurrence
      IN=find(II);
      for l=2:n
        Yn(IN,l+1)=2*(nu+l-1)./x(IN).'.*Yn(IN,l) - Yn(IN,l-1);
      end
    end
    if ~all(II)                       % Iterate using Wronskian relation
      IN=find(~II);
      for l=2:n
        Yn(IN,l+1)=(Jn(IN,l+1).*Yn(IN,l)-2/pi./x(IN).')./Jn(IN,l);
      end
    end
  end
  
  if any(z), Yn(find(z),:)=-inf*ones(size(Yn(find(z),:))); end
end

if nu==0
  if n_is_negative                                 % negative order
    Jn(:,2:2:n+1)=-Jn(:,2:2:n+1); Jn=fliplr(Jn);
    if Return_Y
      Yn(:,2:2:n+1)=-Yn(:,2:2:n+1); Yn=fliplr(Yn);
    end
  end
  if sym & n>0                                     % symmetric output
    Rev=n:-1:1; Rev=2*(floor(Rev/2)*2==Rev)-1;
    Jn=[Jn(:,n+1:-1:2).*(ones(xlen,1)*Rev) Jn];
    if Return_Y
      Yn=[Yn(:,n+1:-1:2).*(ones(xlen,1)*Rev) Yn];
    end
  end
else
  if n_is_negative | sym
    ord=nu+(0:n); co=cos(pi*ord); so=sin(pi*ord);
    Jminus=Jn.*(ones(xlen,1)*co)-Yn.*(ones(xlen,1)*so);
    if Return_Y
      ord=nu+(0:n);co=cos(pi*ord); so=sin(-pi*ord);
      Yminus=(Jminus.*(ones(xlen,1)*co)-Jn)./(ones(xlen,1)*so);
    end
  end
  if sym
    if real(nu)==0
      Jn=[fliplr(Jminus(:,2:n+1)), Jn];
    else
      Jn=[fliplr(Jminus), Jn];
    end
    if Return_Y
      if real(nu)==0
        Yn=[fliplr(Yminus(:,2:n+1)), Yn];
      else
        Yn=[fliplr(Yminus), Yn];
      end
    end
  elseif n_is_negative
    Jn=fliplr(Jminus);
    if Return_Y, Yn=fliplr(Yminus); end
  end
end
if Return_Last
  [r,c]=size(Jn);
  if n_is_negative, c=1; end
  Jn=Jn(:,c);
  if calc_Y, Yn=Yn(:,c); end
end
if nargout<2 & Return_Y, Jn=Yn; end


---------------------------------cgamma.m-------------------------
function Y=cgamma(Z);
% function Y=cgamma(Z) returns the value of the Gamma function for complex
% argument matrix Z.

% Algorithm:
% ----------
% 1) for real Z, use the builtin gamma function.
% 2) If real(Z)<0, use the reflection formula for Z -> -Z.
% 3) For |Z|>6 use an order-15 Stirling's approximation.
% 4) For the rest, use gamma(z+1)=z*gamma(z) and evaluate gamma(z+N) such that
%    |z+N|>6 .

Asymp_coef=[
        +1.0;
        +0.83333333333333333333e-1;
        +0.34722222222222222222e-2;
        -0.26813271604938271605e-2;
        -0.22947209362139917695e-3;
        +0.78403922172006662747e-3;
        +0.69728137583658577743e-4;
        -0.59216643735369388286e-3;
        -0.51717909082605921934e-4;
        +0.83949872067208727999e-3;
        +0.72048954160200105591e-4;
        -0.19144384985654775265e-2;
        -0.16251626278391581690e-3;
        +0.64033628338080697948e-2;
        +0.54016476789260451518e-3;
];

[r,c]=size(Z); Z=reshape(Z,1,prod(size(Z))); Y=zeros(size(Z));

n_real=imag(Z)==0;
n_complex=~n_real;
if any(n_real)
  IN=find(n_real);
  Y(IN)=gamma(real(Z(IN)));
end
if any(n_complex)
  IN=find(n_complex);
  z=Z(IN);
  FN=find(real(z)<0);
  z(FN)=-z(FN);
  N=36-abs(z).^2;
  II=find(N>0); I2=find(N<=0);
  N(II)=ceil(sqrt(N(II)+real(z(II)).^2)-real(z(II)));
  N(I2)=0*N(I2);
  z(II)=z(II)+N(II);
  
y=sqrt(2*pi./z).*z.^z.*exp(-z).*(exp(-(0:length(Asymp_coef)-1)'*log(z)).'*Asymp_coef).';
  for l=II
    y(l)=y(l)/prod(z(l)-(1:N(l)));    % Take back the added N
  end
  z(II)=z(II)-N(II);
  y(FN)=-pi./(z(FN).*sin(pi*z(FN)).*y(FN));
  Y(IN)=y;
end

Y=reshape(Y,r,c);


---------------------------------cgammaln.m-------------------------
function Y=cgammaln(Z);
% function Y=cgammaln(Z) returns the value of the log-gamma function for
% complex argument matrix Z.

% Algorithm:
% ----------
% 1) for real Z, use the builtin log-gamma function.
% 2) If real(Z)<0, use the reflection formula for Z -> -Z.
% 3) For |Z|>6 use an order-15 Stirling's approximation.
% 4) Near the zeros |Z+/-1|<=0.1 and |Z+/-2|<=0.1 use a series expansion.
% 5) For the rest, use gamma(z+1)=z*gamma(z) and evaluate gamma(z+N) such that
%    |z+N|>6 .

One_coef=[
    -0.57721566490153286061;
    +0.82246703342411321826;
    -0.40068563438653142846;
    +0.27058080842778454789;
    -0.20738555102867398526;
    +0.16955717699740818997;
    -0.14404989676884611811;
    +0.12550966952474304243;
    -0.11133426586956469049;
    +0.10009945751278180854;
    -0.90954017145829042236e-1;
    +0.83353840546109004037e-1;
    -0.76932516411352191469e-1;
    +0.71432946295361336071e-1;
%    -0.66668705882420468034e-1;
%    +0.62500955141213040754e-1;
%    -0.58823978658684582341e-1;
%    +0.55555767627403611114e-1;
%    -0.52631679379616660732e-1;
];
Asymp_coef=[
 0.08333333333333333;         %  1/12
-0.002777777777777778;        % -1/360
 0.0007936507936507937;       %  1/1260
-0.0005952380952380953;       % -1/1680
 0.0008417508417508417;       %  1/1188
-0.001917526917526918;        % -691/360360
 0.00641025641025641;         %  1/156
];

[r,c]=size(Z); Z=reshape(Z,1,prod(size(Z))); Y=zeros(size(Z));

n_real=imag(Z)==0;
n_one=~n_real & (abs(Z-1)<=0.1 | abs(Z+1)<=0.1);
n_two=~n_real & (abs(Z-2)<=0.1 | abs(Z+2)<=0.1);
n_complex=~(n_real | n_one | n_two);
if any(n_real)
  IN=find(n_real);
  if exist('OCTAVE_VERSION')
    Y(IN)=lgamma(real(Z(IN)));
  else
    Y(IN)=gammaln(real(Z(IN)));
  end
end
if any(n_one)
  IN=find(n_one);
  z=Z(IN);
  FN=find(abs(z+1)<=0.1);
  z(FN)=-z(FN);
  y=exp((1:length(One_coef))'*log(z-1)).'*One_coef;
  y(FN)=i*pi+log(pi)-log(z)-log(sin(pi*z(FN)))-y(FN);
  Y(IN)=y;  
end
if any(n_two)
  IN=find(n_two);
  z=Z(IN);
  FN=find(abs(z+2)<=0.1);
  z(FN)=-z(FN);
  y=exp((1:length(One_coef))'*log(z-2)).'*One_coef+log(z-1);
  y(FN)=i*pi+log(pi)-log(z)-log(sin(pi*z(FN)))-y(FN);
  Y(IN)=y;  
end
if any(n_complex)
  IN=find(n_complex);
  z=Z(IN);
  FN=find(real(z)<0);
  z(FN)=-z(FN);
  N=36-abs(z).^2;
  II=find(N>0); I2=find(N<=0);
  N(II)=ceil(sqrt(N(II)+real(z(II)).^2)-real(z(II)));
  N(I2)=0*N(I2);
  z(II)=z(II)+N(II);                     % Make it asymptotic
  lz=log(z);
  
y=(log(2*pi)-lz)/2+(lz-1).*z+(exp(-(1:2:2*length(Asymp_coef)-1)'*lz).'*Asymp_coef).';
  for l=II
    y(l)=y(l)-sum(log(z(l)-(1:N(l))));    % Take back the added N
  end
  z(II)=z(II)-N(II);
  y(FN)=log(-pi./(z(FN).*sin(pi*z(FN))))-y(FN);
  Y(IN)=y;
end

Y=reshape(Y,r,c);


---------------------------------cdigamma.m-------------------------
function Y=cdigamma(Z);
% function Y=cdigamma(Z) returns the value of the digamma function for
% complex argument matrix Z.

% Algorithm:
% ----------
% 1) If real(Z)<0, use the reflection formula for Z -> -Z.
% 2) For |Z|>6 use an order-15 Stirling's approximation.
% 3) For the rest, use Psi(z+1)=1/z+Psi(z) and evaluate Psi(z+N) such that
%    |z+N|>6 .

Asymp_coef=[
-0.08333333333333333;    %    -1/12
 0.008333333333333333;   %    +1/120
-0.003968253968253968;   %    -1/252
 0.004166666666666667;   %    +1/240
-0.007575757575757576;   %    -1/132
 0.0210927960927961;     %    +691/32760
-0.08333333333333333;    %    -1/12
];

[r,c]=size(Z); Z=reshape(Z,1,prod(size(Z))); Y=zeros(size(Z));

FN=find(real(Z)<0);
Z(FN)=-Z(FN);
N=36-abs(Z).^2;
II=find(N>0); I2=find(N<=0);
N(II)=ceil(sqrt(N(II)+real(Z(II)).^2)-real(Z(II)));
N(I2)=0*N(I2);
Z(II)=Z(II)+N(II);                     % Make it asymptotic
lZ=log(Z);
Y=lZ-1./Z/2+(exp(-(2:2:2*length(Asymp_coef))'*lZ).'*Asymp_coef).';
for l=II
  Y(l)=Y(l)-sum(1./(Z(l)-(1:N(l))));    % Take back the added N
end
Z(II)=Z(II)-N(II);
II=(imag(Z(FN))==0) & (abs(ceil(Z(FN)-0.5)-(Z(FN)-0.5))<=eps);
F1=FN(find(II)); F2=FN(find(~II));
Y(F1)=Y(F1)+1./Z(F1);
Y(F2)=Y(F2)+1./Z(F2)+pi./tan(pi*Z(F2));

Y=reshape(Y,r,c);