bpqnn(n,a)

function [b,k] = bpqnn(n,a)
%BPQNN  Bessel polynomial linearization coefficients (n,n).
%
%   Special case of expansion in Bessel polynomials 
%
%       q_n(a*u)*q_n((1-a)*u) = sum b_k^{(n,n)}(a)*q_k(u)
%
%   k = n to 2*n. The q_k(u) are Bessel polynomials of order k. 
%
%   [b,k] = bpqnn(n,a)  Column vector b of non-zero Bessel polynomial 
%   linearization coefficients and indexes k the same size, k = n to 2*n.
%
%   [b,k] = bpqnn(n)  Defaults to a = 1/2, an important special case
%   that applies to the expansion of a squared Bessel polynomial.
%
%   Example:
%   format short g; n = 2; a = 0.25; [b,k] = bpqnn(n,a); disp([k b]);
%             2      0.23828
%             3      0.35156
%             4      0.41016
%
%   Algorithm: Equation from Theorem 2.2 of [1]. Also see bpqcc, bpqlc.
% 
%   References:
%   [1] Christian Berg and Christophe Vignat, Linearization coefficients of
%   Bessel polynomials, arXiv:math/0506458v1 [math.PR] 

%   2013-06-01 Robert Dickson

% verify n is a valid Bessel polynomial order
if ~isscalar(n) || rem(n,1) ~= 0 || n < 0
    error('bpqnn:arg','n must be a non-negative integer scalar.');
end

if nargin == 1
    % default to a = 1/2
    vk = (0:n)';
    b = exp(log(2)*(-2*n)+3*gammaln(n+1)-2*gammaln(2*n+1) ...
        +gammaln(2*(n-vk)+1)+gammaln(2*(n+vk)+1) ...
        -gammaln(n-vk+1)-gammaln(n+vk+1) ...
        -gammaln(n-vk+1)-gammaln(vk+1));
    k = n + vk;
    return;
end

% Argument a is on [0,1]
if ~isscalar(a) || a < 0 || a > 1
    error('bpqnn:arg','a must be a scalar on [0,1].');
end

b = zeros(2*n+1,1);
for kk = 0:n
    fac = (4*a*(1-a))^kk*exp(2*(gammaln(n+1)-gammaln(2*n+1)) ...
        +log(2)*(-2*n)+gammaln(2*(n-kk)+1)+gammaln(2*(n+kk)+1) ...
        -gammaln(n-kk+1)-gammaln(n+kk+1));
    v = zeros(n-kk+1,1);
    for jj = 0:n-kk
        v(jj+1) = (2*a-1)^(2*jj) *  exp(...
           gammaln(2*n+2)-gammaln(2*(n-jj)+2)-gammaln(2*jj+1) ... 
        	+gammaln(n-jj+1)-gammaln(n-jj-kk+1)-gammaln(kk+1));
    end
    b(n+kk+1) = fac*sum(v);
end

vk = (0:2*n)';
k = vk(b>0);
b = b(b>0);

end % bpqnn