#################################################################################
## TruncatedRiemannZeta.txt Save this file as TruncatedRiemannZeta.txt         ##
## to use it, stay in the same directory, get into Maple                       ##
## (by typing: maple <Enter> ) and then type:                                  ##
## read TruncatedRiemannZeta.txt <Enter>                                       ##
## Then follow the instructions given there                                    ##
##                                                                             ##
## Written by Edna Jones, Yukun Yao and Doron Zeilberger, Rutgers University , ## 
## elj44@math.rutgers.edu, yao@math.rutgers.edu and DoronZeil@gmail.com        ##
#################################################################################

#Version of Dec. 17, 2018
with(combinat):
with(linalg):
with(plots):
with(LinearAlgebra):
with(numtheory):

print(`First Written: May 2018`):
print(`Version of Dec. 17, 2018:`):
print():
print(`This is TruncatedRiemannZeta.txt, the  Maple packages for`):
print(`the class project of Spring 2018 Experimental course: `):

print():
print(`The most current version is available on WWW at:`):
print(` https://sites.math.rutgers.edu/~yao/TruncatedRiemannZeta.txt .`):
print(`Please report all bugs to: yao at math dot rutgers dot edu .`):
print():
print(`For an overview of all functions type “Help()”`):
print(`For specific help type "ezra(procedure_name);" `):





Help:=proc(): print(`ZNt(N,t) , ZNtR(N,t) , FindIC(f,t,T,res), OneS(f,t,t0),RN( f,t,t0,N,err), GlobalMinimumZN(N,T) , GlobalMaximumZN(N,T) , Z(N,t) , FindAZ(f,t,T,res) , OS(f,t,P) , ZI(f,t,P,eps) , FindZ(f,t,T,res,eps) , MMS(f,t,T,res,eps) , Ze(N,t), AsyTRZ`): end:


ezra:=proc()

if args=NULL then
print(`The main procedures are: ZNt , ZNtR , FindIC, OneS,RN, GlobalMinimumZN, GlobalMaximumZN, Z , FindAZ, OS, ZI , FindZ , MMS , Ze`):

elif nops([args])=1 and op(1,[args])=ZNt then
print(`ZNt:=(N,t): the real part of zeta_N(1/2+I*t)*zeta_N(1/2-I*t) where zeta_N means the truncated Riemann Zeta function up to N.`):

elif nops([args])=1 and op(1,[args])=ZNtR then
print(`ZNtR(N,t): the real part of zeta_N(1/2+I*t)*zeta_N(1/2-I*t) where zeta_N means the truncated Riemann Zeta function up to N. This procedure doesn’t require calculation involving imaginary numbers. `):

elif nops([args])=1 and op(1,[args])=FindIC then
print(`FindIC(f,t,T,res): inputs a non-negative function f of t, and a pos. nunber T and a resolution, and finds appx. minima. `):

elif nops([args])=1 and op(1,[args])=OneS then
print(`OneS(f,t,t0): applies one step of Raphson-Newton. `):

elif nops([args])=1 and op(1,[args])=RN then
print(`RN(f,t,t0,N,err): estimates the closest zero of f of t to t0 as soon as two consecutive iterations are less than err apart, or we reached N iterations and return FAIL. `):

elif nops([args])=1 and op(1,[args])=GlobalMinimumZN then
print(`GlobalMinimumZN(N,T) that inputs a positive integer N and a large positive integer T, and outputs the location and value of the global minimum of Z(N,t) in 0<=t<=T. `):

elif nops([args])=1 and op(1,[args])=GlobalMaximumZN then
print(`GlobalMaximumZN(N,T) that inputs a positive integer N and a large positive integer T, and outputs the location and value of the global maximum of Z(N,t) in 0<=t<=T. `):

elif nops([args])=1 and op(1,[args])=RN then
print(`RN(f,t,t0,N,err): estimates the closest zero of f of t to t0 as soon as two consecutive iterations are less than err apart, or we reached N iterations and return FAIL. `):

elif nops([args])=1 and op(1,[args])=Z then
print(`Z(N,t): add(1/n,n=1..N) + add(n1^(1/2+I*t)*n2^(1/2-I*t)= ADD(sqrt(n1*n2)* 2(cos(log(n1)-log(n2))*t),1<=n1<n2<=N), n^(I*t)=exp(log(n)*I*t), i.e. we want to calculate zeta_N(1/2+I*t)*zeta_N(1/2-I*t): w/o I.`):

elif nops([args])=1 and op(1,[args])=FindAZ then
print(`FindAZ(f,t,T,res): inputs a function f of t, a pos. integer T, and a pixel-size res and finds all the intervals [res*(i-1),res*i] where f changes sign.`):

elif nops([args])=1 and op(1,[args])=OS then
print(`OS(f,t,P): inputs a function f of t and an interval P where the sign of f is opposite and outputs the left or right interval with the same property. :`):

elif nops([args])=1 and op(1,[args])=ZI then
print(`ZI(f,t,P,eps): inputs f,t,P, finds an approximation to the zero of f in P. :`):

elif nops([args])=1 and op(1,[args])=FindZ then
print(`FindZ(f,t,T,res,eps): inputs a function f of t, a pos. integer T, and a pixel-size res and an error eps, and outputs the (hopefully) all the approximate zeros in [0,T].`):

elif nops([args])=1 and op(1,[args])=MMS then
print(`MMS(f,t,T,res,eps): inputs a nice differentiable function f of t and a positive number T and outputs ListMin and ListMax with [location, value]. :`):

elif nops([args])=1 and op(1,[args])=Ze then
print(`Ze(N,t): add(1/n,n=1..N) + add(n1^(1/2+I*t)*n2^(1/2-I*t)= ADD(sqrt(n1*n2)* 2(cos(log(n1)-log(n2))*t),1<=n1<n2<=N)n^(I*t)=exp(log(n)*I*t) EXACT expression of zeta_N(1/2+I*t)*zeta_N(1/2-I*t): w/o I.`):

elif nops([args])=1 and op(1,[args])=AsyTRZ then
print(`AsyTRZ(N,t): inputs large N and t, using integral to estimate truncated Riemann Zeta function.`):







fi:
end:



### Following is from C22.txt ###

#ZNt(N,t): the real part of zeta_N(1/2+I*t)*zeta_N(1/2-I*t) where zeta_N means the truncated Riemann Zeta function up to N.
ZNt:=proc(N,t) local n:
Re(add(1/n^(1/2+I*t),n=1..N)*add(1/n^(1/2-I*t),n=1..N));
end:

#ZNtR(N,t): the real part of zeta_N(1/2+I*t)*zeta_N(1/2-I*t) where zeta_N means the truncated Riemann Zeta function up to N. This procedure doesn’t require calculation involving imaginary numbers. 
#add(1/n,n=1..N) + add(n1^(1/2+I*t)*n2^(1/2-I*t)= ADD(sqrt(n1*n2)* 2(cos(log(n1)-log(n2))*t),1<=n1<n2<=N)
#n^(I*t)=exp(log(n)*I*t)
# zeta_N(1/2+I*t)*zeta_N(1/2-I*t): w/o I
ZNtR:=proc(N,t) local n1,n2:
add(evalf(1/n),n=1..N)+2*add(add( evalf(1/sqrt(n1*n2)*cos((log(n2)-log(n1))*t))    , n2=n1+1..N), n1=1..N):
end:

#FindIC(f,t,T,res): inputs a non-negative function f of t, and a pos. number T
#and a resolution, and finds appx. minima
FindIC:=proc(f,t,T,res) local i,ej,L:
ej:=evalf([seq(subs(t=res*i,f),i=0..trunc(T/res))]):

L:=[]:

for i from 2 to nops(ej)-1 do
  if ej[i]<ej[i-1] and ej[i]<ej[i+1] then
   L:=[op(L),i*res]:
  fi:
od:

L:

end:

#OneS(f,t,t0): applies one step of Raphson-Newton
OneS:=proc(f,t,t0) 
t0- subs(t=t0,f)/subs(t=t0,diff(f,t)):
end:

#RN(f,t,t0,N,err): estimates the closest zero of f of t to t0 as soon as two consecutive
#iterations are less than err apart, or we reached N iterations and return FAIL
RN:=proc(f,t,t0,N,err) local t1,t2,t3,i:

t1:=t0:
t2:=evalf(OneS(f,t,t1)):
for i from 2 to N do
 t3:=evalf(OneS(f,t,t2)):
 t1:=t2:
 t2:=t3:
  if abs(t1-t2)<err then
     RETURN(t2):
  fi:
od:

FAIL:
end:

### End of C22.txt

### Following is from HW23 ###

##Problem 1## Write a program
#GlobalMinimumZN(N,T) that inputs a positive integer N and a large positive integer T, and outputs the location and value of the global minimum of Z(N,t) in 0<=t<=T. What are GlobalMinimumZN(N,1000) for N from 2 to 10?

GlobalMinimumZN:=proc(N,T) local f,t,L,i,m,min:
f:=diff(Z(N,t),t):
L:=MMS(f,t,T,0.01,10^(-9))[1]:
L:=[[0,Z(N,0)],op(L),[T,Z(N,T)]]:
m:=0:
min:=Z(N,0):
for i from 2 to nops(L) do
  if L[i][2]<min then
    min:=L[i][2]:
    m:=L[i][1]:
  fi:
od:
[m,min]:
end:

#seq(GlobalMinimumZN(n, 100), n = 2 .. 10) = [38.5250612, -.98025814346854718452], [1.67718748, -2.3298776161770970476], [1.41808189, -4.0174514965565673205], [1.26701788, -5.9672343446878533521], [1.16592338, -8.1276795873644019889], [1.09248242, -10.462482947038641401], [1.03614092, -12.944965633528193395], [.991204860, -15.554820131834691620], [.954307496, -18.276142769321468855].


##Problem 2## Write a program
#GlobalMaximumZN(N,T) that inputs a positive integer N and a large positive integer T, and outputs the location and value of the global maximum of Z(N,t) in 0<=t<=T. What are GlobalMaximumZN(N,1000) for N from 2 to 10?

GlobalMaximumZN:=proc(N,T) local f,t,L,i,m,max:
f:=diff(Z(N,t),t):
L:=MMS(f,t,T,0.01,10^(-9))[2]:
L:=[[0,Z(N,0)],op(L),[T,Z(N,T)]]:
m:=0:
max:=Z(N,0):
for i from 2 to nops(L) do
  if L[i][2]>max then
    max:=L[i][2]:
    m:=L[i][1]:
  fi:
od:
[m,max]:
end:

#seq(GlobalMaximumZN(n, 100), n = 2 .. 10) = [0, 2.9142135623730950488], [0, 5.2187440150134059439], [0, 7.7532010653895792329], [0, 10.443695163417451168], [0, 13.249009863126475074], [0, 16.143387091123723596], [0, 19.109459695891483099], [0, 22.134862006964257700], [0, 25.210419904701369970].

### End of HW23 ###

### Following is from C23.txt ###

Digits:=20:

#add(1/n,n=1..N) + add(n1^(1/2+I*t)*n2^(1/2-I*t)= ADD(sqrt(n1*n2)* 2(cos(log(n1)-log(n2))*t),1<=n1<n2<=N)
#n^(I*t)=exp(log(n)*I*t)
# zeta_N(1/2+I*t)*zeta_N(1/2-I*t): w/o I
Z:=proc(N,t) local n1,n2:
add(evalf(1/n),n=1..N)+2*add(add( evalf(1/sqrt(n1*n2)*cos((log(n2)-log(n1))*t))    , 
n2=n1+1..N), n1=1..N):
end:

#FindAZ(f,t,T,res): inputs a function f of t, a pos. integer T, and a pixel-size res
#and finds all the intervals [res*(i-1),res*i] where f changes sign
FindAZ:=proc(f,t,T,res) local L,i,P:
L:=[]:
for i from 0 to trunc(T/res)-1 do
 P:=[i*res,(i+1)*res]:
  if evalf(subs(t=P[1],f)*subs(t=P[2],f))<0 then
    L:=[op(L),P]:
  fi:
od:
L:

end:

#OS(f,t,P): inputs a function f t and an interval P where the sign of f is opposite
#and outputs the left or right interval with the same property
OS:=proc(f,t,P) :
 if evalf(subs(t=P[1],f)*subs(t=P[2],f))>=0 then
    FAIL:
  elif evalf(subs(t=P[1],f)*subs(t=(P[1]+P[2])/2,f))<0 then
   [P[1],(P[1]+P[2])/2]:
  else
 
  [(P[1]+P[2])/2,P[2]]:
fi:
end:

#ZI(f,t,P,eps): inputs f,t,P, finds an approximation to the zero of f in P
ZI:=proc(f,t,P,eps) local P1,N,i:
 

N:=ceil(log[2](abs((P[2]-P[1])/eps))):
P1:=P:

for i from 1 to N do
  P1:=OS(f,t,P1):
od:
evalf((P1[1]+P1[2])/2, trunc(log[10](1/eps))):
end:

#FindZ(f,t,T,res,eps): inputs a function f of t, a pos. integer T, and a pixel-size res
#and an error eps,
#and outputs the (hopefully) all the approximate zeros in [0,T]
FindZ:=proc(f,t,T,res,eps) local L:
L:=FindAZ(f,t,T,res):
[seq(ZI(f,t,L[i],eps),i=1..nops(L))]:
end:
    
#MMS(f,t,T,res,eps): inputs a nice differentiable function f of t and a positive number T and outputs ListMin and ListMax with [location, value]
#This version, by Edna Jones, corrects a previous buggy version
MMS := proc(f,t,T,res,eps) local g, L, L1:
    g := diff(f,t);
    L := FindZ(g,t,T,res,eps);
    L1 := [seq([t0, evalf(subs(t=t0, f))], t0 in L)];
    if L1[1][2] < L1[2][2] then
        [seq(L1[2*j+1], j=0..(nops(L)-1)/2)], [seq(L1[2*j], j=1..nops(L)/2)];
    else
        [seq(L1[2*j], j=1..nops(L)/2)], [seq(L1[2*j+1], j=0..(nops(L)-1)/2)];
    fi:
end:




#add(1/n,n=1..N) + add(n1^(1/2+I*t)*n2^(1/2-I*t)= ADD(sqrt(n1*n2)* 2(cos(log(n1)-log(n2))*t),1<=n1<n2<=N)
#n^(I*t)=exp(log(n)*I*t) EXACT expression
# zeta_N(1/2+I*t)*zeta_N(1/2-I*t): w/o I
Ze:=proc(N,t) local n1,n2:
add(1/n,n=1..N)+2*add(add( 1/sqrt(n1*n2)*cos((log(n2)-log(n1))*t)    , 
n2=n1+1..N), n1=1..N):
end:

### End of C23.txt ###


#AsyTRZ(N,t): inputs large N and t, using integral to estimate truncated Riemann Zeta function. 
AsyTRZ:=proc(N,t) local n, n1, n2,ans:
ans:=int(1/n, n=1..N) + 2*int(int(1/sqrt(n1*n2)*cos(t*(ln(n2)-ln(n1))), n1=1..n2), n2=1..N):
Re(evalf(ans)):
end: