Help:=proc(): print(` DioA(a,eps) , LogRat(n,eps), ZeRat(N,t,eps) `): end:

Digits:=20:

###From C23.txt

#C23.txt, April 12, 2018, Expmath, Math640, Spring 2018 (Rutgers University)
Help23:=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,n:
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,i:
  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,j,t0:
    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 n,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 from C23.txt


   #DioA(a,eps): enters a number, a, and error eps, outputs the rational number p/q with LEAST q
    #such that abs(a-p/q)<eps
   DioA:=proc(a,eps) local q,i:
   convert(evalf(a),confrac,30,q):
    for i from 1 to nops(q) while abs(evalf(a)-q[i])>eps do od:
      
     q:=q[i]:

      if abs(evalf(a)-q)>eps then
       RETURN(FAIL):
      else
        RETURN(q):
       fi:
     end:

 #LogRat(n,eps): The rational appx. of log(n) with error <eps
 LogRat:=proc(n,eps) 
   option remember:
   if not type(n,posint) then
      RETURN(FAIL):
    fi:
  if n=1 then
    RETURN(0):
  fi:
  DioA(log(n),eps):
 end:


#ZeRat(N,t,eps): the rational appx. to Ze(N,t) with error (of the log(n)) <eps
ZeRat:=proc(N,t,eps) local n1,n2,n:
add(1/n,n=1..N)+2*add(add( 1/sqrt(n1*n2)*cos((LogRat(n2,eps)-LogRat(n1,eps))*t)    , 
n2=n1+1..N), n1=1..N):
end: