#!/usr/local/bin/maple
# -*- maplev -*-

# Nathaniel Shar
# HW 15
# Experimental Mathematics

#############
# Problem 1 #
#############

# Condition on the location of n. Say it is in the (2k)th place. The
# permutation preceding n is an up-down permutatation. The permutation
# following it is also an up-down permutation. There are binomial(n-1
# choose 2k-1) ways to choose which entries to put on the left.

# Therefore, the number of up-down permutations A(n) of length n is

An := proc(n) local k: option remember:
    if n <= 1 then: return 1: fi:
    return add(binomial(n-1, 2*k-1)*An(2*k-1)*An(n-2*k), k=1..floor(n/2)):
end:

# When I calculate 2/(An(200)/(200*An(199))), I get something within
# 10^(-94) of pi.

#############
# Problem 2 #
#############

##########################
# Borrowed from Class 15:

 #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:
#########################

# Sleazy solution: A word of length nops(w) in base d is just the same
# as a digit in base d^nops(w). However, this will only catch
# occurences if they are properly aligned. So we just have to shift
# the number a few times.

CountDigits := proc(a,d,k,digit) local L, i, counter:
    L := RtoS(a,d,k):
    counter := 0:
    for i from 1 to nops(L) do:
        if L[i] = digit then:
            counter := counter + 1:
        fi:
    od:
    return counter:
end:

# I'm modifying the return value of the function slightly

# This returns the PROPORTION of strings matching w, out of the first N
# strings. It'll be in the form of a fraction.

MaxDev := proc(a,d,w,N) local i, digit, counter:
    # Produce enough digits, please!
    if Digits <= 2*nops(w) + evalf(log(d)*nops(w)/log(10))*N then:
        printf("Not enough digits! Increase Digits.\n"):
    end:
    digit := add(w[i]*d^(nops(w)-i), i=1..nops(w)):
    counter := 0:
    for i from 1 to nops(w) do:
        counter := counter + CountDigits(frac(a*10^(i-1)), d^nops(w), floor(N/nops(w))+1, digit):
    od:
    return counter/N:
end:

# for this version:

# MaxDev(Pi-3,2,[0,1,0],1000) is 119/1000 -- that is, the string 010
# appears 119 times in the first 1002 binary digits of pi
# MaxDev(Pi-3,10,[8,9],10000) is 93/10000

#############
# Problem 3 #
#############

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

# I don't understand what's going on here. This is probably not the
# integral you wanted us to compute -- it's just a polynomial -- but I
# can't figure out what the right one is.