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