# restart; read `\c:/Users/Asus/Google Drive/Aek/
# C2Finite/Ansatz.txt`;

# First version: August 24, 2021

# This program accompanies paper
# of Tipaluck and myself: 
# Ansatz in a Nutshell

# The section numbers go according 
# to those in the paper

with(SolveTools,Linear):
with(LinearAlgebra):

################################
# Section 2: Existing ansatz 
#
# Section 2.1 Polynomial
#
# Functions:
# GuessPol(L,start,mu),
# GenPol(Po,n1,x),
#
################################

# Guess Polynomial
# Input: List of numbers A, starting value,
# variable mu
# Output: The corresponding polynomial in mu
# Try: 
# A:= [seq(add(i^2,i=0..n),n=0..20)];
# GuessPol(A,0,n);
# A:= [seq(n^2+3*n,n=3..20)];
# GuessPol(A,3,n);
GuessPol := proc(L,start,mu) option remember;
local eq,var,i,j,C,deg;

deg := 0;
var := op({});

while var = op({}) do 
deg := deg+1;
eq := {seq( add(C[j]*i^j,j=0..deg)=L[i-start+1]
       ,i=start..nops(L)+start-1)}; 
var := {seq(C[j],j=0..deg)}; 
var := Linear(eq,var); 
od:
#print("Degree is", deg);
if deg = nops(L)-1 then
	return("No Solution");
else 
    return(factor(add(subs(var,C[j])*mu^j,j=0..deg)));
fi:

end:


# Guess generating function of 
# polynomial as a sequence
# Input: Polynomial Po, variable n,x
# Output: The corresponding 
# generating function in x
# according to proposition 3.
# Try: 
# For a_n := add(i^2,i=0..n):
# GenPol(n*(n+1)*(2*n+1)/6,n,x);
# GenPol(n^2+3*n,n,x);
GenPol := proc(Po,n1,x) option remember;
local i,n,a,k;

k := degree(Po,n1);
a := n -> subs(n1=n,Po);

add(add(  binomial(k+1,i)*(-1)^i
*a(n-i)*x^n,i=0..n),n=0..k)/(1-x)^(k+1);
end:



################################
# Section 2.2 C-finite
# 
# Functions:
# CRec(n,C,A), GuessC(A,N), 
# GuessGenRat(A,s,x), 
# GenC(A,x),
# 
# CAddition(A,B,x), CCauchy(A,B,x),
# CParSum(A,x), CSubSeq(m,R,N), 
# CTermWise(R1,R2,N),
# CloseC(A,start,n), Multi(R),
#
################################

# Calcuation of C-finite term
# Input: numeric n, List of numbers C, 
# List of numbers A
# Output: The corresponding term a_n where
# a(n) = C[1]*a(n-1)+C[2]*a(n-2)+...,
# with initial conditions a(0)=A[1], a(1)=A[2],...
# Try: 
# Fibonacci numbers
# [seq(CRec(n,[1,1],[0,1]),n=0..30)];
# Fibonacci_{2n}
# [seq(CRec(n,[3,-1],[0,1]),n=0..30)];
CRec := proc(n,C,A) option remember;
local i,m;

if nops(C) <> nops(A) then
	ERROR("InputSize");
elif n <0 then
	ERROR("n must be positive");
fi:

m := nops(C);

if n < m then
	return(A[n+1]);
fi:

add( C[i]*CRec(n-i,C,A),i=1..m);
end:


# Guess the recurrence of C-finite
# Input: List of numbers A, 
# the left shift operator N
# Output: The corresponding 
# C-finite recurrence in N
# Try: 
# A :=[1,0,3,0,11,0,41,0,153,0,571,0,2131]:
# GuessC(A,N);
# A:=[seq(  floor((n/2)^2),n=0..30)];
# GuessC(A,N);
GuessC := proc(A,N)   option remember;
local i,j,ord,eq,var,C;

ord := 0;
var := op({});

while 2*ord < nops(A)-2 and var = op({}) do

	ord := ord+1;
	eq := {seq( add(C[j]*A[i+j],j=0..ord-1)+A[i+ord]=0
       ,i=1..nops(A)-ord),C[ord]=1};
	var := {seq(C[j],j=0..ord)}; 
	var := Linear(eq,var); 
od:

if ord >= nops(A)/2-1 then
	return("NoSolution");
else
    return(factor(subs(var,add(C[j]*N^j,j=0..ord)=0)));
fi:

end:

############################################

# Guess rational generating function
# for C-finite sequence
# Input: List of numbers A, 
# starting value, variable x
# Output: The corresponding 
# generating function in x
# Try: 
# A := [seq(2^i,i=0..20)];
# GuessGenRat(A,0,x);
# A:=[seq(  floor((n/2)^2),n=0..30)];
# GuessGenRat(A,0,x);

GuessGenRat := proc(A,s,x) option remember;
local i,deg,PolA,P,Q;

Q := lhs(GuessC(A,x));
deg := degree(Q,x);
Q := factor(subs(x=1/x,Q)*x^deg);
 
PolA := add(A[i]*x^(i-1+s),i=1..nops(A));

P := expand(PolA*Q);
P := P-add(coeff(P,x,i)*x^i,i=deg+s+3..degree(P,x));

sort(factor(P/Q));
end:


# Generating function of 
# a C-finite sequence
# Input: List of numbers A, 
# the variable x
# Output: The corresponding 
# generating function in x
# according to proposition 7.
# Try: 
# For A := [1,0,3,0,11,0,41,0,153,0,571,0,2131]:
# GenC(N^4-4*N^2+1,[1,0,3,0],N,x);
# For a_n := floor((n/2)^2):
# GenC(N^4-2*N^3+2*N-1,[0,0,1,2],N,x);
GenC := proc(Rec,A,N,x) option remember;
local i,n,r,C,S;

r := degree(Rec,N); 
C := [seq(coeff(Rec,N,i),i=0..r)];

if nops(A) < r then
	ERROR("Not Enough A");
fi:	

S := add(add( C[i+1]*A[n-r+i+1]*x^n
,i=r-n..r),n=0..r-1)
/factor(x^r*subs(N=1/x,Rec));

factor(S);
end:


#############################
# Closure properties of 
# addition, Cauchy product,
# partial sum, sub-sequence 
# and term-wise multiplication
#
#############################

# The generating of addition of 
# two C-finite sequences
# Input: A,B the generating 
# functions of sequence a_n 
# and b_n, also variable x
# Output: the generating function
# of a_n+b_n
# Try: 
# A := x^2/(1+x)/(1-x)^3;
# B := x/(1-x-x^2);
# CAddition(A,B,x);
CAddition := proc(A,B,x) option remember;
factor(A+B);
end:

# The generating of Cauchy product 
# of two C-finite sequences
# Input: A,B the generating 
# functions of sequence a_n 
# and b_n, also variable x
# Output: the generating function
# of c_n = add(a_i*b_{n-i},i=0..n);
# Try: 
# A := x^2/(1+x)/(1-x)^3;
# B := x/(1-x-x^2);
# CCauchy(A,B,x);
CCauchy := proc(A,B,x) option remember;
factor(A*B);
end:

