# OK to post homework
# Lucy Martinez, 03-05-2026, Assignment 12

with(combinat):


# Question 1:
# Write a procedure A41c(N) that inputs and outputs the same things
#  as A41s(N) and AS41ok(N), but using The recurrence
#  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 that p(n)=0 if n is negative (so the above sums are only for j
#  such that (3*j^2-j)/2 and (3*j^2+j)/2, respectively, are <=n
#  [Hint: first write A41c1(n) to compute p(n), and use option remember]
# How does time(A41c(1000)) compare to time(A41ok(1000))? 
# ANSWER:
#  time(A41c(1000)) returns 1.437 s
#  time(A41ok(1000)) returns 11.453 s

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


A41c1:=proc(n) local i,j:
option remember:

if n<=0 then 
 return(0):
fi:

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


for j from 1 to n do
 if (3*j^2-j)/2<=n and (3*j^2+j)/2<=n then
  return(add((-1)^(j-1)*A41c1(n-(3*j^2-j)/2),j=1..n)+ add((-1)^(j-1)*A41c1(n-(3*j^2+j)/2),j=1..n)):
 fi:
od:

end:



#Question 2:
# Read and understand the proof of the above recurrence
#  (due to Euler, based on the Euler Pentagonal Theorem) given in
#  Bressoud-Zeilberger gem: https://sites.math.rutgers.edu/~zeilberg/mamarimY/Zeilberger_y1985_p54.pdf
# Programs it in Maple


#BZphi(n,lambda): Given a positive integer n and a partition lambda,
# applies the function phi as described on page 55 of Bressoud-Zeilberger
BZphi:=proc(n,lambda) local t,m,par,i,j,k:
t:=nops(lambda):
m:=add(par, par in lambda):
j:=findj(n,m):

if t+3*j>=lambda[1] then
 return([t+3*j-1,seq(lambda[i]-1,i=1..t)]):
elif t+3*j<lambda[1] then
 k:=lambda[1]-3*j-t-1:
 return([seq(lambda[i]+1,i=2..t),1$k]):
fi:

end:


#findj(n,m): given integers n and m, finds the integer j
# such that n-aj(j)=m
findj:=proc(n,m) local j:
for j from 0 while n-aj(j)<>m do
od:
j:
end:


#aj(j): returns the integer (3*j^2+j)/2
aj:=proc(j):
return (3*j^2+j)/2:
end:


#Question 3:
# Write procedures OddToDis(p) and DisToOdd(p) that implement
#  Sylvester's bijection described in this gem:
#  https://sites.math.rutgers.edu/~zeilberg/mamarimY/syl84.pdf
# NOTE: I only did the forward direction:

OddToDis:=proc(p) local m,r,a,p1,p2,a1,d1,d2,i:
if nops(p)=0 then
 return(0):
fi:

m:=0:
while nops(p)-m > 0 and p[nops(p)-m] = 1 do
 m:=m+1:
od:
    
# The parts not equal to 1: (2a_1+1, ..., 2a_r+1)
r:=nops(p)-m:

#base case: there are only 1s in the partition
if r=0 then
 return(m):
fi:

#a is the list of the integers a_1,a_2,...,a_r
#2a_i+1 means a_i = (p[i]-1)/2
a:=[seq((p[i]-1)/2, i=1..r)]:
a1:=a[1]:

p1:=[seq(2*a[i]-1,i=2..r)]:
p2:=OddToDis(p1):

d1:=a1+r+m:
d2:=a1+r-1:

if p2=0 then
  return([d1,d2]):
else
  return([d1,d2,op(p2)]):
fi:
end:




#Question 4: [Optional challenge, 5 dollars, for each pair]
# Using A41c(N), try to find pairs of integers (A,B) such that for all n
#  p(A*n+B) mod A=0
#  How far can you verify it? 
# ANSWER: It looks like the pair [5,5] satisfies the property for 1<=n<=3,000



#finpair(N,K): tries to find pairs of integers (A,B) such that for all n
# p(A*n+B) mod A=0, it first constructs the cross product 
# of {1,...,N} with {1,...,N}, and then it takes all pairs from the cross product
# and tests the desired property for all positive integers from 1 to K
findpair:=proc(N,K) local i,s,n,S1,S2,S,S3,T,allPass:
S1:={seq(i,i=1..N)}:
S2:={seq(i,i=1..N)}:
S:=CP(S1,S2):
S3:={}:

for i from 1 to nops(S) do
  s:=S[i]:
  T:={seq(n, n=1..K)}:
  allPass:=true:
  for n in T do
    if A41c1(s[1]*n + s[2]) mod s[1] <> 0 then
      allPass:=false:
      break:
    fi:
  od:
  if allPass then
    S3:=S3 union {s}:
  fi:
od:
S3:
end:


#From a previous class:
#CP(A,B): the Cartesian product of two sets (possibly given as lists)
CP:=proc(A,B) local a,b:
{seq(seq([op(a),op(b)],a in A),b in B)}:
end:


##################################From previous classes:
# C12.txt, March 02, 2026
with(combinat):

Help12:=proc():
 print(`Park(n,k), Par(n), A41s(N), A41ok(N), EulerM(N,q), ConjP(N), EPTtest(N)`):
end:

#Park(n,k): The set of partitions of n into exactly k parts
Park:=proc(n,k) local S,k1,S1,s1:
option remember:

if n<k then
   RETURN({}):
fi:

if k=0 then
 if n=0 then
   RETURN({[]}):
   else
    RETURN({}):
fi:
fi:

S:={}:
for k1 from 0 to k do
 S1:=Park(n-k,k1):
 S:= S union {seq([op(s1+[1$k1]),1$(k-k1)], s1 in S1)}:
od:
S:
end:


#Par(n): The set of partitions of n, the clever way.
Par:=proc(n) local k:
{seq(op(Park(n,k)), k=1..n)}:
end:


#A41s(N): OEIS sequence A41 done stupidly
# (construct the actual set) and take the cardinality
A41s:=proc(N) local n:
[seq(nops(Par(n)),n=1..N)]:
end:


#A41ok(N): OEIS sequence A41 done much more efficiently
# using Euler's generating function
# Sum(p(n)*q^n,n=0..infinity)=Prod(1/(1-q^k), k=1..infinity)
A41ok:=proc(N) local k,f,q:
f:=mul(1/(1-q^k),k=1..N):
f:=taylor(f,q=0,N+4):
[seq(coeff(f,q,k),k=1..N)]:
end:


#EulerM(N,q): The truncation to 1^N of Prod(1-q^k, k=1..infinity)
EulerM:=proc(N,q) local f,k:
f:=mul(1-q^k,k=1..N):
f:=taylor(f,q=0,N+1):
end:


#ConjP(N): The conjectured powers of the polynomial in EulerM
ConjP:=proc(N) local j,L:
L:=[]:
for j from 1 to N while (3*j^2-j)/2<=N do
 L:=[op(L),(3*j^2-j)/2]:
od:

for j from 1 to N while (3*j^2+j)/2<=N do
 L:=[op(L),(3*j^2+j)/2]:
od:

sort(L):
end:


#EPTtest(N): testing Euler's pentagonal numbers (3*j^2+j)/2
EPTtest:=proc(N) local q,j,f,g,k:
f:=mul(1-q^k,k=1..N):
g:=add((-1)^j*q^((3*j^2+j)/2), j=1..trunc(2*sqrt(N/3)+1))+add((-1)^j*q^((3*j^2-j)/2), j=1..trunc(2*sqrt(N/3)+1)):
taylor(f-g,q=0,N+1):
end:


#NOTES:

#A000041: a(n) is the number of partitions of n (the partition numbers).
#Let p(n) be the number of partitions of the integer n
#Generating function:
# Sum(p(n)*q^n,n=0..infinity)

#Recall that the gen. function of the compositions of n is
# Sum(c(n)*q^n,n=0..infinity)=q/(1-2*q)
# g.f. q+q^2+...=q/(1-q)

#In general, if A=q/(1-q) then
# 1+A+A^2+...=1/(1-A)

#Suppose we want to find a generating function for p(n)
# then we can make a decision for the following list:

# No 1s, One 1, Two 1s, Three 1s, etc
#  {[], 1, 11, 111, 1111, ..., etc}
#  The weight enumerator with partitions that only use 1s:
#   1+q+q^2+...=1/(1-q)
# No 2s, One 2, Two 2s, Three 2s, etc
#  {[], 2, 22, 222, 2222, ..., etc}
#  The weight enumerator with partitions that only use 2s:
#   1+q^2+q^4+...=1/(1-q^2)

# IN GENERAL: let k be an integer then
# No ks, One k, Two ks, Three ks, etc
#  {[], k, kk, ..., etc}
#  The weight enumerator with partitions that only use 2s:
#   1+q^(k)+q^(2k)+...=1/(1-q^k)

#SO: Sum(p(n)*q^n,n=0..infinity)=Prod(1/(1-q^k), k=1..infinity)

#Euler:
#Sum(p(n)*q^n,n=0..infinity)=Prod(1/(1-q^k), k=1..infinity)
# Euler discovered something neat about the following
# Sum(p(n)*q^n,n=0..infinity)*Prod(1-q^k, k=1..infinity)=1

#Prod(1-q^k, k=1..infinity)=Sum((-1)^j*q^((3*j^2+j)/2), j=-infinity..infinity)
#
#Sum(p(n)*q^n,n=0..infinity)*(1 - q - q^2 + q^5 + q^7 - q^12 - q^15 + q^22 + q^26 - ....)=1
# So,
# (p(0)*q^0+p(1)*q^1+p(2)*q^2+...)*(1 - q - q^2 + q^5 + q^7 - q^12 - q^15 + q^22 + q^26 - ....)=1

# p(n)-p(n-1)-p(n-2)+p(n-5)+p(n-7)-...=0 (n>0)
# Hence,
# p(n)=p(n-1)+p(n-2)-p(n-5)-p(n-7)-...
# p(n)=Sum( (-1)^j*p(n-(3*j^2-j)/2), j= )