\\% Fast computations of recurrences via matrix or polynomial powering.
\\ Author: Joerg Arndt
\\ License: GPL version 3 or later
\\ online at http://www.jjj.de/pari/
\\ version: 2011-January-19 (12:51)


vec2charpol(m)=
{ /* return characteristic polynomial for the recursion coefficients in m */
    local(d,p);
    d = length(m);
\\    p = x^d - Pol(m, x);
\\    print(":: m=", m);
    p = 'x^d - Pol(m, 'x);
\\    print(":: p=", p);
    return( p );
} /* ----- */

mrec(v, m, k)=
{
/* compute recurrence using matrix power */
/* Example: mrec([0,1],[1,1],15) ==> [610, 987] */
    local(p,M);
    p = vec2charpol(m);
    M = matcompanion(p);
    M = M^k;
    return ( v * M );
} /* ----- */

recstep(v, m, k)=
{ /* update v by k steps according to the recursion coefficients in m */
    local(n,r);
    if ( k<=0, return(v) );  \\ negative k is forbidden
    n = length(m);
    r = vector(n);
    for (i=1, k,
        for (j=1, n-1, r[j]=v[j+1] ); \\ shift left
        r[n] = sum(j=1,n, m[n+1-j]*v[j]);  \\ new element (convolution)
        v = r;
    );
    return( r );
} /* ----- */

frec(v, m, k)=
{
/* compute recurrence using polynomial power (fast!) */
/* Example: frec([0,1],[1,1],15) ==> [610, 987] */
    local(n, pc, pv, pp, px, r, t);
    n = length(m);
    if ( k<=n, return( recstep(v, m, k) ) );  \\ small indices by definition
    pc = vec2charpol(m);
    pp = Mod( 'x, pc );
    px = pp^(k);
    r = vector(n);
    for (i=1, n,
\\        t = component(px,2);
        t = lift(px);
        r[i] = sum(j=1,n, v[j]*polcoeff(t,j-1,'x));
        px *= pp;
    );
    return( r );
} /* ----- */


\\ ==== end of file ====