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

with(combinat,fibonacci):

# First version: Dec 27, 2016
# Second version: Jan 5, 2020

##########################
# Functions:
#
# Section 1: X-Recursive
# T1(n), Tal(n), 
# M0(n), TestS(n),
# M1(n), M2(n), 
# M3(n), M4(n),
# 
# Section 2: Closure properties
# Add1(ordA,ordB,n), 
# TestAdd1(n), TestAdd2(n),
# Hada(ordA,ordB,n),
# TestHada1(n), TestHada2(n),
# Sum1(ord,n), 
# TestSum1(n), TestSum2(n),
# MyA(i,ord,C,a,n),
#
# Section 3: Guess
# GuessCX(A,r1,r2),
# GSolve(eq,var),
# Yi(n), 
#
# Section 4: Non-holonomic seq.
# Somos4(n), Somos5(n), 
# Bell(n), Tan(N), Tan1(n), 
# C1(n), Ber(N), Ber1(n),
# GraphD(n), 
#
##########################

# Section 1.1: Examples:
# Try: 
# T1 is C-finte, T1(n)=2*T1(n-1).
# A:=[seq(T1(n),n=0..20)];
# GuessC(A,N);
T1 := proc(n) option remember;
local i;
if n = 0 then
	return(1);
fi:	
	
add(T1(i),i=0..n-1);
end:


# Test the sum recurrence
# These sequences 
# (except the last one)
# are holonomic. All has formulas
# Try: 
# A:=[seq(Tal(n),n=1..20)];
# Guess1D(A,1,5,1,n,N);
Tal := proc(n) option remember;
local k;

if n=0 then
	return(1);
fi:	

# add((k+1)*Tal(k),k=0..n-1);
# add((k+1)^2*Tal(k),k=0..n-1);
# n*add(Tal(k),k=0..n-1);
# n^2*add(Tal(k),k=0..n-1);
# add(Tal(k)*binomial(n-1,k),k=0..n-1);
end:


# New sequence: a(n)=F_n*a(n-1);
# Not holonomic, but has non-linear rec.
# Try: 
# A:= [seq(M0(n),n=1..10)];
# Guess1D(A,1,1,1,n,N);
# GuessCX(A,1,2);

M0 := proc(n) option remember;
if n =0 then
	return(1);
fi:

fibonacci(n)*M0(n-1);
end:

# Try: [seq(TestS(n),n=1..10)];
TestS := proc(n) option remember;
if n=0 then
	return(1);
elif n=1 then
	return(1);
elif n=2 then
	return(2);
fi:

TestS(n-1)*(TestS(n-2)/TestS(n-3)+TestS(n-1)/TestS(n-2));	
end:

# New sequence: a(n)= a(n-1)+2^n*a(n-2);
# Not holonomic
# Try:
# A := [seq(M1(n),n=1..20)];
# GuessCX(A,2,1);
M1 := proc(n) option remember;
if n =0 or n = 1 then
	return(1);
fi:

M1(n-1)+2^n*M1(n-2);
end:



# a(n)= sum(F_i*a(i));
# Not holonmic. But a(n)=(F_(n-1)+1)*a(n-1)
# Try: 
# A:= [seq(M2(n),n=2..24)];
# GuessCX(A,1,4);

M2 := proc(n) option remember;
local i;
if n = 1 then
	return(1);
fi:	
	
add(fibonacci(i)*M2(i),i=1..n-1);
end:


# New sequence: a(n)= F(n)*a(n-1)+F(n-1)*a(n-2);
# Not holonomic
# Try: 
# A:=[seq(M3(n),n=1..20)];
# GuessCX(A,2,2);

M3 := proc(n) option remember;
if n =0 or n = 1 then
	return(1);
fi:

fibonacci(n)*M3(n-1)+fibonacci(n-1)*M3(n-2);
end:


M4 := proc(n) option remember;
if n =0 or n = 1 then
	return(1);
fi:

(2^n+1)*M4(n-1);
end:


#################################################
# Section 2: Closure properties

# Closure in the addition of a_n and b_n
# Try: Add1(1,1,n);
Add1 := proc(ordA,ordB,n) option remember;
local i,ord,a,b,c,C,D,alp,S,eq,var;

ord := ordA+ordB:

S := add(alp[i]*(MyA(n+i,ordA,C,a,n)+MyA(n+i,ordB,D,b,n))
,i=0..ord);
S := expand(S);  

eq := { 
seq(coeff(S,a(n-i),1)=0,i=0..ordA-1),
seq(coeff(S,b(n-i),1)=0,i=0..ordB-1),
alp[ord]=1};

var := {seq(alp[i],i=0..ord)};
var := solve(eq,var);

0=subs(var,add(alp[i]*c(n+i),i=0..ord));
end:

# Try: 
# E:= [seq(M0(n)+2^(n*(n+1)/2),n=1..10)];
# [seq(TestAdd1(n),n=1..10)];
TestAdd1 := proc(n) option remember;
local C,D;

if n=1 then
	return(3);
elif n=2 then
	return(9);
fi:

C := n-> fibonacci(n);
D := n-> 2^n;

-(C(n)*C(n-1)-D(n)*D(n-1))
/(-C(n-1)+D(n-1))*TestAdd1(n-1)
+C(n-1)*D(n-1)*(-D(n)+C(n))
/(-C(n-1)+D(n-1))*TestAdd1(n-2);

end:

# Try: 
# [seq(M3(n)+M4(n),n=1..10)];
# [seq(TestAdd2(n),n=1..10)];
TestAdd2 := proc(n) option remember;
local C,D;

if n=1 then
	return(2);
elif n=2 then
	return(7);
elif n=3 then
	return(50);
fi:

C[0] := n-> fibonacci(n);
C[1] := n-> fibonacci(n-1);
D[0] := n-> 2^n+1;


+(D[0](n)*D[0](n-1)*D[0](n-2)-C[0](n)*C[1](n-1)
-D[0](n-2)*C[0](n)*C[0](n-1)-D[0](n-2)*C[1](n))
/(-C[0](n-1)*D[0](n-2)+D[0](n-1)*D[0](n-2)-C[1](n-1))
*TestAdd2(n-1)
-(C[0](n-1)*D[0](n)*D[0](n-1)*D[0](n-2)
-C[0](n)*C[0](n-1)*D[0](n-1)*D[0](n-2)-C[1](n)
*D[0](n-1)*D[0](n-2)+C[1](n)*C[1](n-1))
/(-C[0](n-1)*D[0](n-2)+D[0](n-1)*D[0](n-2)-C[1](n-1))
*TestAdd2(n-2)
-D[0](n-2)*C[1](n-1)*(-C[1](n)-D[0](n-1)
*C[0](n)+D[0](n)*D[0](n-1))
/(-C[0](n-1)*D[0](n-2)+D[0](n-1)*D[0](n-2)-C[1](n-1))
*TestAdd2(n-3);

end:


# Closure in the multiplication of a_n and b_n
# Try: Hada(2,1,n);
Hada := proc(ordA,ordB,n) option remember;
local i,j,ord,a,b,c,C,D,alp,S,eq,var;

ord := ordA*ordB:

S := add(alp[i]*(MyA(n+i,ordA,C,a,n)*MyA(n+i,ordB,D,b,n))
,i=0..ord);
S := expand(S); 

eq := { 
seq(seq(coeff(coeff(S,a(n-i),1),b(n-j),1)=0
,i=0..ordA-1),j=0..ordB-1),
alp[ord]=1};

var := {seq(alp[i],i=0..ord)};
var := solve(eq,var);

0=subs(var,add(alp[i]*c(n+i),i=0..ord));
end:


# Try: 
# [seq(M3(n)*M4(n),n=1..14)];
# [seq( TestHada1(n),n=1..14)];
TestHada1 := proc(n) option remember;
local C,D;

if n=1 then
	return(1);
elif n=2 then
	return(10);
fi:

C[0] := n-> fibonacci(n);
C[1] := n-> fibonacci(n-1);
D[0] := n-> 2^n+1;

D[0](n)*C[0](n)*TestHada1(n-1)+
D[0](n)*D[0](n-1)*C[1](n)*TestHada1(n-2);

end:

# Try: 
# [seq(M3(n)*M3(n-1),n=1..14)];
# [seq( TestHada2(n),n=1..14)];
TestHada2 := proc(n) option remember;
local C,D;

if n=1 then
	return(1);
elif n=2 then
	return(2);
elif n=3 then
	return(10);
elif n=4 then
	return(95);	
fi:

C[0] := n-> fibonacci(n);
C[1] := n-> fibonacci(n-1);
D[0] := n-> fibonacci(n-1);
D[1] := n-> fibonacci(n-2);

