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

# Nathaniel Shar
# HW 21
# Experimental Mathematics

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

Help := proc(): print(`ParBest(N), ParBestSeq(N), ParSeq(N), ConjRC(N, K)`): end:


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

ParBest := proc(N) local i, total: option remember:
    if 0 = N then: return 1: fi:
    total := 0:
    for i from ceil(evalf((-1/2 - sqrt(1/4 + 6*N))/3)) to floor(evalf((-1/2 + sqrt(1/4 + 6*N))/3)) do:
        if i <> 0 then:
            total := total + (-1)^(i-1)*ParBest(N-abs(i*(3*i+1)/2)):
        fi:
    od:
    return total:
end:

ParBestSeq := proc(N) local i:
    return [seq(ParBest(i), i=1..N)]:
end:

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

# From C21:

ParSeq:=proc(N) local F,q,i:

#F:=(1+x+x^2+x^3+x^4+....)*(1+x^2+x^4+x^6+....)*(1+x^3+x^6+...)
#F=1/(1-x)/(1-x^2)/(1-x^3)/......

F:=mul(1/(1-q^i),i=1..N):

F:=taylor(F,q=0,N+1):

[ seq(coeff(F,q,i),i=1..N)]:

end:

# ParBestSeq seems to run in approximately linear time, with even
# ParBestSeq(10000) taking only 9 seconds.

# ParSeq needs at least quadratic time. I found that ParSeq(800) took
# 9 seconds. So ParSeq(10000) would be approximately 100 times
# slower than ParBestSeq(10000).

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

ConjRC := proc(N, K) local a, m, R, L:
    R := {}:
    L := ParBestSeq(N*K + K):
    for a from 1 to K do:
        for m from a to K do:
            if {seq(L[n*m+a] mod m, n=0..N)} = {0} then:
                R := R union {[a,m]}:
            fi:
        od:
    od:
    return R:
end:

# ConjRC(400, 30) = {[1, 1], [4, 5], [5, 7], [6, 11], [24, 25]}