Code covered by the BSD License

# Bessel Polynomial Linearization Coefficients

### Robert Dickson (view profile)

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

```