# The generating of partial sum 
# of a C-finite sequence
# Input: A the generating 
# functions of sequence a_n, 
# also variable x
# Output: the generating function
# of c_n = add(a_i,i=0..n);
# Try: 
# A := x^2/(1+x)/(1-x)^3;
# CParsum(A,x);
CParSum := proc(A,x) option remember;
factor(A/(1-x));
end:


# Compute the recurrence of
# subsequence of a_n i.e. c_n=a_(m*n)
# Input: positive integer m,
# R the recurrence of a_n, symbolic N
# Output: the recurrence of c_n=a_(m*n)
# Try: 
# CSubSeq(2,N^4-2*N^3+2*N-1,N);
# CSubSeq(3,N^2-N-1,N);
CSubSeq := proc(m,R,N) option remember;
local i,j,r,A,M,eq,var,d;

r := degree(R,N);
M := [0$(r+1)];

for i from 0 to r do
	A := HoRec(i*m,R,n,N);
	M[i+1] := [seq(coeff(A,N,j),j=0..r-1)]; 
od:

eq := {seq(add(d[i]*M[i+1][j+1],i=0..r)=0,j=0..r-1),
d[0]=1};
var := {seq(d[i],i=0..r)};
var := solve(eq,var);

factor(subs(var,add(d[i]*N^i,i=0..r)=0));
end:


# Compute the recurrence of
# termwise multiplicaiton c_n=a_n*b_n
# Input: R1,R2 the recurrences of
# a_n and b_n respectively, symbolic N
# Output: the recurrence of c_n=a_n*b_n
# Try: 
# CTermWise(N^4-2*N^3+2*N-1,N^2-N-1,N);
CTermWise := proc(R1,R2,N) option remember;
local i,j,k,r,s,A,B,M,eq,var,d;

r := degree(R1,N);
s := degree(R2,N);

M := [0$(r*s+1)];

for i from 0 to r*s do
	A := HoRec(i,R1,n,N);
	B := HoRec(i,R2,n,N);
	M[i+1] := [seq(seq(coeff(A,N,j)*coeff(B,N,k)
			  ,k=0..s-1),j=0..r-1)]; 
od:

eq := {seq(add(d[i]*M[i+1][j+1]
,i=0..r*s)=0,j=0..r*s-1),d[0]=1};
var := {seq(d[i],i=0..r*s)};
var := solve(eq,var);

if var=Null then
	return("NoSolution");
fi:	

factor(subs(var,add(d[i]*N^i,i=0..r*s)=0));
end:


#######################
# Close form formulas
#######################

# Closed form formula for the 
# C-finite sequence A
# Input: The recurrence R, list of 
# initial values A, starting value
# and symbolic N,n
# Output: its closed form solution in n
# Try:
# Example 1:
# CloseC(N^4-2*N^3+2*N-1,[0,0,1,2],0,N,n);

# Example 2:
# T := CloseC(N^2-N-1,[0,1],0,N,n);
# seq(simplify(subs(n=i,T)),i=0..10);

# Example 3:
# T := CloseC(N^2+1,[0,1],0,N,n);
# seq(simplify(subs(n=i,T)),i=0..10);

CloseC := proc(TN,A,start,N,n) option remember;
local i,j,c,r,T,S,eq,var;

r := degree(TN,N);
if nops(A) < r then
	ERROR("Not enough info.");
fi:	

T := Multi([solve(TN,N)]);
S := add( add( c[i][j]*n^i
,i=0..T[j][1]-1)*T[j][2]^n  
,j=1..nops(T));

eq := {seq(subs(n=start+i-1,S)=A[i],i=1..r)};
var := {seq(seq(c[i][j],i=0..T[j][1]-1)
, j=1..nops(T))};
var := solve(eq,var);

subs(var,S);
end:


# Count mutiplicities of list R
# Input: List R
# Output: List of element 
# and its multiplicities
# Try:
# Multi([-1,1,1,1]);
Multi := proc(R) option remember;
local s,r,A,count;

A := [];
for s in {op(R)} do
	count :=0;
	for r in R do
		if s=r then
			count := count+1;
		fi:	
	od:
	A := [op(A),[count,s]];
od:
A;
end:



################################
# Section 2.3: Holonomic
# 
# Functions:
# GuessHo(A,ord,deg,startN,n,N), 
# SubsConst(var,mv,ORD,deg,d),
# SmallD(S,ORD,deg,d), SmallOrd(ORD),
#
# GuessDiff(A,ord,deg,startN,x,D,c),
# HoToDiff(R,A,n,N,x,D), 
# FallFac(n,j), SolveC(s,t), 
#
# HoToDiffHom(R,A,n,N,x,D),
# DiffToHo(L,x,D,n,N),
#
# HoRec(t,R,n,N),
# DEqRec(t,R,x,D),
# HoAdd(R1,R2,n,N,c),
# HoTermWise(R1,R2,n,N,c), 
# HoSubSeq(m,R,n,N,c),
# HoParSum(R,n,N,c),
# HoCauchy(DA,DB,x,D,c),
#
# WimpExpand(t,n,k,mu),
# binom(a,b), TyP1(pow,K,cen,n),
# HoAsymp21(K,n),
# HoAsymp22(K,n),
#
################################

# Guess recurrence of one variable.
# Since the recurrence can be written 
# in many difference ways, user needs 
# to specify the order and max deg of pol. 
# themselves we should start with small 
# ord and deg then increase.

# Input: 1 dimension list, order of rec. 
# and max deg of pol., start point of seq.,
# symbolic n and N.
# Output: recurrence according 
# to Holonomic Ansatz.

# Try:
# GuessHo([seq(n!,n=0..20)],1,1,0,n,N);
# GuessHo([seq(binomial(2*n,n)/(n+1),n=0..25)]
# ,1,1,0,n,N);
# GuessHo([seq(add(1/i,i=1..n),n=1..35)],2,1,1,n,N);
# GuessHo([seq(floor((n/2)^2),n=0..15)],2,1,0,n,N);
# GuessHo([seq(floor((n/2)^2),n=0..15)],4,0,0,n,N);
# GuessHo([seq(n!+n,n=0..35)],4,1,0,n,N);
# Example in the paper
# A:=[1,1,5,23,135,925,7285,64755,641075,
# 6993545,83339745];
# GuessHo(A,2,1,0,n,N);

GuessHo := proc(A,ord,deg,startN,n,N) option remember;
local i,j,k,T,e,eq,mv,var,f,d,na,s,shift;

if type(A,list) = false then
	return("Input must be List");
elif (ord+1)*(deg+1)+ord-1 >= nops(A)-2 then
	return("Not Enough Info");
fi:

e := add(add(d[i][j]*na^j
,j=0..deg)*f[n+i],i=0..ord)=0;

eq := {seq( subs({na=k-1+startN,n=k},e)
,k=1..nops(A)-ord)};
eq := subs({seq(f[k]=A[k]
,k=1..nops(A))},eq);

mv := {seq(seq(d[i][j],j=0..deg),i=0..ord)}; 
var := Linear(eq,mv);

if {seq(rhs(T), T in var)} = {0} then
	return(NoSolution);
fi:

# (May cause a problem here, not sure)	
# In case there are more than one solution
var := [var][1];

