###################################################################### ## HIT2.txt: Save this file as HIT2.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, # #and Doron Zeilbeerger, Rutgers University # ###################################################################### #VERSION OF Jan. 20, 2023 with(combinat): Digits:=200: print(`First Written: Dec. 2022: tested for Maple 20 `): print(): print(`This is HIT2.txt, One of the Maple packages`): print(`accompanying 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/HIT2.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 "ezra();". For specific help type "ezra(procedure_name);" `): print(`-------------------------------`): print(`For a list of the supporting functions type: ezra1();`): print(`For help with a specific procedure, type "ezra(procedure_name);"`): print(): print(`-------------------------------`): print(`-------------------------------`): print(`For a list of the Story functions type: ezraS();`): print(`For help with a specific procedure, type "ezra(procedure_name);"`): print(): print(`-------------------------------`): ezraS:=proc() if args=NULL then print(` The Story procedures are: GFpGv, Story1, Story2, Story3, Story3a `): print(` `): else ezra(args): fi: end: ezra1:=proc() if args=NULL then print(` The supporting procedures are: IsDP, GFpC, plotDist, OneTimeG, SimuG, Stat `): print(` `): else ezra(args): fi: end: ezra:=proc() if args=NULL then print(` HIT2.txt: A Maple package for studying the statistics of the time it gets to reach your favorite kind of numbers rolling a die, Using a Symbolic-Computational approach `): print(`The main procedures are: BiStat, BiStatLo, EDerr, EstED, SimuGF, GFer, GFlo, GFp, GFpG, StatE, StatS `): elif nargs=1 and args[1]=BiStat then print(`BiStat(n,L,init,P,K): Inputs a positive integer n >=2 (number of faces), a positive integer, L, the number of rounds, a non-neg. integer, init (the starter), a property P, and a positive integer K`): print(`outputs the list`): print(`[surv,[AverageDuration, StandardDeviation, Skewness, Kurtosis],[AverageDestination(compared to init), StandardDeviation, Skewness, Kurtosis],correlation]. Try:`): print(`BiStat(6,100,0,isprime,4);`): elif nargs=1 and args[1]=BiStatLo then print(`BiStatLo(g,x,L,init,P,K): Inputs g, a prob. gen. function of a die, phrased in terms of x,`): print(`a positive integer, L, the number of rounds, a non-neg. integer, init (the starter), a property P, and a positive integer K`): print(`outputs the list`): print(`[surv,[AverageDuration, StandardDeviation, Skewness, Kurtosis],[AverageDestination(compared to init), StandardDeviation, Skewness, Kurtosis],correlation]. Try:`): print(`BiStatLo((x+x^2+x^3+x^4+x^5+x^6)/6,x,100,0,isprime,4);`): elif nargs=1 and args[1]=EDerr then print(`EDerr(n,init,P,er,L): The exact expected duration (for ever) until reaching property P after the minimal number of rounds (<=L), to guarantee termination with prob. 1-er, starting at init. Try:`): print(`EDerr(6,0,isprime,10^(-10),300);`): elif nargs=1 and args[1]=EstED then print(`EstED(n,init,P,er,L): A non-rigorous estimate of the expected duration (for ever) until reaching property P, starting at init, assuming that you only live up to L rounds. Returns L otherwise. Try:`): print(`EstED(6,0,isprime,10^(-10),200);`): elif nargs=1 and args[1]=ExtractP then print(`ExtractP(f,x,P): given a polynomial f in x, and a property P, outputs the part with exponents that satisy P. Try:`): print(`ExtractP(4*x^5+8*x^6,x,isprime);`): elif nargs=1 and args[1]=GFer then print(`GFer(n,init,t,P,er,L): inputs a positive integer n, corresponding to the number of faces in a fair die, `): print(`it also inputs a large integer L for the max. number of rounds that you allow.`): print(`outputs the number of`): print(`rounds that gurantees that you do NOT finish by then is less than er, followed by the cond. `): print(`prob. generating function in t, for number of terms until getting to your goal, conditions on finishing before that number`): print(`If the prob. of NOT finishing in <=L rounds is more than er than it returns FAIL. Try:`): print(`GFer(6,0,t,isprime,10^(-6),200);`): elif nargs=1 and args[1]=GFerG then print(`GFerG(n,init,t,x,P,er,L): inputs a positive integer n, corresponding to the number of faces in a fair die, `): print(`it also inputs a large integer L for the max. number of rounds that you allow.`): print(`outputs the number of`): print(`rounds that gurantees that you do NOT finish by then is less than er, followed by the cond. `): print(`prob. generating function in t and x, for number of terms until getting to your goal, conditions on finishing before that number`): print(`If the prob. of NOT finishing in <=L rounds is more than er than it returns FAIL. `): print(`The coeff. of t^i*x^j is the prob. that you reach integer j with property P after exactly i rounds. Try:`): print(`GFerG(6,0,t,x,isprime,10^(-6),200);`): elif nargs=1 and args[1]=GFlo then print(`GFlo(g,L,init,t,x,P): inputs a prob. generating function of a die expressed in terms of the variable x, outputs the bi-variate prob. generating function in t and x,`): print(`whose coefficient of t^i*x^j is the probability, starting with init, of hitting for the first time an integer with the property P after EXACTLY i times and that happening to be j`): print(`(i.e. the probability that if you roll a die, i times, starting at init, and you add-up the die outcomes, the first i-1 partial sums do not have property P`): print(`but the i-th partial sum does it happens to be j. Try:`): print(`GFlo((x+x^2+x^3+x^4+x^5+x^6)/6,100,0,t,x,isprime);`): elif nargs=1 and args[1]=GFp then print(`GFp(n,L,init,t,P): inputs a positive integer n, corresponding to the number of faces in a fair die, positive integers L and init, variable t and property P.`): print(`Outputs the prob. generating function in t,`): print(`whose coefficient of t^i is the probability of, starting at init, and adding up the die outcomes, of gettingg an integer with property P for the FIRST time after`): print(`exactly i die rolls, by rolling the die at most L times. In other words you are only allowed up to L rolls, or you are out of the game.Try:`): print(`GFp(6,100,0,t,isprime);`): elif nargs=1 and args[1]=GFpG then print(`GFpG(n,L,init,t,x,P): inputs a positive integer n, corresponding to the number of faces in a fair die, outputs the bi-variate prob. generating function in t and x,`): print(`whose coefficient of t^i*x^j is the probability, starting with init, of hitting for the first time an integer with the property P after EXACTLY i times and that happening to be j`): print(`(i.e. the probability that if you roll a die, i times, starting at init, and you add-up the die outcomes, the first i-1 partial sums do not have property P`): print(`but the i-th partial sum does it happens to be j. Try:`): print(`GFpG(6,100,0,t,x,isprime);`): elif nargs=1 and args[1]=GFpGv then print(`GFpGv(n,L,init,t,x,P): Verbose form, and more information of GFpG(n,L,init,t,x,P) (q.v.). Try:`): print(`GFpGv(6,5,0,t,x,isprime);`): elif nargs=1 and args[1]=GFpC then print(`GFpC(n,L,init,t,P): Similar to GFp(n,L,init,P) (q.v.) but conditioned on the life of the game taking <=L rolls, followed by the probability that you will NOT terminate in <=L rolls`): print(`GFpC(6,100,0,t,isprime);`): elif nargs=1 and args[1]=IsDP then print(`IsDP(n,k): Is n a product of k distinct primes? Try:`): print(`IsDP(6,2);`): 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]=plotDist then print(`plotDist(f,x,K): Given a prob. gen. function f(x) that has a `): print(`Taylor series`): print(`for a discrete r.v.`): print(`plots its normalized version (X-mu)/sig between mu-K1*sig andmu+K2*sig. For example, try:`): print(`plotDist((1+x)^40,x,5,5);`): elif nargs=1 and args[1]=Story1 then print(`Story1(P,init,er,L,R,N): Inputs a small er, a maximum number of rounds L and a positive integer N outputs a paper`): print(`with the statistics [Expectation, StandardDeviation,HigherScaledMomentsUpTo the R] for`): print(`reaching an integer with property P, rolling an n-faced die, with probability at least 1-er. It also outputs`): print(`the number of rounds guaranteeing it. Try:`): print(`Story1(isprime,0,10^(-6),200,4,20);`): elif nargs=1 and args[1]=Story2 then print(`Story2(P,init,L,R,N): Inputs a number of rounds L and a positive integer N outputs a paper`): print(`with the statistics `): print(`[error of surviving, [ExpectationDuration, StandardDeviation,HigherScaledMomentsUpTo the R],`): print(`[ExpectatedDisplacement, StandardDeviation,HigherScaledMomentsUpTo the R], condional correlation]`): print(`reaching an integer with property P, rolling an n-faced die, up to L rolls.`): print(`Story2(isprime,0,200,4,10);`): elif nargs=1 and args[1]=Story3 then print(`Story3(P,MaxInit,L,R,n): Fixing a property, P, a postive integer MaxInit, an allowed max. number of rounds, a positive integer R, and a positive integer n, outputs a paper`): print(`about starting with all integers from 0 to MaxInit that do NOT have property P, the prob. of NOT making it in that mumber of rows, followed by`): print(`the statistical information for duration, destination(shifted from the beginning) and correlation. Try:`): print(`Story3(isprime,10,200,4,6);`): elif nargs=1 and args[1]=Story3a then print(`Story3a(P,a,b,L,R,n): Fixing a property, P, a postive integer MaxInit, an allowed max. number of rounds, a positive integer R, and a positive integer n, outputs a paper`): print(`about starting with all integers that are powers of a from , a^1 to a^b, that do NOT have property P (if they do the first one that doesn't), the prob. of NOT making it in that mumber of rows, followed by`): print(`the statistical information for duration, destination(shifted from the beginning) and correlation. Try:`): print(`Story3a(isprime,10,10,200,4,6);`): 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);`): elif nargs=1 and args[1]=SimuGF then print(` SimuGF(F,init,t,P,K) The probability generating function, in t, by doing OneTimeG(P,init,F) K times. Try:`): print(`f:=SimuGF(6,0,t,isprime,10000); Stat(f,t,4);`): elif nargs=1 and args[1]=Stat then print(`Stat(F,t,K): Given a generating function F of t, for a r.v. finds the statistical information. Try:`): print(`Stat((1+t)^10000,4);`): elif nargs=1 and args[1]=StatE then print(`StatE(n,L,init,P,r): If you roll a fair n-faced die up to L times, starting at init, and looking at the cumulative sums and exiting as soon as the running sum`): print(`has property P. Assuming that you finished in <=L rolls, the expected time it took you, followed by the standard deviation, and then the scaled moments up to the`): print(`r-th. It also gives you the probability of NOT finishing in <=L rolls. Try:`): print(`StatE(6,100,0,isprime,6);`): elif nargs=1 and args[1]=StatS then print(`StatS(n,L,init,P,r): Simulating the game L times and getting the statistics (the first r moments) followed by the largest time. Try:`): print(`StatS(6,1000,0,isprime,6);`): else print(`There is no such thing as`, args): fi: end: with(combinat): with(NumberTheory): #plotDist(f,x,K): Given a prob. gen. function f(x) that has a #Taylor series #for a discrete r.v. #plots its normalized version (X-mu)/sig between mu-K1*sig and #mu+K2*sig #For example, try: #plotDist((1+x)^40,x,5,5); plotDist:=proc(f,x,K1,K2) local mu,f1,lu,gu,sig,i,j1: f1:=f/subs(x=1,f): mu:=subs(x=1,diff(f1,x)): gu:=f1/x^mu: sig:=sqrt(subs(x=1,x*diff(x*diff(gu,x),x))): lu:=taylor(f1,x=0,trunc(mu)+K2*trunc(sig)+10): lu:=[seq([i,coeff(lu,x,i)],i=max(0,trunc(mu-K1*sig))..trunc(mu+K2*sig))]: lu:=evalf([seq([(lu[j1][1]-mu)/sig,lu[j1][2]*sig],j1=1..nops(lu))]): plot(lu,scaling=constrained): 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: #Stat(F,t,K): Given a generating function F of t, for a r.v. finds the statistical information. Try: #Stat((1_t)^10000,4); Stat:=proc(F,t,K) local gu,mu,sd,k,f: f:=evalf(F/subs(t=1,F),20): mu:=subs(t=1,diff(f,t)): gu:=[mu]: if K=1 then RETURN(gu): fi: f:=f/t^mu: f:=t*diff(t*diff(f,t),t): sd:=sqrt(subs(t=1,f)): gu:=[op(gu),sd]: for k from 3 to K do f:=t*diff(f,t): gu:=[op(gu),subs(t=1,f)/sd^k]: od: evalf(gu,20): end: #IsDP(n,k): Is n a product of k distinct primes? IsDP:=proc(n,k) local gu,i: if n=0 or n=1 then RETURN(false): fi: gu:=ifactors(n)[2]: if nops(gu)=k and {seq(gu[i][2],i=1..nops(gu))}={1} then RETURN(true): else RETURN(false): fi: end: #ExtractP(f,x,P): given a polynomial f in x, and a property P, outputs the part with exponents that satisy P. Try #ExtractP(4*x^5+8*x^6,x,isprime); ExtractP:=proc(f,x,P) local f1,i: f1:=0: for i from ldegree(f,x) to degree(f,x) do if P(i) then f1:=f1+coeff(f,x,i)*x^i: fi: od: f1: end: #GFp(n,L,init,t,P): inputs a positive integer n, corresponding to the number of faces in a fair die, outputs the prob. generating function in t, #whose coefficient of t^i is the probability, starting with init, of hitting for the first time an integer with the property P after EXACTLY i times #(i.e. the probability that if you roll a die, i times, starting at init, and you add-up the die outcomes, the first i-1 partial sums do not have property P #but the i-th partial sum does. Try: #GFp(6,100,0,t,isprime); GFp:=proc(n,L,init,t,P) local g,i,f,F,f1,x: if P(init) then RETURN(1): fi: g:=add(x^i,i=1..n)/n: f:=x^init: F:=0: for i from 1 to L do f:=expand(f*g): f1:=ExtractP(f,x,P): F:=expand(F+f1*t^i): f:=f-f1: od: F:=subs(x=1,F): collect(F,t): end: # SimuGF(F,init,t,P,K) The probability generating function, in t, by doing OneTimeG(P,init,F) K times. Try: #Stat(SimuGF(isprime,0,6,10000,t),t,4); SimuGF:=proc(F,init,t,P,K) local i: add(t^OneTimeG(P,init,F),i=1..K)/K: end: #GFpC(n,L,init,t,P): Similar to GFp(n,L,init,P) (q.v.) but conditioned on the life of the game taking <=L rolls, followed by the probability that you #would NOT get your goal in <=L rolls. Try: #GFpC(6,100,0,t,isprime); GFpC:=proc(n,L,init,t,P) local f,c: f:=GFp(n,L,init,t,P): c:=subs(t=1,f): [f/c,1-c]: end: #StatE(n,L,init,P,r): If you roll a fair n-faced die up to L times, starting at init, and looking at the cumulative sums and exiting as soon as the running sum #has property P. Assuming that you finished in <=L rolls, the expected time it took you, followed by the standard deviation, and then the scaled moments up to the #r-th. It also gives you the probability of NOT finishing in <=L rolls. Try: #StatE(6,100,0,isprime,6); StatE:=proc(n,L,init,P,r) local gu,c,t: gu:=GFpC(n,L,init,t,P): c:=evalf(gu[2]): gu:=gu[1]: Stat(gu,t,r),c: end: #StatS(n,L,init,P,r): Simulating the game K times and getting the statistics followed by the largest time. Try #StatS(6,1000,0,isprime,4); StatS:=proc(n,L,init,P,r) local f,t,g: f:=SimuGF(n,init,t,P,L): g:=degree(f,t): Stat(f,t,r),g: end: #GFer(n,init,t,P,er,L): inputs a positive integer n, corresponding to the number of faces in a fair die, #it also inputs a large integer L for the max. number of rounds that you allow. #outputs the number of #rounds that gurantees that you do NOT finish by then is less than er, followed by the cond. #prob. generating function in t, for number of terms until getting to your goal, conditions on finishing before that number #If the prob. of NOT finishing in <=L rounds is more than er than it returns FAIL. Try #GFer(6,0,t,isprime,10^(-6),200); GFer:=proc(n,init,t,P,er,L) local g,i,f,F,f1,x: if P(init) then RETURN([0,1]): fi: g:=add(x^i,i=1..n)/n: f:=x^init: F:=0: for i from 1 to L while subs({x=1,t=1},F)<1-er do f:=expand(f*g): f1:=ExtractP(f,x,P): F:=expand(F+f1*t^i): f:=f-f1: od: if i=L+1 then RETURN(FAIL): else F:=evalf(subs(x=1,F)/subs({x=1,t=1},F)): F:=collect(F,t): RETURN([i-1,F]): fi: end: #Story1(P,init,er,L,R,N): Inputs a small er, a maximum number of rounds L and a positive integer N outputs a paper #with the statistics [Expectation, StandardDeviation,HigherScaledMomentsUpTo the R] for #reaching an integer with property P, rolling an n-faced die, with probability at least 1-er. It also outputs #the number of rounds guaranteeing it. Try: #Story1(isprime,0,10^(-6),200,4,20); Story1:=proc(P,init,er,L,R,N) local gu,n,t0,t: if P(init) then print(init, `already has property`, convert(P,string)): RETURN(FAIL): fi: t0:=time(): print(`Statistics of the Number of Rounds for Reaching an integer with property`, convert(P,string), `rolling a fair die with n faces for n from 2 to`, N): print(``): print(`By Shalsoh B. Ekhad `): print(``): print(`Suppose you start at`, init, `and roll a fair die with n faces, with 1 through n dots, and walk forward according to the number of dots, and stop as soon as you hit an integer with property`, convert(P,string)): print(``): print(`You are allowed up to`, L, `die-rolls, but you are interested in the smallest number of rolls that would guarantee that you would reach your goal with prob at least`, evalf(1-er)): print(``): print(`The table below tells you that number of rolls, followed by, conditioned on finishing by that number of rolls, the expected duration of such a game, followed by the standard-deviation and the scaled moments`): print(`up to the`, R): print(``): print(`with Fair Dice with number of faces from 2 to`, N): print(``): for n from 2 to N do gu:=GFer(n,init,t,P,er,L): if gu<>FAIL then print(`Number of faces=`,n, `Number of rounds=`, gu[1]): print(``): print( `[Ave,SD,ThirdScaledMoment,Etc]= `, Stat(gu[2],t,R)): print(``): print(`---------------------`): else print(`With the number of faces=`,n, `You need more than`, L, `rounds to guarantee that you would make it to your goal with probability more than`, evalf(1-er)): fi: od: print(`This ends this paper that took`, time()-t0, `seconds to produce. `): end: #EDerr(n,init,P,er,L): The conditional expected duration until reaching property P, starting at init. Try: #EDerr(6,0,isprime,10^(-10),200); EDerr:=proc(n,init,P,er,L) local f,t: f:=GFer(n,init,t,P,er/10^10,L): if f=FAIL then RETURN(FAIL): fi: f:=f[2]: Stat(f,t,1)[1]: end: #EstED(n,init,P,er,L): A non-rigorous estimate of the expected duration until reaching property P, starting at init. Try: #EstED(6,0,isprime,10^(-10),200); EstED:=proc(n,init,P,er,L) local lu1,lu2,ku: lu1:=EDerr(n,init,P,er,L): if lu1=FAIL then RETURN(FAIL): fi: lu2:=EDerr(n,init,P,er/10^5,L): if lu2=FAIL then RETURN(FAIL): fi: ku:=trunc(-log[10](abs(lu1-lu2))): evalf(lu2,ku-2): end: #GFerG(n,init,t,x,P,er,L): inputs a positive integer n, corresponding to the number of faces in a fair die, #it also inputs a large integer L for the max. number of rounds that you allow. #outputs the number of #rounds that gurantees that you do NOT finish by then is less than er, followed by the cond. #prob. generating function in t and x, for number of terms until getting to your goal, conditions on finishing before that number #If the prob. of NOT finishing in <=L rounds is more than er than it returns FAIL. #The coeff. of t^i*x^j is the prob. that you reach integer j with property P after exactly i rounds. Try: #GFerG(6,0,t,x,isprime,10^(-6),200); GFerG:=proc(n,init,t,x,P,er,L) local g,i,f,F,f1: if P(init) then RETURN([0,1]): fi: g:=add(x^i,i=1..n)/n: f:=x^init: F:=0: for i from 1 to L while subs({x=1,t=1},F)<1-er do f:=expand(f*g): f1:=ExtractP(f,x,P): F:=expand(F+f1*t^i): f:=f-f1: od: if i=L+1 then RETURN(FAIL): else F:=F/subs({x=1,t=1},F): F:=collect(F,t): RETURN([i-1,F]): fi: end: #GFpG(n,L,init,t,x,P): inputs a positive integer n, corresponding to the number of faces in a fair die, outputs the bi-variate prob. generating function in t and x, #whose coefficient of t^i*x^j is the probability, starting with init, of hitting for the first time an integer with the property P after EXACTLY i times and that happening to be j #(i.e. the probability that if you roll a die, i times, starting at init, and you add-up the die outcomes, the first i-1 partial sums do not have property P #but the i-th partial sum does it happens to be j. Try: #GFpG(6,100,0,t,x,isprime); GFpG:=proc(n,L,init,t,x,P) local g,i,f,F,f1: if P(init) then RETURN(1): fi: g:=add(x^i,i=1..n)/n: f:=x^init: F:=0: for i from 1 to L do f:=expand(f*g): f1:=ExtractP(f,x,P): F:=expand(F+f1*t^i): f:=f-f1: od: collect(F,t): end: #GFpGv(n,L,init,t,x,P): A verbose version of GFpG(n,L,init,t,x,P) (q.v.). Try: #GFpGv(6,100,0,t,x,isprime); GFpGv:=proc(n,L,init,t,x,P) local g,i,f,F,f1,err,gu1,gu2,ku: if P(init) then RETURN(1): fi: g:=add(x^i,i=1..n)/n: f:=x^init: print(`We roll a fair die with`, n, `faces, so the probability generating function with respect to the outcome is`): print(``): print(g): print(``): print(`and in Maple notation`): print(``): lprint(g): print(``): print(`We start at`, init, `and we can play at most`, L, `rounds (life is finite). As soon as the running sum has property`, convert(P,string), `quit. this is called the duration of the game `): print(``): print(`We want to construct the bivariate probability generating function in x and t whose coeff. of t^i*x^j is the probability that the game ended after i rounds, hitting integer j`): print(``): F:=0: print(`At round 0 the prob. generating function is`, f): for i from 1 to L do print(`We are at round`, i, `multiplying the previous partial generating function by`,g, ` gives `): print(``): f:=expand(f*g): print(``): lprint(f): print(``): print(`We extract those monomials whose powers have property`, convert(P,string), `The new inductees of round`, i, `are `): f1:=ExtractP(f,x,P): print(``): lprint(f1): print(``): print(`We add it to the current generating function getting that the current generating function, a polynomial in x and t of degree`, i, ` in t is `): F:=expand(F+f1*t^i): print(``): lprint(F): print(``): f:=f-f1: print(`We subtract those monomials that correpond to the new arrivals of round`, i, `and are left with the survivors, those with monomials that have powers that do not have property`, convert(P,string)): print(``): print(`The generating function of the survivors after round`, i, ` is `): print(``): lprint(f): od: F:=expand(F/x^init): F:=collect(F,t): print(``): print(`After we finished, i.e. after `, L, `rounds , the generating function is, where now we look at the dispalment from`, init): print(``): lprint(F): print(``): print(`and in floating point`): print(``): lprint(evalf(F,10)): print(``): err:=1-subs({x=1,t=1},F): print(``): print(`The probability of not finishing the game in <=`, L, `rounds is`, evalf(err,10)): print(``): print(`The bi-variate prob. generating function conditioned on finishing in <=`, L, `rounds is`): print(``): F:=evalf(F/(1-err),10): print(``): lprint(F): print(``): gu1:=Stat(subs(x=1,F),t,2): lprint(``): print(`The conditional expectation of the number of rounds is`, evalf(gu1[1],10), `and the variance is`, evalf(gu1[2]^2,10)): lprint(``): gu2:=Stat(subs(t=1,F),x,2): lprint(``): print(`The conditional expectation of the final destination minus the initial value `,init, ` is `, evalf(gu2[1],10), `and the variance is`, evalf(gu2[2]^2,10)): lprint(``): print(`Finally the conditional correlation between the number of rounds (time) and space (location) is `): lprint(``): ku:=subs({x=1,t=1},diff(diff(F,x),t)): ku:=ku-gu1[1]*gu2[1]: ku:=ku/(gu1[2]*gu2[2]): lprint(evalf(ku,10)): print(``): end: #BiStat(n,L,init,P,K): Inputs a positive integer n >=2 (number of faces), a positive integer, L, the number of rounds, a non-neg. integer, init (the starter), a property P, and a positive integer K #outputs the list #[surv,[AverageDuration, StandardDeviation, Skewness, Kurtosis],[AverageDestination(compared to init), StandardDeviation, Skewness, Kurtosis],correlation]. Try: #BiStat(6,100,0,isprime,4); BiStat:=proc(n,L,init,P,K) local F,t,x,err,gu1,gu2,ku: F:=GFpG(n,L,init,t,x,P): err:=evalf(1-subs({x=1,t=1},F)): F:=evalf(F/(1-err)): F:=expand(F/x^init): gu1:=evalf(Stat(subs(x=1,F),t,K),20): gu2:=evalf(Stat(subs(t=1,F),x,K),20): ku:=subs({x=1,t=1},diff(diff(F,x),t)): ku:=ku-gu1[1]*gu2[1]: ku:=ku/(gu1[2]*gu2[2]): ku:=evalf(ku,10): [evalf(err,10),gu1,gu2,ku]: end: #Story2(P,init,L,R,N): Inputs a number of rounds L and a positive integer N outputs a paper #with the statistics #[error of surviving, [ExpectationDuration, StandardDeviation,HigherScaledMomentsUpTo the R], #[ExpectatedDisplacement, StandardDeviation,HigherScaledMomentsUpTo the R], condional correlation] #reaching an integer with property P, rolling an n-faced die, up to L rolls. #Story2(isprime,0,200,4,10); Story2:=proc(P,init,L,R,N) local gu,n,t0,t: if P(init) then print(init, `already has property`, convert(P,string)): RETURN(FAIL): fi: t0:=time(): print(`Bi Statistics of the Number of Rounds, and Final Destination, for Reaching an integer with property`, convert(P,string), `rolling a fair die with n faces for n from 2 to`, N, `starting at `, init): print(``): print(`By Shalsoh B. Ekhad `): print(``): print(`Suppose you start at`, init, `and roll a fair die with n faces, with 1 through n dots, and walk forward according to the number of dots, and stop as soon as you hit an integer with property`, convert(P,string)): print(``): print(`You are allowed up to`, L, `die-rolls `): print(``): print(`The table below tells you the probability of not finishing in <=`, L, `rounds followed by, conditioned on finishing by that number of rolls, the expected duration of such a game, followed by the standard-deviation and the scaled moments`): print(`followed by the analogous displacement of final destination, followed by the condional correlation between duration and destination`): print(`up to the`, R, `moment `): print(``): print(`with Fair Dice with number of faces from 2 to`, N): print(``): for n from 2 to N do gu:=BiStat(n,L,init,P,R): print(`Number of faces=`,n, ` prob. of not finishing=`,gu[1]): print(``): print( `[AveDuration,SD,ThirdScaledMoment,Etc]= `, gu[2]): print(``): print( `[AveDisplacement,SD,ThirdScaledMoment,Etc]= `, gu[3]): print(``): print(`conditioned correlation=`, gu[4]): print(``): print(`---------------------`): od: print(`This ends this paper that took`, time()-t0, `seconds to produce. `): end: #Story3(P,MaxInit,L,R,n): Fixing a property, P, a postive integer MaxInit, an allowed max. number of rounds, a positive integer R, and a positive integer n, outputs a paper #about starting with all integers from 0 to MaxInit that do NOT have property P, the prob. of NOT making it in that mumber of rows, followed by #the statistical information for duration, destination(shifted from the beginning) and correlation. Try: #Story3(isprime,10,200,4,6); Story3:=proc(P,MaxInit,L,R,n) local init,t0: t0:=time(): print(`Bi Statistics of the Number of Rounds, and Final Destination (from start), for Reaching an integer with property`, convert(P,string), `rolling a fair die with`, n , `faces starting at all integers from 0 to`, MaxInit): print(`that do not have property `, convert(P,string) ): print(``): print(`By Shalsoh B. Ekhad `): print(``): print(`Suppose you start somewhere and roll a fair die with`, n , `faces and walk forward according to the number of dots, and stop as soon as you hit an integer with property`, convert(P,string)): print(`for all starting integers from 1 to`, MaxInit, `that do not have property`, convert(P,string)): print(``): print(`You are allowed up to`, L, `die-rolls `): print(``): print(`The table below tells you the probability of not finishing in <=`, L, `rounds followed by, conditioned on finishing by that number of rolls, the expected duration of such a game, followed by the standard-deviation and the scaled moments`): print(`followed by the analogous displacement of final destination, followed by the condional correlation between duration and destination`): print(`up to the`, R, `moment `): print(``): print(`with Fair Dice with number of faces`, n): print(``): for init from 0 to MaxInit do if not P(init) then print(`If you start at`, init, `then the prob. of not finishing in <=`, L, `rounds followed by the statistical information ([Ave,S.D, skewness, kurtosis,..]) for duration and (relative progress) followed by the correlation are`): print(``): lprint(BiStat(n,L,init,P,R)): print(``): fi: od: print(``): print(`-----------------------------------`): print(``): print(`This took`, time()-t0, `seconds. `): end: #Story3a(P,a,b,L,R,n): Fixing a property, P, a postive integer MaxInit, an allowed max. number of rounds, a positive integer R, and a positive integer n, outputs a paper #about starting with all integers that are powers of a from , a^1 to a^b, that do NOT have property P (if they do the first one that doesn't), the prob. of NOT making it in that mumber of rows, followed by #the statistical information for duration, destination(shifted from the beginning) and correlation. Try: #Story3a(isprime,10,10,200,4,6); Story3a:=proc(P,a,b,L,R,n) local init,t0,b1,init1: t0:=time(): print(`Bi Statistics of the Number of Rounds, and Final Destination (from start), for Reaching an integer with property`, convert(P,string), `rolling a fair die with`, n , `faces starting at all powers of`, a, `up to the `, b): print(`that do not have property `, convert(P,string), `or as close as possible ` ): print(``): print(`By Shalsoh B. Ekhad `): print(``): print(`Suppose you start somewhere and roll a fair die with`, n , `faces and walk forward according to the number of dots, and stop as soon as you hit an integer with property`, convert(P,string)): print(`for all starting integers that are powers of`, a, ` up the `, b, `power that do not have property`, convert(P,string)): print(``): print(`You are allowed up to`, L, `die-rolls `): print(``): print(`The table below tells you the probability of not finishing in <=`, L, `rounds followed by, conditioned on finishing by that number of rolls, the expected duration of such a game, followed by the standard-deviation and the scaled moments`): print(`followed by the analogous displacement of final destination, followed by the condional correlation between duration and destination`): print(`up to the`, R, `moment `): print(``): print(`with Fair Dice with number of faces`, n): print(``): for b1 from 0 to b do init:=a^b1: if not P(init) then init1:=init: else for init1 from init while P(init1) do od: fi: print(`If you start at`, init1, `then the prob. of not finishing in <=`, L, `rounds followed by the statistical information ([Ave,S.D, skewness, kurtosis,..]) for duration and (relative progress) followed by the correlation are`): print(``): lprint(BiStat(n,L,init1,P,R)): print(``): od: print(``): print(`-----------------------------------`): print(``): print(`This took`, time()-t0, `seconds. `): end: #GFlo(g,L,init,t,x,P): inputs a prob. generating function of a die expressed in terms of the variable x, outputs the bi-variate prob. generating function in t and x, #whose coefficient of t^i*x^j is the probability, starting with init, of hitting for the first time an integer with the property P after EXACTLY i times and that happening to be j #(i.e. the probability that if you roll a die, i times, starting at init, and you add-up the die outcomes, the first i-1 partial sums do not have property P #but the i-th partial sum does it happens to be j. Try: #GFlo((x+x^2+x^3+x^4+x^5+x^6)/6,100,0,t,x,isprime); GFlo:=proc(g,L,init,t,x,P) local i,f,F,f1,i1: if min(seq(coeff(g,x,i1),i1=0..degree(g,x)))<0 or max(seq(coeff(g,x,i1),i1=0..degree(g,x)))>1 then print(g, `is not a prob. dist.`): RETURN(FAIL): fi: if subs(x=1,g)<>1 then print(g, ` is not a prob. dist.`): RETURN(FAIL): fi: if P(init) then RETURN(1): fi: f:=x^init: F:=0: for i from 1 to L do f:=expand(f*g): f1:=ExtractP(f,x,P): F:=expand(F+f1*t^i): f:=f-f1: od: collect(F,t): end: #BiStatLo(g,x,L,init,P,K): Inputs g, a prob. gen. function of a die, phrased in terms of x, #a positive integer, L, the number of rounds, a non-neg. integer, init (the starter), a property P, and a positive integer K #outputs the list #[surv,[AverageDuration, StandardDeviation, Skewness, Kurtosis],[AverageDestination(compared to init), StandardDeviation, Skewness, Kurtosis],correlation]. Try: #BiStatLo((x+x^2+x^3+x^4+x^5+x^6)/6,x,100,0,isprime,4); BiStatLo:=proc(g,x,L,init,P,K) local F,t,err,gu1,gu2,ku,i1: if min(seq(coeff(g,x,i1),i1=0..degree(g,x)))<0 or max(seq(coeff(g,x,i1),i1=0..degree(g,x)))>1 then print(g, `is not a prob. dist.`): RETURN(FAIL): fi: if subs(x=1,g)<>1 then print(g, ` is not a prob. dist.`): RETURN(FAIL): fi: F:=GFlo(g,L,init,t,x,P): err:=evalf(1-subs({x=1,t=1},F)): F:=evalf(F/(1-err)): F:=expand(F/x^init): gu1:=evalf(Stat(subs(x=1,F),t,K),20): gu2:=evalf(Stat(subs(t=1,F),x,K),20): ku:=subs({x=1,t=1},diff(diff(F,x),t)): ku:=ku-gu1[1]*gu2[1]: ku:=ku/(gu1[2]*gu2[2]): ku:=evalf(ku,10): [evalf(err,10),gu1,gu2,ku]: end: