#C23.txt, April 12, 2018, Expmath, Math640, Spring 2018 (Rutgers University)
Help:=proc(): print(`Z(N,t), FindAZ(f,t,T,res),OS(f,t,P),ZI(f,t,P,eps), FindZ(f,t,T,res,eps) `): 
 print(` MMS(f,t,T,res,eps), Ze(N,t) `):
end:

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 pos. number T
     #outputs ListMin, ListMax, [location,value]
    MMS:=proc(f,t,T,res,eps) local g,L,L1:
    g:=diff(f,t):
    L:=FindZ(g,t,T,res,eps):
    L1:=[seq([L[i],evalf(subs(t=L[i],f))],i=1..nops(L))]:
    if L1[1][2]<L1[2][2] then
      [seq(L1[2*i-1],i=1..ceil(nops(L)/2))],[seq(L1[2*i],i=1..nops(L)/2)]
     else
    [seq(L1[2*i],i=1..nops(L)/2)], [seq(L1[2*i-1],i=1..ceil(nops(L)/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: