Code covered by the BSD License  

Highlights from
Bessel Polynomial Linearization Coefficients

image thumbnail

Bessel Polynomial Linearization Coefficients

by

 

Linearization coefficients for Bessel polynomials, and certain rational polynomial spectra.

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

Contact us