e := subs(SubsConst(var,mv,ord,deg,d),e);
shift := lhs(subs({na=n,
seq(f[n+i]=N^i,i=0..ord)},e));

k := ilcm(seq(denom(coeff(shift,N,i)),i=0..ord));
T := factor(add(factor(k*coeff(shift,N,i))*N^i,i=0..ord))=0;

# Another choice of output
T := expand(lhs(T)); 
add(factor(coeff(T,N,i))*N^i,i=0..ord)=0;
end:


###################
# Aux. Functions
###################

#substitute constant into the solution
SubsConst := proc(var,mv,ORD,deg,d) option remember;
local S,v,V,T;

V := var;
S := {};

for v in V do
	if type(rhs(v),constant)=false then
		S := S union {lhs(v)};
	fi:
od:

if S = {} then
	return(V);
fi:

V := subs(SmallD(S,ORD,deg,d)=1, V) 
	union{SmallD(S,ORD,deg,d)=1};
V := Linear(V,mv);

T := {};
while V <> {} do 
	for v in V do		
		if rhs(v) = lhs(v) then
			V := V minus {v};
			V := subs(rhs(v)=0,V);
			T := T union {rhs(v)=0};
		elif type(rhs(v),constant) then
			V := V minus {v};
			V := subs(v,V);
			T := T union {v};
		fi:
	od:
od:

return(T);
end:


#Smallest order
SmallD := proc(S,ORD,deg,d) option remember;
local i,j,k,s,T,dim1,Aek;  

dim1 := nops(ORD[1]); 
Aek := {op(ORD)};

if dim1 = 1 then
	for i from 0 to ORD do
	for j from 0 to deg do
		if member(d[i][j],S) then 
			return(d[i][j]);
		fi:
	od:
	od:

elif dim1 = 2 then   
	while Aek <> {} do
		T := SmallOrd(Aek); 
		Aek := Aek minus {T};
		for i from 0 to deg do
		for j from 0 to deg do
			if member(d[T[1]][T[2]][i][j],S) then 
				return(d[T[1]][T[2]][i][j]);
			fi:
	od:od:
	od:

else
	ERROR("NoCaseYet");
fi:

ERROR("NoSmallest");
end:

# Try: SmallOrd([[1,3],[4,2],[1,5]]);
SmallOrd := proc(ORD) option remember;
local s,m,NO;

if nops(ORD) = 1 then
	return(ORD[1]);
fi:

NO := [];
m := min(seq( s[1],s in ORD));

for s in ORD do
	if s[1] = m then
		NO := [op(NO),[op(2..nops(s),s)]];
	fi
od:

[m,op(SmallOrd(NO))];
end:



##############################################

# Find diff eq. by guessing.
# For direct calculation try HoToDiff().

# Input: 1 dimension list, order of rec. 
# and max deg of pol., start point of seq.,
# symbolic n and N.
# Output: the diff. eq. with polynomial 
# coeff. that satisfies the generating 
# function F(x) = sum a_n*x^n.

# Try:
# GuessDiff([seq(n!,n=0..20)],1,1,0,x,D,c);
# GuessDiff([seq(binomial(2*n,n)/(n+1),n=0..25)]
# ,1,1,0,x,D,c);
# GuessDiff([seq(add(1/i,i=1..n),n=1..35)],2,1,1,x,D,c);
# GuessDiff([seq(floor((n/2)^2),n=0..15)],2,1,0,x,D,c);
# GuessDiff([seq(floor((n/2)^2),n=0..15)],4,0,0,x,D,c);
# GuessDiff([seq(n!+n,n=0..35)],4,1,0,x,D,c);
GuessDiff := proc(A,ord,deg,startN,x,D,c) 
option remember;
local i,j,n,N,B,f1,rp,kp,L,eq,var,v,R,S;

B := GuessHo(A,ord,deg,startN,n,N); 
if B = "Input must be List" or 
B = "Not Enough Info" then
	return(B);
fi:

if nops(A) <= (ord+deg+1)*(deg+1)+3 then
	ERROR("Not enough info");
fi:

f1 := add( A[i-startN+1]*x^i
,i=startN..startN+nops(A)-1);
rp := deg;
kp := deg+ord;

L := add(add(c[i,j]*x^i*diff(f1,x$j)
,i=0..kp),j=1..rp)
+add(c[i,0]*x^i*f1,i=0..kp);

eq := {seq( coeff(L,x,i)=0,
i=startN+ord..startN+nops(A)-1 )};
var := {seq(seq(c[i,j],i=0..kp),j=0..rp)};
var := solve(eq,var);
#print(var);

R := add(add(c[i,j]*x^i*D^j,i=0..kp),j=0..rp);
R := subs(var,R);
L := simplify(subs(var,L));

sort(factor(R))*f(x)=add(coeff(L,x,i)*x^i
,i=startN..startN+ord-1);
end:


######################################

# Direct calculaiton of the diff. eq. 
# according to theorem 12 in the paper.

# Input: the holonomic relation of a_n,
# the list of actual sequence A,
# symbolic n,N,x,D
# Output: the diff. eq. with polynomial 
# coeff. that satisfies the generating 
# function F(x) = sum a_n*x^n. 

# Try:
# Example 1:
# A:=[seq(n!,n=0..20)];
# GuessDiff(A,1,1,0,x,D,c);
# HoToDiff(n+1-N,[1],n,N,x,D);
#
# Example 2:
# B := [seq(binomial(2*n,n)/(n+1),n=0..25)]:
# GuessDiff(B,1,1,0,x,D,c);
# HoToDiff(4*n+2-(n+2)*N,[1],n,N,x,D);
#
# Example 3:
# C:=[seq(add(1/i,i=1..n),n=1..35)]:
# GuessDiff(C,2,1,1,x,D,c);
# HoToDiff(n+1-(2*n+3)*N+(2+n)*N^2,[0,1],n,N,x,D);
# 
# Example 4:
# E := [seq(floor((n/2)^2),n=0..15)]:
# GuessDiff(E,2,1,0,x,D,c);
# HoToDiff(n+2+2*N-N^2*n,[0,0],n,N,x,D);
#
# Example 5:
# GuessDiff(E,4,0,0,x,D,c); 
# HoToDiff(1-2*N+2*N^3-N^4,[0,0,1,2],n,N,x,D);
#
# Example 6:
# F:=[seq(n!+n,n=0..75)]:
# T:=GuessDiff(F,3,2,0,x,D,c);
# subs({c[4,2]=4,c[5,2]=0},T);
# R := 3*n+3+(2*n+5)*(1+2*n)*N
# +(-1-15*n-4*n^2)*N^2+(-1+4*n)*N^3;
# HoToDiff(R,[1,2,4],n,N,x,D);

HoToDiff := proc(R,A,n,N,x,D) option remember;
local i,j,t,s,m,r,k,b,T,B;

r := degree(R,N);
k := degree(R,n);

if nops(A) < r then
	ERROR("Not enough Info.");
fi:	

for t from 0 to r do
for s from 0 to k do
	b[s,t] := coeff(coeff(R,N,t),n,s);
od: od:

T := add(add( b[s,t]*add(  
SolveC(s,t)[j+1]*x^(j+r-t)*D^j
,j=0..s),s=0..k),t=0..r);

B := add(add( b[s,t]*add(  
A[m+1]*(m-t)^s*x^(m+r-t)
,m=0..t-1),s=0..k),t=0..r);

add( factor(coeff(T,D,i))*D^i,i=0..k)*f(x)=factor(B);
end:


# Falling factorial
# Input: integers n and j
# Output: Falling factorial (n)_j
# Try: FallFac(n,2);
FallFac := proc(n,j) option remember;
local i;
mul(n-i,i=0..j-1);
end:


# Solve for c_j^(s,t) 
# in HoToDiff of theorem 8
# Input: non-negative 
# integers s and t
# Output: list of c^(s,t) of length s+1

# Try: SolveC(3,2);
SolveC := proc(s,t) option remember;
local j,n,c,eq,var,T;

T := add( c[j]*FallFac(n+t,j),j=0..s);
T := expand(T);

eq := {seq(coeff(T,n,j)=0,j=0..s-1)}
union {coeff(T,n,s)=1};
var := {seq(c[j],j=0..s)};
var := solve(eq,var);

subs(var,[seq(c[j],j=0..s)]);
end:


###############################################

# Calculation of the homogeneous 
# diff. eq. according to corollary 13

# Input: the holonomic relation of a_n,
# the list of actual sequence A,
# symbolic n,N,x,D
# Output: the homogeneous diff. eq. 
# with polynomial coeff. that 
# satisfies the generating function 
# F(x) = sum a_n*x^n. 

# Try: 
# Example 1:
# A:=[seq(n!,n=0..20)];
# HoToDiff(n+1-N,[1],n,N,x,D);
# HoToDiffHom(n+1-N,[1],n,N,x,D);

# Example 2:
# B := [seq(binomial(2*n,n)/(n+1),n=0..25)]:
# HoToDiff(4*n+2-(n+2)*N,[1],n,N,x,D);
# HoToDiffHom(4*n+2-(n+2)*N,[1],n,N,x,D);

# Example 3:
# C := [seq(add(1/i,i=1..n),n=1..35)]:
# HoToDiff(n+1-(2*n+3)*N+(2+n)*N^2,[0,1],n,N,x,D);
# HoToDiffHom(n+1-(2*n+3)*N+(2+n)*N^2,[0,1],n,N,x,D);
# Test:
# F := add(C[i]*x^i,i=1..55):
# expand(-x^2*F-3*x^2*(-1+x)*diff(F,x)
# -x^2*(-1+x)^2*diff(F,x,x));

HoToDiffHom := proc(R,A,n,N,x,D) option remember;
local i,eq,d,L1,R1,L2,R2,T;

eq := HoToDiff(R,A,n,N,x,D);
L1 := expand(lhs(eq)/f(x));
R1 := rhs(eq);

d := degree(L1,D);

if R1=0 then
	return(eq);
fi:
	
L2 := diff(L1,x)+L1*D;
R2 := diff(R1,x);

T := R2*L1-R1*L2;

add( factor(coeff(T,D,i))*D^i,i=0..d+1)*f(x)=0;
end:


# Theorem 14 in the paper
# Convert diff eq. of f(x)
# to holonomic recurrence of a_n

# Input: homogeneous diff. eq. L,
# symbolic x,D,n,a
# Ouput: holonomic recurrence
# in term of a and n. 

# Try: 
# Example 1: n!
# DiffToHo(x^2*D^2+(3*x-1)*D+1,x,D,n,N);

# Example 2: add(1/i,i=1..n)
# L := (x^3-2*x^2+x)*D^3
# +(7*x^2-9*x+2)*D^2+(10*x-6)*D+2:
# DiffToHo(L,x,D,n,N);

DiffToHo := proc(L,x,D,n,N) option remember;
local s,t,r,k,M;

r := degree(L,D);
k := degree(L,x);

M := add(add(
coeff(coeff(L,x,s),D,t)
*FallFac(n+t-s,t)*N^(t-s)   
,s=0..k),t=0..r);

factor(M);
end:



####################################
# Section 2.3: 5 Closure properties: 

# Substitute recurrence to a_{n+t}
# Input: non-negative integer t,
# annihilator relation R, 
# symbolic n,N
# Output: annihilator of a_{n+t} 
# written with the basis 1,N,...,N^r.

# Try: 
# HoRec(1,n+1-N,n,N);
# HoRec(2,n+1-N,n,N);
# HoRec(1,4*n+2-(n+1)*N,n,N);
# HoRec(2,4*n+2-(n+1)*N,n,N);
HoRec := proc(t,R,n,N) option remember;
local i,ord,k,S;

ord := degree(R,N);
if t < ord then
	return(N^t);
fi:

k := t-ord; 
S := expand(subs({n=n+k},R)*N^k);
S := expand(solve(S=0,N^t));

for i from t-1 by -1 to ord do
	S := expand(subs(N^i=HoRec(i,R,n,N),S));
od:
factor(S);
end:



# Substitute relation to diff(f(x),x$t)
# Input: non-negative integer t,
# annihilator relation R, 
# symbolic x,D
# Output: annihilator of diff(f(x),x$t) 
# written with the basis 1,D,...,D^r.

# Try: 
# DEqRec(1,x^2*D^2+(3*x-1)*D+1,x,D);
# DEqRec(3,x^2*D^2+(3*x-1)*D+1,x,D);

DEqRec := proc(t,R,x,D) option remember;
local i,ord,S;

ord := degree(R,D);
if t < ord then
	return(D^t);
fi:

S := R;
for i from 1 to t-ord do 
	S := expand(S*D+diff(S,x));
od:	
S := expand(solve(S=0,D^t));

for i from t-1 by -1 to ord do
	S := expand(subs(D^i=DEqRec(i,R,x,D),S));
od:
factor(S);
end:


# Closure properties
# 1.addition c_n=a_n+b_n
#
# Input: annihilator relation 
# R1 and R2, symbolic n,N,t
# Output: annihilator of a_n+b_n 
# written as combination of n and N.
# 
# Try: 

# Example 1:
# A:=[seq(binomial(2*n,n)/(n+1),n=0..30)]:
# R1:=lhs(GuessHo(A,1,1,0,n,N));
# B:=[seq(sum(1/i,i=1..n),n=1..30)]:  
# R2 :=lhs(GuessHo(B,2,1,1,n,N));
# HoAdd(R1,R2,n,N,c);
# Test:
# C:=[seq(A[n+1]+B[n],n=1..30)]:
# seq(-2*(n+1)*(3*n+7)*(1+2*n)*(2+n)^2*C[n]
# +(5+3*n)*(2+n)*(9*n^3+43*n^2+58*n+20)*C[n+1]
# -(216*n+241*n^2+111*n^3+64+18*n^4)*(3+n)*C[n+2]
# +(3+n)*(4+n)*(3*n+4)*(n+1)^2*C[n+3],n=1..27);

# Example 2:
# A:=[seq(n!,n=0..20)]:
# R1:=lhs(GuessHo(A,1,1,0,n,N));
# B:=[seq(binomial(2*n,n),n=0..20)]:  
# R2 :=lhs(GuessHo(B,1,1,0,n,N));
# HoAdd(n+1-N,4*n+2-(n+1)*N,n,N,c);
# Test:
# C:=[seq(A[n]+B[n],n=1..21)];
# seq(2*(1+2*n)*(n+1)*(n^2-2)*C[n+1]
# -(n^4+6*n^3-3*n^2-20*n-8)*C[n+2]
# +(n+2)*(-1-2*n+n^2)*C[n+3],n=0..18);

HoAdd := proc(R1,R2,n,N,c) option remember;
local i,j,ord,M,eq,var,S,de;

ord := degree(R1,N)+degree(R2,N);
#Remark: rank could even be smaller than ord.

M := Matrix([seq([ 
seq( factor(coeff(HoRec(i,R1,n,N),N,j))
,j=0..degree(R1,N)-1),
seq( factor(coeff(HoRec(i,R2,n,N),N,j))
,j=0..degree(R2,N)-1)
],i=0..ord)]);

eq := {
seq( add(c[i]*M[i+1,j],i=0..ord)=0, 
j=1..degree(R1,N)+degree(R2,N))};
var := {seq(c[i],i=0..ord)};
var := solve(eq,var);

S := [seq(subs(var,c[i]),i=0..ord)]; 
de := lcm( seq(denom(S[i+1]),i=0..ord));
add(factor(S[i+1]*de)*N^i,i=0..ord)=0;
end:


# Closure properties
# 2. termwise multiplicaiton c_n=a_n*b_n

# Input: R1,R2 the recurrences of
# a_n and b_n respectively, symbolic n,N,c
# Output: the recurrence of c_n=a_n*b_n

# Try: 
# Example 1:
# A:=[seq(binomial(2*n,n)/(n+1),n=0..30)]:
# R1:=lhs(GuessHo(A,1,1,0,n,N));
# B:=[seq(sum(1/i,i=1..n),n=1..30)]:  
# R2 :=lhs(GuessHo(B,2,1,1,n,N));
# HoTermWise(R1,R2,n,N,c);
# Test:
# C:=[seq(A[n+1]*B[n],n=1..30)]:
# seq( 4*(3+2*n)*(1+2*n)*(n+1)*C[n]
# -2*(3+2*n)^2*(2+n)*C[n+1]
# +(2+n)^2*(3+n)*C[n+2],n=1..28);

# Example 2:
# A:=[seq(n!,n=0..20)]:
# R1:=lhs(GuessHo(A,1,1,0,n,N));
# B:=[seq(binomial(2*n,n),n=0..20)]:  
# R2 :=lhs(GuessHo(B,1,1,0,n,N));
# HoTermWise(n+1-N,4*n+2-(n+1)*N,n,N,c);
# Test:
# C:=[seq(A[n]*B[n],n=1..21)];
# seq(-2*(2*n+1)*C[n+1]+C[n+2],n=0..19);

HoTermWise := proc(R1,R2,n,N,c) option remember;
local i,j,k,r,s,A,B,M,eq,var,S,de;

r := degree(R1,N);
s := degree(R2,N);

M := [0$(r*s+1)];

for i from 0 to r*s do
	A := HoRec(i,R1,n,N);
	B := HoRec(i,R2,n,N);
	M[i+1] := [seq(seq(coeff(A,N,j)*coeff(B,N,k)
			  ,k=0..s-1),j=0..r-1)]; 
od:

eq := {seq(add(c[i]*M[i+1][j+1]
,i=0..r*s)=0,j=0..r*s-1)};
var := {seq(c[i],i=0..r*s)};
var := solve(eq,var);

if var=Null then
	ERROR("NoSolution");
fi:	

S := [seq(subs(var,c[i]),i=0..r*s)]; 
de := lcm( seq(denom(S[i+1]),i=0..r*s));
add(factor(S[i+1]*de)*N^i,i=0..r*s)=0;
end:



# Closure properties
# 3. sub-sequence of a_n i.e. c_n=a_(m*n)

# Input: positive integer m,
# R the recurrence of a_n, symbolic n,N
# Output: the recurrence of c_n=a_(m*n)

# Try: 
# Example 1:
# A :=[seq( binomial(2*n,n)/(n+1),n=0..31)]:
# R :=lhs(GuessHo(A,1,1,0,n,N));
# HoSubSeq(2,4*n+2-(n+2)*N,n,N,c);
# Test:
# seq(-2*(4*n+1)*(4*n+3)*B[2*n+1]
# +(n+1)*(2*n+3)*B[2*n+3],n=0..14);
#
# Example 2:
# A :=[seq( add(1/i,i=1..n),n=1..31)]:
# R :=lhs(GuessHo(A,2,1,1,n,N));
# HoSubSeq(2,n+1-(2*n+3)*N+(2+n)*N^2,n,N,c);
# Test:
# seq( (4*n+7)*(1+2*n)*(n+1)*A[2*n]
# -(4*n+5)*(4*n^2+10*n+5)*A[2*n+2]
# +(3+4*n)*(3+2*n)*(n+2)*A[2*n+4],n=1..13);

HoSubSeq := proc(m,R,n,N,c) option remember;
local i,j,ord,A,M,eq,var,S,de;

ord := degree(R,N);
M := [0$(ord+1)];

for i from 0 to ord do
	A := HoRec(i*m,R,n,N);
	A := subs(n=m*n,A);
	M[i+1] := [seq(coeff(A,N,j),j=0..ord-1)]; 
od:

eq := {seq(add(c[i]*M[i+1][j+1]
,i=0..ord)=0,j=0..ord-1)};
var := {seq(c[i],i=0..ord)};
var := solve(eq,var);

S := [seq(subs(var,c[i]),i=0..ord)]; 
de := lcm( seq(denom(S[i+1]),i=0..ord));
add(factor(S[i+1]*de)*N^i,i=0..ord)=0;
end:


# Closure properties
# 4. partial sum
#
# Input: annihilator relation 
# R, symbolic n,N,t
# Output: annihilator of sum_{i=0}^n a_i
# written as combination of n and N.
# 
# Try: 
# Example 1:
# A :=[seq( add(1/i,i=1..n),n=1..31)]:  
# R :=lhs(GuessHo(A,2,1,1,n,N));
# HoParSum(n+1+(-3-2*n)*N+(2+n)*N^2,n,N,c);
# Test:
# B :=[seq( add(A[i],i=1..n),n=1..31)]:
# seq( -(2+n)*B[n]+(7+3*n)*B[n+1]
# -(8+3*n)*B[n+2]+(3+n)*B[n+3],n=1..28);

# Example 2:
# A :=[seq(binomial(2*n,n),n=0..31)]:  
# R :=lhs(GuessHo(A,1,1,0,n,N));
# HoParSum(4*n+2-(n+1)*N,n,N,c);
# Test:
# B :=[seq( add(A[i],i=1..n),n=1..31)];
# seq(2*(2*n+3)*B[n+1]-(5*n+8)*B[n+2]
# +(n+2)*B[n+3],n=0..28);
HoParSum := proc(R,n,N,c) option remember;
local i,j,ord,M,eq,var,S,de;

ord := degree(R,N);

M := Matrix([
seq([0$i,1,0$(ord-i-1)],i=0..ord-1),
[seq(-coeff(R,N,j)/coeff(R,N,ord),j=0..ord-1)]
]);

eq := {
seq( add(c[i]*M[i+1,j],i=0..ord)=0, j=1..ord)};
var := {seq(c[i],i=0..ord)};
var := solve(eq,var);

S := [seq(subs(var,c[i]),i=0..ord)];
de := lcm( seq(denom(S[i+1]),i=0..ord));
S := add(factor(S[i+1]*de)*(N^i-N^(i-1)),i=0..ord);
S := expand(subs({n=n+1},S)*N);
add(factor(coeff(S,N,i))*N^i,i=0..ord+1)=0;
end:


# Closure properties
# 5. Cauchy products
#
# Input: DA,DB the homogeneous 
# diff.eq. relation of A(x) and B(x) 
# respectively, symbolic x,D,c
# Output: the recurrence of f(x)=A(x)*B(x)

# Try: 
# Example 1:
# c_n = sum_{i=1}^n i!*(n-i)!
# DA := x^2*D^2+(3*x-1)*D+1; DB:=DA:
# T := HoCauchy(DA,DB,x,D,c); 
# Test:
# sort(simplify(subs({c[4]=0,c[3]=1},T)/x^2));
# C:= [seq(add(i!*(n-i)!,i=0..n),n=0..25)]:
# HoToDiffHom(2*N^2-(3*n+7)*N+(n+2)^2,C,0,n,N,x,D);

# Example 2:
# A:=[seq(binomial(2*n,n)/(n+1),n=0..50)]:
# R1:=lhs(GuessHo(A,1,1,0,n,N));
# DA:=lhs(HoToDiffHom(R1,A,0,n,N,x,D))/f(x);
# B := [seq(n!,n=0..50)]:
# R2:=lhs(GuessHo(B,1,1,0,n,N));
# DB:=lhs(HoToDiffHom(R2,B,0,n,N,x,D))/f(x);
# HoCauchy(DA,DB,x,D,c); 
# Test:
# C := [seq(add(A[i+1]*B[n+1-i]
# ,i=0..n),n=0..50)]:
# S := sort(add(C[i]*x^(i-1),i=1..nops(C)));
# And plug S back to HoCauchy to check

# Example 3:
# MN:=100:A:=[seq(binomial(2*n,n)/(n+1),n=0..MN)]:
# R1:=lhs(GuessHo(A,1,1,0,n,N));
# B:=[seq(sum(1/i,i=1..n),n=0..MN)]:  
# R2 :=lhs(GuessHo(B,2,1,0,n,N));
# DA:=lhs(HoToDiffHom(R1,A,0,n,N,x,D))/f(x);
# DB:=lhs(HoToDiffHom(R2,B,0,n,N,x,D))/f(x);
# T := HoCauchy(DA,DB,x,D,c);
# Test:
# C:= [seq(add(A[i+1]*B[n+1-i]
# ,i=0..n),n=0..MN)]:
# S := sort(add(C[i]*x^(i-1),i=1..nops(C)));
# MT := expand(subs({c[0]=1},T)); 
# subs(D^k,diff(S,x$k)), the result seems fine 

HoCauchy := proc(DA,DB,x,D,c) option remember;
local i,j,k,t,r,s,A,B,M,eq,var,S,de;

r := degree(DA,D);
s := degree(DB,D);

M := [0$(r*s+1)];

for i from 0 to r*s do
	A := DEqRec(i,DA,x,D);
	B := DEqRec(i,DB,x,D); 
	M[i+1] := [seq(seq(add(binomial(i,t)
			*coeff(DEqRec(t,DA,x,D),D,j)
			*coeff(DEqRec(i-t,DB,x,D),D,k),
			t=0..i),k=0..s-1),j=0..r-1)]; 
od:

eq := {seq(add(c[i]*M[i+1][j+1]
,i=0..r*s)=0,j=0..r*s-1)};
var := {seq(c[i],i=0..r*s)};
var := solve(eq,var);

if var=Null then
	ERROR("NoSolution");
fi:	

S := [seq(subs(var,c[i]),i=0..r*s)]; 
de := lcm( seq(denom(S[i+1]),i=0..r*s));
add(factor(S[i+1]*de)*D^i,i=0..r*s)*f(x)=0;
end:


############################
# Section 2.3: 6 Asymptotic 
# approximation solutions 
############################

# Expansion of main formula 
# in Wimp
# Input: number of term t,
# symbolic n,k,mu
# Output: the expansion up
# to t terms
# Try:
# WimpExpand(2,n,k,mu);
# WimpExpand(3,n,k,mu);
WimpExpand := proc(t,n,k,mu) option remember;
local i,j,K,N,T;
K := add(binom(mu*k,i)*(k/n)^i,i=0..t);
K := expand(K);
N := i -> add((mu*k*(-k/n)^i/(i+1))^j/j!,j=0..t/i);

T := K*mul(N(i),i=1..t);
T := expand(T);

[seq(factor(sort(coeff(T,n,-i))),i=0..t)];
end:

# Symbolic binomial
# Input: symbolic a and numeric b
# Output: binomial of a choose b
# Try: binom(n,3);
binom := proc(a,b) option remember;
local i;
mul(a-i,i=0..b-1)/b!;
end:


# Input: power pow, number of 
# terms K, center of expansion of x 
# cen and symbolic n
# Output: Taylor series expansion 
# of (n+x)^pow up to K terms
# (then set x=cen).
# Try:
# TyP1(1/2,3,1,n);
TyP1 := proc(pow,K,cen,n) option remember;
local i,x,A,M;

A := taylor( (n+x)^pow,x=0,K); 
M := add( coeff(A,x,i)*x^i,i=0..K-1);
subs(x=cen,M);
end:


# Asymptotic approximation Example 2
# Method of undertermined coeff. for
# the first solution
# Input: number of term K and symbolic n
# Output: list of parameters c1,c2,...,ck 
# Try: HoAsymp21(3,n);
HoAsymp21 := proc(K,n) option remember;
local i,c,y,L,eq,var;

y := n -> n^2*(1+add(c[i]*n^(-i),i=1..K));
L := (n+2)*y(n)+2*y(n+1)-n*y(n+2);
L := expand(L);

L := subs({seq(1/(n+1)^i=TyP1(-i,K+2,1,n),i=1..K)},L);
L := subs({seq(1/(n+2)^i=TyP1(-i,K+2,2,n),i=1..K)},L);
L := expand(L);

eq := [seq(coeff(L,n,-i),i=0..K)]; 
var := {seq(c[i],i=1..K)};
solve(eq,var);
end:


# Asymptotic approximation Example 2
# Method of undertermined coeff. for
# the second solution
# Input: number of term K and symbolic n
# Output: list of parameters c1,c2,...,ck 
# Try: HoAsymp22(3,n);
HoAsymp22 := proc(K,n) option remember;
local i,c,y,L,eq,var;

y := n -> (-1)^n*(1+add(c[i]*n^(-i),i=1..K));
L := (n+2)*y(n)+2*y(n+1)-n*y(n+2);
L := expand(L);

L := subs({seq(1/(n+1)^i=TyP1(-i,K,1,n),i=1..K)},L);
L := subs({seq(1/(n+2)^i=TyP1(-i,K,2,n),i=1..K)},L);
L := expand(L/(-1)^n);

eq := [seq(coeff(L,n,-i),i=0..K)]; print(eq);
var := {seq(c[i],i=1..K)};
solve(eq,var);
end:



#############################
# Section 3: C^2-finite
#
# Functions:
# C2Rec(n,CC,A),
# GuessC2c1(A,n,N), GuessC2c2(A,n,N),
# 
# C2ToDiff(R,S,A,n,N,x,D),
# C2ToDiffHom(R,S,A,n,N,x,D),
# DiffToC2(E,S,x,D,n,N),
#
# C2Close(),
#
##############################

# Section 3: 1,2 ansatz and examples

# Generate C^2-finite sequences
# Input: non-negative integer n,
# List of C-finites, list of
# initial condition of a_n.
# Output: the n_{th} term 
# of the C^2-finite sequence.

# Try: 
# Example 1: a_{n+1}=F_{n+2}*a_n
# [seq(C2Rec(n,[[[1,1],[1,2]]],[1]),n=0..20)];
# Example 2: a_{n+2}=a_{n+1}+2^n*a_n
# [seq(C2Rec(n,[[[1],[1]],[[2],[1]]]
# ,[1,1]),n=0..20)];
# Example 3: a_{n+2}=F_{n+1}*a_{n+1}+F_n*a_n
# [seq(C2Rec(n,[[[1,1],[1,1]],[[1,1],[0,1]]]
# ,[1,1]),n=0..20)];

C2Rec := proc(n,CC,A) option remember;
local i,m;

if nops(CC) <> nops(A) then
	ERROR("InputSize");
elif n < 0 then
	ERROR("n must be positive");
fi:

m := nops(CC);

if n < m then
	return(A[n+1]);
fi:

add( CRec(n-m,CC[i][1],CC[i][2])
*C2Rec(n-i,CC,A),i=1..m);
end:


############################################
# Section 3: 3. Guessing

# Guess C2 of order 2 with coeff. of order 1
# Input: List of sequence A, symbolic n and N
# Output: the C^2-finite relation of
# order 2 with coeffs as C-finite of order 1

# Try:
# A:=[seq(C2Rec(n,[[[3],[2]]
# ,[[2],[5]]],[1,1]),n=1..20)]:
# GuessC2c1(A,n,N);
GuessC2c1 := proc(A,n,N) option remember;
local i,eq,var,a,b,c1,c2;

eq :={seq(A[i+2]-c1*a^i*A[i+1]-c2*b^i*A[i]=0
,i=1..nops(A)-2)};
var :={a,b,c1,c2};
var := solve(eq,var);

if var = op({}) then
	return("NoSolution");
fi:

subs(var,N^2-c1*a^n*N-c2*b^n=0);
end:


# Guess C2 of order 2 with coeff. of order 2
# Input: List of sequence A, symbolic n and N
# Output: the C^2-finite relation of
# order 2 with coeffs as C-finite of order 2

# Try:
# A:=[seq(C2Rec(n,[[[1,1],[1,1]]
# ,[[1,1],[0,1]]],[1,1]),n=1..20)]:
# GuessC2c2(A,n,N);
GuessC2c2 := proc(A,n,N) option remember;
local i,eq,var,a,b,d,e,c1,c2,c3,c4;

eq :={seq(A[i+2]-(c1*a^i+c2*b^i)*A[i+1]
-(c3*d^i+c4*e^i)*A[i]=0
,i=1..nops(A)-2)};
var :={a,b,d,e,c1,c2,c3,c4};
var := solve(eq,var);

if var = op({}) then
	return("NoSolution");
fi:

subs(var,N^2-(c1*a^n+c2*b^n)*N-(c3*d^n+c4*e^n)=0);
end:


############################################
# Section 3: 4.Generating functions

# Direct calculaiton of the diff. eq. 
# according to theorem 17 in the paper.

# Input: the C^2-finite relation of a_n,
# set of roots alphas, the list of 
# actual sequence A, symbolic n,N,x,D
# Output: the diff. eq. with polynomial 
# coeff. that satisfies the generating 
# function F(x) = sum a_n*x^n. 

# Try:
# Example 1:
# A := [seq(C2Rec(n,[[[1,1],[1,2]]],[1]),n=0..20)]:
# C2ToDiff(N-(c1*a^(n+2)+c2*b^(n+2))
# ,{1,a,b},[1],n,N,x,D);
#
# Example 2:
# A := [seq(C2Rec(n,[[[4,-4],[1,4]]],[1]),n=0..20)]:
# C2ToDiff(N-(n+1)*2^n,{1,2},[1],n,N,x,D);
#
# Example 3:
# A := [seq(C2Rec(n,[[[1],[1]],[[2],[1]]]
# ,[1,1]),n=0..20)];
# C2ToDiff(N^2-N-2^n,{1,2},[1,1],n,N,x,D);

# Example 4:
# A := [seq(C2Rec(n,[[[1,1],[1,1]],[[1,1],[0,1]]]
# ,[1,1]),n=0..20)]:
# C2ToDiff(N^2-(c1*a^(n+1)+c2*b^(n+1))*N
# -(c1*a^n+c2*b^n),{1,a,b},[1,1],n,N,x,D);

# Example 5:
# A := [seq(C2Rec(n,[[[2],[1]],[[1,1],[1,1]]
# ,[[1,1],[0,1]]],[1,1,2]),n=0..30)]:
# C2ToDiff(N^3-2^n*N^2-(c1*a^(n+1)+c2*b^(n+1))*N
# -(c1*a^n+c2*b^n),{1,2,a,b},[1,1,2],n,N,x,D);
#
# Check
# F := x->add(A[i]*x^(i-1),i=1..nops(A)):
# C:= 1/sqrt(5): B:= [(1+sqrt(5))/2,(1-sqrt(5))/2];
# T:=F(x)-1/4*x*F(2*x)-C*x^3*F(B[1]*x)
# -C*x^2*F(B[1]*x)+C*x^3*F(B[2]*x)+C*x^2*F(B[2]*x):
# sort(simplify(T));

C2ToDiff := proc(R,S,A,n,N,x,D) option remember;
local j,t,s,alp,m,r,k,V,W,b,T,B,X;

V := expand(R);
V := add(coeff(V,alp^n)*alp^n,alp in S minus {1});
W := expand(R-V);

r := max(degree(V,N),degree(W,N));
k := max(seq(degree(coeff(V,alp^n),n)
,alp in S minus {1}),degree(W,n));

if nops(A) < r then
	ERROR("Not enough info.");
fi:	

for s from 0 to k do
for t from 0 to r do
for alp in S do
	if alp <> 1 then
		b[s,t,alp] := 
		coeff(coeff(coeff(V,alp^n),N,t),n,s);
	else
		b[s,t,1] := 
		coeff(coeff(W,N,t),n,s);
	fi:
od: od: od:

T := add(add(add( b[s,t,alp]*alp^(-t)*add(  
SolveC(s,t)[j+1]*x^(j+r-t)*D^j*f(alp*x)
,j=0..s),s=0..k),t=0..r),alp in S);

B := add(add(add( b[s,t,alp]*add(  
A[m+1]*(m-t)^s*alp^(m-t)*x^(m+r-t)
,m=0..t-1),s=0..k),t=0..r),alp in S);

print(add(coeff(T,f(alp*x),1)*f(alp*x),alp in S)=B);

X := [op(S)];
[X,[seq(add(add( b[s,t,alp]*alp^(-t)*add(  
SolveC(s,t)[j+1]*x^(j+r-t)*D^j
,j=0..s),s=0..k),t=0..r),alp in X)],B];

end:


###############################################

# Calculation of the homogeneous 
# diff. eq. according to corollary 18

# Input: the C^2-finite relation of a_n,
# set of roots alphas, the list of 
# actual sequence A, symbolic n,N,x,D
# Output: the homogeneous diff. eq. 
# with polynomial coeff. that 
# satisfies the generating function 
# F(x) = sum a_n*x^n. 

# Try: 
# Example 1:
# A := [seq(C2Rec(n,[[[1,1],[1,2]]],[1]),n=0..20)]:
# C2ToDiffHom(N-(c1*a^(n+2)+c2*b^(n+2)),{1,a,b}
# ,[1],n,N,x,D);
#
# Example 2:
# A := [seq(C2Rec(n,[[[4,-4],[1,4]]],[1]),n=0..20)]:
# C2ToDiffHom(N-(n+1)*2^n,{1,2},[1],n,N,x,D);
#
# Example 3:
# A := [seq(C2Rec(n,[[[1],[1]],[[2],[1]]]
# ,[1,1]),n=0..20)]:
# C2ToDiffHom(N^2-N-2^n,{1,2},[1,1],n,N,x,D);
#
# Example 4:
# A := [seq(C2Rec(n,[[[1,1],[1,1]],[[1,1],[0,1]]]
# ,[1,1]),n=0..20)]:
# C2ToDiffHom(N^2-(c1*a^(n+1)+c2*b^(n+1))*N
# -(c1*a^n+c2*b^n),{1,a,b},[1,1],n,N,x,D);

# Example 5:
# A := [seq(C2Rec(n,[[[2],[1]],[[1,1],[1,1]]
# ,[[1,1],[0,1]]],[1,1,2]),n=0..30)]:
# C2ToDiffHom(N^3-2^n*N^2-(c1*a^(n+1)+c2*b^(n+1))*N
# -(c1*a^n+c2*b^n),{1,2,a,b},[1,1,2],n,N,x,D);

C2ToDiffHom := proc(R,S,A,n,N,x,D) option remember;
local i,j,k,M,L,T;

M := C2ToDiff(R,S,A,n,N,x,D);

if M[3]=0 then
	print( add( M[2][i]*f(M[1][i]*x) 
	, i=1..nops(M[2]))=M[3] );
	return(M);
fi:
	
L := M;	
for i from 1 to nops(L[2]) do
	L[2][i] := expand(diff(M[2][i],x)+M[2][i]*D);
od:
L[3] := diff(M[3],x);

T := [M[1], 
[seq( factor(L[3]*M[2][i]-M[3]*L[2][i])
,i=1..nops(M[2]))], 0];

print( add( add(
coeff(T[2][i],D,j)*D^j,
j=0..degree(M[2][i])+1)
*f(T[1][i]*x)
,i=1..nops(T[2]))=0 );

return(T);
end:


# Theorem 19 in the paper
# Convert diff eq. of f(x)
# to C^2-finite recurrence of a_n

# Input: List E for coeff. of 
# each f(alp*x), list of the roots 
# (alphas) S, symbolic x,D,n,a
# Ouput: C^2-finite recurrence
# in term of n and N. 

# Try: 
# Example 1: N-(c1*a^(n+2)+c2*b^(n+2))
# E :=[D,-c1*a^2-c1*a^2*x*D,-c2*b^2-c2*b^2*x*D];
# DiffToC2(E,[1,a,b],x,D,n,N);

# Example 2: N-(n+1)*2^n
# DiffToC2([D,-1-3*x*D-x^2*D^2],[1,2],x,D,n,N);

# Example 3: N^2-N-2^n
# DiffToC2([-1-x*D+D,-2*x-x^2*D],[1,2],x,D,n,N);

# Example 4: N^2-(c1*a^(n+1)+c2*b^(n+1))*N
# -(c1*a^n+c2*b^n)
# E := [D^2,-2*c1-4*D*c1*x-2*D*c1-D^2*c1*x^2-D^2*c1*x,
# -2*c2-4*D*c2*x-2*D*c2-D^2*c2*x^2-D^2*c2*x]:
# DiffToC2(E,[1,a,b],x,D,n,N);

# Example 5: N^3-2^n*N^2-(c1*a^(n+1)+c2*b^(n+1))*N
# -(c1*a^n+c2*b^n)
# E := [D^3,-3/4*D^2-1/4*x*D^3, 
# -6*c1-18*D*c1*x-6*D*c1-9*D^2*c1*x^2
# -6*D^2*c1*x-D^3*c1*x^3-D^3*c1*x^2, 
# -6*c2-18*D*c2*x-6*D*c2-9*D^2*c2*x^2
# -6*D^2*c2*x-D^3*c2*x^3-D^3*c2*x^2]:
# DiffToC2(E,[1,2,a,b],x,D,n,N);

DiffToC2 := proc(E,S,x,D,n,N) option remember;
local j,s,t,r,k;

r := max(seq(degree(E[j],D),j=1..nops(E)));
k := max(seq(degree(E[j],x),j=1..nops(E)));

factor(add(add(add(
coeff(coeff(E[j],x,s),D,t)
*FallFac(n+t-s,t)*S[j]^(n+t-s)
*N^(t-s)  
,s=0..k),t=0..r),j=1..nops(E))) ;

end:



###################################
# Section 3: 5.Closure properties

# Example 2 in the paper
# Test the recurrence of c_n and d_n
# Input: Null
# Output: print the values from 
# recurrences of c_n and d_n.

# Try:
# HoAdd(N-F(n+2),N^2-N-2^n,n,N,c);
# HoTermWise(N-F(n+2),N^2-N-2^n,n,N,c); 
# C2Close();
C2Close := proc() option remember;
local i,n,c,d,F,A,B,C,D;

F := n -> CRec(n,[1,1],[0,1]);
A := n -> C2Rec(n,[[[1,1],[1,2]]],[1]);
B := n -> C2Rec(n,[[[1],[1]],[[2],[1]]],[1,1]);
C := n -> A(n)+B(n);
D := n -> A(n)*B(n);

c[0] := n-> 2^n*F(n+2)*(F(n+4)*F(n+3)-F(n+3)-2^(n+1));
c[1] := n -> F(n+4)*F(n+3)*F(n+2)+2^(2*n+1)
-2^(n+1)*F(n+3)*F(n+2)-F(n+3)*F(n+2);
c[2] := n -> 2^(n+1)*F(n+2)+F(n+2)
-F(n+4)*F(n+3)*F(n+2)+2^n;
c[3] := n -> F(n+3)*F(n+2)-F(n+2)-2^n;

print([seq(add(c[i](n)*C(n+i),i=0..3),n=0..30)]);

d[0] := n -> -F(n+2)*F(n+3)*2^n;
d[1] := n -> -F(n+3);
d[2] := n -> 1;
print([seq(add(d[i](n)*D(n+i),i=0..2),n=0..30)]);

return();
end: