#Yusra Naqvi: HW #15
#OK to post

################################################################################

#1
#An(n): comuptes the number of up-down permutations of {1, ...,n} 

An:=proc(n) local k,i:

option remember:

if n=0 then
 return(1):
elif n=1 then
 return(1):
else
 k:=trunc(n/2):
 return(add((An(2*i-1)*An(n-2*i))*binomial(n-1,2*i-1),i=1..k)):
fi:

end:
################################################################################

#2

###

#RtoS(a,d,k): list of the first k digits in the base d representation of the 
#real number a between 0 and 1

RtoS:=proc(a,d,k) local a1,L,i,kir:

a1:=evalf(a):
if a1<0 or a1>1 then
 return(FAIL):
fi: 
 
L:=[]:

for i from 1 to k do
 a1:=a1*d:
 kir:=trunc(a1):
 L:=[op(L),kir]:
 a1:=a1-kir:
od:

if add(L[i]/d^i,i=1..k)-evalf(a)>1/d^(k-1) then
 return(FAIL):
else:
 return(L):
fi:

end:

###

#MaxDev(a,d,w,N): inputs a real number a between 0 and 1, a pos. integer d ≤ 2,
#a word w in the alphabet {0,...,d-1} (written as a list) and a pos. integer N,
#and outputs the number of occurrences of the word w as consecutive substring
#in the base-d representation of a in the first N*d^(nops(w)) digits

MaxDev:=proc(a,d,w,N) local ad,co,i:

ad:=RtoS(a,d,N*d^(nops(w))):

co:=0:

for i from 1 to nops(ad)-nops(w)+1 do
 if ad[i..i+nops(w)-1]=w then
  co:=co+1:
 fi: 
od:

co:

end:

###

#MaxDev(Pi-3,2,[0,1,0],1000);
#985

#MaxDev(Pi-3,10,[8,9],100);
#93

################################################################################

#3
#J(n): inputs a positive integer n and outputs
#int(t^(4*n)*(1-t)^(4*n),t=0..1)
#in terms of rational numbers and Pi

J:=proc(n):

int(t^(4*n)*(1-t)^(4*n),t=0..1):

end:

###

#Using
#GuessH([seq(J(i),i=1..40)],n,N);
#we get
#(1/32*(4*n+1))*(2*n+1)*(4*n+3)*(n+1)-(1/16*(8*n+5))*(8*n+7)*(8*n+9)*(8*n+3)*N

###

#Using
#ope:=(1/32*(4*n+1))*(2*n+1)*(4*n+3)*(n+1)-(1/16*(8*n+5))*(8*n+7)*(8*n+9)*(8*n+3)*N:
#and 
#HtoSeq(ope,n,N,[J(1)],4000)[4000];
#we get the answer almost instantly.
#On the other hand,
#J(4000);
#took so long that I cancelled it before it finished. 
#We can measure the difference for n=4000 exactly using the following procedure.

TimeCompJ:=proc(n) local t1,t2,t3,ope:

t1:=time():
ope:=(1/32*(4*n+1))*(2*n+1)*(4*n+3)*(n+1)-(1/16*(8*n+5))*(8*n+7)*(8*n+9)*(8*n+3)*N:
HtoSeq(ope,n,N,[J(1)],n)[n];
t2:=time():
J(n):
t3:=time():

(t3-t2)-(t2-t1): 

end:


################################################################################