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

Help:=proc(): print(`hw15ApproxN(N), K(n), Kc(n)`):   end:

### Problem 1
# Another sequence of rational numbers that converges to Pi is
# 16^n/(n*binomial(2*n,n)^2), n=1,2,. ..

# 1. Using Maple, prove it!

# [seq(evalf(16^n/(n*binomial(2*n,n)^2)-Pi), n=1..50)];


# 2. 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]?

hw15ApproxN := proc(N) local n: 
[seq(16^n/(n*binomial(2*n,n)^2), n=1..N)]:
end:

# Digits:=10000:
# -trunc(log[10](abs(evalf(hw15ApproxN(1000)[-1]-Pi)))): 
# output: 3
# -trunc(log[10](abs(evalf(PiAppxVC(1000)[-1]-Pi)))):
# output: 476 



### Problem 2
# Code K(n) according to the definition (using the integral) and also, calling it Kc(n), using the second-order recurrence (don't forget option remember).

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


Kc := proc(n) local m: 
option remember:
if n=0 then 
    RETURN(Pi/4):
elif n=1 then 
    RETURN(22/7-Pi):
fi:

# K(n+2) = K(m) , so we introduce a dummy variable instead of reindexing
m := n-2: 
                                      
(-4*(4*m + 1)*(4*m + 3)*(820*m^3  + 3993*m^2 + 6428*m + 3420)*(m + 1)*(2*m + 1)*K(m)
+ (26843520*m^7 + 211258528*m^6 + 688917624*m^5 + 1202271190*m^4 + 1207256235*m^3 + 693651577*m^2 + 209656416*m + 25472340)*K(m+1) )
/(-2*(8*m + 9)*(8*m + 11)*(8*m + 13)*(8*m + 15)*(820*m^3 + 1533*m^2 + 902*m + 165)):

end:


# What is faster
# [seq(K(n),n=1..200)] OR [seq(Kc(n),n=1..200)]?
# [Use the Maple command package]


# time([seq(K(n),n=1..200)]);
# output: 467.250

# time([seq(Kc(n),n=1..200)]);
# output: 1.906






### copied from C15.txt 

#C15.txt March 12, 2026 Pi  Day special

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

with(combinat):

IsUD:=proc(pi) local n:
n:=nops(pi):
if n=1 then
 RETURN(true):
fi:

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:=proc(pi) local n:
n:=nops(pi):
if n=1 then
 RETURN(true):
fi:

if pi[1]<pi[2] 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,...n} (aka Zigzag permuations, aka Andre permutation)
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 (stupid) way to compute the first n terms of the OEIS sequence A111
 #(nops(UDperms(n)))
 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
 #the stupid way
 PiAppxS:=proc(n) local L,i:
 L:=A111s(n):
 
 [seq(2*i*L[i-1]/L[i],i=2..n)]:
 end:
 
 #UDni(n,i): the NUMBER of up-down permutations of {1, ...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:
 
 #DUni(n,i): the NUMBER of down-up permutations of {1, ...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:=proc(n) local i:
 add(UDni(n,i),i=1..n):
 end:
 
 #A111c(n): The clevered way to compute the first n terms of the OEIS sequence A111
 #(nops(UDperms(n)))
 A111c:=proc(n) local i:
 [seq(UDn(i),i=1..n)]:
 end:
 
 #PiAppx(n): The first n approximations to the combinatorial way to compute pi via up-down permutations
 #the clever way
 PiAppx:=proc(n) local L,i:
 L:=A111c(n):
 
 [seq(2*i*L[i-1]/L[i],i=2..n)]:
 end:
 
 #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 A111
 #(nops(UDperms(n)))
 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
 #the very clever way
 PiAppxVC:=proc(n) local L,i:
 L:=A111vc(n):
 
 [seq(2*i*L[i-1]/L[i],i=2..n)]:
 end: