#John Kim
#Use whatever you like

#read "C:/Users/John Y. Kim/Documents/Maple/hw15.txt":

Help:=proc(): print(`IsUD(pi),AnStupid(n), rf(a,k) ,SR1(M) `): 
print(` RtoS(a,d,k), RtoS(a,d,k) `):  end:


 with(combinat):

  #IsUD(pi): Is the permutation pi Up-Down
  #For example IsUD([1,3,2]); is true 
   #IsUD([1,2,3]); is false
  IsUD:=proc(pi) local i:
  evalb(
  {seq( evalb(pi[2*i-1]<pi[2*i]), i=1..nops(pi)/2),
   seq( evalb(pi[2*i]>pi[2*i+1]), i=1..(nops(pi)-1)/2)}=
     {true}):
   end:

   #AnStupid(n): the number of up-down permutations of {1, ...,n} 
   AnStupid:=proc(n) local co, S, pi:
   S:=convert(permute(n),set):

   co:=0:

   for pi in S do
     if IsUD(pi) then
         co:=co+1:
     fi:
  od:

  co:
  end:



 #rf(a,k):=a(a+1)(a+2)...*(a+k-1)
 rf:=proc(a,k) local i:
  mul(a+i,i=0..k-1):
 end:

 #using M terms in Ramanujan's slow formula for 2/Pi
 SR1:=proc(M) local k:
 2/evalf(add((-1)^k*(4*k+1)*rf(1/2,k)^3/k!^3,k=0..M)):
 end:


 #RtoS(a,d,k): the list of the first k digits in the base-d
 #rep. of the real number a between 0 and 1
 RtoS:=proc(a,d,k) local a1,L,kirsten,i:

  a1:=evalf(a):
  if a1>=1 or a1<=0 then
    RETURN(FAIL):
  fi:

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


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

    RETURN(FAIL):

   else
    RETURN(L):

   fi:

end:


#An(n): the number of up-down permutations on the first n positive integers.
An:=proc(n) option remember: local i:

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

end:

#evalf(An(199)/An(200)*200*2-Pi); gives 1.5770*10^(-95), so it is accurate to 95 digits.
   


#MaxDev(a,d,w,N): inputs a real number a between 0 and 1 (that later we will take to be Pi-3),
#a pos. integer d<=10, a word w in the alphabet {0, ..., d-1} (written as list) and a positive 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. If you believe that a is normal, 
#then it should be (for large N) roughly N times.
#What are:
#MaxDev(Pi-3,2,[0,1,0],1000)? 985
#MaxDev(Pi-3,10,[8,9],100)? 93

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

co:=0:
L:=RtoS(a,d,N*d^nops(w)):

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

co:

end:




#J(n): inputs a positive integer n and outputs
#int(t^(4*n)*(1-t)^(4*n)/(1-t^2)^(1/2),t=0..1) 

J:=proc(n) local t:
 int(t^(4*n)*(1-t)^(4*n)/(1-t^2)^(1/2),t=0..1):
end: