# OK to post homework
# Aurora Hiveley, 2/24/25, Assignment 10

Help:=proc(): print(`CompsMod(n,a,S), OddParStupid(n), DistinctParStupid(n)`):   end:



## Problem 1
# Using procedure CompsG(n,A)

# CompsMod(n,a,S): inputs a positive integer n, a positive integer a and a subset S of {0,1,..., a-1} 
# and outputs the set of compositions of n whose entries mod a belong to S. 
# For example CompsMod(n,1,{0}) is the same as Comps(n), 
# and Comps(n,2,{1}) gives the set of compositions into odd parts.

# [Hint: the set of integers of the form s (mod a) is {seq(a*i+s,i=0..infinity)} 
# but you can safely replace it by {seq(a*i+s,i=0..n)} (and even {seq(a*i+s,i=0..trunc(n/a))} (why?)) ]

CompsMod := proc(n,a,S) local A,i,s:
A := {  seq( seq(a*i+s,i=0..trunc(n/a)), s in S)  } minus {0}:  # 0 causes recursion error 
CompsG(n,A):
end:



## Problem 2
# Which of the following sequences are in the OEIS, if there are what are their A numbers. what is the desription?

# [seq( nops(CompsMod(n,2,{1})),n=1..10)];
# [1, 1, 2, 3, 5, 8, 13, 21, 34, 55] --> A000045, the Fibonacci numbers!
# [seq( nops(CompsMod(n,3,{1})),n=1..10)];
# [1, 1, 1, 2, 3, 4, 6, 9, 13, 19] --> A000930, Narayana's cow sequence. a(n) = a(n-1) + a(n-3) with a(0) = a(1) = a(2) = 1
# [seq( nops(CompsMod(n,5,{1,4})),n=1..10)];
# [1, 1, 1, 2, 3, 5, 7, 10, 15, 23] --> A116975, Number of comps of n into parts of sizes = 1 mod 5 or 4 mod 5



## Problem 3
# A partition is odd if all its members are odd. 
# OddParStupid(n): looks at all the members of Par(n) and selects those that only have odd parts.

OddParStupid := proc(n) local A,B,a,k: 
A := Par(n):    # set of all partitions
B := {}:        # sets with odd parts, which we want to keep

for a in A do 
    if convert(a,set) subset {seq(1+2*k, k=0..trunc((n-1)/2))} then 
        B := B union {a}:
    fi:
od:

B:
end:



## Problem 4
# A partition is distinct if all its members are distinct.
# DistinctParStupid(n): looks at all the members of Par(n) and selects those that have distinct parts.

DistinctParStupid := proc(n) local A,B,a:
A := Par(n):    # set of all partitions
B := {}:        # sets with odd parts, which we want to keep

for a in A do 
    if nops(convert(a,set))=nops(a) then 
        B := B union {a}:
    fi:
od:

B:
end:



## Problem 5
# Verify that for n up to 15, nops(DistinctParStupid(n))=nops(OddParStupid(n))

# [seq(evalb(nops(DistinctParStupid(n))=nops(OddParStupid(n))), n=1..15)];
# [true, true, true, true, true, true, true, true, true, true, true, true, true, true, true]






### copied from C10.txt 


Help10:=proc(): print(`AD(k), Comps(n), CompsG(n,A), ParS(n), Park(n,k), Par(n) `): end:

HelpHW8:=proc(): print(`xnk(n,k), xnkC(n,k), xnkP(n,k,p)`):end:


#xnk(n,k): The solution of the recurrence  xn(n)=(1+add(xn(i)^2,i=0..n-1))/n:
#From Austin DeCicco and others
xnk:=proc(n,k) local i: option remember:
if n=0 then
 1:
 else
 (1+add(xnk(i,k)^k,i=0..n-1))/n:
fi:
end:

#xnkC(n,k): clever version of xn(n) only having to remember what happened yesterday
#From Austin DeCicco and others
# xnk(n,k)= (1+xnk(0,k)^k+...+xnk(n-1,k)^k)/n
#         = (1+xnk(0,k)^k+...+xnk(n-2,k)^k)/n + (xnk(n-1,k)^k)/n
#         = (n-1)/n* (1+xnk(0,k)^k+...+xnk(n-2,k)^k)/(n-1) + (xnk(n-1,k)^k)/n
#         = (n-1)/n* xnk(n-1,k) + (xnk(n-1,k)^k)/n

xnkC:=proc(n,k) option remember:
if n=1 then
2:
else
(xnkC(n-1,k)^(k-1)+n-1)*xnkC(n-1,k)/n:
fi:
end:

#From Austin DeCicco and others
#xnkP(n,k,p): For a prime p, a fast way to compute xnk(n,k) mod p for n<p, as long as it is an integer
xnkP:=proc(n,k,p) option remember:
if n=1 then
2:
else
(xnkP(n-1,k,p)^(k-1)+n-1)*xnkP(n-1,k,p)*n^(-1) mod p:
fi:
end:





#AD(k): Inspired by Austin DeCicco: the first term of the sequence xnk(n,K) that is NOT an integer
AD:=proc(k) local i,ip:
ip:=ithprime:
for i from 1 to 1000 while (xnkP(ip(i) - 1,k,ip(i))^(k-1)+ip(i) - 1)*xnkP(ip(i) - 1,k,ip(i)) mod ip(i)=0 do od:
ip(i):
end:


##Comps(n): The set of compositions of n. For example,
#Comps(3) is {[3],[1,2],[2,1],[1,1,1]}.
Comps:=proc(n) local S,i,S1,s1:
option remember:

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

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

S:={}:

for i from 1 to n do
 S1:=Comps(n-i):
 S:=S union {seq([i,op(s1)],s1 in S1)}:
od:
S:
 end:


#CompsG(n,A): The set of compositions of n whose entries are from the set A/
CompsG:=proc(n,A) local S,i,S1,s1:
option remember:

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

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

S:={}:

for i in A do
 S1:=CompsG(n-i,A):
 S:=S union {seq([i,op(s1)],s1 in S1)}:
od:
S:
 end:

#ParS(n): The set of partitions of n, done stupidly, by first finding Comps(n) and then sorting them
ParS:=proc(n) local S,s:
S:=Comps(n):
{seq(sort(s, `>`), s in S)}:
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: