# OK to post

restart:


# 1. Garvan Maple booklet (pp. 91-end). Small sample only.

# exact vs numeric arithmetic
10/26;
evalf(10/26);

# algebra
factor(x^4-1);
expand((x+y)^5);
simplify((x^2-1)/(x-1));

# lists / sets / seq
L := [seq(i^2, i=1..6)];
S := {seq(i^2, i=1..6)};
nops(L), nops(S);

# loops and if/fi, od
T := 0:
for i from 1 to 10 do
  if i mod 2 = 0 then
    T := T + i;
  fi;
od:
T;   # sum of even numbers <= 10

# a procedure (made-up similar example): trapezoid rule
trap := proc(f,a,b,n)
  local i,h,s;
  h := (b-a)/n;
  s := (f(a)+f(b))/2;
  for i from 1 to n-1 do
    s := s + f(a+i*h);
  od;
  return h*s;
end:

f := x -> 1/(1+x^2):
trap(f,0,1,50);
evalf(%);

# package example
with(combinat):
randperm(6);
permute(4)[1..5];

# saving / reading a procedure (file I/O)
save trap, "trap.m":
read "trap.m":

# writing a short output file
writeto("hw4FirstNameLastName_output.txt"):
printf("trap(f,0,1,100) = %a\n", trap(f,0,1,100));
printf("evalf(trap(f,0,1,100)) = %.12f\n", evalf(trap(f,0,1,100)));
writeto(terminal):

# LaTeX printing
expr := (x^2+1)/(x-1):
latex(expr);


#2. Done

# Now read C4.txt for the remaining questions


read("C4.txt"):


# 3. Run 10 times:
# S:={seq(randperm(3),i=1..3)}; [seq(nops(AvoidPer(n,S)),n=1..8)];
# How many are in OEIS? For those not, conjecture a formula.
#

# My recorded 10 outputs (t, S, L):
# 1, {[3,1,2],[3,2,1]},                 [1,2,4,8,16,32,64,128]
# 2, {[1,2,3],[1,3,2],[2,3,1]},         [1,2,3,4,5,6,7,8]
# 3, {[1,2,3],[2,1,3],[2,3,1]},         [1,2,3,4,5,6,7,8]
# 4, {[1,2,3],[1,3,2],[3,2,1]},         [1,2,3,1,0,0,0,0]
# 5, {[1,2,3],[2,3,1],[3,2,1]},         [1,2,3,1,0,0,0,0]
# 6, {[1,3,2],[2,1,3],[3,2,1]},         [1,2,3,4,5,6,7,8]
# 7, {[1,2,3],[2,3,1],[3,2,1]},         [1,2,3,1,0,0,0,0]
# 8, {[1,3,2],[2,1,3]},                 [1,2,4,8,16,32,64,128]
# 9, {[1,3,2],[2,1,3],[3,1,2]},         [1,2,3,4,5,6,7,8]
#10, {[1,3,2],[3,1,2],[3,2,1]},         [1,2,3,4,5,6,7,8]

# OEIS matches:
# [1,2,3,4,5,6,7,8]       = A000027 (positive integers).  -> 5 times
# [1,2,4,8,16,32,64,128]  = A000079 (powers of 2, shifted). -> 2 times
# So in OEIS: 5 + 2 = 7 out of 10.

# Non-OEIS type observed:
# [1,2,3,1,0,0,0,0]  (3 times)
# Conjecture (for this S):
#   a(1)=1, a(2)=2, a(3)=3, a(4)=1, and a(n)=0 for n>=5.
a3 := n -> piecewise(n<=3, n, n=4, 1, 0);


# 4. Run 10 times:
# S:={seq(randperm(3),i=1..2), seq(randperm(4),i=1..4)};
# [seq(nops(AvoidPer(n,S)),n=1..7)];
# How many in OEIS? For those not, conjecture a formula.


# My recorded 10 outputs (t, S, L):
# 1: [1,2,5,13,34,89,233]
# 2: [1,2,4,7,10,13,16]
# 3: [1,2,4,7,12,20,33]
# 4: [1,2,4,7,12,20,33]
# 5: [1,2,4,7,11,16,22]
# 6: [1,2,4,7,13,24,44]
# 7: [1,2,4,4,1,0,0]
# 8: [1,2,5,12,26,51,92]
# 9: [1,2,4,8,16,32,64]
#10: [1,2,4,8,16,32,64]

# OEIS matches:
# 1: [1,2,5,13,34,89,233] is Fibonacci odd-index subsequence (shift of A001519).
# 2: [1,2,4,7,10,13,16]  is A033627.
# 3/4: [1,2,4,7,12,20,33] is shift of A000071 (= Fibonacci(n)-1).
# 5: [1,2,4,7,11,16,22]  is A000124.
# 6: [1,2,4,7,13,24,44]  is shift of A276661.
# 8: [1,2,5,12,26,51,92] is A027927.
# 9/10: [1,2,4,8,16,32,64] is A000079.

# So in OEIS: 9 out of 10 (everything except run 7).

# For the non-OEIS output (run 7):
# Conjecture:
#   a(1)=1, a(2)=2, a(3)=4, a(4)=4, a(5)=1, a(n)=0 for n>=6.

a4 := n -> piecewise(n=1,1, n=2,2, n=3,4, n=4,4, n=5,1, 0);


# 5. Optional Challenge: Identify b(n) in C4.txt and prove integrality.


# The sequence b(n) is the Somos-4 sequence (OEIS A006720). 
# Recurrence (typical form):
#   b(0)=b(1)=b(2)=b(3)=1,
#   b(n) = (b(n-1)*b(n-3) + b(n-2)^2)/b(n-4)   for n>=4.
#
# Proof of integrality:
#
# Step 1 (symbolic version).  Start with indeterminates x0,x1,x2,x3 and define xn by
#   xn*x(n-4) = x(n-1)*x(n-3) + x(n-2)^2.
#
# Step 2 (Laurent phenomenon).  Fomin–Zelevinsky proved the “Laurent phenomenon”:
# for recurrences of this cluster-algebra type (including Somos-4), every term xn is
# actually a Laurent polynomial in the initial variables:
#   xn ∈ Z[x0^{±1}, x1^{±1}, x2^{±1}, x3^{±1}].
# In particular, denominators are only monomials in x0,x1,x2,x3 (no other factors).
#
# Step 3 (subtract).  Our integer sequence b(n) is obtained by setting x0=x1=x2=x3=1.
# Since any monomial in the denominator becomes 1, we get xn(1,1,1,1) ∈ Z for every n.
# Therefore b(n) is an integer for all n.
#


#6. I can show there will be no integers after a(17)

# THEOREM.
# Let (a_n)_{n≥1} be defined by
#     a_1 = 2,
#     a_n = (1/(n-1)) * Σ_{i=1}^{n-1} a_i^2    for n≥2.
# Then
#     for every n ≥ 18,  a_n ∉ Z.
#
# PROOF.
#
# 1)recurrence.
# From the definition,
#     (n-1) a_n = Σ_{i=1}^{n-1} a_i^2.                           (1)
# Also,
#     n a_{n+1} = Σ_{i=1}^{n} a_i^2
#               = (Σ_{i=1}^{n-1} a_i^2) + a_n^2
#               = (n-1) a_n + a_n^2     (by (1)).
# Hence for every n≥1,
#     a_{n+1} = a_n(a_n + n - 1)/n.                               (2)
#
# 2) finite computation: a_17 ≡ 9 (mod 17), hence a_18 is not an integer.
# Using (2) (and the fact that, by direct computation, a_1,...,a_17 are integers),
# we may reduce them modulo 17. Reducing (2) mod 17 for n=1,2,...,16 yields the
# residue sequence
#     a_1,...,a_17  ≡  2,4,10,6,5,9,4,13,7,6,9,14,15,5,4,15,9   (mod 17),
# so in particular
#     a_17 ≡ 9 (mod 17).                                          (3)
#
# Now apply (2) at n=17:
#     a_18 = a_17(a_17 + 16)/17 = a_17(a_17 - 1)/17.              (4)
# By (3),
#     a_17(a_17 - 1) ≡ 9·8 = 72 ≡ 4 (mod 17),
# so the numerator in (4) is NOT divisible by 17. Therefore the factor 17 in the
# denominator does not cancel, and
#     a_18 ∉ Z, and in lowest terms the denominator of a_18 is divisible by 17.  (5)
#
# 3) Persistence: once 17 divides a denominator, it can never disappear.
# LEMMA.
# Suppose a_n = p/q is in lowest terms (gcd(p,q)=1) and 17 | q. Then a_{n+1} is
# not an integer; in fact, in lowest terms its denominator is still divisible by 17.
#
# Proof of Lemma.
# From (2),
#     a_{n+1} = (p/q) * ( (p/q) + (n-1) ) / n
#             = p(p + (n-1)q) / (n q^2).                          (6)
# Since 17 | q and gcd(p,q)=1, we have 17 ∤ p. Also (n-1)q ≡ 0 (mod 17), hence
#     p + (n-1)q ≡ p (mod 17),
# so 17 ∤ (p + (n-1)q) as well. Thus the numerator p(p + (n-1)q) is NOT divisible
# by 17, while the denominator n q^2 IS divisible by 17. Therefore no factor 17 can
# cancel when reducing (6), so the reduced denominator of a_{n+1} is divisible by 17.
# In particular, a_{n+1} cannot be an integer.  □
#
# 4) Conclusion.
# By (5), a_18 has reduced denominator divisible by 17. Applying the Lemma repeatedly
# for n=18,19,20,... shows that every a_n for n≥18 has reduced denominator divisible
# by 17, hence is not an integer.
#
# Therefore, for every n ≥ 18,  a_n ∉ Z.  □



# 7. 132-avoiding permutations:

# Let a(n) = # permutations of length n avoiding pattern 132.
# Set a(0)=1 for convenience (empty permutation).
#
# Proof of recurrence:
# Take a 132-avoiding permutation π of length n and locate the entry n.
# Suppose n has k entries to its left, so n is in position k+1, with 0<=k<=n-1.
#
# Claim: every entry left of n is larger than every entry right of n.
# Otherwise, pick x on the left and y on the right with x<y. Then (x, n, y)
# forms a 132 pattern (since x<y<n and positions are left < n < right), contradiction.
#
# Therefore the left block consists of the k largest values {n-k,...,n-1}
# and the right block consists of {1,...,n-k-1}.
# Standardizing each block gives two independent 132-avoiding permutations
# of lengths k and n-1-k. Conversely, any such pair gives a unique π.
# Hence, for each k, the number is a(k)*a(n-1-k), and summing gives:
#   a(n) = Sum_{k=0..n-1} a(k)*a(n-1-k).

# Maple code to crank out terms (with remember):
A := proc(n)
option remember;
local k;
if n=0 then
  return 1;
else
  return add( A(k)*A(n-1-k), k=0..n-1 );
fi;
end:

seq(A(n), n=0..15);

# Output:
# 1, 1, 2, 5, 14, 42, 132, 429, 1430, 4862, 16796, 58786, 208012, 742900, 2674440, 9694845, ...

# Conjecture explicit formula (Catalan numbers):
#   a(n) = C_n = binomial(2n,n)/(n+1).
Catalan := n -> binomial(2*n,n)/(n+1):

# quick check for n=0..10:
[seq(A(n)=Catalan(n), n=0..10)];