######################################################################
## HIT2gen.txt: Save this file as  HIT2gen.txt to use it,            #
# stay in the                                                        #
## same directory, get into Maple (by typing: maple <Enter> )        #
#then type: read `HIT2gen.txt`: <Enter>                              #
## Then follow the instructions given there                          #
##                                                                   #
## Written by   Lucy Martinez,                                       #
#and Doron Zeilbeerger, Rutgers University                           #
######################################################################

#VERSION OF Nov. 2024

with(combinat):
Digits:=200:


print(`First Written: Nov. 2024: tested for Maple 20  `):
print():
print(`This is HIT2gen.txt, A Maple package`):
print(`that includes all the functions of HIT2.txt that accompanied`):
print(` Lucy Martinez and Doron Zeilberger's  article: `):
print(` How many  Dice Rolls Would It Take to Reach Your Favorite Kind of Number? `):
print(`This extended package contains analogous functions for the statistics of how long it takes to  visit the primes (or any other property) for the k-th time`):
print(`Answering an interesting question of Yaakov Malinovsky`):
print(`The new functions are listed by typing ezraK();`):

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(`-------------------------------`):


print(`-------------------------------`):
print(`For a list of the Reaching in k functions type:  ezraK();`):
print(`For help with a specific procedure, type "ezra(procedure_name);"`):
print():
 print(`-------------------------------`):




ezraK:=proc()
if args=NULL then

print(` The reaching in k steps procedures are: OneTimeGk`):
print(` GFpCk, GFpK, OneTimeGk, SimuGFk, Story4 `):

else
 ezra(args):
fi:
end:

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]=GFpCk then
print(`GFpCk(n,L,init,t,P,k): Similar to GFpK(n,L,init,P,k) (q.v.) but conditioned on the life of the game (i.e. reacing the property P for the k-th time) taking <=L rolls, followed by the probability that you`):
print(`would NOT get your goal in <=L rolls. Try:`):
print(`GFpCk(6,100,0,t,isprime,2);`):

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]=GFpK then
print(`GFpK(n,L,init,t,P,k): inputs a positive integer n, corresponding to the number of faces in a fair die, outputs the prob. generating function in t,`):
print(`whose coefficient of t^i is the probability, starting with init,  of hitting for the k-th  time an integer with the property P after EXACTLY i times`):
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 only visited`):
print(`the desired property less than k times but the i-th partial sum does. Try:`):
print(`GFpK(6,100,0,t,isprime,3);`):

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]=OneTimeGk then
print(` OneTimeGk(P,init,F,k) starts with sum equal to an initial value (init), adds up `):
print(` the result from rolling a F-faced die, it keeps rolling the die while `):
print(` property P has not been met less than k times. As soon as the "sum" meets P for the k-th time, 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(`OneTimeGk(isprime,0,6,3);`):


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]=Story4 then
print(`Story4(P,init,L,R,F,K): Inputs a property P, an initial location init, a number of rounds L, `):
print(`a pos. integer R, and another pos. integer F and yet another pos. integer K`):
print(`outputs  an article with the statistics [Expectation, StandardDeviation,HigherScaledMomentsUpTo the R] for`):
print(`reaching an integer with property P, rolling an F-faced die at most L times`):
print(`and stopping as soon as visiting the set obeying P k times, for k from 1 to K. `):
print(`It also telling you the prob. of reaching the goal of visiting at least k times in <=L rolls of the dice`):
print(`and the conditional expectation (that is a good estimate for the the actual expectation with no limit to the number of rounds)`):
print(`and the s.d. followed by the moments up to the R-th. Try`):
print(`Story4(isprime,0,200,4,6,3);`):

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]=SimuGFk then
print(`SimuGFk(F,init,t,P,K,k) The probability generating function, in t,  by doing  OneTimeG(P,init,F,k) K times. Try:`):
print(` Try:Stat(SimuGFk(6,0,t,isprime,10000,2),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:



# OneTimeGk(P,init,F,k) 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 less than k times. As soon as the "sum" meets P for the k-th time, 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:
#OneTimeGk(isprime,0,6,3);

OneTimeGk:=proc(P,init,F,k) local s,i,c:
s:=init:
c:=0:
for i from 1 to 1000 while c<k do
 s:=s+rand(1..F)():
if P(s) then
 c:=c+1:
fi:
od:
i-1:
end:


# SimuGFk(F,init,t,P,K,k) The probability generating function, in t,  by doing  OneTimeG(P,init,F,k) K times. Try:
#Stat(SimuGFk(6,0,t,isprime,10000,2),t,4);
SimuGFk:=proc(F,init,t,P,K,k) local i:
 add(t^OneTimeGk(P,init,F,k),i=1..K)/K:
end:





#GFpK(n,L,init,t,P,k): 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 k-th  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 only visited
#the desired property less than k times but the i-th partial sum does. Try:
#GFpK(6,100,0,t,isprime,3);
GFpK:=proc(n,L,init,t,P,k) local g,i,f,F,f1,f2,x,i1,s,f1a:

if P(init) then
 f:=x^init*s:
else
 f:=x^init:
fi:

g:=add(x^i,i=1..n)/n:
F:=0:

for i from 1 to L do

 f:=expand(f*g):
 f1:=0:
 f2:=0:
 for i1 from ldegree(f,x) to degree(f,x) do
  if P(i1) then
   f1:=f1+coeff(f,x,i1)*x^i1:
 else
   f2:=f2+coeff(f,x,i1)*x^i1:
fi:
od:
f:=expand(s*f1+f2):


 f1a:=coeff(f,s,k):

 F:=expand(F+f1a*t^i):
 f:=expand(f-f1a*s^k):
od:
F:=subs(x=1,F):
collect(F,t):

end:



#GFpCk(n,L,init,t,P,k): Similar to GFpK(n,L,init,P,k) (q.v.) but conditioned on the life of the game (i.e. reacing the property P for the k-th time) taking <=L rolls, followed by the probability that you
#would NOT get your goal in <=L rolls. Try:
#GFpCk(6,100,0,t,isprime,2);
GFpCk:=proc(n,L,init,t,P,k) local f,c:

f:=GFpK(n,L,init,t,P,k):
c:=subs(t=1,f):
[f/c,1-c]:

end:





#Story4(P,init,L,R,F,K): Inputs a property P, an initial location init, a number of rounds L, 
#a pos. integer R, and another pos. integer F and yet another pos. integer K
#outputs  an article with the statistics [Expectation, StandardDeviation,HigherScaledMomentsUpTo the R] for
#reaching an integer with property P, rolling an F-faced die at most L times
#and stopping as soon as visiting the set obeying P k times, for k from 1 to K. 
#It also telling you the prob. of reaching the goal of visiting at least k times in <=L rolls of the dice
#and the conditional expectation (that is a good estimate for the the actual expectation with no limit to the number of rounds)
#and the s.d. followed by the moments up to the R-th. Try
#Story4(isprime,0,200,4,6,3);
Story4:=proc(P,init,L,R,F,K) local gu,t0,t,k,lu,LI:
t0:=time():
print(``):
print(`How long should it take until you  visit the integers with property`, convert(P,string), `for the k-th time if you roll a fair die with`, F, `faces for k from 1 to`, K):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
LI:=[]:
for k from 1 to K do
print(``):
print(`If you roll a fair die with`, F, `faces  starting at`, init, `and look at the running total and keep going until that running total has the property`, convert(isprime,string), `for `, k, `times `):
print(``):
gu:=GFpCk(F,L,init,t,P,k):
print(`The prob. of reaching that goal in <=`, L, `rounds is`, evalf(1-gu[2],20)):
print(``):
lu:=evalf(Stat(gu[1],t,R),20):
print(``):
LI:=[op(LI),lu[1]]:
print(`Conditioned that you are succesful the expected number of rounds for that goal is`, lu[1], ` (and this is a good estimate for unconditional expectation) and the standard deviation is`, lu[2]):
print(``):
if R>=3 then
print(`The skewness is`, lu[3]):
fi:
if R>=4 then
print(`The kurtosis is`, lu[4]):
fi:
if R>4 then
print(`The next scaled moments up to the`,R, `are `):
print([op(5..R,lu)]):
fi:
print(``):
print(`-----------------------`):
print(``):
od:
print(``):
print(`To sum up the expected time it takes to visit`, convert(P,string), `k times for k from 1 to `, K, `are `):
print(``):
lprint(evalf(LI,10)):
print(``):
print(`----------------`):
print(`This ends this paper that took`, time()-t0, `seconds to produce.`):
print(``):
end: