bpqcc(n,a)

function [c,k] = bpqcc(n,a)
%BPQCC   Bessel polynomial connection coefficients (n,0).
%
%   Special case of expansion in Bessel polynomials ([1], eq. 7)
% 
%       q_n(a u) = SUM_k c_k^{(n)}(a) q_k(u)
%
%   k = 0 to n where the q_k(u) are Bessel polynomials of order k.
%   Connection coefficients c_k^{(n)}(a) are equal to the special case of
%   linearization coefficients beta_k^{(n,0)}(a).
%
%   [c,n] = bpqcc(n,a)  Column vector c of n+1 Bessel polynomial connection
%   coefficients, and column vector k of indexes 0 to n, describing an
%   expansion of q_n(a u) into a weighted sum of q_k(u).
%
%   Example:
%    n = 3; a = 0.3; [c,k] = bpqcc(n,a);
%    bp = bpqlc(n,a); [b,kk] = bp.bkvec(n,0);
%    format short g; disp([k c kk b]);
%         0          0.7            0          0.7
%         1       0.1974            1       0.1974
%         2       0.0756            2       0.0756
%         3        0.027            3        0.027
%
% Algorithm: Equation from Theorem 2.1 of [1].
%
% References:
% [1] Christian Berg and Christophe Vignat, Linearization coefficients of
% Bessel polynomials, arXiv:math/0506458v1 [math.PR] 

% 2013-06-01 Robert Dickson

ik = 0:n-1;
fac = a.^ik.*exp(gammaln(n+1)-gammaln(n-ik+1)-gammaln(ik+1)- ...
    gammaln(2*n+1)+gammaln(2*(n-ik)+1)+gammaln(2*ik+1));

c = [zeros(n,1); a^n];

for kk = ik
    r = 1:min(n-kk,kk+1);
    v = (1-a).^r.*exp(gammaln(n+2)-gammaln(n-kk+r+1)-gammaln(kk+2-r)+ ... 
    	gammaln(n-kk)-gammaln(n-kk-r+1)-gammaln(r)); 
    c(kk+1) = fac(kk+1)*sum(v);
end

k = (0:n)';

end