tocol(A)=if(type(A)=="t_COL",A,A~);
zip(f,a,b)=vector(min(#a,#b),i,f(a[i],b[i]));

/* A linear recurrence relation R with polynomial coefficients is given by a
 * matrix whose first column contains the shifts and the second columns
 * contains the polynomials.  the sequence V start
 * at s, which by default is set to 0.
 * that is [-1, x; 0, x^2; 1, x^3] means n*u_{n-1} + n^2*u_n + n^3*u_{n+1} = 0
 */

addhelp("reccoefs","reccoefs(R,V,n,s=0): extends the sequence V up to V_n while satisfying the linear recurrence relation R");

reccoefs(R,V,n,s=0)=
{
  my(A=R[,1],U=R[,2]);
  my(x=variable(U),a,b,ia,ib);
  a=vecmin(A,&ia);
  b=vecmax(A,&ib);
  my(W=concat(tocol(V),vectorv(1+n-#V)));
  for(i=-a,n-b,
    my(r=subst(U,x,s+i));
    W[1+i+b]=sum(j=1,#A,if(j!=ib,r[j]*W[1+i+A[j]]))/-r[ib];
  );
  W;
}

addhelp("reccheck","reccheck(R,V,s=0): return 0 if the sequence V satisfies the linear recurrence relation R");

reccheck(R,V,s=0)=
{
  my(A=R[,1],U=R[,2]);
  my(x=variable(U),a,b,ia,ib);
  a=vecmin(A,&ia);
  b=vecmax(A,&ib);
  for(i=-a,#V-b-1,
    my(r=subst(U,x,s+i));
    if (sum(j=1,#A,r[j]*V[1+i+A[j]]),return(1+i+b))
  );
  0;
}

addhelp("recfind","recfind(V,A,D,s=0): given a vector of shifts, search for polynomials of degree at most D such that V satisfies a linear recurrence relation");

recfind(V,A,D,s=0)=
{
  my(a,b);
  a=vecmin(A); b=vecmax(A);
  my(N=#V+a-b);D++;
  my(M=Mat([concat(vector(#A,k,vector(D,d,(n+s)^(d-1)*V[1+n+A[k]])))| n<-[-a..#V-b-1]]~));
  [Mat([tocol(A),vectorv(#A,k,Polrev(x[1+(k-1)*D..k*D]))]) | x<-matker(M)];
}

addhelp("recshift","recshift(R,d): shift R by d without changing the solutions");

recshift(R,d)=
{
  my(A=R[,1],U=R[,2]);
  my(x=variable(U),a);
  a=vecmin(A);
  Mat([[a+d|a<-A]~,subst(U,x,x+d)]);
}

test()=
{
  my(n=x);
  my(R = [-2,(n-1)^3;-1,-(34*n^3-51*n^2+27*n-5);0,n^3]);
  my(L1=reccoefs(R,[0,6],30),L2=reccoefs(R,[1,5],30));
  my(L3=L2); L3[15]++;
  if([reccheck(R,L1),reccheck(R,L2),reccheck(R,L3)]!=[0,0,15],error("reccheck"));
  my(R2 = recshift(R,2));
  if([reccheck(R2,L1),reccheck(R2,L2),reccheck(R2,L3)]!=[0,0,15],error("reccheck2"));
  my(S1=recfind(L1,[-2..0],3)[1],S2=recfind(L2,[-2..0]~,3)[1]);
  if(S1!=R || S2!=R,error("recfind"));
  my(S13=recfind(L1[4..-1],[-2..0],3,3)[1],S23=recfind(L2[4..-1],[-2..0]~,3,3)[1]);
  if(S13!=R || S23!=R,error("recfind3"));
  my(L4=reccoefs(R,L1[4..5],27,3),L5=reccoefs(R,L2[4..5],27,3));
  if(L4!=L1[4..-1] || L5!=L2[4..-1],error("reccoefs2"));
  zip((x,y)->x/y,L1,L2)/zeta(3);
}