+(C[0](n)*C[1](n-1)*D[1](n)*D[0](n-2)+C[0](n)*C[1](n-1)
*D[0](n)*D[0](n-1)*D[0](n-2)-C[1](n)*C[0](n-2)*D[0](n)
*D[1](n-1)-C[0](n)*C[0](n-1)*C[0](n-2)*D[0](n)*D[1](n-1))
/(C[1](n-1)*D[0](n-1)*D[0](n-2)-C[0](n-1)*C[0](n-2)*D[1](n-1))
*TestHada2(n-1)
-(-D[0](n-1)*C[1](n-1)*C[1](n)
*D[0](n)*D[1](n-1)+C[0](n)*C[0](n-1)^2*C[0](n-2)*D[1](n)
*D[1](n-1)-C[1](n)*D[1](n)*C[1](n-1)*D[0](n-1)*D[0](n-2)
+C[1](n)*C[0](n-2)*D[1](n)*C[0](n-1)*D[1](n-1)+C[0](n)
*C[1](n-1)*D[1](n)*C[0](n-1)*D[1](n-1)-C[1](n)*D[0](n)
*D[0](n-1)^2*C[1](n-1)*D[0](n-2))
/(C[1](n-1)*D[0](n-1)*D[0](n-2)-C[0](n-1)*C[0](n-2)*D[1](n-1))
*TestHada2(n-2)
-C[1](n-1)*D[1](n-1)*(C[0](n)*C[1](n-1)*D[1](n)*D[0](n-2)
-C[1](n)*C[0](n-2)*D[0](n)*D[1](n-1)-C[1](n)*C[0](n-2)
*D[0](n)*D[0](n-1)*D[0](n-2)+C[0](n)*C[0](n-1)*C[0](n-2)
*D[1](n)*D[0](n-2))
/(C[1](n-1)*D[0](n-1)*D[0](n-2)-C[0](n-1)*C[0](n-2)*D[1](n-1))
*TestHada2(n-3)
+C[1](n-2)*C[1](n-1)*D[1](n-2)*D[1](n-1)
*(-D[0](n-1)*C[1](n)*D[0](n)+C[0](n)*D[1](n)*C[0](n-1))
/(C[1](n-1)*D[0](n-1)*D[0](n-2)-C[0](n-1)*C[0](n-2)*D[1](n-1))
*TestHada2(n-4);

end:



# Closure in the sum a_i 
# Try: Sum1(1,n);
Sum1 := proc(ord,n) option remember;
local i,a,C,S,c;

S := {seq(a(n-i)=c(n-i)-c(n-i-1),i=0..ord)};
S := subs( S ,a(n)-MyA(n,ord,C,a,n-1));
expand(S);
end:


# Try: 
# [seq(add(M0(i),i=1..n),n=2..24)];
# [seq( TestSum1(n),n=2..24)];
TestSum1 := proc(n) option remember;
local C;

if n=1 then
	return(1);
elif n=2 then
	return(2);
fi:

C[0] := n-> fibonacci(n);

(1+C[0](n))*TestSum1(n-1)
-C[0](n)*TestSum1(n-2);
end:


# Try: 
# [seq(add(M1(i),i=1..n),n=1..14)];
# [seq( TestSum2(n),n=1..14)];
TestSum2 := proc(n) option remember;
local C;

if n=1 then
	return(1);
elif n=2 then
	return(6);
elif n=3 then
	return(19);	
fi:

C[0] := n-> 1;
C[1] := n-> 2^n;

(1+C[0](n))*TestSum2(n-1)
-(C[0](n)-C[1](n))*TestSum2(n-2)
-C[1](n)*TestSum2(n-3);
end:

# Try: MyA(n+2,2,C,a,n);
MyA := proc(i,ord,C,a,n) option remember;
local j;

if i-n <= 0 then
	return(a(i));
fi:

add(C[j](i)*MyA(i-j-1,ord,C,a,n),j=0..ord-1);
end:


#################################################
# Section 3: Guess function

# Try: 
# A:=[seq(M0(n),n=1..10)];
# GuessCX(A,1,2);

# Closure
# 1. (add)
# A:= [seq(1+M4(n),n=1..18)]:
# GuessCX(A,2,3);

# 2. (product)
# A:= [seq(n*M4(n),n=1..18)]:
# GuessCX(A,2,1);

# 3. (summation)
# A:=[seq( add(M4(j),j=1..n),n=2..20)];
# t:=time():GuessCX(A,2,2); time()-t;
GuessCX := proc(A,ordA,ordC) option remember;
local m,i,j,c,C,eq,var,top;

top := (ordA+1)*(ordC+1)+3;

# Step 1: Check Input
if type(A,list) = false then
	return("Input must be List");
elif nops(A) <= top   then
	ERROR("NeedMoreInputs");	
fi:


# Step 2: Make eqs and vars
eq := {seq(A[m]=add(C[i,m]*A[m-i]
,i=1..ordA),m=ordA+1..top)}
union {seq(seq(C[i,m]=add(c[i,j]*C[i,m-j],j=1..ordC)
,i=1..ordA),m=ordA+1..top)};

