Understanding an algorithm for computing a matrix polynomial I'm trying to understand this algorithm by Charles Van Loan for evaluating a matrix polynomial $p(\mathbf A)=\sum\limits_{k=0}^q b_k \mathbf A^k=b_0\mathbf I+b_1\mathbf A+\cdots$ (where $\mathbf A$ is $n\times n$):

$\displaystyle s=\lfloor q\rfloor, \quad r=\lfloor q/s\rfloor$
  $\displaystyle \mathbf Y=\mathbf A^s$
  $\displaystyle \text{for}\quad j=1,\dots,n$
  $\displaystyle\quad\quad \mathbf y_0^{(j)}=\mathbf e_j$ ($\mathbf e_j$ is the $j$th column of the identity matrix)
  $\displaystyle\quad\quad \text{for}\quad k=1,\dots,s-1$
  $\displaystyle\quad\quad\quad \mathbf y_k^{(j)}=\mathbf A\mathbf y_{k-1}^{(j)}$
  $\displaystyle\quad\quad \mathbf f_0^{(j)}=\sum_{k=q}^{rs}b_k\mathbf y_{k-rs}^{(j)}$
  $\displaystyle\quad\quad \text{for}\quad k=1,\dots,r$
  $\displaystyle\quad\quad\quad \mathbf f_k^{(j)}=\mathbf Y\mathbf p_{k-1}^{(j)}+\sum_{i=0}^{s-1}b_{s(r-k)+i}\mathbf y_{i}^{(j)}$  

I've tried going through it a number of times, but it seems that the quantity (vector?) $\mathbf p_{k-1}^{(j)}$ was never defined anywhere in the paper. What might it be? Also, the paper claims (or more probably it's how I (mis)understood the paper) that only three matrices need to be stored for the evaluation: $\mathbf A$, $\mathbf Y$, and the final $p(\mathbf A)$. However, I can't see how to properly organize things so that an extra array for storing the $\mathbf y_k^{(j)}$ is not needed (otherwise, four matrices are now needed instead of the three claimed in the paper).
Could anyone enlighten me on how to implement this algorithm properly? Better yet, is there any code available anywhere that implements this algorithm (I'm language-agnostic, so whatever code you post, I can likely translate into the language I'm currently using).
Thanks for any help!
 A: The storage required for the $y^{(j)}_k$ is less than a "full" matrix in proportion as $s \ll n$.  That is, $s$ such vectors are required for each step $j$, and the storage for one step may be reused in the next step.
The Question states $s$ to be the floor of $q$, the degree of the polynomial, but this is mistaken.  Per the discussion of the paper preceding Algorithm I, $s$ can be any integer strictly between 1 and $q$.  While this does not guarantee $s \ll n$, it allows for it.  See the discussion of applying Algorithm II to the Taylor polynomial for $e^A$, where $s=4$ and $q=16$ were chosen.  [Note also the discussion of minimizing "work count" in Algorithm II by choosing $s$ around $\sqrt{q/2}$.]
That leaves still the mysterious $p^{(j)}_k$.  My guess is it is a typo, that $f^{(j)}_k$ is being computed by adjusting $Yf^{(j)}_{k-1}$, after initializing $f^{(j)}_0$, but I need to do a bit of scribbling to confirm that works out.
Added:  My guess above is correct, that we should read $f^{(j)}_{k-1}$ instead of $p^{(j)}_{k-1}$ at one point in Algorithm II where the latter appears.
The algorithm may be understood as combining two ideas.  One is that the matrix polynomial $p(A)$ will be evaluated column by column, i.e. by applying it successively to unit vectors $e_j$, $j=1,..,n$.
The other idea is simple though a bit cumbersome to spell out.   Recall that $1 \lt s \lt q$ where $q$ is the degree of $p$.  Let the polynomial:
$$ p(x) = \sum^q_{i=0} b_i x^i$$
be grouped according to the largest power of $x^s$ that divides each term and rewrite $p(x)$ in an expansion involving powers of $x^s$ and "coefficients" that are polynomials of degree at most $s-1$.  Taking $r = \lfloor q/s \rfloor$, and assuming any indicated coefficients $b_m$ where $m \gt q$ are zero:
$$ p(x) = \sum^r_{k=0} x^{sk} \beta_k(x)$$
where $\beta_k(x) = \sum^{s-1}_{i=0} b_{sk+i} x^i$ are the polynomial "coefficients" as above.
Assume that $Y = A^s$ has already been computed, e.g. by binary exponentiation or other scheme.
The outermost loop of Algorithm II is then iteratively applying $p(A)$ to $n$ unit vectors $v = e_j$ as explained above.  For simplicity we just describe computing $p(A)v$ for an arbitrary vector $v$ of compatible dimension, a column of length $n$ where $A$ is $n \times n$.  In effect we thus omit superscript $(j)$ from the remainder of the discussion.
First is a loop to initialize $y_0 = v$ and $y_i = A y_{i-1}$ for $i = 1,..,s-1$.
Second is a loop that essentially applies Horner's method to get $p(A)v$ using the expansion in powers of $A^s$ and "coefficients" $\beta_k(A)v$.  Note that the latter are column vectors, and for $k = 0,..,r$ we evaluate $\beta_k(A)v$ by taking the appropriate linear combination of previously initialized $y_i$.  That is:
$$ f_0 = \beta_r(A) v $$
$$ f_k = Y f_{k-1} + \beta_{r-k}(A) v $$
for $k = 1,..,r$.  This shows us that the term $p^{(j)}_{k-1}$ is a misprint for $f^{(j)}_{k-1}$ in Algorithm II.
A: Here is some octave code that implements "Algorithm II" quite faithfully. The arrays named with "_shifted" are indexed using $i+1$ in the last component wherever the paper uses $i,$ because arrays in matlab are indexed from 1 instead of 0. All the arithmetic is done mod 29 to make it easy to check against a simpler reference implementation. I haven't tried to make it particularly fast; it's really just meant to demonstrate that the algorithm works.
q=20;
s=floor(sqrt(2*q));
r=floor(q/s);
n=100;
modulo=29;
A=randi(modulo,n,n);
b_shifted=randi(modulo,q+1,1);

pA_van_loan = zeros(n,n);

# Compute A^s
As=eye(n);
for k = 1:s
  As=A*As;
  As=mod(As, modulo);
end

# Van Loan algorithm II
for j = 1:n
  y_shifted = zeros(n,s);
  y_shifted(j,1) = 1;
  for k = 1:s-1
    y_shifted(:,k+1) = A*y_shifted(:,k-1+1);
    y_shifted(:,k+1) = mod(y_shifted(:,k+1), modulo);
  end

  f = zeros(n,1);
  for m = 0:(q-s*r)
    f = f + b_shifted(s*r+m+1) * y_shifted(:,m+1);
  end
  f = mod(f, modulo);
  for k = 1:r
    f = As * f;
    for m = 0:s-1
      f = f + b_shifted(s*(r-k)+m+1) * y_shifted(:,m+1);
    end
    f = mod(f, modulo);
  end

  pA_van_loan(:,j) = f;
end

# Compute polynomial naively for reference
pA_reference = zeros(n,n);
Aj = eye(n);
for j = 0:q
  pA_reference = pA_reference + b_shifted(j+1)*Aj;
  Aj=A*Aj;
  Aj=mod(Aj,modulo);
end
pA_reference=mod(pA_reference,modulo);

assert(pA_van_loan == pA_reference);

