# hw11JikeLiu.txt
# Math 640 — Homework 11
# Jike Liu
#

read "C10.txt":
read "C11.txt":

# ------------------------------------------------------------
# 1) CountNuCompsClever(a,b,A,B,N)
# ------------------------------------------------------------

# CountNuCompsClever(a,b,A,B,N):
# returns [c1,...,cN] where cn = number of qualifying compositions of n.
CountNuCompsClever := proc(a,b,A,B,N) local x,f;
  f := GenCompsGF(a,b,A,B,x):
  [seq( coeff(taylor(f, x=0, N+3), x, i), i=1..N )]:
end:


# ------------------------------------------------------------
# 2) CountNuCompsStupid1 and CountNuCompsStupid
# ------------------------------------------------------------

# CountNuCompsStupid1(a,b,A,B,n):
# brute force: loop over all compositions L of n and check
# (i) member(nops(L) mod b, B)
# (ii) for all parts, member(L[i] mod a, A)
CountNuCompsStupid1 := proc(a,b,A,B,n) local L, S, ok, i;
  S := {}:

  for L in Comps(n) do

    # length condition
    if not member(nops(L) mod b, B) then
      next:
    fi:

    # residue condition on every part
    ok := true:
    for i from 1 to nops(L) do
      if not member(L[i] mod a, A) then
        ok := false: break:
      fi:
    od:

    if ok then
      S := S union {L}:
    fi:
  od:

  nops(S):
end:

# CountNuCompsStupid(a,b,A,B,N): first N terms
CountNuCompsStupid := proc(a,b,A,B,N) local n;
  [seq(CountNuCompsStupid1(a,b,A,B,n), n=1..N)]:
end:


# Check the two given parameter sets:
# (a,b,A,B) = (3,4,{1,2},{0,3})
evalb( CountNuCompsClever(3,4,{1,2},{0,3},15)
     = CountNuCompsStupid(3,4,{1,2},{0,3},15) );

# (a,b,A,B) = (2,3,{1},{2})
evalb( CountNuCompsClever(2,3,{1},{2},15)
     = CountNuCompsStupid(2,3,{1},{2},15) );


# ------------------------------------------------------------
# 3) Find many (a,b,A,B) whose sequences are in OEIS
# ------------------------------------------------------------

# Below are choices of (a,b,A,B) for CountNuCompsClever(a,b,A,B,20),
# together with the OEIS sequence they match (up to the usual shift since we list n=1..).

# 1) Powers of 2 (all compositions, no restrictions)
# (a,b,A,B) = (1,1,{1},{0})
# Output starts: 1,2,4,8,16,32,64,128,256,512,...
# OEIS: A000079 (powers of 2), shifted (our term is 2^(n-1)).

# 2) All ones (only part size 1 effectively allowed for n<=20)
# (a,b,A,B) = (1000,1,{1},{0})   # any a>20 works for first 20 terms
# Output starts: 1,1,1,1,1,1,1,1,1,1,...
# OEIS: A000012 (all 1's).

# 3) Fibonacci via odd parts
# (a,b,A,B) = (2,1,{1},{0})      # parts ≡1 (mod 2), i.e. odd parts
# Output starts: 1,1,2,3,5,8,13,21,34,55,89,144,...
# OEIS: A000045 (Fibonacci), shifted.

# 4) Fibonacci via parts {1,2}
# (a,b,A,B) = (21,1,{1,2},{0})   # any a>20 works for first 20 terms
# Output starts: 1,2,3,5,8,13,21,34,55,89,144,233,...
# OEIS: A000045 (Fibonacci), shifted.

# 5) Tribonacci via parts {1,2,3}
# (a,b,A,B) = (21,1,{1,2,3},{0}) # any a>20 works for first 20 terms
# Output starts: 1,2,4,7,13,24,44,81,149,274,504,927,...
# OEIS: A000073 (Tribonacci), shifted.

# 6) Tetranacci via parts {1,2,3,4}
# (a,b,A,B) = (21,1,{1,2,3,4},{0}) # any a>20 works for first 20 terms
# Output starts: 1,2,4,8,15,29,56,108,208,401,773,1490,...
# OEIS: A000078 (Tetranacci), shifted.

# 7) Parts not divisible by 3 (residues 1 or 2 mod 3)
# (a,b,A,B) = (3,1,{1,2},{0})
# Output starts: 1,2,3,6,11,20,37,68,125,230,423,778,...
# OEIS: A001590 (a tribonacci-type sequence), shifted.

# 8) Parts ≡ 1 (mod 3)  (Narayana's cows)
# (a,b,A,B) = (3,1,{1},{0})
# Output starts: 1,1,1,2,3,4,6,9,13,19,28,41,...
# OEIS: A000930 (Narayana's cows).

# ------------------------------------------------------------
# 4) Dice pips with same sum distribution as standard dice (Sicherman dice)
# ------------------------------------------------------------

# NOTE: CP in C11 returns a SET, so it loses multiplicities.
# For dice distributions we need multiplicities, so we use a list-based CPm.

CPm := proc(A,B) local a,b;
  [seq(seq([a,b], a in A), b in B)]:
end:

SumDist := proc(D1,D2,x);
  expand(WtE(CPm(D1,D2), x)):
end:

# Standard dice and Sicherman dice (nonstandard)
Dstd := [1,2,3,4,5,6]:
D1   := [1,2,2,3,3,4]:
D2   := [1,3,4,5,6,8]:

Pstd := SumDist(Dstd, Dstd, x):
Psic := SumDist(D1, D2, x):

# Check they match
expand(Psic - Pstd);   # output 0

# Print the counts for sums 2..12 (output: 1,2,3,4,5,6,5,4,3,2,1)
[seq(coeff(Psic,x,s), s=2..12)];


# End