var := {seq(seq(c[i,j],j=1..ordC),i=1..ordA),
seq(seq(C[i,m],m=ordA+1-ordC..top),i=1..ordA)};

# Step 3: Solve for solutions
var := GSolve(eq,var);
if var = op({}) then
	return(NoSolution);
fi:

# Step 4: Output
print(a(n)=add(C[i,n]*a(n-i),i=1..ordA));
print("where");
for i from 1 to ordA do
	print(C[i,n]=  simplify(
	add(subs(var,c[i,j])*C[i,n-j],j=1..ordC)));
od:
print("with initial conditions");
for i from 1 to ordA do
	print(seq(C[i,m]=subs(var,C[i,m])
	,m=ordA+1-ordC..ordA));
od:
return();
end:




# Groebner solve
# Try: 
# eq :={x^2+y^2+z^2=1,x^2+y^2+z^2=2*x,2*x-3*y-z=0}:
# var:={x,y,z}: GSolve(eq,var);
GSolve := proc(eq,var) option remember;
local t,neq,G;

neq := [seq(rhs(t)-lhs(t), t in eq)];
G := Groebner[Basis](neq,plex(op(var)));
solve(G,var);
end:

# Closure of the sum
# Try: 
# A:= [seq(M0(n)+1,n=1..13)];
# [seq(Yi(n),n=1..13)];
Yi := proc(n) option remember;
local M; 
if n = 1 then
	return(2);
elif n=2 then
	return(2);
elif n=3 then
	return(3);
fi:

M := (fibonacci(n)*fibonacci(n-1)-1)*Yi(n-1)
+fibonacci(n-1)*(1-fibonacci(n))*Yi(n-2);	

M/(fibonacci(n-1)-1);
end:


##############################################
# Section 4: non-holonomic seq.

# Try: [seq(Somos4(n),n=1..20)];
Somos4 := proc(n) option remember;
if n=1 or n=2 or n=3 or n=4 then
	return(1);
fi:	

(Somos4(n-1)*Somos4(n-3)
+Somos4(n-2)^2)/Somos4(n-4);
end:

# Try: [seq(Somos5(n),n=1..20)];
Somos5 := proc(n) option remember;
if n=1 or n=2 or n=3 or n=4 or n=5 then
	return(1);
fi:	

(Somos5(n-1)*Somos5(n-4)
+Somos5(n-2)*Somos5(n-3))/Somos5(n-5);
end:



# Bell numbers, most likely not holonomic.
# the exponential generating function
# B'(x) = exp(x)*B(x)
# Alos Bell(n)= sum(k^n/k!,k=0..infinity)/exp(1);
# Try: 
# A := [seq(Bell(n),n=1..10)];
Bell := proc(n) option remember;
local k;

if n=0 then
	return(1);
fi:	

add(binomial(n-1,k)*Bell(k),k=0..n-1);
end:


# Try: Tan(14);
Tan := proc(N) option remember;
local i,x,A;
A:=taylor(tan(x),x=0,N+1);
[seq(coeff(A,x,i),i=0..N)];
end:

# Coeffi of power series of tan(x).
# Try: 
# [seq(Tan1(n),n=0..14)];
# Tan(14);
Tan1 := proc(n) option remember;
local k;

if n=0 then
	return(0);
elif n=1 then
	return(1);
fi:	

add( k*C1(n-k+1)*Tan1(k),k=1..n-1)/(1-n);
end:


# Try: [seq(C1(n),n=0..10)];
C1 := proc(n) option remember;
local l;

if n = 0 then
	0;
elif n = 1 then
	1;
else
	-4/n/(n-1)*C1(n-2);
fi:	
end:

# Try: Ber(14);
Ber := proc(N) option remember;
local i,x,A;
A:=taylor(x/(exp(x)-1),x=0,N+2);
[seq(coeff(A,x,i),i=0..N)];
end:


# Try: 
# [seq(Ber1(n),n=0..14)];
# Ber(14);
Ber1 := proc(n) option remember;
local k;

if n=0 then
	return(1);
fi:	

-add(  Ber1(k)/(n+1-k)!,k=0..n-1);
end:

# From Wilf generating functions.
# It does not seem to have pattern.
# page 87
# Try: A:=[seq(GraphD(n),n=1..10)];
GraphD := proc(n) option remember;
local k;

2^binomial(n,2)
-add( binomial(n-1,k-1)*GraphD(k)*2^binomial(n-k,2) ,k=1..n-1);
end: