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

#VERSION OF July 2023

with(Statistics):


read `StPeteData.txt`:

print(`First Written: July 2023: tested for Maple 20  `):
print():
print(`This is StPete.txt, a Maple package `):
print(`accompanying Lucy Martinez and Doron Zeilberger's  article: `):
print(`A Guide to the Risk-Averse Gambler and Resolving the St. Petersburg Paradox Once and For All`):

print(`-------------------------------`):

print(`-------------------------------`):
print():
print(`The most current version is available at the url:`):
print(` https://sites.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/stpete.html`):
print(`Please report all bugs to: DoronZeil@gmail.com  .`):
print():
print(`-------------------------------`):
print(`For general help, and a list of the MAIN functions,`):
print(` type "Help();". For specific help type "Help(procedure_name);" `):

print(`-------------------------------`):
print(`For a list of the supporting functions type:  Help1();`):
print(`For help with a specific procedure, type "Help(procedure_name);"`):
print():
 print(`-------------------------------`):




print(`For a list of the  general functions (not described in the paper) where cutoff of loss is general, type:  HelpG();`):
print(`For help with a specific procedure, type "Help(procedure_name);"`):
print():



print(`-------------------------------`):
print(`For a list of the STORY functions type:  HelpS();`):
print(`For help with a specific procedure, type "Help(procedure_name);"`):
print():
 print(`-------------------------------`):




HelpS:=proc()
if args=NULL then

print(` The STORY procedures are: Paper1, Paper1Ab, Paper1a, Paper1aAb, Paper2, Paper2a, Paper3, Paper3A, Paper4  `):
print(` `):

else
 Help(args):
fi:
end:

HelpG:=proc()
if args=NULL then

print(` The Generalized procedures are:  OpeProbMore, ProbMore, ProbMoreSeq, ProbMoreSeqF`):
print(` `):

else
 Help(args):
fi:
end:

Help1:=proc()
if args=NULL then

print(` The supporting procedures are:  CollectOpe, ExpVal, MyPT, OpeK, OpeKpc, OpeProbPos, PGF, plotDist, ProbPos, ProbPosA,`):
print(` `):

else
 Help(args):
fi:
end:

Help:=proc()

if args=NULL then

print(` StPete.txt: A Maple package for  helping risk averse gamblers and experimenting with the St. Petesburg paradox `):
print(`The main procedures are:  CutOff, CutOffF, CutOffA,  CutOffF, CutOffK, ProbPosSeq, ProbPosSeqF, Simu, Simu1, StPetePT, `):




elif nargs=1 and args[1]=CollectOpe then 
print(`CollectOpe(K,n,Sn): inputs a positive integer K>=2 and outputs a list of operators in (n,Sn) annihilating ProbPos([[-1, (i-1)/i],[i,1/i]]   ],    n) for i from N. Try:`):
print(`CollectOpe(3,n,Sn):`):

elif nargs=1 and args[1]=CutOff then 
print(`CutOff(P,RA) inputs prob. table, and a pos. number RA between 0 and 1`):
print(`and outputs the smallest integer n such that if you play P n times, `):
print(`your probability of being positive is larger than 1-RA. Try:`):
print(`CutOff([[-1,2/3],[3,1/3]],0.01);`):


elif nargs=1 and args[1]=CutOffA then 
print(`CutOffA(P,RA): an approximation of CutOff(P,RA) using the central limit theorem approximation. Only valid for high rish aversion, i.e. small RA. Try:`):
print(`CutOffA([[-1,2/3],[3,1/3]],0.01);`):

elif nargs=1 and args[1]=CutOffF then 
print(` CutOffF(P,RA) same as CutOff(P,RA), but (hopefully) faster, using the operator. Try:`):
print(`CutOff([[-1,4/5],[5,1/5]],0.01);`):

elif nargs=1 and args[1]=CutOffK then 
print(` CutOffK(K,RA) same as CutOfF([[-1,(K-1)/K],[K,1/K]],RA).  It also returns the approximate value using the Central Limit Theorem. Try:`):
print(`CutOffK(5,0.01);`):

elif nargs=1 and args[1]=ExpVal then 
print(`ExpVal(T) computes the expectation value for the probability table T, followed by the standard-deviation`):


elif nargs=1 and args[1]=HowManyRounds then
print(`HowManyRounds(P,k) inputs a probability table P and a positive integer k,`):
print(`and outputs a list of integers of length k that would instruct a risk-averse gambler`):
print(`how many rounds to play with risk-averseness 1/10^i for i=1,...,k and plots the data. Try:`):
print(`HowManyRounds([[1,3/5],[-1,2/5]],5)`):


#elif nargs=1 and args[1]=MetaSimu then
#print(`MetaSimu(T,N) first finds the expected gain of T (if you don't have to pay anything), called K, and runs Simu(P,i,N) for 1<=i<=K. Try`):
#print(`MetaSimu(StPetePT(7,3),1000);`):


elif nargs=1 and args[1]=MyPT then
print(`MyPT(P) inputs P=[[a1,p1],[a2,p2],...,[ar,pr]] where p1+p2+...+pr=1, and uses the Maple Sample function to output ai with probability pi. Try`):
print(`MyPT([[-1,1/3],[1,1/6],[5,1/2]]);`):


elif nargs=1 and args[1]=OpeK then
print(`OpeK(K,n,Sn): The operator annihilating ProbPos([[-1,(K-1)/K],[K,1/K]]) together with the initial condtions. In other words OpeProbPos([[-1,(K-1)/K],[K,1/K]]). Try:`):
print(`OpeK(5,n,Sn);`):


elif nargs=1 and args[1]=OpeKpc then
print(`OpeKpc(K,n,Sn): The PRE-COMPUTED OpeK(K,n,Sn) for K from 2 to 10. Try:`):
print(`OpeKpc(10,n,Sn);`):

elif nargs=1 and args[1]=OpeProbMore then
print(`OpeProbMore(P,n,Sn,L): inputs a probabiity table P and symbols n and Sn (the shift in n), as well as an integer (positive or negative) L, and outputs the linear recurrence operator in n, Sn, annihilating`):
print(`ProbMore(P,n,L), along with the initial conditions. Try:`):
print(`OpeProbMore([[-1,4/5],[5,1/5]],n,Sn,1);`):



elif nargs=1 and args[1]=OpeProbPos then
print(`OpeProbPos(P,n,Sn): inputs a probabiity table P and symbols n and Sn (the shift in n) and outputs the linear recurrence operator in n, Sn, annihilating`):
print(`ProbPos(P,n), along with the initial conditions. Try:`):
print(`OpeProbPos([[-1,4/5],[5,1/5]],n,Sn);`):

elif nargs=1 and args[1]=Paper1 then
print(`Paper1(k,K1,K2): Advice to the risk-averse gambler if the probability table is [[-1,(1-1)/i],[i,1/i]] for i from 2 to k. Do:`):
print(`Paper1(5,100,2000):`):


elif nargs=1 and args[1]=Paper1Ab then
print(`Paper1Ab(k,K1,K2): abbreviated version of Paper1(k,K1,K2) without displaying the recurrences. Try: `):
print(`Paper1Ab(5,100,2000):`):


elif nargs=1 and args[1]=Paper1a then
print(`Paper1a(P,K1,K2): inputs a probability table P with positive Expectation output a paper with advise to the risk-averse gambler. K2 is the maximum allowed number of rounds and K1 is how often you are told.Try:`):
print(`Paper1a([[-1,4/5],[5,1/5]],100,1000):`):


elif nargs=1 and args[1]=Paper1aAb then
print(`Paper1aAb(P,K1,K2): Abbreviated version of Paper1a(P,K1,K2) not displaying the operators. Try`):
print(`Paper1aAb([[-1,4/5],[5,1/5]],100,1000):`):


elif nargs=1 and args[1]=Paper2 then
print(`Paper2(k,K1,K2): Advice to the risk-averse gambler if you are playing the St. Petersburg casino with i rounds, entrance fee i-1, from  i to k.It starts at i=4 and goes to i=k. It tells youthe prob. of making money from K1 times to K2 times`):
print(`with intervals of K1. Try:`):
print(`Paper2(6,10,50):`):


elif nargs=1 and args[1]=Paper2a then
print(`Paper2a(k,K1,K2): inputs a positive integer k and outputs advice to the risk-averse gambler entering the St. Petersburg casino with k rounds (whose expected gain is k),`):
print(`paying k-1 dollars to enter the casino, and hence the expectation is 1 dollar K1 and K2 are such that you are given`):
print(`the prob. of exiting with positive money from K1 to K2 with increments of K1. Try:`):
print(`Paper2a(5,100,1000):`):


elif nargs=1 and args[1]=Paper3 then
print(`Paper3(K,N): A paper about the cutoffs for risk-aversness for playing gambling games [[-1,(i-1)/i],[i,1/i]] for i from 2 to 10, and risk-averseness 10^(-j) for j from 1 to  N. `):
print(`It also gives the Central Limit approximations`):
print(`Paper3(6,3);`):


elif nargs=1 and args[1]=Paper3A then
print(`Paper3A(K,N): A paper about the cutoffs for risk-aversness for playing gambling games [[-1,(i-1)/i],[i,1/i]] for i from 2 to 10, and risk-averseness 10^(-j) for j from 1 to  N. `):
print(` using the Central Limit approximations. Try: `):
print(`Paper3A(10,4);`):

elif nargs=1 and args[1]=Paper4 then
print(`Paper4(K,N): A paper about the cutoffs, using the CLT approximations, for risk-aversness for playing The Finite St. Petersburg game with i rounds from i=3 to i=K with`):
print(`risk-averseness 1/10^j, for j from 1 to N. Try:`):
print(`Paper4(7,4);`):


elif nargs=1 and args[1]=PGF then
print(`PGF(P,t) inputs a probability table P and a variable t and outputs the probability generating function Sum(P[i][2]*t^P[i][1],i=1..nops(P))`):



elif nargs=1 and args[1]=ProbMore then
print(` ProbMore(P,N,L) outputs the exact prob. of exiting with >=L dollars. ProbMore(P,N,0) is the same as ProbPos(P,N). Try:`):
print(`ProbMore([[-1,2/3],[3,1/3]],10,-2);`):


elif nargs=1 and args[1]=ProbMoreSeq then
print(`ProbMoreSeq(P,N,L): The first N terms of ProbMoreSeq(P,n,L), done directly. Try:`):
print(`ProbMoreSeq([[-1,4/5],[5,1/5]],10,0);`):

elif nargs=1 and args[1]=ProbMoreSeqF then
print(`ProbMoreSeqF(P,N,L): A fast version of ProbMoreSeq(P,N,L) using the recurrence gotten by the Almkvist-Zeilberger algorithm. Try:`):
print(`ProbMoreSeqF([[-1,4/5],[5,1/5]],100,0);`):

elif nargs=1 and args[1]=ProbPos then
print(` ProbPos(P,N) outputs the EXACT (as a rational number) prob. of winning a positive amount of money if you play the "die" P, N times. Try`):
print(`ProbPos([[-1,9/10],[10,1/10]],10);`):



elif nargs=1 and args[1]=ProbPosA then
print(`ProbPosA(P,N) outputs the approximates prob. of not losing money if you play the "die" P, N times`):
print(`using the approximation by the Central Limit Theorem, pretending that N is large. Try:`):
print(`ProbPosA([[-1,4/5],[5,1/5]],10);`):

elif nargs=1 and args[1]=ProbPosSeq then
print(`ProbPosSeq(P,N): The first N terms of ProbPos(P,n), done directly. Try:`):
print(`ProbPosSeq([[-1,4/5],[5,1/5]],10);`):


elif nargs=1 and args[1]=ProbPosSeqF then
print(`ProbPosSeqF(P,N): A fast vesion of ProbPosSeq(P,N), using the recurrence gotten from the Almkvist-Zeilberger algorithm. Try:`):
print(`ProbPosSeqF([[-1,4/5],[5,1/5]],40);`):



elif nargs=1 and args[1]=plotDist then
print(`plotDist(f,x): Given a prob. gen. polynnmial f(x), plots it` ):

elif nargs=1 and args[1]=Simu then
print(`Simu(P,n,N) runs Simu1(P,n) N times, and outputs the average of your gains and the fraction of times that it was positive. Try `):
print(`Simu([[-1,1/3],[1,1/6],[5,1/2]],100,1000);`):

elif nargs=1 and args[1]=Simu1 then
print(`Simu1(P,n) adds MyPT(P) n times. Try `):
print(`Simu1([[-1,1/3],[1,1/6],[5,1/2]],100);`):



elif nargs=1 and args[1]=StPetePT then
print(`StPetePT(k,M) describes the outcome of the (finite version of) the Saint Petersburg Paradox where you pay M dollars to enter and the game lasts up to k rounds. Hence, your initial capital is -M. Try`):
print(`StPetePT(3,2);`):



else
 print(`There is no such thing as`, args):

fi:

end:




###############################Start from EKHAD.txt
 

pashetd:=proc(p,N)
local gu,p1,mekh,i:
p1:=normal(p):
mekh:=denom(p1):
p1:=numer(p1):
p1:=expand(p1):
gu:=0:
for i from 0 to degree(p1,N) do
gu:=gu+factor(coeff(p1,N,i))*N^i:
od:
RETURN(gu,mekh):
end:
 
goremd:=proc(N,ORDER,bpc)
local i, gu:
gu:=0:
for i from 0 to ORDER do
gu:=gu+(bpc[i])*N^i:
od:
RETURN(gu):
end:

duisd:= proc(A,ORDER,y,n,N)
local gorem, conj, yakhas,lu1a,LU1,P,Q,R,S,j1,g,Q1,Q2,l1,l2,mumu,
mekb1,fu,meka1,k1,gugu,eq,ia1,va1,va2,degg,shad,KAMA,i1,bpc,apc:
KAMA:=20:
gorem:=goremd(N,ORDER,bpc):
 
 
conj:=gorem:
yakhas:=0:
 
for i1 from 0 to degree(conj,N) do
 
yakhas:=yakhas+coeff(conj,N,i1)*simplify(subs(n=n+i1,A)/A):
od:
 
lu1a:=yakhas:
LU1:=numer(yakhas):
yakhas:=1/denom(yakhas):
yakhas:=normal(diff(yakhas,y)/yakhas+diff(A,y)/A):
P:=LU1:
Q:=numer(yakhas):
R:=denom(yakhas):
j1:=0:
while j1 <= KAMA do
g:=gcd(R,Q-j1*diff(R,y)):
if g <> 1 then
Q2:=(Q-j1*diff(R,y))/g:
R:=normal(R/g):
Q:=normal(Q2+j1*diff(R,y)):
P:=P*g^j1:
j1:=-1:
fi:
j1:=j1+1:
od:
P:=expand(P):
R:=expand(R):
Q:=expand(Q):
Q1:=Q+diff(R,y):
Q1:=expand(Q1):
l1:=degree(R,y):
l2:=degree(Q1,y):
meka1:=coeff(R,y,l1):
mekb1:=coeff(Q1,y,l2):
if l1-1 <>l2 
then
k1:=degree(P,y)-max(l1-1,l2):
else
 mumu:= -mekb1/meka1:
  if type(mumu,integer) and mumu > 0
  then
   k1:=max(mumu, degree(P,y)-l2):
  else
   k1:=degree(P,y)-l2:
  fi:
fi:
fu:=0:
if k1 < 0 then
RETURN(0):
fi:
for ia1 from 0 to k1 do
fu:=fu+apc[ia1]*y^ia1:
od:
gugu:=Q1*fu+R*diff(fu,y)-P:
gugu:=expand(gugu):
degg:=degree(gugu,y):
for ia1 from 0 to degg do
eq[ia1+1]:=coeff(gugu,y,ia1)=0:
od:
for ia1 from 0 to k1 do
va1[ia1+1]:=apc[ia1]:
od:
for ia1 from 0 to ORDER do
va2[ia1+1]:=bpc[ia1]:
od:
eq:=convert(eq,set):
va1:=convert(va1,set):
va2:=convert(va2,set):
va1:=va1 union va2 :
va1:=solve(eq,va1):
fu:=subs(va1,fu):
gorem:=subs(va1,gorem):
if fu=0 and gorem=0 then
 RETURN(0):
fi:
for ia1 from 0 to k1 do
gorem:=subs(apc[ia1]=1,gorem):
od:
for ia1 from 0 to ORDER do
gorem:=subs(bpc[ia1]=1,gorem):
od:
fu:=normal(fu):


shad:=pashetd(gorem,N):
S:=lu1a*R*fu*shad[2]/P:
S:=subs(va1,S):
S:=normal(S):
S:=factor(S):
for ia1 from 0 to k1 do
S:=subs(apc[ia1]=1,S):
od:
for ia1 from 0 to ORDER do
S:=subs(bpc[ia1]=1,S):
od:

RETURN(shad[1],S):
end:

AZd:= proc(A,y,n,N)
local ORDER,gu,KAMA:
KAMA:=40:

for ORDER from 0 to KAMA do
gu:=duisd(A,ORDER,y,n,N):
if gu<>0 then
 RETURN(gu):
fi:
od:
0:
end:

Yafe:=proc(ope,N) local i,ope1,coe1,L: 
if ope=0 then
 RETURN(1,0):
fi:
ope1:=expand(ope):
L:=degree(ope1,N):
coe1:=coeff(ope1,N,L):
ope1:=normal(ope1/coe1):
ope1:=normal(ope1):
ope1:=
convert(
[seq(factor(coeff(ope1,N,i))*N^i,i=ldegree(ope1,N)..degree(ope1,N))],`+`):
factor(coe1),ope1:
end:

#SeqFromRec(ope,n,N,Ini,K): Given the first L-1
#terms of the sequence Ini=[f(1), ..., f(L-1)]
#satisfied by the recurrence ope(n,N)f(n)=0
#extends it to the first K values
SeqFromRec:=proc(ope,n,N,Ini,K)
local ope1,gu,L,n1,j1:
ope1:=Yafe(ope,N)[2]:
L:=degree(ope1,N):
if nops(Ini)<>L then
 ERROR(`Ini should be of length`, L):
fi:

ope1:=expand(subs(n=n-L,ope1)/N^L):

gu:=Ini:

for n1 from nops(Ini)+1 to K do
gu:=[op(gu), -add(gu[nops(gu)+1-j1]*subs(n=n1,coeff(ope1,N,-j1)),
j1=1..L)]:
od:

gu:

end:

########################33

with(Statistics):

# MyPT(P) inputs P=[[a1,p1],[a2,p2],...,[ar,pr]] where p1+p2+...+pr=1, and uses the Maple Sample function to
#  output ai with probability pi
MyPT:=proc(P) local n,newP,i,X,j,k:
n:=nops(P):
newP:=[]:
for i from 1 to n do
 newP:=[op(newP),P[i,2]]:
od:
if add(j, j in newP)=1 then
 X:=RandomVariable(ProbabilityTable(newP)):
 k:=convert(Sample(X,1)[1],integer):
else
 print(`Sum of all the probabilities (elements in the second entry of the given list) should add up to 1`):
fi:
P[k,1]:
end:


# Simu1(P,n) adds MyPT(P) n times
Simu1:=proc(P,n) local i:
add(MyPT(P),i=1..n):
end:


# Simu(P,n,N) runs Simu1(P,n) N times, and outputs the average of your gains and the fraction of times
#  that it was positive
Simu:=proc(P,n,N) local i,s,t,co:
s:=0:
co:=0:
for i from 1 to N do
 t:=Simu1(P,n):
 s:=s+t:
 if t>0 then
  co:=co+1:
 fi:
od:
evalf(s/N),evalf(co/N):
end:


# StPetePT(k,M) describes the outcome of the (finite version of) the Saint Petersburg Paradox where 
# you pay M dollars to enter the Casino (your initial capital is -M) and the game lasts up to k rounds
# StPetePT(3,0) should output [[1,1/2],[2,1/4],[4,1/8],[4,1/8]]
# StPetePT(3,2) should output [[-1,1/2],[0,1/4],[2,1/8],[2,1/8]]
StPetePT:=proc(k,M) local L,i:
L:=[]:

for i from 0 to k-2 do
 L:=[op(L),[-M+2^(i+1),1/2^(i+1)]]:
od:

[op(L),[-M+2^(k-1),1/2^(k-1)]]:



end:


# MetaSimu(T,N) first finds the expected gain of T (if you don't have to pay anything), call it K,
#  and then for i=1,2,...,K runs Simu(P,i,N)
#MetaSimu:=proc(T,N) local K,i,L:
#K:=trunc(add(T[i][1]*T[i][2],i=1..nops(T))):
#if K>=1 then
# L:=[]:
 # Here we might have to change since Simu outputs two numbers
# for i from 1 to K do
#  L:=[op(L),[Simu(T,i,N)]]:
# od:
#else
# return `the expected gain is less than 0`:
#fi:
#L:
#end:

# PGF(P,t) inputs a probability table P and a variable t and outputs the probability generating function
#  Sum(P[i][1]*t^P[i][2],i=1..nops(P))
PGF:=proc(P,t) local i:
sum(P[i][2]*t^P[i][1],i=1..nops(P)):
end:


# ProbPosOld(P,N) outputs the exact prob. of not losing money if you play the "die" P, N times
ProbPosOld:=proc(P,N) local f,l,m,s,i,t :
f:=expand(PGF(P,t)^N):
l:=ldegree(f):
m:=degree(f):
s:=0:
for i from l to m do
 if i>0 then
  s:=s+coeff(f,t,i):
 fi:
od:
s:
end:




# ProbPos(P,N) outputs the exact prob. of not losing money if you play the "die" P, N times
ProbPos:=proc(P,N) local t,f,f1,i:

if not (type(P,list) and  add(P[i][2],i=1..nops(P))=1 and  min(seq(P[i][2],i=1..nops(P)))>=0 and max(seq(P[i][2],i=1..nops(P)))<=1)  then
 print(P, `is not a probability table`):
 RETURN(FAIL):
fi:

for i from 1 to nops(P) do
if not type(P[i][1],integer) then
print(`All stakes must be integers`):
 RETURN(FAIL):
fi:
od:

f:=expand(PGF(P,t)^N):

if not type(f,`+`) then
RETURN(FAIL):
fi:

f1:=0:

for i from 1 to nops(f) do
  if degree(op(i,f),t)>0 then
   f1:=f1+op(i,f):
  fi:
od:
subs(t=1,f1):

end:



# ExpVal computes the expectation value for the probability table T, followed by the standard-deviation
ExpVal:=proc(T) local i,mu,sig:

if not (type(T,list)  and min(seq(T[i][2],i=1..nops(T)))>=0 and max(seq(T[i][2],i=1..nops(T)))<=1 and add(T[i][2],i=1..nops(T))=1) then
 print(T , ` is bad `):
 RETURN(FAIL):
fi:

mu:=add(T[i][1]*T[i][2],i=1..nops(T)):
sig:=sqrt(add((T[i][1]-mu)^2*T[i][2],i=1..nops(T))):
[mu,sig]:
end:


# CutOff(P,RA) inputs prob. table, and a pos. number RA between 0 and 1
#  and outputs the smallest integer n such that if you play P n times your probability
#  of being pos. is larger than 1-RA
CutOff:=proc(P,RA) local d,co,i:

if not (type(P,list) and  add(P[i][2],i=1..nops(P))=1 and  min(seq(P[i][2],i=1..nops(P)))>=0 and max(seq(P[i][2],i=1..nops(P)))<=1)  then
 print(P, `is not a probability table`):
 RETURN(FAIL):
fi:

d:=1-RA:
co:=1:
while ProbPos(P,co)<=d do
 co:=co+1:
od:
co:
end:


# CutOffF(P,RA) same as CutOff(P,RA), but (hopefully) faster, using the operator
#Try: CutOff([[-1,4/5],[5,1/5]],0.01);
CutOffF:=proc(P,RA) local N,d,co,i,L:

if not (type(P,list) and  add(P[i][2],i=1..nops(P))=1 and  min(seq(P[i][2],i=1..nops(P)))>=0 and max(seq(P[i][2],i=1..nops(P)))<=1)  then
 print(P, `is not a probability table`):
 RETURN(FAIL):
fi:

N:=CutOffA(P,RA)+100:
L:=ProbPosSeqF(P,N):

d:=1-RA:

for co from 1 to nops(L) while L[co]<=d do od:


if co>=nops(L) then
 RETURN(FAIL):
else
 RETURN(co):
fi:

end:




# CutOffK(K,RA) same as CutOfF([[-1,(K-1)/K],[K,1/K]],RA). It also returns the approximate value using the Central Limit Theorem.Try:
#CutOffK(5,0.01);
CutOffK:=proc(K,RA) local P,ope,L,N,d,co:
 P:=[[-1,(K-1)/K],[K,1/K]]:

if K<=10 then
 ope:=OpeKpc(K,n,Sn):
else
 ope:=OpeK(K,n,Sn):
fi:



N:=CutOffA(P,RA)+100:

L:=SeqFromRec(ope[1],n,Sn,ope[2],N):


d:=1-RA:

for co from 1 to nops(L) while L[co]<=d do od:

if co>=nops(L) then
 RETURN(FAIL):
else
 RETURN([co, CutOffA(P,RA)]):
fi:

end:





# HowManyRounds(P,k) inputs a probability table P and a positive integer k,
#  and outputs a list of integers of length k that would instruct a risk-averse gambler
#  how many rounds to play with risk-averseness 1/10^i for i=1,...,k and plots the data
HowManyRounds:=proc(P,k) local L,i:
L:=[seq(CutOff(P,1/10^i),i=1..k)]:
plot([seq([i,L[i]],i=1..k)]):
end:



#ProbPosSeq(P,N): The first N terms of ProbPos(P,n), done directly. Try:
#ProbPosSeq([[-1,4/5],[5,1/5]],10);
ProbPosSeq:=proc(P,N) local n,i:

for i from 1 to nops(P) do
if not type(P[i][1],integer) then
print(`All stakes must be integers`):
 RETURN(FAIL):
fi:
od:

[seq(ProbPos(P,n),n=1..N)]:
end:


#OpeProbPos(P,n,Sn): inputs a probabiity table P and symbols n and Sn (the shift in n) and outputs the linear recurrence operator in n, Sn, annihilating
#ProbPos(P,n), along with the initial conditions. Try:
#OpeProbPos([[-1,4/5],[5,1/5]],n,Sn);
OpeProbPos:=proc(P,n,Sn) local f,t,ope:
option remember:
f:=PGF(P,t):
ope:=AZd(f^n/(1-t)/t,t,n,Sn):
if ope=FAIL then
 RETURN(FAIL):
fi:
ope:=ope[1]:
[ope,ProbPosSeq(P,degree(ope,Sn))]:
end:

#OpeK(K,n,Sn): The operator annihilating ProbPos([[-1,(K-1)/K],[K,1/K]]) together with the initial condtions. In other words OpeProbPos([[-1,(K-1)/K],[K,1/K]]). Try:
#OpeK(5,n,Sn);
OpeK:=proc(K,n,Sn): 
if not (type(K,integer) and K>=2) then
 RETURN(FAIL):
fi:
OpeProbPos([[-1,(K-1)/K],[K,1/K]],n,Sn):
end:

#ProbPosSeqF(P,N): A fast version of ProbPosSeq(P,N) using the recurrence gotten by the Almkvist-Zeilberger algorithm. Try:
#ProbPosSeqF([[-1,4/5],[5,1/5]],100);
ProbPosSeqF:=proc(P,N) local gu,n,Sn,i:

for i from 1 to nops(P) do
if not type(P[i][1],integer) then
print(`All stakes must be integers`):
 RETURN(FAIL):
fi:
od:

gu:=OpeProbPos(P,n,Sn):
if gu=FAIL then
 RETURN(FAIL):
fi:
SeqFromRec(gu[1],n,Sn,gu[2],N):
end:

##start generalized procedures



# ProbMore(P,N,L) outputs the exact prob. of exiting with >=L dollars. ProbMore(P,N,0) is the same as ProbPos(P,N). Try:
#ProbMore([[-1,2/3],[3,1/3]],10,-2);
ProbMore:=proc(P,N,L) local t,f,f1,i:


if not (type(P,list) and  add(P[i][2],i=1..nops(P))=1 and  min(seq(P[i][2],i=1..nops(P)))>=0 and max(seq(P[i][2],i=1..nops(P)))<=1)  then
 print(P, `is not a probability table`):
 RETURN(FAIL):
fi:

for i from 1 to nops(P) do
if not type(P[i][1],integer) then
print(`All stakes must be integers`):
 RETURN(FAIL):
fi:
od:


f:=expand(PGF(P,t)^N):

if not type(f,`+`) then
RETURN(FAIL):
fi:

f1:=0:

for i from 1 to nops(f) do
  if degree(op(i,f),t)>=L then
   f1:=f1+op(i,f):
  fi:
od:
subs(t=1,f1):

end:




#ProbMoresSeq(P,N,L): The first N terms of ProbMore(P,n,L), done directly. Try:
#ProbMoreSeq([[-1,4/5],[5,1/5]],10,1);
ProbMoreSeq:=proc(P,N,L) local n,i:

for i from 1 to nops(P) do
if not type(P[i][1],integer) then
print(`All stakes must be integers`):
 RETURN(FAIL):
fi:
od:

[seq(ProbMore(P,n,L),n=1..N)]:
end:

#OpeProbMore(P,n,Sn,L): inputs a probabiity table P and symbols n and Sn (the shift in n), as well as an integer (positive or negative) L, and outputs the linear recurrence operator in n, Sn, annihilating
#ProbMore(P,n,L), along with the initial conditions. Try:
#OpeProbMore([[-1,4/5],[5,1/5]],n,Sn,1);
OpeProbMore:=proc(P,n,Sn,L) local f,t,ope:
option remember:
f:=PGF(P,t):
ope:=AZd(f^n/(1-t)/t^L,t,n,Sn):
if ope=FAIL then
 RETURN(FAIL):
fi:
ope:=ope[1]:
[ope,ProbMoreSeq(P,degree(ope,Sn),L)]:
end:




#ProbMoreSeqF(P,N,L): A fast version of ProbMoreSeq(P,N,L) using the recurrence gotten by the Almkvist-Zeilberger algorithm. Try:
#ProbMoreSeqF([[-1,4/5],[5,1/5]],100,0);
ProbMoreSeqF:=proc(P,N,L) local gu,n,Sn,i:

for i from 1 to nops(P) do
if not type(P[i][1],integer) then
print(`All stakes must be integers`):
 RETURN(FAIL):
fi:
od:

gu:=OpeProbMore(P,n,Sn,L):
if gu=FAIL then
 RETURN(FAIL):
fi:
SeqFromRec(gu[1],n,Sn,gu[2],N):
end:




#plotDist(f,x): Given a prob. gen. polynmial f(x) that has a 
plotDist:=proc(f,x) local i,f1,L:
f1:=expand(f):
L:=[]:

for i from ldegree(f1,x) to degree(f1,x) do
if coeff(f1,x,i)<>0 then
 L:=[op(L),[i,coeff(f1,x,i)]]:
fi:
od:

plot(L)


end:




#ProbPosA(P,N) outputs the approximates prob. of not losing money if you play the "die" P, N times
#using the approximation by the Central Limit Theorem, pretending that N is large. Try
#ProbPosA([[-1,4/5],[5,1/5]],10);
ProbPosA:=proc(P,N) local gu,i,a:

if not (type(P,list) and  add(P[i][2],i=1..nops(P))=1 and  min(seq(P[i][2],i=1..nops(P)))>=0 and max(seq(P[i][2],i=1..nops(P)))<=1)  then
 print(P, `is not a probability table`):
 RETURN(FAIL):
fi:

for i from 1 to nops(P) do
if not type(P[i][1],integer) then
print(`All stakes must be integers`):
 RETURN(FAIL):
fi:
od:


gu:=ExpVal(P):

a:=-gu[1]*sqrt(N)/gu[2]:

(1-erf(a/sqrt(2)))/2:
end:




# CutOffA(P,RA) inputs prob. table, and a pos. number RA between 0 and 1
#  and outputs the approximation, only valid for high risk-aversion of CutOff(P,RA).
CutOffA:=proc(P,RA) local d,co,i:

if not (type(P,list) and  add(P[i][2],i=1..nops(P))=1 and  min(seq(P[i][2],i=1..nops(P)))>=0 and max(seq(P[i][2],i=1..nops(P)))<=1)  then
 print(P, `is not a probability table`):
 RETURN(FAIL):
fi:

d:=1-RA:

for co from 1 by 1000 while evalf(ProbPosA(P,co))<=d do od:

for i from co-1000+1 to co while evalf(ProbPosA(P,i))<=d do od:
i:

end:

# CutOffAoriginal(P,RA) inputs prob. table, and a pos. number RA between 0 and 1
#  and outputs the approximation, only valid for high risk-aversion of CutOff(P,RA).
CutOffAoriginal:=proc(P,RA) local d,co,i:

if not (type(P,list) and  add(P[i][2],i=1..nops(P))=1 and  min(seq(P[i][2],i=1..nops(P)))>=0 and max(seq(P[i][2],i=1..nops(P)))<=1)  then
 print(P, `is not a probability table`):
 RETURN(FAIL):
fi:

d:=1-RA:
co:=1:
while evalf(ProbPosA(P,co))<=d do
 co:=co+1:
od:
end:

#CollectOpe(K,n,Sn): inputs a positive integer K>=2 and outputs a list of operators in (n,Sn) annihilating ProbPos([[-1, (i-1)/i],[i,1/i]]   ],    n) for i from N. Try:
#CollectOpe(3,n,Sn):
CollectOpe:=proc(N,n,Sn) local gu,P,i,ope:
gu:=[]:
for i from 1 to N do
P:=[[-1,(i-1)/i],[i,1/i]]:
ope:=OpeProbPos(P,n,Sn):
gu:=[op(gu),ope]:
od:
gu:
end:



#Paper1a(P,K1,K2): inputs a probability table P with positive Expectation output a paper with advise to the risk-averse gambler. K2 is the maximum allowed number of rounds and it tells you every K1 rounds .Try:
#Paper1a([[-1,4/5],[5,1/5]],100,1000):
Paper1a:=proc(P,K1,K2) local t0, ope,lu,ProbL,i,INI,Sn,n,LI:

if not (type(P,list) and  add(P[i][2],i=1..nops(P))=1 and  min(seq(P[i][2],i=1..nops(P)))>=0 and max(seq(P[i][2],i=1..nops(P)))<=1)  then
 print(P, `is not a probability table`):
 RETURN(FAIL):
fi:

t0:=time():
ProbL:=0:
for i from 1 to nops(P) do
 if P[i][1]<0 then
 ProbL:=ProbL+P[i][2]:
 fi:
od:


print(`Suppose you are offered a bet where`):
ptint(``):

for i from 1 to nops(P) do
 if P[i][1]<0 then
  print(`With probability`, P[i][2], `you lose `, -P[i][1] , `dollars . `):
 else
   print(`With probability`, P[i][2], `you win `, P[i][1] , `dollars . `):
 fi:
od:

lu:=ExpVal(P):

if lu[1]<=0 then
 print(` The expectation `, lu[1], ` is not positive, so it is a bad deal, refuse!`):
 RETURN(FAIL):
fi:

print(` The expectation `, lu[1], ` is positive, so conventional wisdom tells you that it is a good deal, and you should accept it, but beware,`):
print(``):
print(`the standard-deviation is`, evalf(lu[2]), `that should raise a red flag.`):
print(``):
print(`Indeed, it you are only allowed to play once, then the probability that you would lose money is`, ProbL ):
print(``):
print(`It would be stupid to take this bet.`):
print(``):
print(`However, by the law of large numbers, and more  quantitatively by the Central Limit Theorem, if you are allowed to play many times,  you can make the probability of`):
print(`losing money as small as you wish (of course, every gamble carries some rish), the question is how many times should you insist on being able to play? `):
print(``):


ope:=OpeProbPos(P,n,Sn):
INI:=ope[2]:
ope:=ope[1]:

print(``):

print(`Let p(n) be the probability that if you play n times, you would gain (at least some) money`):
print(``):
print(`Thanks to the amazing Almkvist-Zeilberger algorithm, there is a fast way to compute this sequence, via the following recurrence `):

print(``):
print(add(coeff(ope,Sn,i)*p(n+i),i=0..degree(ope,Sn))=0):
print(``):
print(`subject to the initial conditions `):
print(``):
print(seq(p(i)=INI[i],i=1..nops(INI))):
print(``):

print(`and in Maple notation`):

print(``):
lprint(add(coeff(ope,Sn,i)*p(n+i),i=0..degree(ope,Sn))=0):
print(``):
print(`subject to the initial conditions `):
print(``):
lprint(seq(p(i)=INI[i],i=1..nops(INI))):
print(``):

LI:=SeqFromRec(ope,n,Sn,evalf(INI),K2):


print(`Using this recurrence we can compute the following probabilities`):
print(``):
for i from K1  by K1 to K2 do

print(`If you play`, i, `times then your probability of winning a positive amount of money is`, LI[i]):

od:

print(`-------------------------------`):
print(``):
print(`Here is a plot of the number of rounds vs. the probability of winding up with a positive amount of money`):
print(``):
print(plot([seq([i,LI[i]],i=1..nops(LI))]));
print(``):
print(`------------------------------`):
print(``):
print(`This ends this chapter that took`, time()-t0, `seconds to generate.`):
print(``):
end:


#Paper1aAb(P,K1,K2): Abbreviated version of Paper1(P,K1,K2). Try:
#Paper1aAb([[-1,4/5],[5,1/5]],100,1000):
Paper1aAb:=proc(P,K1,K2) local t0, ope,lu,ProbL,i,INI,Sn,n,LI:

if not (type(P,list) and  add(P[i][2],i=1..nops(P))=1 and  min(seq(P[i][2],i=1..nops(P)))>=0 and max(seq(P[i][2],i=1..nops(P)))<=1)  then
 print(P, `is not a probability table`):
 RETURN(FAIL):
fi:

t0:=time():
ProbL:=0:
for i from 1 to nops(P) do
 if P[i][1]<0 then
 ProbL:=ProbL+P[i][2]:
 fi:
od:


print(`Suppose you are offered a bet where`):
ptint(``):

for i from 1 to nops(P) do
 if P[i][1]<0 then
  print(`With probability`, P[i][2], `you lose `, -P[i][1] , `dollars . `):
 else
   print(`With probability`, P[i][2], `you win `, P[i][1] , `dollars . `):
 fi:
od:

lu:=ExpVal(P):

if lu[1]<=0 then
 print(` The expectation `, lu[1], ` is not positive, so it is a bad deal, refuse!`):
 RETURN(FAIL):
fi:

print(` The expectation `, lu[1], ` is positive, so conventional wisdom tells you that it is a good deal, and you should accept it, but beware,`):
print(``):
print(`the standard-deviation is`, evalf(lu[2]), `that should raise a red flag.`):
print(``):
print(`Indeed, it you are only allowed to play once, then the probability that you would lose money is`, ProbL ):
print(``):
print(`It would be stupid to take this bet.`):
print(``):
print(`However, by the law of large numbers, and more  quantitatively by the Central Limit Theorem, if you are allowed to play many times,  you can make the probability of`):
print(`losing money as small as you wish (of course, every gamble carries some rish), the question is how many times should you insist on being able to play? `):
print(``):


ope:=OpeProbPos(P,n,Sn):
INI:=ope[2]:
ope:=ope[1]:

print(``):

print(`Let p(n) be the probability that if you play n times, you would gain (at least some) money`):
print(``):
print(`Thanks to the amazing Almkvist-Zeilberger algorithm, there is a fast way to compute this sequence, via a certain recurrence that we do not display `):
print(`of order`, degree(ope,Sn), `and degree of coefficients `, degree(ope,n) ):

LI:=SeqFromRec(ope,n,Sn,evalf(INI),K2):


print(`Using this recurrence we can compute the following probabilities`):
print(``):
for i from K1  by K1 to K2 do

print(`If you play`, i, `times then your probability of winning a positive amount of money is`, LI[i]):

od:

print(`-------------------------------`):
print(``):
print(`Here is a plot of the number of rounds vs. the probability of winding up with a positive amount of money`):
print(``):
print(plot([seq([i,LI[i]],i=1..nops(LI))]));
print(``):
print(`------------------------------`):
print(``):
print(`This ends this chapter that took`, time()-t0, `seconds to generate.`):
print(``):
end:









#Paper2a(k,K1,K2): inputs a positive integer k and outputs advice to the risk-averse gambler entering the St. Petersburg casino with k rounds (whose expected gain is k),
#paying k-1 dollars to enter the casino, and hence the expectation is 1 dollar K1 and K2 are such that you are given
#the prob. of exiting with positive money from K1 to K2 with increments of K1. Try
#Paper2a(5,100,1000):
Paper2a:=proc(k,K1,K2) local P,LI,t0,lu,i, ProbL:

t0:=time():
print(`Suppose you are at St. Petersburg and are ready to enter its famous casino, where there are`, k , ` rounds `):
print(`With prob. 1/2 you win two ducats, and must leave`):
print(`If you still here, With prob. 1/4 you four two ducats, and must leave`):
print(`etc. If you  survived`, k, ` rounds, you get`, 2^k, `ducats `):
print(``):
print(`Note that the expected gain is`, k, `ducats, hence if you pay an entrance fee of`, k-1, `ducats you still expect to win one ducat, so if you are rational you should play,  or shold you?`):
print(``):
print(`Depending on your risk-averseness, you should insist on the permission to play it multiple times`):
print(``):

P:=StPetePT(k,k-1):


ProbL:=0:
for i from 1 to nops(P) do
 if P[i][1]<0 then
 ProbL:=ProbL+P[i][2]:
 fi:
od:


print(`So here is the probability distribution`):
ptint(``):

for i from 1 to nops(P) do
 if P[i][1]<0 then
  print(`With probability`, P[i][2], `you lose `, -P[i][1] , `dollars . `):
 else
   print(`With probability`, P[i][2], `you win `, P[i][1] , `dollars . `):
 fi:
od:


lu:=ExpVal(P):

print(` The expectation  is indeed 1,  so conventional wisdom tells you that it is a good deal, and you should accept it, but beware,`):
print(``):
print(`the standard-deviation is`, evalf(lu[2]), `that should raise a red flag.`):
print(``):
print(`Indeed, it you are only allowed to play once, then the probability that you would lose money is`, ProbL ):
print(``):
print(`It would be stupid to take this bet.`):
print(``):
print(`However, by the law of large numbers, and more  quantitatively by the Central Limit Theorem, if you are allowed to play many times,  you can make the probability of`):
print(`losing money as small as you wish (of course, every gamble carries some rish), the question is how many times should you insist on being able to play? `):
print(``):


LI:=evalf(ProbPosSeq(P,K2)):


print(``):
for i from K1  by K1 to K2 do

print(`If you play`, i, `times then your probability of winning a positive amount of money is`, LI[i]):

od:

print(`-------------------------------`):
print(``):
print(`Here is a plot of the number of rounds vs. the probability of winding up with a positive amount of money`):
print(``):
print(plot([seq([i,LI[i]],i=1..nops(LI))]));
print(``):
print(`------------------------------`):
print(``):
print(`This ends this chapter that took`, time()-t0, `seconds to generate.`):
print(``):
end:



#Paper1(k,K1,K2): Advice to the risk-averse gambler if the probability table is [[-1,(1-1)/i],[i,1/i]] for i from 2 to k. Do:
#Paper1(5,100,2000):
Paper1:=proc(k,K1,K2) local i,t0:
t0:=time():
print(``):
print(`Advice to the risk averse gambler if the Probality Table is [[-1,(i-1)/i],[i,1/i]] for i from 2 to`, k):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`----------------------------------------------------`):
for i from 2 to k do
print(``):
Paper1a([[-1,(i-1)/i],[i,1/i]],K1,K2):
print(``):
od:
print(`--------------------------------------`):
print(`This ends this paper that took`, time()-t0, `seconds to produce. `):
print(``):
end:


#Paper1Ab(k,K1,K2): Abbreviated version of Paper1(k,K1,K2), w/o printing the operators. Try:
#Paper1Ab(5,100,2000):
Paper1Ab:=proc(k,K1,K2) local i,t0:
t0:=time():
print(``):
print(`Advice to the risk averse gambler if the Probality Table is [[-1,(i-1)/i],[i,1/i]] for i from 2 to`, k):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`----------------------------------------------------`):
for i from 2 to k do
print(``):
Paper1aAb([[-1,(i-1)/i],[i,1/i]],K1,K2):
print(``):
od:
print(`--------------------------------------`):
print(`This ends this paper that took`, time()-t0, `seconds to produce. `):
print(``):
end:


#Paper2(k,K1,K2): Advice to the risk-averse gambler for playing the St. Peterburg with i rounds (and entrance fee i-1) for i from 4 to k. Do
#Paper2(5,10,50):
Paper2:=proc(k,K1,K2) local i,t0:
t0:=time():
print(``):
print(`Advice to the risk averse gambler about playing at the St. Petersburg casino with i rounds from i=4 to i=`, k):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`----------------------------------------------------`):
for i from 4 to k do
print(``):
Paper2a(i,K1,K2):
print(``):
od:
print(`--------------------------------------`):
print(`This ends this paper that took`, time()-t0, `seconds to produce. `):
print(``):
end:



#Paper3(K,N): A paper about the cutoffs for risk-aversness for playing gambling games [[-1,(i-1)/i],[i,1/i]] for i from 2 to K, and risk-averseness 10^(-j) for j from 1 to  N. Try:
#Paper3(4);
#also gives the Central Limit approximations
Paper3:=proc(K,N) local j,lu,i,t0:

t0:=time():
print(`How many times should a risk-averse gambler inisist on  playing with Probability table [[-1,(i-1)/i],[i,1/i]] for i from 2 to 10 with risk-averseness 10^(-j) for j from 1 to `, N):
print(``):
print(`By Shalosh B. Ekhad`):
print(``):

for i from  2 to K do
print(``):
print(`------------------------------------`):
print(``):
print(`Suppose you are offered a bet such that `):
print(`You lose`, i , `dollar with probability`, (i-1)/i, `and you win`, i^2, `dollars with probability`, 1/i ):
print(``):
print(`The expected gain is one dollar, so if you are rational you should accept this bet. Alas you hate to lose money.`):
print(``):

for j from 1 to N do
lu:=CutOffK(i,10^(-j)):
print(`If you  are willing to take a chance of `,1/10^j, `of losing some money then you should insist on playing this game`, lu[1], `times. BTW, the approx. using CLT gives`, lu[2]):
od:

od:

print(`--------------------------------------`):
print(`This ends this paper that took`, time()-t0, `seconds to produce. `):
print(``):

end:




#Paper3A(K,N): A paper about the cutoffs, using the CLT approximations, for risk-aversness for playing gambling games [[-1,(i-1)/i],[i,1/i]] for i from 2 to K, and risk-averseness 10^(-j) for j from 1 to  N. Try:
#Paper3A(4);
Paper3A:=proc(K,N) local j,lu,i,t0:
t0:=time():

print(`How many times should a risk-averse gambler inisist on  playing with Probability table [[-1,(i-1)/i],[i,1/i]] for i from 2 to 10 with risk-averseness 10^(-j) for j from 1 to `, N):
print(``):
print(`By Shalosh B. Ekhad`):
print(``):

print(`The above values are based on the Central Limit theorem approximations, and are slightlt higher than necessary for the stated risk-averseness`):
print(``):
for i from  2 to K do
print(``):
print(`------------------------------------`):
print(``):
print(`Suppose you are offered a bet such that `):
print(`You lose`, i , `dollar with probability`, (i-1)/i, `and you win`, i^2, `dollars with probability`, 1/i ):
print(``):
print(`The expected gain is one dollar, so if you are rational you should accept this bet. Alas you hate to lose money.`):
print(``):

for j from 1 to N do
lu:=CutOffA([[-1,(i-1)/i],[i,1/i]],10^(-j)):
print(`If you  are willing to take a chance of `,1/10^j, `of losing some money then you should insist on playing this game`, lu, `times `):
od:

od:

print(`--------------------------------------`):
print(`This ends this paper that took`, time()-t0, `seconds to produce. `):
print(``):

end:



#Paper4(K,N): A paper about the cutoffs, using the CLT approximations, for risk-aversness for playing The Finite St. Petersburg game with i rounds from i=3 to i=K with
#risk-averseness 1/10^j, for j from 1 to N. Try:
#Paper4(7,4);
Paper4:=proc(K,N) local j,lu,i,P,t0:
t0:=time():

print(`How many times should a risk-averse gambler inisist on  playing with Probability table [[-1,(i-1)/i],[i,1/i]] for i from 2 to 10 with risk-averseness 10^(-j) for j from 1 to `, N):
print(``):
print(`By Shalosh B. Ekhad`):
print(``):

print(`The above values are based on the Central Limit theorem approximations, and are slightlt higher than necessary for the stated risk-averseness`):
print(``):
for i from  3 to K do
print(``):
print(`------------------------------------`):
print(``):
P:=StPetePT(i,i-1):
print(`For the St. Petersburg finite version with`, i, `rounds and entrance fee`, i-1, `so the expected gain of one shot is 1, the probability table is`):
print(``):
lprint(P):
print(``):
for j from 1 to N do
lu:=CutOffA(P,10^(-j)):
print(`If you  are willing to take a chance of `,1/10^j, `of losing some money then you should insist on playing this game`, lu, `times `):
od:

od:

print(`--------------------------------------`):
print(`This ends this paper that took`, time()-t0, `seconds to produce. `):
print(``):

end: