#!/usr/local/bin/maple
#-*- maplev -*-

# Nathaniel Shar
# HW 19
# Experimental Mathematics

# It is okay to link to this assignment on the course webpage.

Help := proc(): print(`Adn(d,n), IsBad(d,S), adnS(d,n), NdnS(d,n), Wn(A,FP,B,n), AdnM(d,n), adnM(d,n), NdnM(d,n), WGn(A,FP,B,n,BFP), AdjustBFP(a,B,FP,BFP), SFP(d,n, B), NdnC(d,n), WnX(A, FP, B, n, x), WGnX(A,FP,B,n,BFP, x), adnC(d,n), GuessRec1(L,r), GuessRec(L), WnR(A, FP, B, N)`): end:

# From C19.txt:

#Dedicated to Endre Szemeredi on the occaison of the Abel Prize 2012

with(combinat):
#Let a_d(n) be the largest size of a subset of {1, ...,n}
#avoiding arthmetical pogressions of length d
#(Szemeredi's famous theorem says that a_d(n)/n->0)
#d=3 is it the already deep Roth theorem who proved that
#a_3(n)/n<=C/log(log(n))




#Adn(d,n): Inputs positive integers d and n and outputs
#the set of all subsets of {1, ...,n} that DO NOT have
#an arithmetical progression of length d
Adn:=proc(d,n) local S,s1,i,T:

    S:=powerset(n):
    T:={}:

    for s1 in S do

        if not IsBad(d,s1) then

#T:=T union {s1}:

            T:={op(T), s1}:

        fi:

    od:

    T:


end:




#IsBad(d,S): Does the set of pos. integers  S contain
#an AP of length d?
IsBad:=proc(d,S) local m,S1,dif,PTM,j:
option remember:

    if nops(S)<d then
        RETURN(false):
    fi:

    m:=max(S):

    S1:=S minus {m}:

    if IsBad(d,S1) then
        RETURN(true):
    fi:

    for dif from 1 to trunc(m/(d-1)) do
        PTM:={seq(m-j*dif, j=1..d-1)}:

        if PTM minus S1={} then
            RETURN(true):
        fi:

    od:


    false:
end:

#adnS(d,n): a stupid way of computing a_d(n) above
adnS:=proc(d,n) local Y,y:
    Y:=Adn(d,n):

    max(seq(nops(y),y in Y)):

end:


#NdnS(d,n): a stupid way of computing the number of
#subsets of {1, ...,n} that do not contain an AP
#of length d
NdnS:=proc(d,n) local Y,y:
    Y:=Adn(d,n):

    nops(Y):

end:




#Wn(A,FP,B,n): inputs
#(i) A: a set of symbols (the alphabet)
#(ii) FP: a set of forbidden patterns, given as a list
#where B denotes Blank. For example, if FP contains
#[1,B,1,B,1] then our words may not contain
#10101, 10111, 11101, 11111
#(iii) A letter B
#n a pos. integer
#Outputs:
#The number of words in the alphabet A of length n that
#AVOID all the patterns of FP.
#For example,
#Wn({0,1},{[1,B,1]},B,3); should return
#6
Wn:=proc(A,FP,B,n)
option remember:
    WGn(A,FP,B,n,{});
end:


#




#AdnM(d,n): mediocre version of Adn(d,n)
#Inputs positive integers d and n and outputs
#the set of all subsets of {1, ...,n} that DO NOT have
#an arithmetical progression of length d
AdnM:=proc(d,n) local S,T,t1,s1:
option remember:

    if n=1 then
        RETURN({{},{1}}):
    fi:

    S:=AdnM(d,n-1):

    T:=S:

    for t1 in T do

        s1:=t1 union {n}:

        if not IsBad(d,s1) then

#S:=S union {s1}:

            S:={op(S), s1}:

        fi:

    od:

    S:


end:

#adnM(d,n): a mediocre way of computing a_d(n) above
adnM:=proc(d,n) local Y,y:
    Y:=AdnM(d,n):

    max(seq(nops(y),y in Y)):

end:


#NdnM(d,n): a stupid way of computing the number of
#subsets of {1, ...,n} that do not contain an AP
#of length d
NdnM:=proc(d,n) local Y,y:
    Y:=AdnM(d,n):

    nops(Y):

end:

#WGn(A,FP,B,n,BFP): inputs
#(i) A: a set of symbols (the alphabet)
#(ii) FP: a set of forbidden patterns, given as a list
#where B denotes Blank. For example, if FP contains
#[1,B,1,B,1] then our words may not contain
#10101, 10111, 11101, 11111
#(iii) A letter B
#(iv) n a pos. integer
#(v): Another set of patterns (same format at for FP)
#Outputs:
#The number of words in the alphabet A of length n that
#AVOID all the patterns of FP and IN ADDITION AVOID
#the set of patterns BFP at the beginning
#For example,
#Wn({0,1},{[1,B,1]},B,3); should return
#6
WGn:=proc(A,FP,B,n,BFP) local co,a,BFP1:
option remember:
    if member([],BFP union  FP) then
        RETURN(0):
    fi:

    if n=0 then
        RETURN(1):
    fi:


    co:=0:

    for a in A do

        BFP1:=AdjustBFP(a,B,FP,BFP):

        co:=co+WGn(A,FP,B,n-1,BFP1):

    od:

    co:

end:


#AdjustBFP(a,B,FP,BFP): inputs a letter a, a symbol B for
#blank, a set of forbidden patterns FP, and another set
#of forbidden patterns at the beginning
#outputs the new BFP for the chopped word assuming
#it starts with a
AdjustBFP:=proc(a,B,FP,BFP)  local  NBFP,fp:

    NBFP:={}:


    if member([],BFP) then
        RETURN(FAIL):
    fi:

    for fp in FP union BFP do

        if fp[1]=a or fp[1]=B then
            NBFP:=NBFP union {fp[2..nops(fp)]}:
        fi:

    od:



    NBFP:

end:

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

#############
# Problem 1 #
#############

with(ListTools):

SFP := proc(d,n, B) local i, j, FP:
    if d = 1 then:
        return [1]:
    fi:
    FP := {}:
    for i from 0 to floor((n-d)/(d-1)) do:
        FP := FP union {Flatten([[1, B$i]$(d-1), 1])}:
    od:
    return FP:
end:

NdnC := proc(d,n) local B:
    return Wn({0,1}, SFP(d,n, B), B, n):
end:

#############
# Problem 2 #
#############

WnX := proc(A, FP, B, n, x):
    WGnX(A, FP, B, n, {}, x):
end:

WGnX := proc(A,FP,B,n,BFP, x) local co,a,BFP1:
option remember:
    if member([],BFP union  FP) then
        RETURN(0):
    fi:

    if n=0 then
        RETURN(1):
    fi:


    co:=0:

    for a in A do

        BFP1:=AdjustBFP(a,B,FP,BFP):

        co:=co+x[a]*WGnX(A,FP,B,n-1,BFP1, x):

    od:

    co:

end:

#############
# Problem 3 #
#############

adnC := proc(d,n) local x:
    return degree(expand(WnX({0,1}, SFP(d,n,B), B, n, x)), x[1]):
end:

#############
# Problem 4 #
#############

# Let S = {(T,U)}, where T and U are sets such that every element of T
# and U is a suffix of a pattern in S.  Note that S is a finite
# set. Let a(T, U, n) be the number of words of length n that match a
# pattern in T or match a pattern in U at the beginning.  Observe that
# a(T, U, n) can be written in terms of the {a(T, U, n-1)}. (In fact,
# this is how we wrote the functions above.) Therefore, the a(T, U, n)
# are described by a system of first-order linear difference
# equations. We conclude that they are all C-finite. In particular,
# A(S, {}, n) is C-finite.  The quantity we want is just |A|^n - A(S,
# {}, n), so it is also C-finite since |A|^n is C-finite. So we are
# done.

#############
# Problem 5 #
#############

#############
# From class 6:

#GuessRec1(L,r): inputs a sequence of numbers, L, and a pos.
#integer r, and tries to guess a linear recurrence equation
#of order r satisfied by it
#L(n)+c1*L(n-1)+...+cr*L(n-r)=0 for all n within the range
#The  output would be given  as a list of r numbers
#[[c1,c2,...,cr],InitialConditions]
#For example, the Fibonacci sequence is
#[[-1,-1], [1,1]]
GuessRec1:=proc(L,r) local c,i,var,eq:

    if nops(L)<=2*r+6 then
        RETURN(FAIL):
    fi:


    var:={seq(c[i],i=1..r)}:


    eq:={seq( L[n]+add(c[i]*L[n-i],i=1..r)=0 , n=r+1..nops(L))}:

    var:=solve(eq,var):

    if var=NULL then
        RETURN(FAIL):
    fi:

    [[seq(subs(var, c[i]),i=1..r)], L[1..r]   ] ;

end:


#GuessRec(L): inputs a sequence of numbers, L,
# and tries to guess a linear recurrence equation
#of some order, r, satisfied by it
#L(n)+c1*L(n-1)+...+cr*L(n-r)=0 for all n within the range
#The  output would be given  as a list of r numbers
#[[c1,c2,...,cr],InitialConditions]
#For example, the Fibonacci sequence is
#[[-1,-1], [1,1]]
GuessRec:=proc(L) local r,ans:

    for r from 1 to (nops(L)-6)/2 do
        ans:=GuessRec1(L,r):
        if ans<>FAIL then
            RETURN(ans):
        fi:
    od:

    FAIL:
end:

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


WnR := proc(A, FP, B, N):
    GuessRec([seq(Wn(A, FP, B, n), n=1..N)]):
end:


# Problem 1: [[-1, -1], [2, 3]]

# Problem 2: [[-1, -1], [2, 4]]

# Problem 3: [[-1, -1, -1], [2, 4, 8]]

# Problem 4: [[-26, 0, 0, 0, 2, -26, 676, -17576, 456976, -1, 26, 0,
# 0, 0, -3, 26, -1352, 0, 0, 0, 0, 0, 0, 0, 1], [26, 676, 17576,
# 456976, 11881375, 308915724, 8031808148, 208826994272,
# 5429500937120, 141167000602370, 3670341397830146, 95428860279966824,
# 2481149949625113728, 64509887831250168736, 1677256801278467538269,
# 43608669492556433587452, 1133825215948714586103456,
# 29479450652365844680379856, 766465587941714581179849616,
# 19928101931970411885560383365, 518130563013877042540242072648,
# 13471392370709989493023397440900, 350256142679548497482735104595720,
# 9106658176736827011456605811559488,
# 236773072738946929326737191050813688]]