# restart; read `G:/My Drive/MyNote/Talks
# /25.5 DRZ75/Primes.txt`;

with(numtheory):
with(Statistics):
with(LinearAlgebra):

# First version: November 2, 2025

################################
# Section 1: Compute expectation
# and second moment
#
# Functions:
# EPrime(s,N,BV), ProbPrime(s,N),
# EPrimeNoMem(s,N,BV),
#
################################

# Input: s=value to be start, N=cutoff value, 
# BV=Boundary value
# Output: Expected number of rolls starting at s.
# Try: 
# Digits:=1003:
# A:= EPrime(0,70010,0.);
# B := EPrime(0,70010,1000.);
# B-A;
EPrime := proc(s,N,BV) option remember;
local i;

if s > N then
	return(BV);
elif isprime(s) then
	return(0);
fi:	

1+1/6*add(EPrime(s+i,N,BV),i=1..6);
end:

# Input: non-negative integer s and N
# Output: Probability that start at s
# and survive (avoid all primes 
# in between s and N) beyond N.
# Try:
# Digits:=170:
# A:=evalf(coeff(EPrime(0,1000,a),a,1));
# B:=ProbPrime(0,1000);
# A-B;
# 
# ProbPrime(0,72000);
ProbPrime := proc(s,N) option remember;
local i,A;

if s > N then
	return(1);
elif isprime(s) then
	return(0);
fi:	

A := 1/6*add(ProbPrime(s+i,N),i=1..6);
evalf(A);
end:

##############################################

# Input: s=value to be start, N=cutoff value, 
# BV=Boundary value
# Output: Expected number of rolls starting at s
# with no memory used
# Try: 
# Digits:=3000:
# A := EPrimeNoMem(0,200000,0.);
# B := EPrimeNoMem(0,200000,1000.);
# A-B;
EPrimeNoMem := proc(s,N,BV) option remember;
local i,Ec,L,Avg;

L := [BV$6];
Avg := BV;

for i from N by -1 to s do
	if isprime(i) then
		Ec := 0;
    else
		Ec := 1+Avg;
	fi:
	Avg := Avg+(-L[6]+Ec)/6;
	L := [Ec,op(1..5,L)]; 
od;

Ec;
end:

################################
# Section 2: Upper Bound U
#
# Functions:
# ProbQ(s,n), SubProbQ(s,n),
# ProbR(s,n), SubProbR(s,n),
#
################################

# Input: non-negative integers s and n
# Output: Probability of starting at s
# and land on n while avoiding primes.
# Try:
# ProbQ(0,10000);
ProbQ := proc(s,n) option remember;
if s=n then
	return(1);
fi:

add(SubProbQ(s+i,n),i=1..6)/6;
end:

# Input: non-negative integers s and n
# Output: Probability of starting at s
# and land on n while avoiding primes.
# Try:
# [seq(SubProbQ(0,n),n=0..20)];
# [seq(SubProbQ(14,n),n=14..50)]; 
SubProbQ := proc(s,n) option remember;
local i;

if s=n then
	return(1);
elif isprime(s) or s>n then
	return(0);
fi:

add(SubProbQ(s+i,n),i=1..6)/6;	
end:

# Landing Probability
# Input: non-negative integers s and n
# Output: Same output as ProbQ but this 
# time we do forward calculation
# Try:
# ProbQ(0,10000);
# ProbR(0,10000);
# 
# [seq(ProbQ(0,n),n=0..20)];
# A:=[seq(ProbR(0,n),n=0..40)];
# plot( [seq(n,n=0..40)],A); 
ProbR := proc(s,n) option remember;
if s=n then
	return(1);
fi:

add(SubProbR(s,n-i),i=1..6)/6;
end:

# Input: non-negative integers s and n
# Output: Probability of starting at s
# and land on n while avoiding primes.
# Try:
# [seq(SubProbR(14,n),n=14..50)]; 
SubProbR := proc(s,n) option remember;
local i;

if s=n then
	return(1);
elif isprime(n) or s>n then
	return(0);
fi:
add(SubProbR(s,n-i),i=1..6)/6;	
end:



#############################################
# Test Quantiative Prime Number Theorem
# N:=2000:
# A:=[seq(pi(n),n=17..N)]:
# R1:=plot([seq(n,n=17..N)],A):
# B:=[seq(x/log(x),x=17..N)]:
# R2:=plot([seq(n,n=17..N)],B,color=blue):
# C:=[seq(1.25506*x/log(x),x=17..N)]:
# R3:=plot([seq(n,n=17..N)],C,color=green):
# plots[display]([R1,R2,R3]);

# Test Proposition 3: bound on R_N(n).
# L:=200:P:=300:
# A:=[seq(evalf(ProbR(80000,80000+i)),i=L..P)]:
# R1 := plot([seq(i,i=L..P)],A):
# B:=[seq(evalf(1/32*(5/6)^(pi(80000+i-1)-pi(80000))),i=L..P)]:
# R2 := plot([seq(i,i=L..P)],B,color=blue):
# plots[display]([R1,R2]);

# The upper bound U computation with N=70010
# Digits:=30: N:=70010:
# add( evalf((ithprime(i)-N)*(5/6)^(i-6939)),i=6939..9593);
# 6/5*evalf(int((x-N)*(5/6)^(x/log(x)-1.25506*N/ln(N)  )  
# ,x=100003..infinity ));



################################
# Section 3: Higher Moments
#
# Functions:
# Prob2Prime(s,N),
# E2NoMem(s,N,BV), VarPrime(s,N,L,U),
# E3NoMem(s,N,BV), MoMean3(s,N,L,U),
# E4NoMem(s,N,BV), MoMean4(s,N,L,U),
#  
################################

# Prob surivive for second moment
# Try:
# Digits:=1203:
# A:=E2Prime(0,7000,1000.):
# B:=E2Prime(0,7000,0.)+Prob2Prime(0,7000)*1000;
# A-B;
Prob2Prime := proc(s,N) option remember;
local i;

if s > N then
	return(1);
elif isprime(s) then
	return(0);
fi:	

1/3*add(ProbPrime(s+i,N),i=1..6)
+1/6*add(Prob2Prime(s+i,N),i=1..6);
end:


# Input: s=value to be start, 
# N=cutoff value, BV=Boundary value
# Output: Second moment of number of 
# rolls starting at s with no memory used
# Try: 
# Digits:=1100:
# E2NoMem(0,72000,[0.,0.]);
E2NoMem := proc(s,N,BV) option remember;
local i,Ec1,Ec2,B1,B2,Avg1,Avg2;

B1 := [BV[1]$6];
B2 := [BV[2]$6];
Avg1 := BV[1];
Avg2 := BV[2];

for i from N by -1 to s do
	if isprime(i) then
		Ec1 := 0;
		Ec2 := 0;
    else
		Ec1 := 1+Avg1;
		Ec2 := 1+2*Avg1+Avg2;
	fi:
	Avg1 := Avg1+(-B1[6]+Ec1)/6;
	B1 := [Ec1,op(1..5,B1)]; 
	Avg2 := Avg2+(-B2[6]+Ec2)/6;
	B2 := [Ec2,op(1..5,B2)]; 
od;

Ec2;
end:


# Variance of X
# Input: s=value to be start, N=cutoff value, 
# L= lower bound of E[X^2] at cutoff N,
# U= upper bound of E[X^2] at cutoff N. 
# Try:
# Digits:=1200:
# VarPrime(0,72000,[0.,0.],[412.,47004.]);
# CU and CL agrees with Malinovsky's previous 
# calculation of 6.242778668279075.
VarPrime := proc(s,N,L,U) option remember;
local L1,L2,U1,U2,CL,CU;

L1 := EPrimeNoMem(s,N,L[1]); 
L2 := E2NoMem(s,N,L); 
U1 := EPrimeNoMem(s,N,U[1]);  
U2 := E2NoMem(s,N,U);   

CU := U2-L1^2;  
print("Upper bound of variance", CU);
CL := L2-U1^2; 
print("Lower bound of variance", CL);
print("Their difference is");
CU-CL;
end:

######################################################
# Input: s=value to be start, 
# N=cutoff value, BV=Boundary value
# Output: Third moment of number of 
# rolls starting at s with no memory used
# Try: 
# Digits:=1200:
# E3NoMem(0,72000,[0.,0.,0.]);
E3NoMem := proc(s,N,BV) option remember;
local i,Ec1,Ec2,Ec3,B1,B2,B3,Avg1,Avg2,Avg3;

B1 := [BV[1]$6];
B2 := [BV[2]$6];
B3 := [BV[3]$6];
Avg1 := BV[1];
Avg2 := BV[2];
Avg3 := BV[3];

for i from N by -1 to s do
	if isprime(i) then
		Ec1 := 0;
		Ec2 := 0;
		Ec3 := 0;
    else
		Ec1 := 1+Avg1;
		Ec2 := 1+2*Avg1+Avg2;
		Ec3 := 1+3*Avg1+3*Avg2+Avg3;
	fi:
	Avg1 := Avg1+(-B1[6]+Ec1)/6;
	Avg2 := Avg2+(-B2[6]+Ec2)/6;
	Avg3 := Avg3+(-B3[6]+Ec3)/6;
	B1 := [Ec1,op(1..5,B1)]; 
	B2 := [Ec2,op(1..5,B2)];
	B3 := [Ec3,op(1..5,B3)];	
od;

Ec3;
end:


# Third Moment about the mean
# Input: s=value to be start, N=cutoff value, 
# L= lower bound of E[X^i], i=1,2,3, at cutoff N,
# U= upper bound of E[X^i], i=1,2,3, at cutoff N. 
# Try:
# Digits:=1200:
# MoMean3(0,72000,0.,8277786.);
# MoMean3(0,72000,[0.,0.,0.],[412.,47004.,8277786.]);
# My answer is 52.8836004...
MoMean3 := proc(s,N,L,U) option remember;
local L1,L2,L3,U1,U2,U3,CL,CU;

L1 := EPrimeNoMem(s,N,L[1]);
L2 := E2NoMem(s,N,L);
L3 := E3NoMem(s,N,L);
U1 := EPrimeNoMem(s,N,U[1]); 
U2 := E2NoMem(s,N,U);  
U3 := E3NoMem(s,N,U);

CU := U3-3*L2*L1+2*U1^3; 
print("Upper bound of E[(X-mu)^3]", CU);
CL := L3-3*U2*U1+2*L1^3;  
print("Lower bound of E[(X-mu)^3]", CL);
print("Their difference is");
CU-CL;
end:


# Input: s=value to be start, 
# N=cutoff value, BV=Boundary value
# Output: Fourth moment of number of 
# rolls starting at s with no memory used
# Try: 
# Digits:=1200:
# E4NoMem(0,72000,[0.,0.,0.,0.]);
E4NoMem := proc(s,N,BV) option remember;
local i,Ec1,Ec2,Ec3,Ec4,B1,B2,B3,B4,Avg1,Avg2,Avg3,Avg4;

B1 := [BV[1]$6];
B2 := [BV[2]$6];
B3 := [BV[3]$6];
B4 := [BV[4]$6];
Avg1 := BV[1];
Avg2 := BV[2];
Avg3 := BV[3];
Avg4 := BV[4];

for i from N by -1 to s do
	if isprime(i) then
		Ec1 := 0;
		Ec2 := 0;
		Ec3 := 0;
		Ec4 := 0;
    else
		Ec1 := 1+Avg1;
		Ec2 := 1+2*Avg1+Avg2;
		Ec3 := 1+3*Avg1+3*Avg2+Avg3;
		Ec4 := 1+4*Avg1+6*Avg2+4*Avg3+Avg4;
	fi:
	Avg1 := Avg1+(-B1[6]+Ec1)/6;
	Avg2 := Avg2+(-B2[6]+Ec2)/6;
	Avg3 := Avg3+(-B3[6]+Ec3)/6;
	Avg4 := Avg4+(-B4[6]+Ec4)/6;
	B1 := [Ec1,op(1..5,B1)]; 
	B2 := [Ec2,op(1..5,B2)];
	B3 := [Ec3,op(1..5,B3)];
	B4 := [Ec4,op(1..5,B4)];	
od;

Ec4;
end:


# Fourth Moment about the mean
# Input: s=value to be start, N=cutoff value, 
# L= lower bound of E[X^i], i=1,2,3,4, at cutoff N,
# U= upper bound of E[X^i], i=1,2,3,4, at cutoff N. 
# Try:
# Digits:=1200:
# MoMean4(0,72000,[0.,0.,0.,0.],[412.,47004.,
# 8277786.,2024915563.]);
# My answer is 803.66497701864...
MoMean4 := proc(s,N,L,U) option remember;
local L1,L2,L3,L4,U1,U2,U3,U4,CL,CU;

L1 := EPrimeNoMem(s,N,L[1]);
L2 := E2NoMem(s,N,L);
L3 := E3NoMem(s,N,L);
L4 := E4NoMem(s,N,L);
U1 := EPrimeNoMem(s,N,U[1]); 
U2 := E2NoMem(s,N,U);  
U3 := E3NoMem(s,N,U);
U4 := E4NoMem(s,N,U);

CU := U4-4*L3*L1+6*U2*U1^2-3*L1^4; 
print("Upper bound of E[(X-mu)^4]", CU);
CL := L4-4*U3*U1+6*L2*L1^2-3*U1^4;  
print("Lower bound of E[(X-mu)^4]", CL);
print("Their difference is");
CU-CL;
end: