###################################################################### ## HIT1.txt: Save this file as HIT1.txt to use it, # # stay in the # ## same directory, get into Maple (by typing: maple ) # #then type: read `HIT1.txt`: # ## Then follow the instructions given there # ## # ## Written by Lucy Martinez, Rutgers University # ###################################################################### #VERSION OF Dec. 06, 2022 with(combinat): print(`First Written: Dec. 2022: tested for Maple 20 `): print(): print(`This is HIT1.txt, One of the Maple packages`): print(`accompanying Robert Dougherty-Bliss, Lucy Martinez and Doron Zeilberger's article: `): print(` How many Dice Rolls Would It Take to Reach Your Favorite Kind of Number? `): print(`-------------------------------`): print(`-------------------------------`): print(): print(`The most current version is available on WWW at:`): print(` http://www.math.rutgers.edu/~zeilberg/tokhniot/HIT.txt .`): print(`Please report all bugs to: zeilberg at math dot rutgers dot edu .`): print(): print(`-------------------------------`): print(`For general help, and a list of the MAIN functions,`): print(` type "Help();". For specific help type "Help(procedure_name);" `): print(`-------------------------------`): print(`For a list of the supporting functions type: Help1();`): print(`For help with a specific procedure, type "Help1(procedure_name);"`): print(): print(`-------------------------------`): Help1:=proc() if args=NULL then print(` The supporting procedures are: IsEven, IsPP, IsPPP, OneTime, OneTimeG, OneTime3, p, pGen, pGenSum, pGenProp, Pow2 `): print(` `): else ezra(args): fi: end: Help:=proc() if args=NULL then print(` HIT1.txt: A Maple package for studying the statistics of the time it gets to reach your favorite kind of numbers rolling a die, Using the Alon-Malinovsky approach `): print(`The main procedures are: Paper, pk, pkG, Simu, SimuG `): elif nargs=1 and args[1]=IsEven then print(`IsEven(n) checks whether n is even, returns true/false`): elif nargs=1 and args[1]=IsPP then print(`IsPP(n) checks whether n is a product of two distinct primes, returns true/false`): elif nargs=1 and args[1]=IsPPP then print(`IsPPP(n) checks whether n is a product of three distinct primes, returns true/false`): elif nargs=1 and args[1]=OneTime then print(` OneTime() rolls a 6-faced die, each time it adds up the result after rolling the die,`): print(` and stops until the sum is equal to a prime number, it then returns the number of iterations`): print(` it took to get a prime sum`): elif nargs=1 and args[1]=OneTimeG then print(` OneTimeG(G,init,K) starts with sum equal to an initial value (init), adds up`): print(` the result from rolling a 6-faced die, it keeps rolling the die while`): print(` property P has not been met. As soon as the "sum" meets property P, it stops and`): print(` returns the number of iterations`): print(` This is for any property on a 6-faced die, we run this procedure until`): print(` we get what we want. Try:`): print(`OneTimeG(isprime,0,6);`): elif nargs=1 and args[1]=OneTime3 then print(` OneTime3(N) rolls a N-faced die, each time it adds up the result after rolling the die`): print(` and stops until the sum is equal to a prime number, it then returns the number of iterations`): print(` it took to get a prime sum. Try OneTime3(7);`): elif nargs=1 and args[1]=p then print(` p(k,n): Inputs a non-prime pos. integer n and a pos. integer k, The probability that when rolling a fair die k times the sum is n and all its partial sums are also non-prime`): print(`in other words the probability of not yet-arriving at primes after k rolls, and the sum is n. Of course k<=n<=6*k for it to be non-zero. Try:`): print(`p(10,30);`): elif nargs=1 and args[1]=Paper then print(`Paper(P,k,N): Given a property P, a positive integer k, and a positive integer N outputs`): print(`a paper with estimates for the expected duration, starting at 0, of hitting a number with property P`): print(`with a fair die with n faces, from n=2 to n=N. It assumes that you roll it up to k times. Try:`): print(`Paper(isprime,100,10):`): elif nargs=1 and args[1]=pGen then print(` pGen(k,n,N) returns the approximation to expectation to get a prime sum with a N-faced die using up to k rolls. `): print(`pGen(k,n,6) is the same as p(k,n);`): elif nargs=1 and args[1]=pGenProp then print(` pGenProp(P,k,n,N) returns the approximation to expectation to get a sum with property P with a N-faced die using up to k rolls and`): print(` it starts with an initial value depending on the property P`): print(` pGenProp(isprime,k,n,6) is the same as p(k,n)`): elif nargs=1 and args[1]=pkG then print(` pkG(P,k,N) adds up pGenProp(P,i,j,N) for 1<=i<=k and k<=j<=kN`): print(` pkG(isprime,k,6) is the same as pk(k); Try:`): print(`pkG(isprime,100,6);`): print(`pkG(isprime,100,9);`): print(`pkG(IsPP,100,9);`): elif nargs=1 and args[1]=pGenSum then print(` pGenSum(k,N) adds up pGen(i,j,N) for 1<=i<=k and k<=j<=kN. pGenSum(k,6) is the same as pk(k). Try:`): print(`pGenSum(10,6);`): elif nargs=1 and args[1]=pk then print(` pk(k) adds up p(i,j) for 1<=i<=k and k<=j<=6*k plus 1. This is the estimated expected number of rolls to hit a prime if k is large. Try:`): print(`pk(100);`): elif nargs=1 and args[1]=Pow2 then print(` Pow2(n) checks whether n is a power of 2, returns true/false `): elif nargs=1 and args[1]=Simu then print(`Simu(K) adds up OneTime() K times and returns its average. Try: `): print(`Simu(10000);`): elif nargs=1 and args[1]=SimuG then print(`SimuG(P,init,F,K) runs OneTimeG(P,init) K times, adds it up and returns the average.Try`): print(`SimuG(isprime,0,6,1000);`): else print(`There is no such thing as`, args): fi: end: with(combinat): with(NumberTheory): # OneTime() rolls a 6-faced die, each time it adds up the result after rolling the die, # and stops until the sum is equal to a prime number, it then returns the number of iterations # it took to get a prime sum OneTime:=proc() local i,s: s:=0: for i from 1 to 1000 while not isprime(s) do s:=s+rand(1..6)(): od: i-1: end: # Simu(K) adds up OneTime() K times and returns its average sum Simu:=proc(K) local s,i: s:=0: for i from 1 to K do s:=s+OneTime(): od: return evalf(s/K): end: # IsEven(n) checks whether n is even, returns true/false IsEven:=proc(n): evalb(n mod 2=0): end: # Pow2(n) checks whether n is a power of 2, returns true/false Pow2:=proc(n) local x: if n<1 then RETURN(false): fi: x:=log(n)/log(2): if type(x,integer) and n>1 then return true: else return false: fi: end: # IsPP(n) checks whether n is a product of two distinct primes, returns true/false IsPP:=proc(n) local S, i,pr: if n<=0 then return false: else S:=PrimeFactors(n): if nops(S)=2 then pr:=1: for i from 1 to nops(S) do pr:=pr*S[i]: od: if pr=n then return true: else return false: fi: else return false: fi: fi: end: # IsPPP(n) checks whether n is a product of three primes, returns true/false IsPPP:=proc(n) local S, i,pr: if n<=0 then return false: else S:=PrimeFactors(n): if nops(S)=3 then pr:=1: for i from 1 to nops(S) do pr:=pr*S[i]: od: if pr=n then return true: else return false: fi: else return false: fi: fi: end: # OneTimeG(P,init,F) starts with sum equal to an initial value (init), adds up # the result from rolling a F-faced die, it keeps rolling the die while # property P has not been met. As soon as the "sum" meets property P, it stops and # returns the number of iterations # This is for any property on a 6-faced die, we run this procedure until # we get what we want. Try: #OneTimeG(isprime,0,6); OneTimeG:=proc(P,init,F) local s,i: s:=init: for i from 1 to 1000 while not P(s) do s:=s+rand(1..F)(): od: i-1: end: # SimuG(P,init,F,K) runs OneTimeG(P,init,F) K times, adds it up and returns the average. SimuG:=proc(P,init,F,K) local i: evalf(add(OneTimeG(P,init,F),i=1..K)/K): end: # p(k,n): Inputs a non-prime pos. integer n and a pos. integer k, The probability that when rolling a fair die k times the sum is n and all its partial sums are also non-prime #in other words the probability of not yet-arriving at primes after k rolls, and the sum is n. Of course k<=n<=6*k for it to be non-zero. Try: #p(10,30); p:=proc(k,n) local i,S: option remember: S:={}: if isprime(n) or n<=0 then return 0: fi: if k=1 and n=1 then return 1/6: elif k=1 and n=4 then return 1/6: elif k=1 and n=6 then return 1/6: else for i from 1 to 6 do if not isprime(n-i) then S:=S union {i}: fi: od: 1/6*add(p(k-1,n-i),i in S): fi: end: # pk(k) adds up p(i,j) for 1<=i<=k and k<=j<=6*k plus 1. This is the estimated expected number of rolls to hit a prime if k is large. Try: #pk(200); pk:=proc(k) local i,j,s: s:=0: for i from 1 to k do for j from i to 6*i do s:=s+p(i,j): od: od: 1+s: end: # pGen(k,n,N) returns the approximation to expectation to get a prime sum with a N-faced die using up to k rolls. #pGen(k,n,6) is the same as p(k,n) pGen:=proc(k,n,N) local i,S: option remember: S:={}: if isprime(n) or n<=0 then return 0: fi: # we need n<=N, otherwise, pGen(1,8,6) is gonna return 1/6 if k=1 and isprime(n)=false and n<=N then return 1/N: else for i from 1 to N do if not isprime(n-i) then S:=S union {i}: fi: od: 1/N*add(pGen(k-1,n-i,N),i in S): fi: end: # pGenSum(k,N) adds up pGen(i,j,N) for 1<=i<=k and k<=j<=kN. pGenSum(k,6) is the same as pk(k). Try: #pGenSum(10,6); pGenSum:=proc(k,N) local s,i,j: option remember: s:=0: for i from 1 to k do for j from i to N*i do s:=s+pGen(i,j,N): od: od: 1+s: end: # pGenProp(P,k,n,N) returns the expectation to get a "P" sum with a N-faced die pGenProp:=proc(P,k,n,N) local i,S: option remember: S:={}: if P(n) or n<=0 then return 0: fi: if k=1 and P(n)=false and n<=N then return 1/N: else for i from 1 to N do if not P(n-i) then S:=S union {i}: fi: od: 1/N*(add(pGenProp(P,k-1,n-i,N),i in S)): fi: end: # pkG(P,k,N) adds up pGenProp(P,init,i,j,N) for 1<=i<=k and k<=j<=kN pkG:=proc(P,k,N) local s,i,j: s:=0: for i from 1 to k do for j from i to N*i do s:=s+pGenProp(P,i,j,N): od: od: 1+s: end: #Paper(P,k,N): Given a property P, a positive integer k, and a positive integer N outputs #a paper with estimates for the expected duration, starting at 0, of hitting a number with property P #with a fair die with n faces, from n=2 to n=N. It assumes that you roll it up to k times. Try: #Paper(isprime,100,10): Paper:=proc(P,k,N) local n,t0: t0:=time(): print(`Suppose that you roll a fair die with n faces, starting at 0, and keep adding up the outcomes, until you reach an integer that has property`, convert(P,string) ): print(``): print(`approximating it with the assumtion that you do it up to`, k, `times , here are the averages from n=2 to n=N`): print(``): for n from 2 to N do print(`If the fair die has`, n, `faces , then the estimated expected duration is`, evalf(pkG(P,k,n))): print(``): od: print(``): print(`------------------------`): print(``): print(`This ends this paper that took`, time()-t0, `seconds to prepare`): print(``): end: