# restart; read `c:/Users/Asus/Google Drive/Guess.txt`; with(SolveTools,Linear): with(combinat,powerset): # Report problem on Guess: ##################################### # Guessing Functions: # GuessPol, GuessPol2D, # GuessRat, # # GuessC, GuessGenRat, # SolveCFinite, # # Guess1D, Guess2D, # # AlgPoly, AlgHyper ##################################### # Try: A:= [seq(n^2+3*n,n=3..20)]; # GuessPol(A,3,x); 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 polynomial of 2 variables. # Try: # A := [seq([seq(i+j^2-2,i=2..5)],j=1..6)]; # GuessPol2D(A,[2,1],[x,y]); GuessPol2D := proc(L,start,mu) option remember; local bound,s,l1,i1,j1,i2,j2,eq,var,C,deg; deg := -1; var := op({}); l1 := nops(L); bound := floor(sqrt(add(nops(s),s in L)))-2; while var = op({}) and deg <= bound do deg := deg+1; eq := {seq(seq( add(add(C[j1,j2]*(i1+start[2]-1)^j1 *(i2+start[1]-1)^j2 ,j2=0..deg),j1=0..deg)=L[i1,i2] ,i2=1..nops(L[i1])),i1=1..nops(L))}; var := {seq(seq(C[j1,j2] ,j2=0..deg),j1=0..deg)}; var := Linear(eq,var); od: if deg > bound then return("No Solution"); else return(factor( add(add( subs(var,C[j1,j2])*mu[2]^j1*mu[1]^j2 ,j2=0..deg),j1=0..deg) )); fi: end: ################################################# # Guess rational function, f(x)/g(x) # Try: (from Regevid) # A:=[seq(XiTwo(n)/Cata(n),n=3..30)]; # GuessRat(A,2,2,3,n); GuessRat := proc(L,top,bot,start,mu) option remember; local eq,var,i,j,k,T,C,B; if top+bot >= nops(L)-2 then return("NeedMoreInput"); fi: eq := {seq( add(C[j]*i^j,j=0..top) =L[i-start+1]*add(B[j]*i^j,j=0..bot) ,i=start..nops(L)+start-1)}; #,C[top]=1}; var := {seq(C[j],j=0..top),seq(B[j],j=0..bot)}; var := Linear(eq,var); T := {seq(subs(var,C[j]),j=0..top) ,seq(subs(var,B[j]),j=0..bot)}; if T ={0} then return("No Solution"); else return(factor(add(subs(var,C[j])*mu^j,j=0..top)) /factor(add(subs(var,B[j])*mu^j,j=0..bot))); fi: end: ################################################## # Guess C-finite. # 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=1..30)]; # GuessC(A,N); GuessC := proc(A,N) option remember; local 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 # Try: A := [seq(2^i,i=0..20)]; # 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: # Solve the closed form of C-finite # Try: A:=[seq(FastSumMeOdd(2,n),n=1..120)]: # M:=SolveCFinite(A,1,n); # [seq(subs(n=i,M)-A[i],i=1..200)]; SolveCFinite := proc(A,s,n) option rememeber; local i,j,k,N,a,S,T,eq,var,exp; S := [solve(GuessC(A,N))]; if member(0,{op(S)}) then return("0 is part of solution"); fi: T := {seq( [a,Count1(a,S)],a in {op(S)})}; k := 0; exp := 0; for i from 1 to nops(T) do exp := exp + add(a[k+j]*n^(j-1),j=1..T[i][2])*T[i][1]^n; k := k+T[i][2]; od: eq := {seq(expand(subs(n=i-1+s,exp))=A[i],i=1..k+10)}; #print(eq); var := {seq(a[i],i=1..k)}; var:=solve(eq,var); subs(var,exp); end: # Count multiplicity of the entry Count1 := proc(a,S) option remember; local i,count; count := 0; for i from 1 to nops(S) do if a= S[i] then count := count+1; fi: od: count; end: #################################################### # 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 of seq. # Output: recurrence according to Holonomic Ansatz. # Try: # Guess1D([seq(sum(binomial(i,k),k=0..i),i=1..45)] # ,1,1,1,n,N); # Guess1D([seq(floor((n/2)^2),n=1..15)],2,1,1,n,N); # Guess1D([seq(n!+n,n=1..35)],4,1,1,n,N); Guess1D := proc(C,ord,deg,startN,n,x) option remember; local i,j,k,T,e,eq,mv,var,f,d,na,s,shift; if type(C,list) = false then return("Input must be List"); elif (ord+1)*(deg+1)+ord-1 >= nops(C)-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(C)-ord)}; eq := subs({seq(f[k]=C[k] ,k=1..nops(C))},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]=x^i,i=0..ord)},e)); k := ilcm(seq(denom(coeff(shift,x,i)),i=0..ord)); T := factor(add(factor(k*coeff(shift,x,i))*x^i,i=0..ord))=0; # Another choice of output T := expand(lhs(T)); add(factor(coeff(T,x,i))*x^i,i=0..ord)=0; end: ###################################################### # Guess recurrences of two variables. # Input: 2 dimensions list, extra variables, # order of rec.n and k, # max deg of pol, starting points. # Output: recurrence according to Holonomic Ansatz. # Try: # A:=[seq([seq(n+2*j,j=0..10)],n=0..10)]: # Guess2D(A,[[0,0],[1,0]],1,0,0,n,j,N,J); # A:=[seq([seq(floor((n+2*j)/3)^2,j=0..20)],n=0..20)]: # Guess2D(A,[[0,0],[1,0],[0,1],[1,1]],4,0,0,n,j,N,J); Guess2D := proc(C,ORD,deg,startn,startk,n,k,N,K) local i,j,e,eq,mv,var,f,d,maxordn,maxordk, s,t,shift; # Step 1: Check Input maxordn := max(seq(ORD[i][1],i=1..nops(ORD))); maxordk := max(seq(ORD[i][2],i=1..nops(ORD))); #Check input if type(C,list) = false then return("Input must be a List"); elif (nops(ORD)+1)*binomial(deg+2,2) > (nops(C[1])-maxordk)*(nops(C)-maxordn)-2 then return("Not Enough Info", (nops(C[1])-maxordk)*(nops(C)-maxordn)-2); fi: #Step2: Make eqs and variables #Starting solving by making eqs and variables e := add(add(add( d[ORD[i][1]][ORD[i][2]][s][t]*na^s*ka^t ,t=0..deg-s),s=0..deg) *f[n+ORD[i][1]][k+ORD[i][2]] ,i=1..nops(ORD))=0; eq := {seq(seq( subs({na=i-1+startn,ka=j-1+startk,n=i,k=j},e) ,j=1..nops(C[i])-maxordk) ,i=1..nops(C)-maxordn)}; eq := subs({seq(seq(f[i][j]=C[i][j] ,j=1..nops(C[i])),i=1..nops(C))},eq); mv := {seq(seq(seq( d[ORD[i][1]][ORD[i][2]][s][t] ,t=0..deg-s),s=0..deg),i=1..nops(ORD))}; #Step3: Solve for relation # Check input first if nops(eq) <= nops(mv) then ERROR("NeedMoreInputs"); fi: var := Linear(eq,mv); if {seq(rhs(T), T in var)} = {0} then return(NoSolution); fi: #in case there are more than one solution var := [var][1]; e := subs(SubsConst(var,mv,ORD,deg,d),e); #Step4: Put things back for output shift := lhs( subs({ka=k,na=n, seq(seq(f[n+s][k+t]=N^s*K^t,s=0..maxordn), t=0..maxordk)},e)); j := ilcm(seq(seq( denom(coeff(coeff(shift,N,s),K,t)) ,t=0..maxordk),s=0..maxordn)); factor(add(add( factor(j*coeff(coeff(shift,N,s),K,t)) *N^s*K^t,t=0..maxordk),s=0..maxordn))=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: ##################################### # Solve for polynomial/ hypergeometric # solution of the holonomic recurrence # according to the book A=B. ##################################### # Solve for the polynomial solution # with degree less than deg # of the given recurrence # Try: AlgPoly(3*y[2]-n*y[1]+(n-1)*y[0],2,3,y,n); AlgPoly := proc(rec,ord,deg,y,n) option remember; local i,a,M,eq,var; y := n->add(a[i]*n^i,i=0..deg); M := subs({seq(y[i]=y(n+i),i=0..ord)},rec); M := expand(M); eq := {seq(coeff(M,n,i)=0,i=0..deg)}; var := {seq(a[i],i=0..deg)}; var := solve(eq,var); factor(subs(var,y(n))); end: # Solve for the Hyper geormetric term # solution with degree less than deg # of the given recurrence # Output: the possible S(n); # Try: # 1. Example 8.4.1 # AlgHyper((n-1)*y[2]-(n^2+3*n-2)*y[1] # +2*n*(n+1)*y[0],2,3,y,n); # 2. My Example # AlgHyper(y[2]-y[1]-y[0],2,3,y,n); # 3. Example 8.4.2 # AlgHyper((n+2)^3*y[2]-(2*n+3)*(17*n^2+51*n+39)*y[1] # +(n+1)^3*y[0],2,3,y,n); # 4. Example 8.4.3 # AlgHyper(y[2]-y[1]-(n+1)*y[0],2,3,y,n); # 5,6. Mathematica page 158,159 # AlgHyper(4*(n+2)*(2*n+3)*(2*n+5)*y[2] # -12*(2*n+3)*(9*n^2+27*n+22)*y[1] # +81*(n+1)*(3*n+2)*(3*n+4)*y[0],2,3,y,n); # AlgHyper(y[2]-(2*n+1)*y[1]+(n^2-2)*y[0],2,3,y,n); # y(n) = add(binomial(k,n-k)^2,k=0..n); # AlgHyper((n+4)*y[4]-(2*n+7)*y[3]-(n+3)*y[2] # -(2*n+5)*y[1]+(n+2)*y[0],4,3,y,n); AlgHyper := proc(rec,ord,deg,y,n) option remember; local a,b,c,cc,z,zz,A,B,Tc,Z,T,Ans; Ans := []; A := SetFac(coeff(rec,y[0]),n); B := SetFac(coeff(rec,y[ord]),n); for a in A do for b in B do b:= subs(n=n-ord+1,b); T := add( z^i*coeff(rec,y[i])* mul(subs(n=n+j,a),j=0..i-1)* mul(subs(n=n+j,b),j=i..ord-1)* c[i],i=0..ord); T := factor(T/subs(n=n,a)/subs(n=n+ord-1,b)); Tc := subs({seq(c[i]=1,i=0..ord)},T); Z := {solve(coeff(Tc,n,degree(Tc,n))=0,z)} minus {0} ; if Z <> {} then for zz in Z do T := subs(z=zz,T); cc := AlgPoly(T,ord,deg,c,n); if cc <> 0 then Ans := [op(Ans), factor(zz*a/b*subs(n=n+1,cc)/cc)]; fi: od: fi: od: od: Ans: end: # Output: all the linear monic # factors of polynomial P # Try: SetFac( 2*n*(n+1),n); # SetFac(n^2-2,n); SetFac := proc(P,n) option remember; local i,a,b,A; A := factors(P,sqrt(2))[2]; A := [seq( A[i][1]$(A[i][2]),i=1..nops(A))]; {seq( mul(b,b in a) ,a in powerset(A))}; end: #Try: [seq(MS(n),n=1..15)]; MS := proc(n) option remember; if n=0 then return(1); elif n=1 then return(1); elif n=2 then return(2); elif n=3 then return(5); fi: ((2*n-1)*MS(n-1)+(n-1)*MS(n-2)+(2*n-3)*MS(n-3) -(n-2)*MS(n-4))/n; end: