# OK to post homework
# Aurora Hiveley, 3/3/26, Assignment 12

Help:=proc(): print(`A41c(N), A41c1(n), OddToDis(p), DisToOdd(p)`):   end:


### Problem 1
# A41c(N): 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)

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


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

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

add((-1)^(j-1)*A41c1(n-(3*j^2-j)/2),j=1..trunc(2*sqrt(n/3)+1) )+ add((-1)^(j-1)*A41c1(n-(3*j^2+j)/2),j=1..trunc(2*sqrt(n/3)+1)):
end:


# How does time(A41c(1000)) compare to time(A41ok(1000))?
# 1.062 versus 3.453 (yikes!!)



### Problem 2
# Read and understand the proof of the above recurrence (due to Euler, based on the Euler Pentagonal Theorem) given in
# Bressoud-Zeilberger gem. Program it in Maple

Phi := proc(n,lambda) local i,j,t,m,S:
t := nops(lambda):
m := add(i, i in lambda):     # m = n - aj

## brute force this 
S := [seq((3*i^2 + i)/2, i=1..trunc(2*sqrt(n/3)+1))]:
member(n-m, S, 'j'):      # aj = the jth pentagonal number

if n-m < 0 then 
    j := 0:
fi:

if j = 'j' then 
    RETURN(FAIL):
fi:

## do the mapping 
if lambda[1] <= (t+3*j) then 
    [t + 3*j -1, op(lambda - [1$t]) ]:
else
    [op( lambda[2..-1] + [1$(t-1)] ), 1$(lambda[1]-3*j-t-1) ]:
fi:

end:


### Problem 3
# OddToDis(p) and DisToOdd(p): implement Sylvester's bijection 

OddToDis := proc(p) local a1,d1,d2,m,r,i:
option remember:

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

# initialize T(1^m) = m
member(1,p,'i'):

if i = 'i' then                 # no 1's appear
    m := 0:
    r := nops(p):

else 
    m := nops(p) - i + 1:
    r := i-1:

    if m = nops(p) then         # p = [1,1,...,1]
        RETURN([m]):
    fi:
fi:

a1 := (p[1]-1)/2:
d1 := a1 + r + m:     
d2 := a1 + r - 1:     # a2 + r - 1

[d1, d2, op( OddToDis(p[2..r]) )]:
end:


DisToOdd := proc(p) local r,m,a1:
option remember:

if p = [] then 
    RETURN([]):
elif nops(p) = 1 then 
    RETURN([1$op(p)]):
fi:

# r-1 = s-3 (s = nops(p))
r := nops(p) - 2:

a1 := p[2] - r + 1:
m := p[1] - p[2] - 1:

[2*a1-1, op(DisToOdd(p[3..-1])),1$m]:
end:


## for testing, read homework 10, which has functions OddParStupid(n) and DistinctParStupid(n)
# read `hw10AuroraHiveley.txt`:


### Problem 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?


# A41c(N) = [seq( p(n) , n=1..N)], so we are looking for a subsequence of L := A41c which satisfies L[A*n + B] mod A = 0 for all n
# can check with seq(L[A*n + B] mod A, n=1..N); (for some sufficiently large N)

# technically, setting A = 1 and B = Z for any integer Z is a solution, since p(n) is always an integer and any integer mod 1 is 0  
# but i will spare dr. z from owing me infinity many dollars 

# L := A41c(1000): # just to be safe, do a lot of terms 

# for A from 2 to 15 do 
#     for B from 1 to 15 do 
#         if [seq(L[A*n + B] mod A, n=1..50)] = [0$50] then 
#             print([A,B]);
#         fi:
#     od: 
# od: 

# results: [5, 4], [5, 9], [5, 14], [7, 5], [7, 12], [11, 6]







### copied from C12.txt 

# 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= )