######################################################################
##Bitcoin.txt: Save this file as  Bitcoin.txt                        #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read Bitcoin.txt<Enter>                            #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#DoronZeil at gmail dot com                                          #
######################################################################
 
#Created: July-Aug. 2018
#Version of  Sept. 5, 2018
 
print(`Created:  July-Aug.  12018. First posted: Aug. 15, 2018`):
print(` This is Bitcoin.txt `):
print(`It the package that accompanies the article `):
print(` A Combinatorial-Probabilstic Analysis of BitCoin Attacks`):
print(`by Evangelos Georgiadis and Doron Zeilberger`):

print(`The most current version of this package  is available from the following url`):
print(` http://sites.math.rutgers.edu/~zeilberg/tokhniot/BitCoin.txt `):
print(`and the most current version of the  article is available from the following url`):
print(`http://sites.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/bitcoin.html`):
print(``):
print(`Please report bugs to DoronZeil at gmail dot com `):
print(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://sites.math.rutgers.edu/~zeilberg/  .`):


print(`---------------------------------------`):
 print(`For a list of the Gambler's Ruin procedures type ezraGR();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the Supporting procedures type ezra1();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the GENERALIZED (with finite credit) procedures type ezraG();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the Simulations procedures type ezraS();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the MAIN procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

with(combinat):


ezraGR:=proc()

if args=NULL then
 print(` The Gambler;s ruin procedures are:  FnXq, fzq,  GRmoms, GRmomsC,GRmomsCs, GRp, GRpN,ZXq `):
 print(``):

else
ezra(args):
fi:

end:

ezraS:=proc()

if args=NULL then
 print(` The Simulations procedures are: BAsim, BAsimS, BAsim1, coin, GRsim, GRsim1, Phase1,  `):
 print(``):

else
ezra(args):
fi:

end:

ezraG:=proc()

if args=NULL then
 print(` The generalized procedures are: rMeniG `):
 print(``):

else
ezra(args):
fi:

end:

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: Alpha, AsyExpEmone ,  AsyExpM2mone, CheckAsyExpE, CheckAsyExpM2, CheckAsyExpVar, EssayLaTex,`):
 print(`EstC1,  Ope, OpeE, OpeGF,  OpeM2, `):
 print(` OpeM21, rMeniEmone, rMeniMOMmone, rSBEEmone, rSBEMOM2, rSBEMOM2mone, rMeniGF1 `):
 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: AsyExp, AsyExpE, AsyExpM2, AsyExpVar,CutOff, CutOffStory, Essay, rGP, rMeni, rMeniE, `):
 print(` rMeniEseq, rMeniGF, rMeniMOM, rMeniSeq, rSBE, rSBEE, rSBEEseq, rSBEseq `):
 print(` `):



elif nops([args])=1 and op(1,[args])=Alpha then
print(`Alpha(f,x,N): Given a probability generating function`):
print(`f(x) (a polynomial, rational function and possibly beyond) `):
print(`returns a list, of length N, whose `):
print(` (i) First entry is the average `):
print(` (ii): Second entry is the variance `):
print(` for i=3...N, the i-th entry is the so-called alpha-coefficients `):
print(` that is the i-th moment about the mean divided by the `):
print(`variance to the power i/2 (For example, i=3 is the skewness`):
print(`and i=4 is the Kurtosis) `):
print(` If f(1) is not 1, than it first normalizes it by dividing `):
print(`by f(1) (if it is not 0) .`):
print(`For example, try:`):
print(`Alpha(((1+x)/2)^100,x,4);`):

elif nops([args])=1 and op(1,[args])=AsyExp then
print(`AsyExp(q,n,m): the asymptotic expansion of the Rosenfeld polynomial when q<1/2`):

elif nops([args])=1 and op(1,[args])=AsyExpE then
print(`AsyExpE(q,n,m): the asymptotic expansion of rMeniE(n,q) (alias rSBEE(n,q) to order m when q<1/2. Try:`):
print(`AsyExpE(q,n,4);`):


elif nops([args])=1 and op(1,[args])=AsyExpEmone then
print(`AsyExpEmone(q,n,m): the asymptotic expression for the Numerator of the expected duration, up to a constant`):

elif nops([args])=1 and op(1,[args])=AsyExpM2 then
print(`AsyExpM2(q,n,m): the asymptotic expansion of rMeniMOM(n,q,2) (alias rSBEMOM2(n,q) to order m when q<1/2. Try:`):
print(`AsyExpM2(q,n,4);`):

elif nops([args])=1 and op(1,[args])=AsyExpM2mone then
print(`AsyExpM2mone(q,n,m): the asymptotic expansion of the numerator of rMeniMOMmone(n,q,2) to order m when q<1/2`):

elif nops([args])=1 and op(1,[args])=AsyExpVar then
print(`AsyExpVar(q,n,m): the asymptotic expansion of the variance of conditional duration. Try:`):
print(`AsyExpVar(q,n,4);`):

elif nops([args])=1 and op(1,[args])=BAsim then
print(`BAsim(n,q,N,M,x): simulates M times BAsim1(n,q,N). q is the prob. of an attacker succedding, z is his current debt, N is his maximuam credit`):
print(`x is a variable for the p.g.f for duration of phase 2, in case the attacker succeeds. The output is a triple`):
print(` [p,av,pgf] where p is the ratio (out of the M tries) where the attacker succeeded, av is the average duration it took him  `):
print(` (when he did succeed) and pgf is the polynomial. Try: `):
print(` BAsim(5,2/5,100,100,x); `):

elif nops([args])=1 and op(1,[args])=BAsim1 then
print(`BAsim1(n,q,N): simulates ONE bitcoin attack with prob. q (q<1/2) of success when the attacker started out with one block`):
print(`and the honest player continues until he gathers n good blcoks N (large) is the credit limit of the dishonest`):
print(`attacker. The output is the pair`):
print(`[1,LengthOfSecondPhase,m], in case of success [0,LengthOfSecondPhase,m], in case of failure. `):
print(`m is the number of blocks the attacker got. Try: `):
print(`BAsim1(10,2/5,200);`):

elif nops([args])=1 and op(1,[args])=BAsimS then
print(`BAsimS(n,q,N,M):  Similiar to BAsim(n,q,N,M,x) (q.v.), but only returns`):
print(`the average in case of success, and the average time it took there`):
print(`it also returns the true prob. of success and true expectation`):
print(`The output is`):
print(`[[SampleProbabiltyOfSuccess,TrueProbability], [SampleAverage, TrueExpectation],sd].`):
print(`where sd is the theoretical standard deviation . `):
print(`Try: `):
print(`BAsimS(5,2/5,100,100); `):

elif nops([args])=1 and op(1,[args])=CheckAsyExpE then
print(`CheckAsyExpE(q,m,M): Checks the asymptotic expansion to order m of the (conditional) expected duration with q<1/2`):
print(`in other words AsyExpE(q,n,m) up to the M-th term. Try:`):
print(`CheckAsyExpE(2/5,4,1000);`):

elif nops([args])=1 and op(1,[args])=CheckAsyExpM2 then
print(`CheckAsyExpM2(q,m,M): Checks the asymptotic expansion to order m of the (conditional) second moment of duration with q<1/2`):
print(`in other words AsyExpM2(q,n,m) up to the M-th term. Try:`):
print(`CheckAsyExpM2(2/5,4,1000);`):

elif nops([args])=1 and op(1,[args])=CheckAsyExpVar then
print(`CheckAsyExpVar(q,m,M): Checks the asymptotic expansion to order m of the (conditional) Variance of duration with q<1/2`):
print(`in other words AsyExpM2(q,n,m) up to the M-th term. Try:`):
print(`CheckAsyExpVar(2/5,4,1000);`):


elif nops([args])=1 and op(1,[args])=CheckAsyExpE then
print(`CheckAsyExpE(q,m,M): Checks the asymptotic expansion to order m of the (conditional) expected duration with q<1/2`):
print(`in other words AsyExpE(q,n,m) up to the M-th term. Try:`):
print(`CheckAsyExpE(2/5,4,1000);`):

elif nops([args])=1 and op(1,[args])=coin then
print(`coin(q): inputs a rational mumber q between 0 and 1, outputs 1 with prob. q and -1 with prob. 1-q. Try:`):
print(`coin(1/3);`):

elif nops([args])=1 and op(1,[args])=CutOff then
print(`CutOff(q,T): the smallest n such that rMeni(n,q)<T. It returns the actual probability, and the expected duration. Try:`):
print(`CutOff(0.45, 10^(-3));`):

elif nops([args])=1 and op(1,[args])=CutOffStory then
print(`CutOffStory(qs,qe,N,K): the cut-offs staring with qs and ending with qe, with N intermediate stops with`):
print(`confidence level 10^(-k), for k from 1 to K. Try:`):
print(`CutOffStory(0.3,0.4,2,3);`):

elif nops([args])=1 and op(1,[args])=Essay then
print(`Essay(q,n,M): inputs symbols q (for probabality of success of an attacker), n (for the number of good blocks generated by the honest party),`):
print(`and a positive integer M (for the order of the asymptotic expansions) and outputs an article about bitcoin atacks. Try:`):
print(`Essay(q,n,4):`):

elif nops([args])=1 and op(1,[args])=EssayLaTex then
print(`EssayLaTex(q,n,M): like Essay(q,n,M) (q.v.) but with latex instead of lprint. Try: `):
print(`EssayLaTex(q,n,4):`):

elif nops([args])=1 and op(1,[args])=EstC1 then
print(`EstC1(q,M,m): estimates the leading constant in the Expected duration, if successful, using`):
print(`the asymptotic expansion to order m, and the M-th term of rSBEE(M,q). M must be at least 1000. Try:`):
print(`EstC1(1/10,1000,4);`):

elif nops([args])=1 and op(1,[args])=FnXq then
print(`FnXq(n,X,q): the conditional prob. gen. function, in X,  of "duration of game" if succeeded in a gambler's ruin`):
print(`game starting at n with prob. q of going down and prob. 1-q going up, conditioned on the fact that it made it to 0.`):
print(`Try : `):
print(`FnXq(n,X,q); `):

elif nops([args])=1 and op(1,[args])=fzq then
print(`fzq(z,q): In a gambler's ruin, `):
print(` the expected time it takes to 0 starting at z in [0,infinity] with probability q of going left, conditioned on it winding at 0.`):
print(`rather than infinity.  Try: `):
print(`fzq(z,q); `):





elif nops([args])=1 and op(1,[args])=GRmoms then
print(`GRmoms(q,n,K): exact expressions for the (uncentraluized) moments, up to the the K-th in terms of q and n`):
print(`for the random variable "duration"  for a gambler's ruin game starting at n, with infinite credit,`):
print(`with probability q of going left (q<1/2) conditioned on it making it to 0. Try:`):
print(`GRmoms(q,n,4);`):
   

elif nops([args])=1 and op(1,[args])=GRmomsC then
print(`GRmomsC(q,n,K): exact expressions for the Centralized moments (i.e. moments about the mean), up to the the K-th in terms of q and n`):
print(`for the random variable "duration"  for a gambler's ruin game starting at n, with infinite credit,`):
print(`with probability q of going left (q<1/2) conditioned on it making it to 0. The first entry is the average, the second the variance. Try:`):
print(`GRmomsC(q,n,4);`):

elif nops([args])=1 and op(1,[args])=GRmomsCs then
print(`GRmomsCs(q,n,K): exact expressions for the Scaled Centralized moments (i.e. moments about the mean), up to the the K-th in terms of q and n`):
print(`for the random variable "duration"  for a gambler's ruin game starting at n, with infinite credit,`):
print(`with probability q of going left (q<1/2) conditioned on it making it to 0. The first entry is the average, the second the variance. `):
print(` The third is the skewness, the fourth, the kurtosis. Try:`):
print(`GRmomsCs(q,n,4);`):

elif nops([args])=1 and op(1,[args])=GRp then
print(`GRp(q,z): the probability of the dishonet attacker who currently has a deficit of z of catching up if`):
print(`his probability of winning a round is q (q<1/2) with unlimited credit. Try:`):
print(`GRp(.4,5);`):

elif nops([args])=1 and op(1,[args])=GRpN then
print(`GRpN(q,z,N): the probability of the dishonet attacker who currently has a deficit of z of catching up if`):
print(`his probability of winning a round is q ( if q<1/2) and his credit is limited to N. Try:`):
print(`GRpN(.4,5,100);`):
print(`GRpN(q,2,10);`):

elif nops([args])=1 and op(1,[args])=GRsim then
print(`GRsim(q,z,N,M,x): simulates M times GRsim1(q,z,N). q is the prob. of an attacker succedding, z is his current debt, N is his maximuam credit`):
print(` x is a variable for the p.g.f for duration in case the attacker succeeds. The output is a pair of triples followed by a number`):
print(` [ [p,av1,pgf1] , [1-p,av2,pgf2],av] `):
print(`where p is the ratio (out of the M tries) where the attacker succeeded, av is the average duration it took him  `):
print(` (when he did succeed), av2 is the average time it took until he failed,  and pgf1, pgf2  are the respective p.g.f. for duration. Try: `):
print(`av is the average time until the end, regardliess. Try: `):
print(`GRsim(11/25,5,20,100,x);`):

elif nops([args])=1 and op(1,[args])=GRsim1 then
print(`GRsim1(q,z,N): simulates ONE gambler's ruin with current debt of z prob. q of gaining a dollar and prob. 1-q of losing one`):
print(`quits when there is a debt of N. Outputs the pair [1,Duration] or [0,Duration] in case of success and failure respectively.`):
print(`Try:`):
print(`GRsim1(2/5,5,100);`):

elif nops([args])=1 and op(1,[args])=Ope then
print(`Ope(q,n,N): the second order recurrence operator annihilating the sequence of Rosenfeld polynomials (prob. of success)`):
print(`Try: Ope(q,n,N);`):

elif nops([args])=1 and op(1,[args])=OpeE then
print(`OpeE(q,n,N): the second order recurrence operator annihilating the sequence of Numetors of the `):
print(`expected time it takes to achieve success of attack (if realized)`):
print(`with prob. q (<1/2) of winning and  n good blocks. `):
print(`Try: OpeE(q,n,N);`):

elif nops([args])=1 and op(1,[args])=OpeGF then
print(`OpeGF(q,z,n,N): the second order recurrence operator annihilating the sequence of prob. gen, function for duration of`):
print(`a succesful attack, if the prob. of success at a single block is q, and the honest person gathered n good blocks. Try:`):
print(`OpeGF(q,z,n,N);`):

elif nops([args])=1 and op(1,[args])=OpeM21 then
print(`OpeM21(q,n,N): the second order recurrence operator (but in N^(-1)) annihilating the sequence of Numerators of the second moment`):

elif nops([args])=1 and op(1,[args])=Phase1 then
print(`Phase1(n,q): inputs a positive integer n and a prob. q ( a rational number) tosses a loaded coin with prob. 1-q of getting a 1`):
print(` and prob. q of getting a -1 `):
print(` keeps playing until there are  n 1's and outputs the the number of -1's. Try `):
print(` Phase1(10,2/5); `):

elif nops([args])=1 and op(1,[args])=rMeni then
print(`rMeni(n,q): the probability of succesful attack  by a dishonet person`):
print(`if the honest party mines n and the probabity of its sucess  is q (and the prob. of failure (i.e. success for the honest person) is 1-q.`):
print(`Using Eq (1) in Merni Rosenfeld's paper "Analysis of hashrate-based double-spending". Try: `):
print(`rMeni(5,1/3); `):
print(` rMeni(5,q); `):


elif nops([args])=1 and op(1,[args])=rMeniE then
print(`rMeniE(n,q): the expected time it takes for the attack by a dishonet person provided he is succesful`):
print(`if the honest party mines n and the probabity of its sucess  is q (and the prob. of failure (i.e. success for the honest person) is 1-q.`):
print(`Using the analog (for expected duration) of Eq (1) in Merni Rosenfeld's paper "Analysis of hashrate-based double-spending". Try: `):
print(`rMeniE(5,1/3); `):
print(` rMeniE(5,q); `):

elif nops([args])=1 and op(1,[args])=rMeniEmone then
print(`rMeniEmone(n,q): the NUMERATOR of the rMeniE(n,q) (q.v.) . Try: `):
print(` rMeniEmone(5,q); `):

elif nops([args])=1 and op(1,[args])=rMeniEseq then
print(`rMeniEseq(N,q): the first N terms of the expected duration, rMeniE(n,q) (q.v.) , done directly. Try:`):
print(` rMeniEseq(5,q); `):

elif nops([args])=1 and op(1,[args])=rMeniG then
print(`rMeniG(n,q,N): the probability of succesful attack  by a dishonet person with FINITE credit (it can be at most N behind)`):
print(`if the honest party mines n and the probabity of its sucess  is q (and the prob. of failure (i.e. success for the honest person) is 1-q.`):
print(`Using a generalization of Eq (1) in Merni Rosenfeld's paper "Analysis of hashrate-based double-spending". `):
print(` arXiv:1402.2009v1 [cs.CR] `):
print(` Try: `):
print(` rMeniG(5,1/3,10); `):
print(`rMeniG(5,q,10);`):

elif nops([args])=1 and op(1,[args])=rMeniGF then
print(`rMeniGF(n,q,z): the prob. generation function for duration (of second phase) time for the attacker, in case he succeeds`):
print(`It uses an adaptation of the Meni Rosenfeld formula. Try:`):
print(` rMeniGF(5,q,z); `):
print(` rMeniGF(5,q,z); `):

elif nops([args])=1 and op(1,[args])=rMeniGF1 then
print(`rMeniGF1(n,q,Z): the prob. generation function, in terms of z, and Z=1-sqrt(1-4*q*(1-q)*z^2)`):
print(` for duration (of second phase) time for the attacker, in case he succeeds`):
print(`It uses an adaptation of the Meni Rosenfeld formula. Try : `):
print(` rMeniGF1(5,q,Z); `):

elif nops([args])=1 and op(1,[args])=rMeniMOM then
print(`rMeniMOM(n,q,k): The k-th moment of "duration" of a bitcoin attack. Try:`):
print(`rMeniMOM(5,q,2);`):

elif nops([args])=1 and op(1,[args])=rMeniMOMmone then
print(`rMeniMOMmone(n,q,k): the numerator of the k-th moment of "duration" of a bitcoin attack. Try:`):
print(`rMeniMOMmone(5,q,2);`):

elif nops([args])=1 and op(1,[args])=rMeniSeq then
print(`rMeniSeq(N,q): the first N terms of probability of succesful attack  by a dishonet person`):
print(`if the honest party mines n and the probabity of its sucess  is q (and the prob. of failure (i.e. success for the honest person) is 1-q.`):
print(`Using Eq (1) in Merni Rosenfeld's paper "Analysis of hashrate-based double-spending". Try: `):
print(` rMeniSeq(5,q); `):


elif nops([args])=1 and op(1,[args])=rGP then
print(`rGP(n,q): the probability of succesful attack  by a dishonet person`):
print(` if the honest party mines n and the probabity of its sucess  is q (and the prob. of failure (i.e. success for the honest person) is 1-q. `):
print(` Using Theorem 1 in   Cyril Grunspan and Richard Perez-Marco's  paper "Double Spend Racing" `):
print(` arXiv:1702.02867v2 [cs.CR]. It should give the same answer as rMeni(n,q); `):
print(` Try: `):
print(`rGP(5,1/3);`):

elif nops([args])=1 and op(1,[args])=rSBE then
print(`rSBE(n,q): the probability of succesful attack  by a dishonet person`):
print(` if the honest party mines n and the probabity of its sucess  is q (and the prob. of failure (i.e. success for the honest person) is 1-q. `):
print(` Using Shalosh B. Ekhad's recurrence.  It should give the same answer as rMeni(n,q);`):
print(` Try: `):
print(`rSBE(5,q);`):

elif nops([args])=1 and op(1,[args])=rSBEE then
print(`rSBEE(n,q): the expected time it takes for the attack by a dishonet person provided he is succesful`):
print(`if the honest party mines n and the probabity of its sucess  is q (and the prob. of failure (i.e. success for the honest person) is 1-q.`):
print(`Using Shalosh B. Ekhad's recurrence. Try:`):
print(`rSBEE(5,1/3); `):
print(` rSBEE(5,q); `):

elif nops([args])=1 and op(1,[args])=rSBEEmone then
print(`rSBEEmone(n,q): rMeniEmone(n,q)  via the recurrence. Try`):
print(` rSBEEmone(5,q); `):

elif nops([args])=1 and op(1,[args])=rSBEEseq then
print(`rSBEEseq(N,q): Like rMeniEseq(N,q) (q.v.) , but much faster, using rSBEE(n,q) that uses the recurrence. Try`):
print(` rSBEEseq(10,q); `):

elif nops([args])=1 and op(1,[args])=rSBEMOM2 then
print(`rSBEMOM2(n1,q): The rMeniMOM(n,q,2) via a recurrence. Try:`):
print(`rSBEMOM2(10,q);`):

elif nops([args])=1 and op(1,[args])=rSBEMOM2mone then
print(`rSBEMOM2mone(n1,q): The rMeniMOMmone(n,q,2) via a recurrence. Try:`):
print(`rSBEMOM2mone(10,q);`):

elif nops([args])=1 and op(1,[args])=rSBEseq then
print(`rSBEseq(N,q): the first N terms of the probability of succesful attack  by a dishonet person`):
print(` if the honest party mines n and the probabity of its sucess  is q (and the prob. of failure (i.e. success for the honest person) is 1-q. `):
print(` Using Shalosh B. Ekhad's recurrence. q can be left symbolic (or numeric <1/2 and >0)`):
print(` Try: `):
print(`rSBEseq(10,q);`):


elif nops([args])=1 and op(1,[args])=ZXq then
print(`ZXq(X,q): both versions of 1\pm sqrt(4*q*(1-q)*X^2). Try`):
print(`ZXq(X,q);`):

else
print(`There is no ezra for`,args):
fi:
 
end:


###From Histabrut
#AveAndMoms(f,x,N): Given a probability generating function
#f(x) (a polynomial, rational function and possibly beyond)
#returns a list whose first entry is the average 
#(alias expectation, alias mean)
#followed by the variance, the third moment (about the mean) and
#so on, until the N-th moment (about the mean).
#If f(1) is not 1, than it first normalizes it by dividing
#by f(1) (if it is not 0) .
#For example, try:
#AveAndMoms(((1+x)/2)^100,x,4);

AveAndMoms:=proc(f,x,N) local mu,F,memu1,gu,i:
mu:=simplify(subs(x=1,f)):

if mu=0 then
print(f, `is neither a prob. generating function nor can it be made so`):
RETURN(FAIL):
fi:

F:=f/mu:


memu1:=simplify(subs(x=1,diff(F,x))):

gu:=[memu1]:

F:=F/x^memu1:

F:=x*diff(F,x):
for i from 2 to N do
 F:=x*diff(F,x):
 gu:=[op(gu),simplify(subs(x=1,F))]:
od:

gu:

end:


#Alpha(f,x,N): Given a probability generating function
#f(x) (a polynomial, rational function and possibly beyond)
#returns a list, of length N, whose 
#(i) First entry is the average
#(ii): Second entry is the variance
#for i=3...N, the i-th entry is the so-called alpha-coefficients
#that is the i-th moment about the mean divided by the
#variance to the power i/2 (For example, i=3 is the skewness
#and i=4 is the Kurtosis)
#If f(1) is not 1, than it first normalizes it by dividing
#by f(1) (if it is not 0) .
#For example, try:
#Alpha(((1+x)/2)^100,x,4);
Alpha:=proc(f,x,N) local gu,i:
 gu:=AveAndMoms(f,x,N):
 if gu=FAIL then
  RETURN(gu):
 fi:

if gu[2]=0 then
print(`The variance is 0`):
  RETURN(FAIL):
fi:

[gu[1],gu[2],seq(gu[i]/gu[2]^(i/2),i=3..N)]:

end:

###End From Histabrut

#rMeni(n,q): the probability of succesful attack  by a dishonet person
#if the honest party mines n and the probabity of its sucess  is q (and the prob. of failure (i.e. success for the honest person) is 1-q.
#Using Eq (1) in Merni Rosenfeld's paper "Analysis of hashrate-based double-spending". 
#arXiv:1402.2009v1 [cs.CR]
#Try
#rMeni(5,1/3);
#rMeni(5,q);
rMeni:=proc(n,q) local m:

if not (type(n, integer) and n>=0 and (type(q,numeric) and q>=0 and q<=1) or type(q,symbol)) then
 print(`Bad input`):
RETURN(FAIL):
fi:

if type(q,numeric) and q>=1/2 then
 RETURN(1):
fi:

expand(1- (1-q)^n*add(binomial(n+m-1,m)*q^m,m=0..n) +q^n*add(binomial(n+m-1,m)*(1-q)^m,m=0..n)):
end:



#rGP(n,q): the probability of succesful attack  by a dishonet person
#if the honest party mines n and the probabity of its sucess  is q (and the prob. of failure (i.e. success for the honest person) is 1-q.
#Using Theorem 1 in   Cyril Grunspan and Richard Perez-Marco's  paper "Double Spend Racing"
#arXiv:1702.02867v2 [cs.CR]
#Try
#rGP(5,1/3);

rGP:=proc(n,q) local t:
if not (type(n, integer) and n>=0 and type(q,numeric) and q>=0 and q<=1) then
 print(`Bad input`):
RETURN(FAIL):
fi:

evalf(GAMMA(n+1/2)/GAMMA(n)/GAMMA(1/2))*evalf(Int(t^(n-1)*(1-t)^(-1/2),t=0..4*q*(1-q))):

end:


#2*q*(q-1)*(-3+2*n)/(-1+n)/N^2-(-6*q^2+6*q+1+4*n*q^2-4*n*q-n)/(-1+n)/N
#rSBE(n,q): the probability of succesful attack  by a dishonet person
#if the honest party mines n and the probabity of its sucess  is q (and the prob. of failure (i.e. success for the honest person) is 1-q.
#Using S.B. Ekhad's recurrence
#rSBE(5,q);
rSBE:=proc(n,q)
option remember:
if n=0 then
1:
elif n=1 then
2*q:
else
expand(-(-6*q^2+6*q+1+4*n*q^2-4*n*q-n)/(n-1)*rSBE(n-1,q)+2*q*(q-1)*(-3+2*n)/(-1+n)*rSBE(n-2,q)):
fi:
end:




#rMeniSeq(N,q): the first N terms of the Meni polynomials
rMeniSeq:=proc(N,q) local n:
[seq(rMeni(n,q),n=1..N)]:
end:

rSBEseq:=proc(N,q) local n:
 [seq(rSBE(n,q),n=1..N)]:
end:

rSBEEseq:=proc(N,q) local n:
 normal([seq(rSBEE(n,q),n=1..N)]):
end:




#Ope(q,n,N): the second order recurrence operator annihilating the sequence of probability of success sequence 
#Try: Ope(q,n,N);
Ope:=proc(q,n,N):
-2*q*(q-1)*(1+2*n)/(1+n)+(2*q^2-2*q-1+4*n*q^2-4*n*q-n)/(1+n)*N+N^2:
end:

#OpeE(q,n,N): the second order recurrence operator annihilating the sequence of Numerators of the expected duration
OpeE:=proc(q,n,N)
-2*q*(q-1)*(1+2*n)/n+(4*n*q^2+2*q^2-4*q*n-2*q-n-2)/(1+n)*N+N^2:
end:

#OpeE1(q,n,N): the second order recurrence operator (but in N^(-1)) annihilating the sequence of Numerators of the expected duration
OpeE1:=proc(q,n,N):1-2*(q-1)*(2*n-3)*q/(n-2)/N^2+(4*q^2*n-4*q*n-n+6*q-6*q^2)/(n-1)/N:end:

#AsyExp(q,n,m): the asymptotic expansion of rMeni(n,q) to order m when q<1/2
AsyExp:=proc(q,n,m) local ope,a,i,mu,gu,N,g,x,var,eq:
ope:=Ope(q,n,N):
g:=4*q*(1-q):
mu:=1/sqrt(n)+add(a[i]/n^(i+1/2),i=1..m):
gu:=add(coeff(ope,N,i)*subs(n=n+i,mu)*g^i,i=0..degree(ope,N)):
gu:=normal(subs(n=1/x^2,gu)):
gu:=simplify(gu,symbolic):
gu:=taylor(gu,x=0,2*m+5);
eq:={seq(expand(coeff(gu,x,i)),i=1..2*m+4)}:
var:={seq(a[i],i=1..m)}:
var:=solve(eq,var):
mu:=1/sqrt(n)+add(factor(subs(var,a[i]))/n^(i+1/2),i=1..m):
mu:=1/(1-2*q)/sqrt(Pi)*g^n*mu:
mu:

end:

#AsyExpEmone(q,n,m): the asymptotic expansion of the numerator of rMeniE(n,q) to order m when q<1/2
AsyExpEmone:=proc(q,n,m) local ope,a,i,mu,gu,N,x,var,eq,g:
ope:=OpeE(q,n,N):
g:=4*q*(1-q):
mu:=1/sqrt(n)+add(a[i]/n^(i+1/2),i=1..m):
gu:=add(coeff(ope,N,i)*subs(n=n+i,mu)*g^i,i=0..degree(ope,N)):
gu:=normal(subs(n=1/x^2,gu)):
gu:=simplify(gu,symbolic):
gu:=taylor(gu,x=0,2*m+5);
eq:={seq(expand(coeff(gu,x,i)),i=1..2*m+4)}:
var:={seq(a[i],i=1..m)}:
var:=solve(eq,var):
mu:=1/sqrt(n)+add(factor(subs(var,a[i]))/n^(i+1/2),i=1..m):
mu:=g^n*mu/sqrt(Pi):
mu:

end:


#rMeniG(n,q,N): the probability of succesful attack  by a dishonet person with FINITE credit (it can be at most N behind)
#if the honest party mines n and the probabity of its sucess  is q (and the prob. of failure (i.e. success for the honest person) is 1-q.
#Using a generalization of Eq (1) in Merni Rosenfeld's paper "Analysis of hashrate-based double-spending". 
#arXiv:1402.2009v1 [cs.CR]
#Try
#rMeniG(5,1/3,10);
#rMeniG(5,q,10);
rMeniG:=proc(n,q,N) local m:

if not (type(n, integer) and n>=0 and (type(q,numeric) and q>=0 and q<=1) or type(q,symbol)) then
 print(`Bad input`):
RETURN(FAIL):
fi:

if type(q,numeric) and q>=1/2 then
 RETURN(1):
fi:

1- (1-q)^n*add(binomial(n+m-1,m)*q^m,m=0..n-1) 
+(1-q)^n*add(binomial(m+n-1,m)*q^m*( (q/(1-q))^(n-m)-(q/(1-q))^(N+1) )/(1-(q/(1-q))^(N+1)) ,m=0..n-1):
end:


#rMeniEseq(N,q): the first N terms of the Meni Expected time polynomials
rMeniEseq:=proc(N,q) local n:
normal([seq(rMeniE(n,q),n=1..N)]):
end:

#rMeniE(n,q): the expected time for the attacker, in case he succeeds
#probability of succesful attack  by a dishonet person
#if the honest party mines n and the probabity of its sucess  is q (and the prob. of failure (i.e. success for the honest person) is 1-q.
#Using Eq (1) in Merni Rosenfeld's paper "Analysis of hashrate-based double-spending". 
#arXiv:1402.2009v1 [cs.CR]
#Try
#rMeniE(5,1/3);
#rMeniE(5,q);
rMeniEold:=proc(n,q) local m:

if not (type(n, integer) and n>=0 and (type(q,numeric) and q>=0 and q<=1) or type(q,symbol)) then
 print(`Bad input`):
RETURN(FAIL):
fi:

if type(q,numeric) and q>=1/2 then
 RETURN(FAIL):
fi:

factor(normal(((1-q)^n*add(binomial(n+m-1,m)*q^m*(n-m),m=0..n))/(1-2*q))):
end:


rSBEE:=proc(n1,q) 
option remember:

normal(rSBEEmone(n1,q)/rSBE(n1,q)):
end:


#The numerator of rSBEE(n1,q)
rSBEEmone:=proc(n1,q) local ope,N,n:
option remember:
ope:=OpeE1(q,n,N):
if n1=0 then
0:
elif n1=1 then
 q/(1-2*q):
elif n1=2 then
 2*q^2*(2-q)/(1-2*q):
else
 normal(-rSBEEmone(n1-1,q)*subs(n=n1,coeff(ope,N,-1))-rSBEEmone(n1-2,q)*subs(n=n1,coeff(ope,N,-2))):
fi:

end:







#coin(q): inputs a rational mumber q between 0 and 1, outputs 1 with prob. q and -1 with prob. 1-q. Try:
#coin(1/3);
coin:=proc(q1) local q,m,n,ra:
if not (type(q1,numeric) and q1>0 and q1<1) then
 RETURN(FAIL):
fi:

if not type(q1, fraction) then
 q:=trunc(q1*10^5)/10^5:
 else
  q:=q1:
fi:



m:=numer(q):
n:=denom(q):
ra:=rand(1..n)():

if ra<=m then
 -1:
else
 1:
fi:

end:


#Phase1(n,q): inputs a positive integer n and a prob. q ( a rational number) tosses a loaded coin with prob. q of getting a 1
#and 1-q of getting a -1
#keeps playing until there are  n 1's and outputs the the number of -1's. Try
#Phase1(10,2/5);
Phase1:=proc(n,q) local gu,m,lu:
gu:=0:
m:=0:

while gu<n do
 lu:=coin(q):
 if lu=1 then
   gu:=gu+1:
 else
 m:=m+1:
 fi:
od:

m:

end;



#GRpN(q,z,N): the probability of the dishonet attacker who currently has a deficit of z of catching up if
#his probability of winning a round is q (q<1/2) and his credit is limited to N. Try
#GRp(.4,5,100);
GRpN:=proc(q,z,N) local p:
p:=1-q:

normal(((q/p)^z-(q/p)^N)/(1-(q/p)^N)):
end:

#GRp(q,z): the probability of the dishonet attacker who currently has a deficit of z of catching up if
#his probability of winning a round is q (q<1/2) with unlimited credit. Try
#GRp(.4,5,100);
GRp:=proc(q,z,N) local p:

p:=1-q:

(q/p)^z:
end:



   
#GRsim1(q,z,N): simulates gambler's ruin with current debt of z prob. q of gaining a dollar and prob. 1-q of losing one
#quits when there is a debt of N. outputs the pair [1,Duration] or [0,Duration] in case of success and failure respectively.
#Try:
#GRsim1(2/5,5,100);
GRsim1:=proc(q,z,N) local gu,lu,c:
gu:=z:
lu:=0:

while gu>0 and gu<N do
 c:=coin(q):
 lu:=lu+1:
 if c=-1 then
   gu:=gu-1
 else
   gu:=gu+1:
 fi:
od:

if gu=0 then
 RETURN([1,lu]):
else
 RETURN([0,lu]):
fi:

end:

#GRsim(q,z,N,M,x): simulates M times GRsim1(q,z,N). q is the prob. of an attacker succedding, z is his current debt, N is his maximuam credit
#x is a variable for the p.g.f for duration in case the attacker succeeds. The output is a triple
#[p,av,pgf] where p is the ratio (out of the M tries) where the attacker succeeded, av is the average duration it took him 
#(when he did succeed) and pgf is the polynomial. Try:
#GRsim(.4,5,20,100,x);
GRsim:=proc(q,z,N,M,x) local mu1,mu2,gu1,gu2,lu,i,mu,av1,av2,av,p:
gu1:=0:
gu2:=0:
mu1:=0:
mu2:=0:

for i from 1 to M do
lu:=GRsim1(q,z,N):


 if lu[1]=1 then
  gu1:=gu1+1:
   mu1:=mu1+x^lu[2]:
 else
  gu2:=gu2+1:
  mu2:=mu2+x^lu[2]:
fi:

od:


mu:=mu1+mu2:

if gu1<>0 then
 av1:=subs(x=1,diff(mu1,x))/gu1:
else
av1:=FAIL:
fi:

if gu2<>0 then
 av2:=subs(x=1,diff(mu2,x))/gu2:
else
av2:=FAIL:
fi:


 av:=subs(x=1,diff(mu,x))/M:



p:=gu1/M:

if av1<>FAIL and av2<>FAIL then
 if p*av1+(1-p)*av2<>av then
  print(`Something is wrong `):
 fi:
fi:

if mu1<>0 then
 RETURN([[p,av1,mu1/subs(x=1,mu1)], [1-p,av2,mu2/subs(x=1,mu2)],av]):
else
 RETURN(0):
fi:

end:

#fzqOld(z,q):  In a gambler's ruin,
#the expected time it takes to 0 starting at z in [0,infinity] with probability q of going left, conditioned on it getting to 0. Try
#fzqOld(z,q)
fzqOld:=proc(z,q)
z/(1-2*q)-q*(1+2*q)/(1-2*q):
end:


#fzq(z,q):  In a gambler's ruin,
#the expected time it takes to 0 starting at z in [0,infinity] with probability q of going left, conditioned on it getting to 0. Try
#fzq(z,q)
fzq:=proc(z,q)
z/(1-2*q):
end:


#ZXq(X,q): both versions of 1\pm sqrt(4*q*(1-q)*X^2). Try
#ZXq(X,q);
ZXq:=proc(X,q) local Z1,Z2:
Z1:=(1-sqrt(1-4*q*(1-q)*X^2))/(2*q*X):
Z2:=(1+sqrt(1-4*q*(1-q)*X^2))/(2*q*X):
[Z1,Z2]:
end:

#FnXq(n,X,q): the conditional prob. gen. function, in X,  of "duration of game" if succeeded in a gambler's ruin
#game starting at n with prob. q of going down and prob. 1-q going up, conditioned on the fact that it made it to 0.
#Try 
#FnXq(n,X,q);
FnXq:=proc(n,X,q) local Z,Z1,Z2:
Z:=ZXq(X,q):
if simplify(subs(X=1,Z[1]),symbolic)=1 then
 RETURN(Z[1]^n):
elif simplify(subs(X=1,Z[2]),symbolic)=1 then
 RETURN(Z[2]^n):
else
 RETURN(FAIL):
fi:

end:

#GRmoms(q,n,K): exact expressions for the (uncentraluized) moments, up to the the K-th in terms of q and n
#for the random variable "duration"  for a gambler's ruin game starting at n, with infinite credit,
#with probability q of going left (q<1/2) conditioned on it making it to 0. Try:
#GRmoms(q,n,4);
GRmoms:=proc(q,n,K) local gu,mu,X,gu1,i,lu:

lu:=ZXq(X,q):

if simplify(subs(X=1,lu[1]),symbolic)=1 then
 mu:=lu[1]:
elif simplify(subs(X=1,lu[2]),symbolic)=1 then
 mu:=lu[2]:
 else
 RETURN(FAIL):
fi:


gu1:=simplify(subs(X=1,mu),symbolic):

if gu1<>1 then
 RETURN(FAIL):
fi:

mu:=mu^n:

mu:=X*diff(mu,X):

gu1:=simplify(subs(X=1,mu),symbolic):
gu:=[gu1]:

for i from 2 to K do
 mu:=X*diff(mu,X):
 gu1:=simplify(subs(X=1,mu),symbolic):
 gu:=[op(gu),gu1]:
od:

factor(gu):

end:



#GRmomsC(q,n,K): exact expressions for the CENTRALIZED moments, up to the the K-th in terms of q and n
#for the random variable "duration"  for a gambler's ruin game starting at n, with infinite credit,
#with probability q of going left (q<1/2) conditioned on it making it to 0. In particular the
#first term gives you the expectation, and the second term gives you the variance
#Try:
#GRmomsC(q,n,4);
GRmomsC:=proc(q,n,K) local gu,av,lu,lu1,r,i:
 gu:=GRmoms(q,n,K):
av:=gu[1]:
lu:=[av]:

for i from 2 to nops(gu) do
 lu1:=expand(add(binomial(i,r)*(-av)^r*gu[i-r],r=0..i-1)+(-av)^i):
 lu:=[op(lu),lu1]:
od:

factor(lu):

end:




#GRmomsCs(q,n,K): exact expressions for the SCALED CENTRALIZED moments, up to the the K-th in terms of q and n
#for the random variable "duration"  for a gambler's ruin game starting at n, with infinite credit,
#with probability q of going left (q<1/2) conditioned on it making it to 0. In particular the
#first term gives you the expectation, and the second term gives you the variance, and the remaining ones are
#the skenewness,
#Try:
#GRmomsCs(q,n,4);
GRmomsCs:=proc(q,n,K) local gu,lu,i:
gu:=GRmomsC(q,n,K):

if K<=2 then
 RETURN(gu):
fi:

lu:=[op(1..2,gu)]:

for i from 3 to nops(gu) do
 lu:=[op(lu),gu[i]/gu[2]^(i/2)]:
od:

lu:

end:



#rMeniGF1(n,q,Z): the prob. generation function, in terms of z, and Z=1-sqrt(1-4*q*(1-q)*z^2)
# for duration (of second phase) time for the attacker, in case he succeeds
#It uses an adaptation of the Meni Rosenfeld formula. Try
#rMeniGF1(5,q,Z);
rMeniGF1:=proc(n,q,Z) local m,gu:

if not (type(n, integer) and n>=0 and (type(q,numeric) and q>=0 and q<=1) or type(q,symbol)) then
 print(`Bad input`):
RETURN(FAIL):
fi:

if type(q,numeric) and q>=1/2 then
 RETURN(FAIL):
fi:

gu:=(q/(1-q))^n*normal((1-q)^n*add(binomial(n+m-1,m)*q^m*Z^(n-m)   , m=0..n-1))  + (1-(1-q)^n*add(binomial(m+n-1,m)*q^m,m=0..n-1)):

gu:=gu/rMeni(n,q):

end:




#rMeniGF(n,q,z): the prob. generation function, in terms of z
# for duration (of second phase) time for the attacker, in case he succeeds
#It uses an adaptation of the Meni Rosenfeld formula. Try
#rMeniGF(5,q,z);
#rMeniGF(5,1/3);
#rMeniGF(5,q,z);
rMeniGF:=proc(n,q,z) local m,gu:

if not (type(n, integer) and n>=0 and (type(q,numeric) and q>=0 and q<=1) or type(q,symbol)) then
 print(`Bad input`):
RETURN(FAIL):
fi:

if type(q,numeric) and q>=1/2 then
 RETURN(1):
fi:

gu:=(1- (1-q)^n*add(binomial(n+m-1,m)*q^m,m=0..n-1)) +q^n*add(binomial(n+m-1,m)*(1-q)^m*FnXq(n-m,z,q),m=0..n-1):

gu/rMeni(n,q):

end:


#rMeniE(n,q): the prob. generation function, in terms of z
# for duration (of second phase) time for the attacker, in case he succeeds
#It uses an adaptation of the Meni Rosenfeld formula. Try
#rMeniE(5,q,z);
#rMeniE(5,1/3);
#rMeniE(5,q,z);
rMeniE:=proc(n,q) local m,gu:

if not (type(n, integer) and n>=0 and (type(q,numeric) and q>=0 and q<=1) or type(q,symbol)) then
 print(`Bad input`):
RETURN(FAIL):
fi:

if type(q,numeric) and q>=1/2 then
 RETURN(1):
fi:

gu:=q^n*add(binomial(n+m-1,m)*(1-q)^m*(n-m)/(1-2*q),m=0..n-1):

normal(gu/rMeni(n,q)):

end:


#BAsim1(n,q,N): simulates ONE bitcoin attack with prob. q (q<1/2) of success when the attacker started out with one block
#and the honest player continues until he gathers n good blcoks N (large) is the credit limit of the dishonet
#attacker. The output is the pair
#[1,LengthOfSecondPhase,m], in case of success [0,m,LengthOfSecondPhase], in case of failure. 
#m is the number of blocks the attacker got.
#Try
#BAsim1(10,2/5,200);
BAsim1:=proc(n,q,N) local m,mu:
if not (type(n,integer) and n>0 and  N>0 and type(N,integer) and type(q,numeric) and q>0 and q<1/2 ) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

m:=Phase1(n,q):

if m>=n then
 RETURN([1,0,m]):
fi:

mu:=GRsim1(q,n-m,N):
[op(mu),m]:
end:



#BAsim(n,q,N,M,x): simulates M times BAsim1(n,q,N). q is the prob. of an attacker succedding, z is his current debt, N is his maximuam credit
#x is a variable for the p.g.f for duration of phase 2, in case the attacker succeeds. The output is a triple
#[p,av,pgf] where p is the ratio (out of the M tries) where the attacker succeeded, av is the average duration it took him 
#(when he did succeed) and pgf is the polynomial. Try:
#BAsim(5,2/5,100,100,x);
BAsim:=proc(n,q,N,M,x) local mu1,mu2,gu1,gu2,lu,i,mu,av1,av2,av,p:

if not (type(n,integer) and n>0 and type(N,integer) and N>0 and type(M,integer) and M>0 and type(x,symbol)) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

gu1:=0:
gu2:=0:
mu1:=0:
mu2:=0:

for i from 1 to M do
lu:=BAsim1(n,q,N):


 if lu[1]=1 then
  gu1:=gu1+1:
   mu1:=mu1+x^lu[2]:
 else
  gu2:=gu2+1:
  mu2:=mu2+x^lu[2]:
fi:

od:


mu:=mu1+mu2:

if gu1<>0 then
 av1:=subs(x=1,diff(mu1,x))/gu1:
else
av1:=FAIL:
fi:

if gu2<>0 then
 av2:=subs(x=1,diff(mu2,x))/gu2:
else
av2:=FAIL:
fi:


 av:=subs(x=1,diff(mu,x))/M:



p:=gu1/M:

if av1<>FAIL and av2<>FAIL then
 if p*av1+(1-p)*av2<>av then
  print(`Something is wrong `):
 fi:
fi:

if mu1<>0 then
 RETURN([[p,av1,mu1/subs(x=1,mu1)], [1-p,av2,mu2/subs(x=1,mu2)],av]):
else
 RETURN(0):
fi:

end:



rMeniEmone:=proc(n,q) local m,gu:

if not (type(n, integer) and n>=0 and (type(q,numeric) and q>=0 and q<=1) or type(q,symbol)) then
 print(`Bad input`):
RETURN(FAIL):
fi:

if type(q,numeric) and q>=1/2 then
 RETURN(1):
fi:

gu:=q^n*add(binomial(n+m-1,m)*(1-q)^m*(n-m)/(1-2*q),m=0..n-1):

normal(gu):

end:


#BAsimS(n,q,N,M):  Similiar to BAsim(n,q,N,M,x) (q.v.), but only returns
#the average in case of success, and the average time it took there
#it also returns the true prob. of success and true expectation
#The output is
#[[SampleProbabiltyOfSuccess,TrueProbability], [SampleAverage, TrueExpectation],sd].
#where sd is the theoretical standard-deviation
#Try:
#BAsimS(5,2/5,100,100);
BAsimS:=proc(n,q,N,M) local gu1,lu,i,mu1,ku,sd,z:
if not (type(n,integer) and n>0 and type(N,integer) and N>0 and type(M,integer) and M>0 ) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

sd:=evalf(sqrt(Alpha(rMeniGF(n,q,z),z,2)[2])):

gu1:=0:

mu1:=0:


for i from 1 to M do
lu:=BAsim1(n,q,N):


 if lu[1]=1 then
  gu1:=gu1+1:
   mu1:=mu1+lu[2]:
 fi:

od:

if gu1<>0 then
[[evalf(gu1/M),evalf(rSBE(n,q))], [evalf(mu1/gu1),evalf(rSBEE(n,q)) ],sd]:
else
 FAIL:
fi:


end:



#EstC1(q,M,m): estimates the leading constant in the Expected duration, if successful, using
#the asymptotic expansion to order m, and the M-th term of rSBEE(M,q). M must be at least 1000. Try:
#EstC1(1/10,1000,4);

EstC1:=proc(q,M,m) local gu,gu1,gu2,mu,lu,n:
Digits:=20:
if M<1000 then
 print(`M must be at least 1000.`):
RETURN(FAIL):
fi:
 
mu:=evalf(rSBEEseq(2*M,q), Digits):

gu1:=AsyExp(q,n,m):
gu2:=AsyExpEmone(q,n,m):
gu:=gu2/gu1:

lu:=[mu[2*M]/subs(n=2*M,gu),mu[M]/subs(n=M,gu)]:

if abs(lu[1]-lu[2])>10^(-10) then

print(lu):

RETURN(FAIL):
else
 RETURN(evalf(lu[1],10)):
fi:

end:






#AsyExpE(q,n,m): the asymptotic expansion of rMeniE(n,q) (alias rSBEE(n,q) to order m when q<1/2. Try:
#AsyExpE(q,n,4);
AsyExpE:=proc(q,n,m) local x,i,gu,gu1,gu2:

gu1:=AsyExp(q,n,m+1):
gu2:=AsyExpEmone(q,n,m+1):
gu:=normal(gu2/gu1):
gu:=normal(subs(n=1/x,gu)):

gu:=taylor(gu,x=0,m+1):

gu:=add(factor(coeff(gu,x,i))/n^i,i=0..m):

(1-q)/(1-2*q)^3*gu:

end:

#CheckAsyExpE(q,m,M): Checks the asymptotic expansion to order m of the (conditional) expected duration with q<1/2
#in other words AsyExpE(q,n,m) up to the M-th term. Try:
#CheckAsyExpE(2/5,4,1000);
CheckAsyExpE:=proc(q,m,M) local n,gu,mu,i:
 gu:=AsyExpE(q,n,m):
 mu:=rSBEEseq(M,q):

[seq(evalf(mu[i]/subs(n=i,gu)),i=1..M)]:

end:




#rMeniMOMmone(n,q,k): the numerator of the k-th moment of "duration" of a bitcoin attack. Try:
#rMeniMOMmone(5,q,2);
rMeniMOMmone:=proc(n,q,k) local m,gu,mu,n1:


if not (type(n, integer) and n>=0 and (type(q,numeric) and q>=0 and q<=1) or type(q,symbol)) then
 print(`Bad input`):
RETURN(FAIL):
fi:

if type(q,numeric) and q>=1/2 then
 RETURN(1):
fi:

mu:=GRmoms(q,n1,k)[k]:

gu:=q^n*add(binomial(n+m-1,m)*(1-q)^m*subs(n1=n-m,mu),m=0..n-1):

normal(gu):

end:

#OpeM21(q,n,N): the second order recurrence operator (but in N^(-1)) annihilating the sequence of Numerators of the second moment
OpeM21:=proc(q,n,N):
1-2*(q-1)*(2*n-3)*(8*n*q^3+8*q^4-8*n*q^2-16*q^3-4*n*q+3*n+11*q-3)*q/(n-2)/(8*n*q^3+8*q^4-8*n*q^2-24*q^3-4*n*q+8*q^2+3*n+15*q-6)/N^2+(32*n^2*q^5+32*n*q^6-64*n^2*q^4-176*n*q^5-48*q^6+8*n^2*q^3+216*n*q^4+192*q^5+
36*n^2*q^2+20*n*q^3-192*q^4-8*n^2*q-126*n*q^2-42*q^3-3*n^2+31*n*q+138*q^2+3*n-54*q+6)/(n-1)/(8*n*q^3+8*q^4-8*n*q^2-24*q^3-4*n*q+8*q^2+3*n+15*q-6)/N
:
end:








#rSBEMOM2mone(n1,q): The rMeniMOMmone(n,q,2) via a recurrence. Try:
#rSBEMOM2mone(10,q);
rSBEMOM2mone:=proc(n1,q) local ope,N,n,lu:
option remember:
lu:=[q*(4*q^2-2*q-1)/(-1+2*q)^3, -2*q^2*(4*q^3-10*q^2+q+3)/(-1+2*q)^3]:
ope:=OpeM21(q,n,N):
if n1=0 then
0:
elif n1=1 then
 lu[1];
elif n1=2 then
 lu[2]:
else
 normal(-rSBEMOM2mone(n1-1,q)*subs(n=n1,coeff(ope,N,-1))-rSBEMOM2mone(n1-2,q)*subs(n=n1,coeff(ope,N,-2))):
fi:

end:


#rSBEMOM2(n1,q): The rMeniMOMmone(n,q,2) via a recurrence. Try:
#rSBEMOM2(10,q);
rSBEMOM2:=proc(n1,q):
rSBEMOM2mone(n1,q)/rSBE(n1,q):
end:

#rMeniMOM(n1,q,k): The k-th moment of "duration" in a succesful bitcoin  attack with prob. 1 and n1 good blocks
#rMeniMOM(10,q,2);
rMeniMOM:=proc(n1,q,k):
rMeniMOMmone(n1,q,k)/rSBE(n1,q):
end:

#OpeM2(q,n,N): the second order recurrence operator  annihilating the sequence of Numerators of the second moment
OpeM2:=proc(q,n,N):
-2*q*(-1+q)*(2*n+1)*(8*n*q^3+8*q^4-8*n*q^2-16*q^2-4*q*n+3*q+3*n+3)/n/(8*n*q^3-8*q^3+8*q^4-8*n*q^2-8*q^2-4*q*n+7*q+3*n)+(-24*q-9*n+30*q^2+30*q^3-16*q^4-32*q^5-q*n+16
*q^6+52*n*q^3-3*n^2-8*q*n^2+32*n^2*q^5+8*n^2*q^3+36*n^2*q^2+32*n*q^6-64*n^2*q^4-40*n*q^4-48*n*q^5+18*n*q^2)/(n+1)/(8*n*q^3-8*q^3+8*q^4-8*n*q^2-8*q^2-4*q*n+7*q+3*n)*
N+N^2:
end:



#AsyExpM2mone(q,n,m): the asymptotic expansion of the numerator of rMeniMOMmone(n,q,2) to order m when q<1/2
AsyExpM2mone:=proc(q,n,m) local ope,a,i,mu,gu,N,x,var,eq,g:
ope:=OpeM2(q,n,N):
g:=4*q*(1-q):
mu:=1/sqrt(n)+add(a[i]/n^(i+1/2),i=1..m):
gu:=add(coeff(ope,N,i)*subs(n=n+i,mu)*g^i,i=0..degree(ope,N)):
gu:=normal(subs(n=1/x^2,gu)):
gu:=simplify(gu,symbolic):
gu:=taylor(gu,x=0,2*m+5);
eq:={seq(expand(coeff(gu,x,i)),i=1..2*m+4)}:
var:={seq(a[i],i=1..m)}:
var:=solve(eq,var):
mu:=1/sqrt(n)+add(factor(subs(var,a[i]))/n^(i+1/2),i=1..m):
mu:=g^n*mu/sqrt(Pi):
mu:

end:


#AsyExpM2(q,n,m): the asymptotic expansion of rMeniMOM(n,q,2) (alias rSBEMOM2(n,q) to order m when q<1/2. Try:
#AsyExpM2(q,n,4);
AsyExpM2:=proc(q,n,m) local x,i,gu,gu1,gu2:

gu1:=AsyExp(q,n,m+1):
gu2:=AsyExpM2mone(q,n,m+1):
gu:=normal(gu2/gu1):
gu:=normal(subs(n=1/x,gu)):

gu:=taylor(gu,x=0,m+1):

gu:=add(factor(coeff(gu,x,i))/n^i,i=0..m):

-(-1+q)*(-3-2*q+4*q^2)/(-1+2*q)^5*gu:

end:


EvalCF:=proc(L) local L1:
if nops(L)=1 then
 RETURN(L[1]):
else
L[1]+1/EvalCF([op(2..nops(L),L)]):
fi:
end:

Kav2:=proc(q,K) local lu,gu,n,n1,ka,mu,i:
Digits:=40:
gu:=AsyExpM2(q,n,6):
lu:=[seq(rSBEMOM2(n1,q),n1=1..2000)]:
mu:=evalf([lu[1000]/subs(n=1000,gu), lu[1500]/subs(n=1500,gu), lu[2000]/subs(n=2000,gu)]):
if abs(mu[3]-mu[2])>1/10^K then
 RETURN(FAIL, abs(mu[3]-mu[2])):
fi:

ka:=convert(mu[3],confrac):

for i from 1 to nops(ka) while  ka[i]<10000 do od:

ka:=[op(1..i-1,ka)]:

EvalCF(ka):

end:




#CheckAsyExpM2(q,m,M): Checks the asymptotic expansion to order m of the SECOND moment (conditional)  duration with q<1/2
#in other words AsyExpM2(q,n,m) up to the M-th term. Try:
#CheckAsyExpM2(2/5,4,1000);
CheckAsyExpM2:=proc(q,m,M) local n,gu,mu,i,n1:
 gu:=AsyExpM2(q,n,m):
 mu:=[seq(rSBEMOM2(n1,q),n1=1..M)]:

[seq(evalf(mu[i]/subs(n=i,gu)),i=1..M)]:

end:




#AsyExpVar(q,n,m): the asymptotic expansion of the variance of conditional duration. Try
#AsyExpVar(q,n,4);
AsyExpVar:=proc(q,n,m) local gu,i:

#asympt(AsyExpM2(q,n,m)-AsyExpE(q,n,m)^2,n,m):

gu:=expand(AsyExpM2(q,n,m)-AsyExpE(q,n,m)^2):

add(factor(coeff(gu,n,-i))/n^i,i=0..m):

end:



#CheckAsyExpVar(q,m,M): Checks the asymptotic expansion to order m of the  (conditional)  Variance duration with q<1/2
#in other words AsyExpM2(q,n,m) up to the M-th term. Try:
#CheckAsyExpVar(2/5,4,1000);
CheckAsyExpVar:=proc(q,m,M) local n,gu,mu,i,n1:
 gu:=AsyExpVar(q,n,m):
 mu:=[seq(rSBEMOM2(n1,q)-rSBEE(n1,q)^2,n1=1..M)]:

[seq(evalf(mu[i]/subs(n=i,gu)),i=1..M)]:

end:

#Essay(q,n,M): inputs symbols q (for probabality of success of an attacker), n (for the number of good blocks generated by the honest party),
#and a positive integer M (for the order of the asymptotic expansions) and outputs an article about bitcoin atacks. Try
#Essay(q,n,4):
Essay:=proc(q,n,M) local ope,gu,m,P,E,N,L,lu,n1,mu,a,Var,lu1,lu2,i,av,sd:
Digits:=30:
print(`An Essay on Bitcoin attacks (Phrased in terms of a Two-Phase Soccer Match) `):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Assume that  two Soccer teams, the Good Guys Team and the Bad Guys Team play each other. We assume that `):
print(`the probability of the Bad Guys Team scoring the next goal is q, `):
print(`and hence that the probability of the Good Guys Team  scoring the next goal is 1-q. `):
print(`We also assume that the outcome of each goal is INDEPENDENT of the outcomes  of the other goals. `):
print(``):
print(`We assume throughout this essay that q is strictly positive but less than 1/2, so the Good Guys have an advantage. `):
print(``):
print(`There are two phases.`):
print(``):
print(`Phase I: The teams play until the Good Guys Team scores n goals. By this time, the Bad Guys Team scored a certain number of goals`):
print(`let's call it m. In the unlikely event that m is larger or equal to n, the game ends immediately, and the Bad Guys Team is declared the winner. `):
print(`Otherwise, the game proceeds to Phase II.`):
print(``):
print(`Phase II: Keep playing until the Bad Guys Team forces a tie, i.e. they both have the same number of goals, and then`):
print(`the Bad Guys Team is declared the winner, or the Good Guys Team is so much ahead (say by a Zillion goals), to make`):
print(`it hopeless for the Bad Guys to catch up, in which case the Good Guys Team won.`):
print(``):
print(`Definition: Let`, P[n](q), ` that we will call the Rosenfeld polynomials (in honor of Meni Rosenfeld, who wrote the great article`):
print(`Analysis of hashrate-based double spending, ArXiv 1402.2009v1 [cs. CR], 9 Feb. 2014. )`):
print(`be the probability of the Bad Guys Team winning. `):
print(``):
print(`Theorem 0 (Meni Rosenfeld, ArXiv 1402.2009v1, Eq. (1), p. 7)`):
print(``):
print(P[n](q)=1- (1-q)^n*Sum(binomial(n+m-1,m)*q^m,m=0..n) +q^n*Sum(binomial(n+m-1,m)*(1-q)^m,m=0..n)):
print(``):
print(`Using the so-called Zeilberger algorithm, we find a second-order linear recurrence that enables the fast`):
print(`computation of the first 10000, say, Rosenfeld polynomials, much faster than using the definition.`):
print(``):
print(``):
print(`Theorem 1: The Rosenfeld polynomials`, P[n](q), `satisfy the following linear recurrence with polynomial coefficients. `):
print( P[n](q)= -(-6*q^2+6*q+1+4*n*q^2-4*n*q-n)/(n-1)*P[n-1](q)+2*q*(q-1)*(-3+2*n)/(-1+n)*P[n-2](q)):
print(``):
print(`with initial conditions: `):
print(``):
print(P[0](q)=1 ,  P[1](q)=2*q ):
print(``):
print(`and in Maple format`):
print(``):
lprint( P[n](q)= -(-6*q^2+6*q+1+4*n*q^2-4*n*q-n)/(n-1)*P[n-1](q)+2*q*(q-1)*(-3+2*n)/(-1+n)*P[n-2](q)):
print(``):
lprint(P[0](q)=1 ,  P[1](q)=2*q ):
print(``):
print(`In another interesting paper,  by Cyril Grunspan and Ricardo Perez-Marco, entitled `):
print(``):
print(`Double Speed Races `):
print(`arXiv: 1702.02867v2 [cs.CR] 17 Feb 2017`):
print(``):
print(`the authors stated and proved the following asymptotic formula for P[n](q) (valid, of course, for 0<q<1/2)  `):
print(``):
print((4*q*(1-q))^n/sqrt(Pi)/(1-2*q)/sqrt(n) ) :
print(``):
print(`They used their representation formula in terms of the incomplete Beta function.`):
print(``):
print(`However, using the Birkhoff-Trijinski method, beautifully exposited and applied in the article `):
print(`"Resurrecting the Asymptotics of linear recurrences", by J.Wimp and D. Zeilberger, J. Math. Anal. Appl. 111 (1985), 162-177`):
print(`available from http://sites.math.rutgers.edu/~zeilberg/mamarimY/Zeilberger_y1985_p162.pdf , `):
print(`and implemented in the Maple package`):
print(` http://sites.math.rutgers.edu/~zeilberg/tokhniot/AsyRec.txt`):
print(`(See http://sites.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/asy.html for instructions)`):
print(``):
print(` we can get a much more precise asymptotics, to order`, M, `in the following theorem. `):
print(``):

gu:=AsyExp(q,n,M):
print(``):
print(`Theorem 2: The Rosenfeld polynomials, P[n](q), i.e. the probability of the Bad Guys Team winning`):
print(`have the following asymptotic expansion to order`, M):
print(``):
print(P[n](q)= gu + O(1/n^(M+1) ) ):
print(``):
print(`and in Maple format ,without the`, O(1/n^(M+1)), `we have`):
print(``):
lprint(P[n](q)= gu ):
print(``):

print(`Suppose that the Bad Guys Team is successful, how long should Phase II last until it finishes?`):
print(``):
print(`Meni Rosenfeld, in deriving Theorem 0, used the fact`):
print(`(e.g. Feller's classic) that if a random walker starts at L, and at each atomic step moves one unit left with probability q`):
print(`and  one unit to the right with probability 1-q, assuming q<1/2, the probability of making it to 0 is`):
print(``):
print((q/(1-q))^L ):
print(``):
print(`Also in Feller (but see the Georgiadis-Zeilberger paper for another proof) `):
print(`is  the famous fact that the (conditional) expected duration, in case of success is`):
print(L/(1-2*q)):
print(``):
print(`By an analogous argument to the one used  by Meni Rosenfeld we have `):
print(``):
print(`Theorem 3: Let  `):
print(``):
print(A[n](q)=q^n*Sum(binomial(n+m-1,m)*(1-q)^m*(n-m)/(1-2*q),m=0..n-1)):
print(``):
print(`then the (conditional) expectation of the duration of Phase II (in case of success to the Bad Guys Team), let's call it`,  E[n](q), `is given by `):
print(``):
print(E[n](q) = A[n](q)/P[n](q) ):
print(``):
ope:=1-OpeE1(q,n,N):
print(`Like with`, P[n](q), ` we can compute many terms of the sequence`, A[n](q) , `and hence `, E[n](q), `using the second-order recurrence given`):
print(`by Theorem 4`):
print(``):
print(`Theorem 4:`, A[n](q) , `satisfies the following second-order linear recurrence equation with polynomial coefficients: `):
print(``):
print(A[n](q)=coeff(ope,N,-1)*A[n-1](q)+coeff(ope,N,-2)*A[n-2](q) ):
print(``):
print(`subject to the initial conditions `):
print(``):
print(A[1](q)= q/(1-2*q), A[2](q)= 2*q^2*(2-q)/(1-2*q) ):
print(``):

print(`and in Maple format `):
print(``):
lprint(A[n](q)=coeff(ope,N,-1)*A[n-1](q)+coeff(ope,N,-2)*A[n-2](q) ):
print(``):
print(`subject to the initial conditions `):
print(``):
lprint(A[1](q)= q/(1-2*q), A[2](q)= 2*q^2*(2-q)/(1-2*q) ):
print(``):
print(`Using AsyRec.txt again, we get an asymptotic formula for `. A[n](q), `that we will not state, that implies an asymptotic formula`):
print(` for the (conditional) expected duration of Phase II (if the Bad Guys Team won) `, E[n](q)):
print(``):
gu:=AsyExpE(q,n,M):
print(``):
print(`Theorem 5: The (conditional) expected duration of Phase II (in case the Bad Guys Team made it)`, E[n](q)):
print(`has the following asymptotic expansion to order`, M):
print(``):
print(E[n](q)= gu + O(1/n^(M+1) ) ):
print(``):
print(`and in Maple format ,without the`, O(1/n^(M+1)), `we have`):
print(``):
lprint(E[n](q)= gu ):
print(``):

lu1:=(1-q)/(1-2*q)^2:
print(`Note that the leading term is`):
print(``):
print(lu1):
print(``):
print(`and in Maple format `):
print(``):
lprint(lu1):
print(``):

print(`and that it is INDEPENDENT of`, n, `the number of goals that the Good Guys Team scored in Phase II`):
print(`it follows that while the probability of success of the Bad Guys Team strongly depends on`, n, `and decays exponentially `):
print(`in the unlikely event that they do succeed, for large n, it takes roughly the same time in Phase II, namely`, lu1):
print(``):
print(`Here is a plot of`, lu1, `for q from 0.1 to 0.45 `):
print(``):
print(plot(lu1,q=0.1..0.45)):
print(``):
print(`For example, if q=0.45, the probability of the Bad Guys Team succeeding, followed by the expected duration of Phase II (if they do succeed) is`):
print(`for various n are `):
 for n1 from 10 by 20 to 400 do
 print(`n1=`,n1, `Prob=`, rSBE(n1,0.45), `Expected Duration=`, rSBEE(n1,0.45)):
 od:

print(`and for q=0.49, the probability of the Bad Guys Team succeeding, followed by the expected duration of Phase II (if they do succeed) is`):
print(`various n are `):
 for n1 from 1000 by 1000 to 20000 do
 print(`n1=`,n1, `Prob=`, rSBE(n1,0.49), `Expected Duration=`, rSBEE(n1,0.49)):
 od:
print(``):
print(`But what about the variance of the random variable (conditional) "duration of Phase II". We have the following theorem`):
print(``):
mu:=GRmoms(q,n,2)[2]:
print(``):
print(`Theorem 6:,Let`, B[n](q), `be defined as follows. `):
print(``):
print(B[n](q)=q^n*Sum(binomial(n+a-1,a)*(1-q)^a*subs(n1=n-a,mu),a=0..n-1)):
print(``):
print(`Then the second moment of the (conditional) random variable "duration of Phase II (in case of the Bad Guys Team winning)", let,s call it`, M2[n](q), `is `):
print(``):
print(M2[n](q)=B[n](q)/P[n](q)):
print(``):
ope:=1-OpeM21(q,n,N):
lu:=[q*(4*q^2-2*q-1)/(-1+2*q)^3, -2*q^2*(4*q^3-10*q^2+q+3)/(-1+2*q)^3]:
print(`Like with`, P[n](q),` and `, A[n](q),  `we can compute many terms of the sequence `, B[n](q)):
print(`and hence `, M2[n](q), `using the second-order recurrence given by the next theorem.`):
print(``):
print(`Theorem 7:`, B[n](q) , `satisfies the following second-order linear recurrence equation with polynomial coefficients: `):
print(``):
print(B[n](q)=coeff(ope,N,-1)*B[n-1](q)+coeff(ope,N,-2)*B[n-2](q) ):
print(``):
print(`subject to the initial conditions `):
print(``):
print(B[1](q)= lu[1], B[2](q)= lu[2] ):
print(``):
print(`and in Maple format`):
print(``):
lprint(B[n](q)=coeff(ope,N,-1)*B[n-1](q)+coeff(ope,N,-2)*B[n-2](q) ):
print(``):
print(`subject to the initial conditions `):
print(``):
lprint(B[1](q)= lu[1], B[2](q)= lu[2] ):
print(``):
gu:=AsyExpVar(q,n,M):

print(`Using AsyRec.txt yet one more time, we get an asymptotic formula for `. B[n](q), `that we will not state, and hence for M2[n](q)`):
print(`and hence for the (conditional) variance, let us call it `, Var[n](q) ):
print(``):
print(`Theorem 8: The (conditional) variance of the duration of Phase II (in case the Bad Guys Team made it)`, Var[n](q) ):
print(`has the following asymptotic expansion to order`, M):
print(``):
print(Var[n](q)= gu + O(1/n^(M+1) ) ):
print(``):
print(`and in Maple format ,without the`, O(1/n^(M+1)), `we have`):
print(``):
lprint(Var[n](q)= gu ):
print(``):
print(`Once again, like the expectation the leading term of the variance does not depend on n`):
lu2:= (q-1)*(4*q^2-3*q-2)/(2*q-1)^4;
print(`and is`, lu2):
print(`and hence the standard deviation is, for large n, tends to a constant (that only depends on q)`):
lu2:=sqrt(lu2):
print(``):
print(lu2):
print(``):
print(`and in Maple format `):
print(``):
lprint(lu2):
print(``):

print(``):
print(`note that there is no concentration about the mean. Not only that, the standard deviation is larger than the expectation`):
print(``):
print(`the duration of Phase II, in case of success, is very unpredictable. `):
print(``):
print(`Let's conclude with a table of the limiting expected duration and standard deviations for q from 0.1 to 0.49`):
print(``):
for i  from 0 to 39 do
print(q=0.1+i*0.01, av=subs(q=0.1+i*0.01,lu1), sd=subs(q=0.1+i*0.01,lu2)):
od:
print(``):
print(`----------------------------------------------------`):
print(``):
print(`This concludes this essay. Thanks for reading it. `):

end:




#EssayLaTex(q,n,M):  Like Essay(q,n,N) (q.v.) but instead of lprint it uses latex . Try:
#EssayLaTex(q,n,4):
EssayLaTex:=proc(q,n,M) local ope,gu,m,P,E,N,L,lu,n1,mu,a,Var,lu1,lu2,i,av,sd:
Digits:=30:
print(`An Essay on Bitcoin attacks (Phrased in terms of a Two-Phase Soccer Match) `):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Assume that  two Soccer teams, the Good Guys Team and the Bad Guys Team play each other. We assume that `):
print(`the probability of the Bad Guys Team scoring the next goal is q, `):
print(`and hence that the probability of the Good Guys Team  scoring the next goal is 1-q. `):
print(`We also assume that the outcome of each goal is INDEPENDENT of the outcomes  of the other goals. `):
print(``):
print(`We assume throughout this essay that q is strictly positive but less than 1/2, so the Good Guys have an advantage. `):
print(``):
print(`There are two phases.`):
print(``):
print(`Phase I: The teams play until the Good Guys Team scores n goals. By this time, the Bad Guys Team scored a certain number of goals`):
print(`let's call it m. In the unlikely event that m is larger or equal to n, the game ends immediately, and the Bad Guys Team is declared the winner. `):
print(`Otherwise, the game proceeds to Phase II.`):
print(``):
print(`Phase II: Keep playing until the Bad Guys Team forces a tie, i.e. they both have the same number of goals, and then`):
print(`the Bad Guys Team is declared the winner, or the Good Guys Team is so much ahead (say by a Zillion goals), to make`):
print(`it hopeless for the Bad Guys to catch up, in which case the Good Guys Team won.`):
print(``):
print(`Definition: Let`, P[n](q), ` that we will call the Rosenfeld polynomials (in honor of Meni Rosenfeld, who wrote the great article`):
print(`Analysis of hashrate-based double spending, ArXiv 1402.2009v1 [cs. CR], 9 Feb. 2014. )`):
print(`be the probability of the Bad Guys Team winning. `):
print(``):
print(`Theorem 0 (Meni Rosenfeld, ArXiv 1402.2009v1, Eq. (1), p. 7)`):
print(``):
print(P[n](q)=1- (1-q)^n*Sum(binomial(n+m-1,m)*q^m,m=0..n) +q^n*Sum(binomial(n+m-1,m)*(1-q)^m,m=0..n)):
print(``):
print(`Using the so-called Zeilberger algorithm, we find a second-order linear recurrence that enables the fast`):
print(`computation of the first 10000, say, Rosenfeld polynomials, much faster than using the definition.`):
print(``):
print(``):
print(`Theorem 1: The Rosenfeld polynomials`, P[n](q), `satisfy the following linear recurrence with polynomial coefficients. `):
print( P[n](q)= -(-6*q^2+6*q+1+4*n*q^2-4*n*q-n)/(n-1)*P[n-1](q)+2*q*(q-1)*(-3+2*n)/(-1+n)*P[n-2](q)):
print(``):
print(`with initial conditions: `):
print(``):
print(P[0](q)=1 ,  P[1](q)=2*q ):
print(``):
print(`and in latex`):
print(``):
lprint( latex(P[n](q)= -(-6*q^2+6*q+1+4*n*q^2-4*n*q-n)/(n-1)*P[n-1](q)+2*q*(q-1)*(-3+2*n)/(-1+n)*P[n-2](q))):
print(``):
lprint(latex(P[0](q)=1) ,  latex(P[1](q)=2*q) ):
print(``):
print(`In another interesting paper,  by Cyril Grunspan and Ricardo Perez-Marco, entitled `):
print(``):
print(`Double Speed Races `):
print(`arXiv: 1702.02867v2 [cs.CR] 17 Feb 2017`):
print(``):
print(`the authors stated and proved the following asymptotic formula for P[n](q) (valid, of course, for 0<q<1/2)  `):
print(``):
print((4*q*(1-q))^n/sqrt(Pi)/(1-2*q)/sqrt(n) ) :
print(``):
print(`They used their representation formula in terms of the incomplete Beta function.`):
print(``):
print(`However, using the Birkhoff-Trijinski method, beautifully exposited and applied in the article `):
print(`"Resurrecting the Asymptotics of linear recurrences", by J.Wimp and D. Zeilberger, J. Math. Anal. Appl. 111 (1985), 162-177`):
print(`available from http://sites.math.rutgers.edu/~zeilberg/mamarimY/Zeilberger_y1985_p162.pdf , `):
print(`and implemented in the Maple package`):
print(` http://sites.math.rutgers.edu/~zeilberg/tokhniot/AsyRec.txt`):
print(`(See http://sites.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/asy.html for instructions)`):
print(``):
print(` we can get a much more precise asymptotics, to order`, M, `in the following theorem. `):
print(``):

gu:=AsyExp(q,n,M):
print(``):
print(`Theorem 2: The Rosenfeld polynomials, P[n](q), i.e. the probability of the Bad Guys Team winning`):
print(`have the following asymptotic expansion to order`, M):
print(``):
print(P[n](q)= gu + O(1/n^(M+1) ) ):
print(``):
print(`and in Maple format ,without the`, O(1/n^(M+1)), `we have`):
print(``):
lprint(latex(P[n](q)= gu) ):
print(``):

print(`Suppose that the Bad Guys Team is successful, how long should Phase II last until it finishes?`):
print(``):
print(`Meni Rosenfeld, in deriving Theorem 0, used the fact`):
print(`(e.g. Feller's classic) that if a random walker starts at L, and at each atomic step moves one unit left with probability q`):
print(`and  one unit to the right with probability 1-q, assuming q<1/2, the probability of making it to 0 is`):
print(``):
print((q/(1-q))^L ):
print(``):
print(`Also in Feller (but see the Georgiadis-Zeilberger paper for another proof) `):
print(`is  the famous fact that the (conditional) expected duration, in case of success is`):
print(L/(1-2*q)):
print(``):
print(`By an analogous argument to the one used  by Meni Rosenfeld we have `):
print(``):
print(`Theorem 3: Let  `):
print(``):
print(A[n](q)=q^n*Sum(binomial(n+m-1,m)*(1-q)^m*(n-m)/(1-2*q),m=0..n-1)):
print(``):
print(`then the (conditional) expectation of the duration of Phase II (in case of success to the Bad Guys Team), let's call it`,  E[n](q), `is given by `):
print(``):
print(E[n](q) = A[n](q)/P[n](q) ):
print(``):
ope:=1-OpeE1(q,n,N):
print(`Like with`, P[n](q), ` we can compute many terms of the sequence`, A[n](q) , `and hence `, E[n](q), `using the second-order recurrence given`):
print(`by Theorem 4`):
print(``):
print(`Theorem 4:`, A[n](q) , `satisfies the following second-order linear recurrence equation with polynomial coefficients: `):
print(``):
print(A[n](q)=coeff(ope,N,-1)*A[n-1](q)+coeff(ope,N,-2)*A[n-2](q) ):
print(``):
print(`subject to the initial conditions `):
print(``):
print(A[1](q)= q/(1-2*q), A[2](q)= 2*q^2*(2-q)/(1-2*q) ):
print(``):

print(`and in Maple format `):
print(``):
lprint(latex(A[n](q)=coeff(ope,N,-1)*A[n-1](q)+coeff(ope,N,-2)*A[n-2](q) )):
print(``):
print(`subject to the initial conditions `):
print(``):
lprint(latex(A[1](q)= q/(1-2*q)), latex(A[2](q)= 2*q^2*(2-q)/(1-2*q))  ):
print(``):
print(`Using AsyRec.txt again, we get an asymptotic formula for `. A[n](q), `that we will not state, that implies an asymptotic formula`):
print(` for the (conditional) expected duration of Phase II (if the Bad Guys Team won) `, E[n](q)):
print(``):
gu:=AsyExpE(q,n,M):
print(``):
print(`Theorem 5: The (conditional) expected duration of Phase II (in case the Bad Guys Team made it)`, E[n](q)):
print(`has the following asymptotic expansion to order`, M):
print(``):
print(E[n](q)= gu + O(1/n^(M+1) ) ):
print(``):
print(`and in Maple format ,without the`, O(1/n^(M+1)), `we have`):
print(``):
lprint(latex(E[n](q)= gu )):
print(``):

lu1:=(1-q)/(1-2*q)^2:
print(`Note that the leading term is`):
print(``):
print(lu1):
print(``):
print(`and in Maple format `):
print(``):
lprint(latex(lu1)):
print(``):

print(`and that it is INDEPENDENT of`, n, `the number of goals that the Good Guys Team scored in Phase II`):
print(`it follows that while the probability of success of the Bad Guys Team strongly depends on`, n, `and decays exponentially `):
print(`in the unlikely event that they do succeed, for large n, it takes roughly the same time in Phase II, namely`, lu1):
print(``):
print(`Here is a plot of`, lu1, `for q from 0.1 to 0.45 `):
print(``):
print(plot(lu1,q=0.1..0.45)):
print(``):
print(`For example, if q=0.45, the probability of the Bad Guys Team succeeding, followed by the expected duration of Phase II (if they do succeed) is`):
print(`for various n are `):
 for n1 from 10 by 20 to 400 do
 print(`n1=`,n1, `Prob=`, rSBE(n1,0.45), `Expected Duration=`, rSBEE(n1,0.45)):
 od:

print(`and for q=0.49, the probability of the Bad Guys Team succeeding, followed by the expected duration of Phase II (if they do succeed) is`):
print(`various n are `):
 for n1 from 1000 by 1000 to 20000 do
 print(`n1=`,n1, `Prob=`, rSBE(n1,0.49), `Expected Duration=`, rSBEE(n1,0.49)):
 od:
print(``):
print(`But what about the variance of the random variable (conditional) "duration of Phase II". We have the following theorem`):
print(``):
mu:=GRmoms(q,n,2)[2]:
print(``):
print(`Theorem 6:,Let`, B[n](q), `be defined as follows. `):
print(``):
print(B[n](q)=q^n*Sum(binomial(n+a-1,a)*(1-q)^a*subs(n1=n-a,mu),a=0..n-1)):
print(``):
print(`Then the second moment of the (conditional) random variable "duration of Phase II (in case of the Bad Guys Team winning)", let,s call it`, M2[n](q), `is `):
print(``):
print(M2[n](q)=B[n](q)/P[n](q)):
print(``):
ope:=1-OpeM21(q,n,N):
lu:=[q*(4*q^2-2*q-1)/(-1+2*q)^3, -2*q^2*(4*q^3-10*q^2+q+3)/(-1+2*q)^3]:
print(`Like with`, P[n](q),` and `, A[n](q),  `we can compute many terms of the sequence `, B[n](q)):
print(`and hence `, M2[n](q), `using the second-order recurrence given by the next theorem.`):
print(``):
print(`Theorem 7:`, B[n](q) , `satisfies the following second-order linear recurrence equation with polynomial coefficients: `):
print(``):
print(B[n](q)=coeff(ope,N,-1)*B[n-1](q)+coeff(ope,N,-2)*B[n-2](q) ):
print(``):
print(`subject to the initial conditions `):
print(``):
print(B[1](q)= lu[1], B[2](q)= lu[2] ):
print(``):
print(`and in Maple format`):
print(``):
lprint(latex(B[n](q)=coeff(ope,N,-1)*B[n-1](q)+coeff(ope,N,-2)*B[n-2](q) )):
print(``):
print(`subject to the initial conditions `):
print(``):
lprint(latex(B[1](q)= lu[1]), latex(B[2](q)= lu[2] )) :
print(``):
gu:=AsyExpVar(q,n,M):

print(`Using AsyRec.txt yet one more time, we get an asymptotic formula for `. B[n](q), `that we will not state, and hence for M2[n](q)`):
print(`and hence for the (conditional) variance, let us call it `, Var[n](q) ):
print(``):
print(`Theorem 8: The (conditional) variance of the duration of Phase II (in case the Bad Guys Team made it)`, Var[n](q) ):
print(`has the following asymptotic expansion to order`, M):
print(``):
print(Var[n](q)= gu + O(1/n^(M+1) ) ):
print(``):
print(`and in Maple format ,without the`, O(1/n^(M+1)), `we have`):
print(``):
lprint(latex(Var[n](q)= gu )):
print(``):
print(`Once again, like the expectation the leading term of the variance does not depend on n`):
lu2:= (q-1)*(4*q^2-3*q-2)/(2*q-1)^4;
print(`and is`, lu2):
print(`and hence the standard deviation is, for large n, tends to a constant (that only depends on q)`):
lu2:=sqrt(lu2):
print(``):
print(lu2):
print(``):
print(`and in Maple format `):
print(``):
lprint(latex(lu2)):
print(``):

print(``):
print(`note that there is no concentration about the mean. Not only that, the standard deviation is larger than the expectation`):
print(``):
print(`the duration of Phase II, in case of success, is very unpredictable. `):
print(``):
print(`Let's conclude with a table of the limiting expected duration and standard deviations for q from 0.1 to 0.49`):
print(``):
for i  from 0 to 39 do
print(q=0.1+i*0.01, av=subs(q=0.1+i*0.01,lu1), sd=subs(q=0.1+i*0.01,lu2)):
od:
print(``):
print(`----------------------------------------------------`):
print(``):
print(`This concludes this essay. Thanks for reading it. `):

end:



#CutOff(q,T): the smallest n such that rMeni(n,q)<T. It returns the actual probability, and the expected duration. Try:
#CutOff(0.45, 10^(-3));
CutOff:=proc(q,T) local n:

for n from 1 while rSBE(n,q)>T do od:

[n,rSBE(n,q),rMeniE(n,q)]  :

end:

#CutOffStory(qs,qe,N,K): the cut-offs staring with qs and ending with qe, with N intermediate stops with
#confidence level 10^(-k), for k from 1 to K. Try:
#CutOffStory(0.3,0.4,2,3);
CutOffStory:=proc(qs,qe,N,K) local del,q,gu,i,k:
del:=(qe-qs)/N:

for i from 0 to N do
q:=qs+i*del:

 print(`If the probability of a bad guy getting the next block is`, q):
print(``):

for k from 1 to K do
 gu:=CutOff(q,10^(-k)):
print(`in order to guarantee prob. of success less than`, 10^(-k), `you have to mint`, gu[1], `good coins `):
print(`and if the dishonest attacker succeeds, he should expect to take`, gu[3], `steps in phase II`):
print(``):

od:
od:
end: