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



with(plots):

print(`First Written: Sept. 2025: tested for Maple 2020 `):
print(`This Version: Sept. 15, 2025  `):
print():
print(`This is  Hodes.txt, a Maple package`):
print(`accompanying  Doron Zeilberger's  article: `):
print(`Every Fifth  Real Number is Evil `):

print():
print(`The most current version is available on WWW at:`):
print(` http://sites.math.rutgers.edu/~zeilberg/tokhniot/Hodes.txt .`):
print(`Please report all bugs to: DoronZeil at gmail dot com .`):

print():
print(`For general help, and a list of the MAIN functions,`):
print(` type "ezra();". For specific help type "ezra(procedure_name);" `):
print(`For a list of the Story functions type: ezraSt();`):
print():





ezraSt:=proc()
if args=NULL then

print(`The Story procedures are: EvilStory `):

else
ezra(args):

fi:


end:

ezra1:=proc()
if args=NULL then

print(`The SUPPORTING procedures are:  Di, EL, GR, MomE, PF, PFg, RL, Stat `):

else
ezra(args):

fi:


end:


ezra:=proc()
if args=NULL then

print(` Hodes.txt: A Maple package for  looking studying evil real numbers`):

print(`The MAIN procedures are:   EP, EvilF, Gxr, Gxtr, MomElist, MomElistC, PiPrimes, Simu `):
print(``):


elif nargs=1 and args[1]=Di then
print(`The list of the first N base-b digits of x (ignoring initial zeroes). Try:`):
print(` Di(Pi,10,20);  `):


elif nargs=1 and args[1]=EL then
print(`EL(L): inputs a list of integers that start with a non-zero integer outputs the first time the partial sum is K or returns FAIL. Try:`):
print(`EL(Di((sqrt(5)-1)/2,1000),666);`):


elif nargs=1 and args[1]=EPb then
print(`EP(x,b,K): the smallest integer k such that the first k base-b  digits  of x add up to K. Try:`):
print(`EP((sqrt(5)-1)/2,10,666);`):


elif nargs=1 and args[1]=EvilF then
print(`EvilF(f,x,b,N,K): inputs an expression f in the variable x finds all integer i<=N such that the base-b expansion of subs(x=i,f) (if irrational)  has`):
print(`a partial sum equal K, followed by the place. Also returns the fraction of those that have it. Try:`):
print(`EvilF(sqrt(x),x,10,10,666);`):

elif nargs=1 and args[1]=EvilStory then
print(`EvilStory(n,b,M): an article about the expectation, variance, and central moments about the mean up to M and more`):
print(`Try: `):
print(`EvilStory(n,b,4);`):

elif nargs=1 and args[1]=GR then
print(`GR(DS,n): guesses a rational function`):
print(`for the data set DS {[a,b]}, Try:`):
print(`GR({seq([i,(i+1)/(i+3)],i=1..10)},n);`):

elif nargs=1 and args[1]=Gxr then
print(`Gxr(x,r): The rational function whose coefficient of x^i in its Taylor expansion aronund x=0 is the probability that a random real number, expressed in base r,`):
print(`has a partial sum that equals i. Try:`):
print(`Gxr(x,10);`):


elif nargs=1 and args[1]=Gxtr then
print(`Gxtr(x,t,r): The rational function in x and t whose coefficient of x^i in its Taylor expansion aronund x=0 is the probability generating function, in t, whose`):
print(`coefficient of t^j is the probability that partial sum up to the j-th digit,  equals i . Try:`):
print(`Gxtr(x,t,10);`):


elif nargs=1 and args[1]=MomE then
print(`MomE(n,b,m,: The leading asymptotic (mod. exponentially small corrections) for the m-th moment of the r.v. "location of evil spot" for a random evil real number. Try:`):
print(`MomE(n,b,2,10); `):

elif nargs=1 and args[1]=MomElist then
print(`MomElist(n,b,M): The expectation, and i-th  moment up to the M of the r.v. location of a random evil real number in base b (symbolic b). Try:`):
print(`MomElist(n,b,4);`):

elif nargs=1 and args[1]=MomElistC then
print(`MomElistC(n,b,M): The expectation, and i-th CENTRAL moment up to the M of the r.v. location of a random evil real number in base b (symbolic b). Try:`):
print(`MomElistC(n,b,4);`):

elif nargs=1 and args[1]=PiPrimes then
print(`PiPrimes(b,K,N): The list of all primes p among the first N primes, such that p*Pi is evil, followed by their frequency, and the statistics of the evil-location in base b. Try:`):
print(`PiPrimes(10,666,100);`):

elif nargs=1 and args[1]=PF then
print(`PF(f,x) :`):
print(`inputs a rational function whose denominator has a power of (1-x), say (1-x)^r, outputs the vector [A1,A2,..Ar] such that`):
print(`the partial fraction expansion of f is A1/(1-x)+ ...+`):

elif nargs=1 and args[1]=PFg then
print(`inputs a rational function that depends on a parameter b, whose denominator has a power of (1-x), say (1-x)^r, outputs the vector [A1,A2,..Ar] such that`):
print(`the partial fraction expansion of f is A1/(1-x)+ ...+ with the symbolic b. K is a guessing parameter. Try:`):
print(`PFg(normal(subs(t=1,diff(-t^2*(-1+x^(r-1))^2*x^2*(t*x^r-2*r*x+2*r-t)/(t*x^r-r*x+r-t)^2/(r-1)/(-1+x),t))),x,r,10);`):

elif nargs=1 and args[1]=RL then
print(`RL(r,N): a random list of integers in {0,1,..,r-1} that start with a non-zero entry. Try:`):
print(`RL(10,1000);`):

elif nargs=1 and args[1]=Simu then
print(`Simu(r,N,K,x): generates K random lists in {0,...,r-1} that start with non-zero, of length N, that outputs`):
print(`the ratio of those that have partial sum K, and  for those the prob. gen. function in x`):
print(`of the location where K is achieved. Try:`):
print(`Simu(10,2000,666,x);`):

elif nargs=1 and args[1]=Stat then
print(`Stat(F,t,K): Given a generating function F of t, for a r.v. `):
print(` finds the statistical information: av,sd,skewness, kurtosis, .... Try:`):
print(`Stat(((1+x)/2)^1000,x,6);`):

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

fi:


end:


###From Stat
#Stat(F,t,K): Given a generating function F of t, for a r.v. 
# finds the statistical information
Stat:=proc(F,t,K) local gu,mu,sd,k,f:
f:=evalf(F/subs(t=1,F),20):
mu:=subs(t=1,diff(f,t)):
gu:=[mu]:
if K=1 then
 RETURN(gu):
fi:
f:=f/t^mu:
f:=t*diff(t*diff(f,t),t):
sd:=sqrt(subs(t=1,f)):
gu:=[op(gu),sd]:
for k from 3 to K do
 f:=t*diff(f,t):
gu:=[op(gu),subs(t=1,f)/sd^k]:
od:

evalf(gu,20):
end:


####
#The list of the first N decimal digits of x
Di:=proc(x,b,N) local p,L,y,i,d:
p:=trunc(log[b](evalf(x)))+1:
if p<0 then
p:=p-1:
fi:

L:=[]:
y:=x/b^p:
for i from 1 to N do
d:=trunc(b*y):
L:=[op(L),d]:
y:=(b*y-d):
od:
L:

if L[1]=0 then
L:=[op(2..nops(L),L)]:
 fi:
 L:
end:




#EL(L): inputs a list of integers that start with a non-zero integer outputs the first time the partial sum is K or returns FAIL.
#Try:
#EL(Di((sqrt(5)-1)/2,1000),666);
EL:=proc(L,K) local i:

if convert(L,`+`)<K then
print(`Sum of list is less than`, K):
print(L):
 RETURN(FAIL):
 fi:

for i from 1 to nops(L) do
if convert([op(1..i,L)],`+`)=K then
 RETURN(i):
 fi:
 od:
 FAIL:
 end:




#EP(x,b,K): the smallest integer k such that the first k base-b digits of x (ignoring leading 0) add up to K. Try:
#EP((sqrt(5)-1)/2,10,666);
EP:=proc(x,b,K) local L,R:

R:=K*3/(b-1):
L:=Di(x,b,R):
if convert(L,`+`)<K then
L:=Di(x,b,3*R/2):
fi:
EL(L,K):
end:





#Gxr(x,r): The rational function whose coefficient of x^i in its MacLaurin expansion is the probability that a random real number, expressed in base r,
#has a partial sum that equals i. Try:
#Gxr(x,10);
Gxr:=proc(x,r):
normal(1+x*(1-x^(r-1))/(r-1)/(1-x)+
x*(1-x^(r-1))/(1-x)/(r-1)* (1/(1-(1-x^r)/(1-x)/r)) * (x*(1-x^(r-1))/(1-x)/r)):
end:


#Gxr(x,r): The rational function in x and t, whose coefficient of x^i in its Maclaurin expansion is the probability generating function
#has a partial sum that equals i. Try:
#Gxr(x,10);

Gxtr:=proc(x,t,r):
normal(
1+x*(1-x^(r-1))/(r-1)/(1-x)*t+
t*x*(1-x^(r-1))/(1-x)/(r-1)* (1/(1-t*(1-x^r)/(1-x)/r)) * (x*t*(1-x^(r-1))/(1-x)/r)):
end:



#RL(b,N): a random list of integers in {0,1,..,r-1} that start with a non-zero entry. Try:
#RL(10,1000);
RL:=proc(b,N) local i,ra:
ra:=rand(0..b-1):
[rand(1..b-1)(), seq(ra(),i=1..N-1)]:
end:



#Simu(b,N,K,x): generates K random lists in {0,...,b-1} that start with non-zero, of length N, that outputs
#the ratio of those that have partial sum K, and  for those the prob. gen. function in x
#of the location where K is achieved
Simu:=proc(b,N,K,x) local f,L,lu,i:

f:=0:

for i from 1 to K do
L:=RL(b,N):
lu:= EL(L,N):
if lu<>FAIL then
 f:=f+x^lu:
 fi:
 od:
 [evalf(subs(x=1,f)/K,10),f]:
 end:
 





#EvilF(f,x,b,N,K): inputs an expression f in the variable x finds all integer 1<=i<=N such that the base-b expansion of subs(x=i,f) (if irrational)  has
#a partial sum equal K, followed by the place. Also returns the fraction of those that have it. Try:
#EvilFb(sqrt(x),x,10,10,666);
EvilF:=proc(f,x,b,N,K) local i,M,ku,ha,tot:
tot:=0:
M:=[]:
for i from 1 to N do
ku:=simplify(subs(x=i,f)):

 if not (type(ku,integer) or type(ku,fraction)) then
  tot:=tot+1:
  ha:=EP(ku,b,K):
  if ha<>FAIL then
    M:=[op(M),[i,ha]]:
  fi:
 fi:
 od:
 [M, evalf(nops(M)/tot)]:
 end:





###Start Guess a rational function
#GR1(DS,d1,d2,n): guesses a rational function in n with top degree d1 and bottom degree d2
#for the data set DS {[a,b]}
GR1:=proc(DS,d1,d2,n) local a,b,T,B,R,eq,var,i:

if nops(DS)<=d1+d2+3 then
  RETURN(FAIL):
  fi:

T:=add(a[i]*n^i,i=0..d1):
B:=add(b[i]*n^i,i=0..d2-1)+n^d2:
R:=T/B:
var:={seq(a[i],i=0..d1),seq(b[i],i=0..d2-1)}:
eq:={seq(numer( subs(n=DS[i][1],R)-DS[i][2]),i=1..nops(DS))}:
var:=solve(eq,var):
if var=NULL then
 RETURN(FAIL):
 else
  factor(subs(var,R)):
  fi:
  end:


#GR(DS,n): guesses a rational function
#for the data set DS {[a,b]}, Try:
#GR({seq([i,(i+1)/(i+3)],i=1..10)},n);
GR:=proc(DS,n) local d1,d2,hope:

for d1 from 0 to nops(DS)-4 do
 for d2 from 0 to nops(DS)-d1-4 do
   hope:=GR1(DS,d1,d2,n):
     if hope<>FAIL then
         RETURN(hope):
	   fi:
	    od:
	    od:
	    FAIL:
	    end:


GP:=proc(DS,n) local k,f,a,eq,var,i,i1:
k:=nops(DS)-4:
f:=add(a[i]*n^i,i=0..k):
var:={seq(a[i],i=0..k)}:
eq:={seq(subs(n=DS[i1][1],f)-DS[i1][2],i1=1..nops(DS))}:
var:=solve(eq,var):
if var=NULL then
 RETURN(FAIL):
fi:

subs(var,f):

end:

##end Guess a rational function

Ord1:=proc(f,x) local f1,i:
f1:=denom(f):

for i from 0 do
if subs(x=1,f1)<>0 then
   RETURN(i):
fi:

f1:=expand(diff(f1,x)):

od:
end:

#inputs a rational function whose denominator has a power of (1-x), say (1-x)^r, outputs the vector [A1,A2,..Ar] such that
#the partial fraction expansion of f is A1/(1-x)+ ...+
PF:=proc(f,x) local r,i,A,L,f1:

r:=Ord1(f,x):

L:=[]:

f1:=f:

for i from r by -1 to 1 do
A:=limit(normal((1-x)^i*f1),x=1):
L:=[op(L),A]:
f1:=normal(f1-A/(1-x)^i):
od:
L:

end:



#inputs a rational function that depends on a parameter b, whose denominator has a power of (1-x), say (1-x)^r, outputs the vector [A1,A2,..Ar] such that
#the partial fraction expansion of f is A1/(1-x)+ ...+ with the symbolic b. Try:
#PFg(normal(subs(t=1,diff(-t^2*(-1+x^(r-1))^2*x^2*(t*x^r-2*r*x+2*r-t)/(t*x^r-r*x+r-t)^2/(r-1)/(-1+x),t))),x,r,10);
PFg:=proc(f,x,r,K) local M,i,R,k,M1,lu,i1:

M:=[seq([i,PF(normal(subs(r=i,f)),x)],i=6..6+K)]:

if nops({seq(nops(M[i][2]),i=1..nops(M))})<>1 then
 RETURN(FAIL):
 fi:

k:=nops(M[1][2]):

R:=[]:

for i from 1  to k do
M1:=[seq([M[i1][1],M[i1][2][i]],i1=1..nops(M))]:
lu:=GR(M1,r):

if lu=FAIL then
 print(`Make `, K, `bigger `):
 RETURN(FAIL):
 fi:
R:=[op(R),lu]:
od:

R:

end:


#MomE(n,b,m,: The leading asymptotic (mod. exponentially small corrections) for the m-th moment of the r.v. "location of evil spot" for a random evil real number. Try:
#MomE(n,b,2,10); 
MomE:=proc(n,b,m) local lu,i,t,x,gu,ku,i1:
option remember:
lu:=Gxtr(x,t,b)/(2/b):
for i from 1 to m do
lu:=normal(t*diff(lu,t)):
od:
lu:=normal(subs(t=1,lu)):

gu:=PFg(lu,x,b,3+3*m):

if gu=FAIL then
 RETURN(FAIL):
 fi:

if nops(gu)<>m+1 then
  RETURN(FAIL):
 fi:

ku:=expand(add(gu[i]*expand(binomial(n+m+1-i,m+1-i)),i=1..m+1)):

add(factor(coeff(ku,n,i1))*n^i1,i1=0..degree(ku,n)):

end:

#MomElist(n,b,M): The expectation, and i-th moment up to the M of the r.v. location of a random evil real number in base b (symbolic b). Try:
#MomElist(n,b,4);
MomElist:=proc(n,b,M) local m:

[seq(MomE(n,b,m),m=1..M)]:

end:


#MomElistC(n,b,M): The expectation, and i-th CENTRAL moment up to the M of the r.v. location of a random evil real number in base b (symbolic b). Try:
#MomElistC(n,b,4);

MomElistC:=proc(n,b,M) local gu,mu,lu,r,i1,lu1:
gu:=MomElist(n,b,M):

mu:=gu[1]:

lu:=[mu]:
#va:=expand(M[2]-M[1]^2):
#va:=add(factor(coeff(va,n,i1))*n^i1,i1=0..degree(va,n)):
for r from 2 to M do
lu1:=expand(add(gu[i1]*(-mu)^(r-i1)*binomial(r,i1),i1=1..r)+(-mu)^r):
lu1:=add(factor(coeff(lu1,n,i1))*n^i1,i1=0..degree(lu1,n)):
lu:=[op(lu),lu1]:
od:
lu:
end:



#EvilStory(n,b,M): an article about the expectation, variance, and central moments about the mean up to M
EvilStory:=proc(n,b,M) local f,gu,x,ru,i,t0:

t0:=time():

print(` The Statistical Distribution of the Evil Location of Evil Real Numbers for any Base`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(` Definition (Isaac Hodes): A real number is K-evil in base b if its base-b representation has the property that one of the partial sums is K`):

print(``):
print(`Theorem 1: Let a[b](n) be the probability that a random real number in base b has one of its partial sums equal to n, then`):
print(``):

f:=Gxr(x,b):
print(Sum(a[b](n)*x^n,n=0..infinity)=f):
print(``):
print(`and in Maple notation `):
print(``):
lprint(Sum(a[b](n)*x^n,n=0..infinity)=f):
print(``):
print(`it follows that as n goes to infinity, the prob. that one of the partial sums of a random real number equals n, goes to 2/b`):
print(``):


print(``):
print(`Theorem 2: Let A[b](n,k) be the probability that a random real number in base b has its k-th partial sum equal to n , then`):
print(``):

f:=Gxtr(x,t,b):
print(Sum(Sum(A[b](n,k)*x^n*t^k,k=0..infinity),n=0..infinity)=f):
print(``):
print(`and in Maple notation `):
print(``):
lprint(Sum(Sum(A[b](n,k)*x^n*t^k,k=0..infinity),n=0..infinity)=f):
print(``):
gu:=MomElistC(n,b,M):
print(``):
print(`it follows that as n goes to infinity, the expected evil location of an n-evil number is asymptotic to`, gu[1]):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu[1]):

print(`In fact we have explicit expressions for all central moments up to the `, M):

print(`The variance is asymptotic (with exponentially decaying error) to`, gu[2]):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu[2]):
print(``):
ru:=[0,1]:

if nops(gu)>=3 then
print(``):
print(`The skewness is`, gu[3]):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu[3]):
print(``):
if degree(gu[3],n)<3/2*degree(gu[2],n) then
 ru:=[op(ru),0]:
fi:
fi:


if nops(gu)>=4 then
print(``):
print(`The kurtosis is`, gu[4]):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu[4]):
print(``):
if degree(gu[4],n)=2*degree(gu[2],n) then
 ru:=[op(ru),lcoeff(gu[4],n)/lcoeff(gu[2],n)^2]:
fi:
fi:

for i from 5 to M do

print(`The `, i, `-th central moment of an n-evil location for random evil real number in base b is`):
print(``):
print(gu[i]):
print(``):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu[i]):
print(``):
if degree(gu[i],n)<(i/2)*degree(gu[2],n) then
ru:=[op(ru),0]:
elif degree(gu[i],n)=(i/2)*degree(gu[2],n) then
ru:=[op(ru),lcoeff(gu[i],n)/lcoeff(gu[2],n)^(i/2)]:
fi:
od:

print(`To sum up the asymptotics, up to exponentially decaying terms for the first `, M, `central moments are`):
print(``):
print(gu):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu):
print(``):
print(`The limits of the scaled moments up to the `, M, `-th are `):
print(``):
print(ru):
print(``):
print(`This confirms that our random variable is asymptotically normal `):
print(``):
print(`----------------------------------`):
print(``):
print(`This ends this article that took`, time()-t0, `seconds to generate. `):


end:



with(numtheory):

#PiPrimes(b,K,N): The list of all primes p among the first N primes, such that p*Pi is evil, followed by their frequency, and the statistics of the evil-location in base b. Try:
#PiPrimes(10,666,100);
PiPrimes:=proc(b,K,N) local L, i,t,gu,f,p:

L:=[]:
f:=0:
for i from 1 to N do
p:=ithprime(i):
gu:=EP(p*Pi,b,K):
if gu<>FAIL then
 L:=[op(L),[p,gu]]:
 f:=f+t^gu:
fi:
od:
[L,evalf(nops(L)/N),Stat(f,t,4)]:
end: