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

with(combinat):


# Question 1:
# Another sequence of rational numbers that converges to Pi is
#  16^n/(n*binomial(2*n,n)^2) for n=1,2,...
#  (a) Using Maple, prove it!
#  (b) Write a little procedure that outputs the first N terms of
#      this sequence, how many digits of PI do you get from the 1000-th term?
#      How does it compare with PiAppxVC(1000)[1000] ?
# ANSWER:
# (a) limit(16^n/(n*binomial(2*n,n)^2),n=infinity); returns Pi
# (b) From the 1000-th term, I only get that it agrees to the first two digits 
#  PiAppxVC(1000)[1000] does much better. PiAppxVC(1000)[1000] agrees
#  with every digit when you set Digits:=120 and compare the first 1,000 terms
PiB:=proc(N) local n:
[seq(16^n/(n*binomial(2*n,n)^2),n=1..N)]:
end:


#Question 2:
# Read, understand and implement:
#  https://sites.math.rutgers.edu/~zeilberg/mamarim/mamarimPDF/pi14.txt
#  Code K(n) according to the definition (using the integral) and also,
#  calling it Kc(n), using the second-order recurrence
# What is faster between the following:
#  [seq(K(n),n=1..200)] 
#  [seq(Kc(n),n=1..200)]
# ANSWER:
#  time([seq(K(n),n=1..200)]); returns 245.421 in my laptop
#  time([seq(Kc(n),n=1..200)]); returns 0.046 in my laptop

K:=proc(n) local x:
int((x^(4*n)*(1-x)^(4*n))/(x^2+1),x=0..1):
end:

Kc:=proc(n) local d:
option remember:

if n=0 then
 return(Pi/4):
fi:

if n=1 then
 return(22/7-Pi):
fi:


(-4*(4*(n-2) + 1)*(4*(n-2) + 3)*(820*(n-2)^3  + 3993*(n-2)^2  + 6428*(n-2) + 3420)*(n-1)*(2*(n-2) + 1)*Kc(n-2)+
(26843520*(n-2)^7  + 211258528*(n-2)^6  + 688917624*(n-2)^5  + 1202271190*(n-2)^4 + 1207256235*(n-2)^3  + 
693651577*(n-2)^2  + 209656416*(n-2) + 25472340)*Kc(n-1))/(-2*((8*(n-2) + 9)*(8*(n-2) + 11)*(8*(n-2) + 13)*(8*(n-2) + 15)*(820*(n-2)^3  + 1533*(n-2)^2  + 902*(n-2) + 165))):

end:


#Question 3:
# Prove that indeed the limit of PiAppx(n) is Pi
# Hint: You must use calculus and the ratio test



#########################From previous class:
# C15.txt, March 12, 2026
# Almost PI Day :)

with(combinat):

Help15:=proc():
 print(`IsUD(pi), IsDU(pi), UDperms(n), A111s(n)`):
 print(`PiAppxS(n), UDni(n,i), DUni(n,i), UDn(n), A111c(n)`):
 print(`PiAppx(n), UDc(n), A111vc(n), PiAppxVC(n)`):
end:


#IsUD(pi): Given a permutation pi, is pi an Up-down permutation?
IsUD:=proc(pi) local n:
n:=nops(pi):
if n=1 then
 return(true):
fi:

#in other words, does the permutation start with a down?
if pi[1]>pi[2] then
 return(false):
fi:

if IsDU([op(2..n,pi)]) then
 return(true):
else
 return(false):
fi:
FAIL:
end:


#IsDU(pi): Given a permutation pi, is pi a down-up permutation?
IsDU:=proc(pi) local n:
n:=nops(pi):
if n=1 then
 return(true):
fi:

#in other words, does the permutation start with a down?
if pi[2]>pi[1] then
 return(false):
fi:

if IsUD([op(2..n,pi)]) then
 return(true):
else
 return(false):
fi:
FAIL:
end:


#UDperms(n): The set of up-down permutations of {1,2,...,n}
# also known as zigzag permutations or Andre permutations
UDperms:=proc(n) local A,G,pi:
A:=permute(n):
G:={}:
for pi in A do
 if IsUD(pi) then
  G:=G union {pi}:
 fi:
od:
G:
end:


#A111s(n): The brute force way to compute the first n terms
# of the sequence of the up-down permutations
# OEIS sequence A000111
A111s:=proc(n) local i:
[seq(nops(UDperms(i)),i=1..n)]:
end:



#PiAppxS(n): The first n approximations to the combinatorial way to compute
# PI via up-down permutations
PiAppxS:=proc(n) local L,i:
L:=A111s(n):
[seq(2*i*L[i-1]/L[i],i=2..n)]
end:


#How can we compute the number of up-down perms ?
#Let UDn(n): number of up-down permutations of {1,2,...,n}
#UDni(n,i): Number of up-down permutions of {1,2,...,n} that start with i 
# where 1<=i<=n
# [i,....]
#if pi=[i,j,...] with i<j then what happens if you chop i?
# you get down-up permutation of {1,2,...,i-1,i+1,...n} of lenght n-1
# that start with j and since i<jj then we need to rename j
#So, 
# UDni(n,i)=add(DUni(n-1,j), j=i..n-1)


#UDni(n,i): The number of up-down permutions of {1,2,...,n} that start with i
UDni:=proc(n,i) local j:
option remember:
if n=1 then
 if i=1 then
  return(1):
 else
  return(0):
 fi:
fi:
add(DUni(n-1,j), j=i..n-1):
end:


#How to count the following set?
#Let pi be a down-up permutation that starts with i 
# so pi=[i,j,...] with i>j and j=1,2,...,i-1 and if you chop i then you get
# and up-down permutation of {1,2,...,i-1,i+1,...n}
# that starts with j

#DUni(n,i): The number of down-up permutions of {1,2,...,n} that start with i
DUni:=proc(n,i) local j:
option remember:
if n=1 then
 if i=1 then
  return(1):
 else
  return(0):
 fi:
fi:
add(UDni(n-1,j), j=1..i-1):
end:


#UDn(n): The number of up-down permutations of fixed length n
UDn:=proc(n) local i:
option remember:
add(UDni(n,i),i=1..n):
end:


#A111c(n): The clever way to compute the first n terms
# of the sequence of the up-down permutations
# OEIS sequence A000111
A111c:=proc(n) local k:
[seq(UDn(k),k=1..n)]:
end:


#PiAppx(n): The first n approximations to the combinatorial way to compute
# PI via up-down permutations
PiAppx:=proc(n) local L,i:
L:=A111c(n):
[seq(2*i*L[i-1]/L[i],i=2..n)]
end:


#A non-linear recurrence to compute UDn(n)?
# Suppose the permutation is as follows:
# ...n...=pi1 n pi2 ; where ... on the left is renamed pi1
# Suppose that n is in location k but it has to be in the EVEN position
#  pi1 is a permutation of some subset of size k-1
#  of {1,2,...,n-1} that is up-down
#  pi2 is a permutation of the complement of the above set is ALSO
#  up-down
#  UDn(n)=Sum_(k=1..n) binomial(n-1,k-1)*UDn(k-1)*UDn(n-k)
#  binomial(n-1,k-1) choose what numbers to place in the perm pi1
#  UDn(k-1) counts the number of perms of length k-1 for pi1
#  UDn(n-k) counts the number of perms of length n-k for pi2


#UDc(n): The number of up-down permutations done even more cleverly
UDc:=proc(n) local k:
option remember:
if n=0 or n=1 then
 return(1):
fi:
add(binomial(n-1,2*k-1)*UDc(2*k-1)*UDc(n-2*k),k=1..trunc(n/2)):
end:



#A111vc(n): The very clever way to compute the first n terms of the OEIS sequence A000111
A111vc:=proc(n) local i:
[seq(UDc(i),i=1..n)]:
end:


#PiAppxVC(n): The first n approximations to the combinatorial way to compute PI via
# up-down permutations
PiAppxVC:=proc(n) local L,i:
L:=A111vc(n):
[seq(2*i*L[i-1]/L[i],i=2..n)]:
end:




###Notes:
#Due to Stirling:
# The asymptotics of n! is (n/e)^n*sqrt(2*PI*n) as n goes to infinity
#UpDown Permutations:
# One way to compute PI is related to counting Up-down permutations
# 
#Let L := [1, 1, 2, 5, 16, 61, 272, 1385, 7936] with n=9
# evalf(2*L[8]*9/L[9])=3.141381048
#Above, we are evaluating the end behavior
#So:
# Limit of 2*n*A111(n-1)/A111(n) as n goes to infinity is the number PI