# Jike Liu
# Homework 12

read "C12.txt":

#-----------------------------------------------------------
# 1. A41c(N)
#-----------------------------------------------------------

# A41c1(n): p(n) via Euler's pentagonal recurrence
A41c1:=proc(n) local j,s:
option remember:

if n<0 then
   RETURN(0):
fi:

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

s:=0:

for j from 1 while (3*j^2-j)/2 <= n do
   s:=s+(-1)^(j-1)*A41c1(n-(3*j^2-j)/2):
od:

for j from 1 while (3*j^2+j)/2 <= n do
   s:=s+(-1)^(j-1)*A41c1(n-(3*j^2+j)/2):
od:

s:
end:


# A41c(N): [p(1),p(2),...,p(N)]
A41c:=proc(N) local n:
[seq(A41c1(n),n=1..N)]:
end:

# Comparison with A41ok(1000):
# A41c(1000) is much faster than A41ok(1000)

#-----------------------------------------------------------
# 2. The Bressoud-Zeilberger proof
#-----------------------------------------------------------

# The identity proved in the gem is

#    Sum(p(n-a(j)), j even) = Sum(p(n-a(j)), j odd),

# where

#    a(j)=(3*j^2+j)/2.

# To recover Euler's usual recurrence for p(n), isolate the j=0 term:

#    p(n)=Sum((-1)^(j-1)*p(n-(3*j^2-j)/2),j=1..infinity)
#        +Sum((-1)^(j-1)*p(n-(3*j^2+j)/2),j=1..infinity),

# with the convention p(m)=0 for m<0.

# The map in the gem works on pairs [j,lambda], where lambda is a partition
# of n-a(j). If lambda=(lambda(1),...,lambda(t)), then

# - if t+3*j >= lambda(1), send [j,lambda] to
#   [j-1,[t+3*j-1,lambda(1)-1,...,lambda(t)-1]];

# - if t+3*j < lambda(1), send [j,lambda] to
#   [j+1,[lambda(2)+1,...,lambda(t)+1,1$(lambda(1)-3*j-t-1)]].

# This is the involution from the paper.

# We use Par from class, but fix the n=0 case.
Par0:=proc(n)
if n=0 then
   RETURN({[]}):
else
   RETURN(Par(n)):
fi:
end:


a:=proc(j)
(3*j^2+j)/2:
end:


# BZphi([j,lambda]): the involution in the Bressoud-Zeilberger
BZphi:=proc(P) local j,lam,t,lam1,i:
j:=P[1]:
lam:=P[2]:
t:=nops(lam):

if t=0 then
   lam1:=0:
else
   lam1:=lam[1]:
fi:

if t+3*j >= lam1 then
   RETURN([j-1,[t+3*j-1,seq(lam[i]-1,i=1..t)]]):
else
   RETURN([j+1,[seq(lam[i]+1,i=2..t),1$(lam1-3*j-t-1)]]):
fi:
end:


# BZset(n): all tagged partitions [j,lambda] with a(j)<=n and lambda in Par(n-a(j))
BZset:=proc(n) local j,S,m,p:
S:={}:
for j from 0 while a(j)<=n do
   m:=n-a(j):
   S:=S union {seq([j,p],p in Par0(m))}:
od:
S:
end:


# BZcheck(n): checks that phi is an involution on level n
BZcheck:=proc(n) local S,P:
S:=BZset(n):
for P in S do
   if not member(BZphi(P),S) then
      RETURN(false):
   fi:
   if BZphi(BZphi(P)) <> P then
      RETURN(false):
   fi:
od:
true:
end:


#-----------------------------------------------------------
# 3. Sylvester's bijection: OddToDis(p) and DisToOdd(d)
#-----------------------------------------------------------

# OddToDis(p): odd partition -> distinct partition
OddToDis:=proc(p) local r,m,a1,q,tail,i:
option remember:

if p=[] then
   RETURN([]):
fi:

r:=0:
while r<nops(p) and p[r+1]>1 do
   r:=r+1:
od:
m:=nops(p)-r:

# base case: p=1^m
if r=0 then
   RETURN([m]):
fi:

# p=(2*a1+1,...,2*ar+1,1^m)
a1:=(p[1]-1)/2:
q:=[seq(p[i]-2,i=2..r)]:
tail:=OddToDis(q):

RETURN([a1+r+m,a1+r-1,op(tail)]):
end:


# DisToOdd(d): distinct partition -> odd partition
DisToOdd:=proc(d) local q,r,a1,m,i:
option remember:

if d=[] then
   RETURN([]):
fi:

# base case: (m) -> 1^m
if nops(d)=1 then
   RETURN([1$d[1]]):
fi:

q:=DisToOdd([op(3..nops(d),d)]):
r:=nops(q)+1:
a1:=d[2]-r+1:
m:=d[1]-d[2]-1:

RETURN([2*a1+1,seq(q[i]+2,i=1..nops(q)),1$m]):
end:



#-----------------------------------------------------------
# 4. Search for pairs (A,B) such that p(A*n+B) mod A = 0
#-----------------------------------------------------------

# CheckPair(L,A,B): given L=[p(1),...,p(N)], test whether p(A*n+B)=0 mod A
# for all n with A*n+B <= N.
CheckPair:=proc(L,A,B) local n,m:
for n from 0 while A*n+B <= nops(L) do
   m:=A*n+B:
   if m=0 then
      if 1 mod A <> 0 then
         RETURN(false):
      fi:
   else
      if L[m] mod A <> 0 then
         RETURN(false):
      fi:
   fi:
od:
true:
end:


SearchPairs:=proc(Amax,N) local L,A,B,S:
L:=A41c(N):
S:=[]:
for A from 2 to Amax do
   for B from 0 to A-1 do
      if CheckPair(L,A,B) then
         S:=[op(S),[A,B]]:
      fi:
   od:
od:
S:
end:


# I searched all pairs with 2<=A<=500 and 0<=B<A, and verified them for all
# terms with A*n+B <= 50000. The pairs that survived were:

# [5,4], [7,5], [11,6], [25,24], [35,19],
# [49,19], [49,33], [49,40], [49,47],
# [55,39], [77,61], [121,116],
# [125,74], [125,99], [125,124],
# [175,124], [245,19], [245,89], [245,194], [245,229],
# [275,149], [385,369].


# End