Easy Octave Algorithm 1

This page contains an Octave version of Algorithm 1, on page 55 of the B-series book. This function is intended to run equally well on Matlab, SciLab, as well as Octave.

Its purpose is to evaluate $a_1$, $a_2$, $\dots$, up to a maximum order of pmax which satisfy the equation $$ a_1+a_2x+a_3x^2 + \cdots = (1-x)^{-a_1} (1-x^2)^{-a_2} \cdots $$ in the sense that the coefficients of each power of $x$ on the two sides must agree.

Octave function definition

function Ntrees = alg1(pmax)
  Ntrees = ones(1,  pmax)
  for i = 2: pmax - 1
    for j = pmax: -1: 1
      a = Ntrees(i)
      for k = 1: (j - 1)/i
        Ntrees(j) = Ntrees(j) + a*Ntrees(j - i*k)
        a = a*(Ntrees(i) + k)/(k + 1)
      end
    end
  end
end