(**************************************************************)
(*
    ARIBAS code for the calculation of the Riemann zeta function
    and its zeros
    File author:
        (C) Otto Forster 2000/2018 <forster AT math.lmu.de>
    Date of last change: 
        2018-01-24

    This code is placed under the GNU General Public License
*)
(**************************************************************)
(*
Functions:

zeta(s: complex): complex;
    Riemann zeta function of a complex argument.
    Complex numbers are represented as pairs (x,y) of real numbers
    Accuracy about 8 decimal places (although more are shown)

Zet(t: real): real;
    Real function such that |Zet(t)| = |zeta((0.5,t))|
    Therefore zeros (of odd multiplicity) of the zeta function 
    on the critical line Re(s) = 0.5 can be located by Zet
    using the mean value theorm

zetazero(t1,t2: real): real;
    If Zet(t1) and Zet(t2) have different sign, then
    this function returns the imaginary part of a zero
    of zeta on the critical line

NTapprox(T: real): real;
    Approximate number of zeros of the zeta function with
    imaginary part 0 < t <= T

gram(n: integer): real;
    n-th Gram point on the critical line (n >= 0). 
    In general (but there are exceptions), 
    between two consecutive Gram points gram(n-1) and gram(n) 
    there is exactly one zero of zeta.
    Therefore  NTapprox(gram(n)) ~= n+1.
*)
(**********************************************************)
(*
** Examples:

==> zeta((2,0)).
-: (1.64493406684749307, 0.000000000000000000)

==> pi**2/6.
-: 1.64493406684822644

==> zeta((-1,0)).
-: (-0.0833333333332963906, 0.000000000000000000)

==> -1/12.
-: -0.0833333333333333333

==> zeta((0.5,100)).
-: (2.69261988568128231, -0.0203860296026174821)

==> T := 1000.
-: 1000

==> zeta((0.5,T)).
-: (0.356334367194383065, 0.931997831233114642)

==> Zet(T).
-: 0.997794637521694977

==> Zet(T-1).
-: -0.711847138112776643

==> t := zetazero(T-1,T).
-: 999.791571553796530

==> zeta((0.5,t)).
-: (2.69435146903045255e-9, -1.63677023707042960e-8)

==> NTapprox(T).
-: 648.616235312967350

==> t1 := gram(648).
-: 1000.47557551225284

==> t2 := gram(647).
-: 999.236223445226699

==> zetazero(t1,t2).
-: 999.791571553043058

*)
(**********************************************************)
var
    BMax: integer;
    Bern: array of array[2];
    Ber: array of real;
end;
(**********************************************************)
function bern_ini(anz: integer): integer;
external
    BMax: integer;
    Ber, Stir: array of real;
    Bern: array of array[2];
var
    i, len, bnum, bden: integer;
begin
    Bern := bernoulli(anz);
    BMax := anz;
    len := anz+1;
    Ber := alloc(array,len,0.0);
    for i := 0 to BMax do
        (bnum, bden) := Bern[i]
        Ber[i] := bnum/bden;
    end;
    return anz;
end;
(**********************************************************)
(*
** Calculates the Bernoulli numbers B(2*k), k = 1,2,...,n
** and stores them in the external array Bern
*)
function bernoulli(n): array of array[2]
var
    Bern: array of array;
    k: integer;
begin
    n := max(n,1);
    Bern := alloc(array,n+1,());
    Bern[0] := (1,1);
    Bern[1] := (1,6);
    for k := 2 to n do
        nextbern(Bern,k);
    end;
    return Bern;
end;
(*-----------------------------------------------------------*)
(*
** Auxiliary function for bernoulli(n).
** Calculates B(2*k) from the previous ones.
** Supposes k >= 2.
**
** The following recursion formula is used:
**
** B(2*k) = (2*k - 1)/(2 * (2*k+1))
**          - binomial(2*k,1) * B(2)/2
**          - binomial(2*k,3) * B(4)/4
**          - ...
**          - binomial(2*k,2*k-3) * B(2*k-2)/(2*k-2).
*)
function nextbern(var Bern: array of array[2]; k: integer): array[2];
var
    xy, xy1: array[2];
    m, m1, m2, cc, bnum, bden: integer;
begin
    xy := (2*k - 1, 4*k + 2);
    cc := 2*k;
    xy1 := (k,6);
    xy := subrat(xy,xy1);
    for m := 2 to k-1 do
        m1 := 2*k - 2*(m-2) - 1;
        m2 := 2*m - 2;
        cc := (cc * m1) div m2;
        cc := (cc * (m1-1)) div (m2+1);
        (bnum, bden) := Bern[m];
        xy1[0] := bnum * cc;
        xy1[1] := bden * 2*m;
        xy := subrat(xy,xy1);
    end;
    return (Bern[k] := xy);
end;
(**********************************************************)
(*
** Addition of rational numbers
*)
function addrat(u,v: array[2]): array[2];
var
    x,y,d: integer;
begin
    x := u[0]*v[1] + u[1]*v[0];
    y := u[1]*v[1];
    d := gcd(x,y);
    return (x, y) div d;
end;
(*--------------------------------------------------------*)
(*
** Subtraction of rational numbers
*)
function subrat(u,v: array[2]): array[2];
var
    x,y,d: integer;
begin
    x := u[0]*v[1] - u[1]*v[0];
    y := u[1]*v[1];
    d := gcd(x,y);
    return (x, y) div d;
end;
(**********************************************************)
(*
** Arithmetic of complex numbers
*)
(*--------------------------------------------------------*)
type
    complex = array[2] of real;
end;
(*--------------------------------------------------------*)
(*
** multiplication of complex numbers
*)
function cmult(x,y: complex): complex;
begin
    return (x[0]*y[0]-x[1]*y[1], x[0]*y[1]+x[1]*y[0]);
end;
(*--------------------------------------------------------*)
(*
** division of complex numbers
*)
function cdiv(x,y: complex): complex;
var
    y0,y1,r2: real;
begin
    (y0,y1) := y;
    r2 := y0*y0 + y1*y1;
    return cmult(x,(y0,-y1)/r2);
end;
(*--------------------------------------------------------*)
(*
** inverse of a complex number
*)
function cinv(x: complex): complex;
var
    y0,y1,r2: real;
begin
    (y0,y1) := x;
    r2 := y0*y0 + y1*y1;
    return (y0, -y1)/r2;
end;
(*--------------------------------------------------------*)
(*
** absolute value of a complex number
*)
function cnorm(s: complex): real;
begin
    return sqrt(s[0]*s[0] + s[1]*s[1]);
end;
(*--------------------------------------------------------*)
(*
** Calculates a**s, where a is real > 0
*)
function cpow(a: real; s: complex): complex;
var
    lna, r, phi: real;
begin
    lna := log(a);
    r := exp(lna * s[0]);
    phi := lna * s[1];
    return (r * cos(phi), r * sin(phi));
end;
(*--------------------------------------------------------*)
(*
** Exponential function with complex argument
*)
function cexp(s: complex): complex;
var
    x,y: real;
begin
    (x,y) := s;
    return exp(x) * (cos(y), sin(y));
end;
(*--------------------------------------------------------*)
(*
** Function cosine with complex argument
*)
function ccos(s: complex): complex;
var
    x, y, t, t1: real;
begin
    (x,y) := s;
    t := exp(y);
    t1 := 1/t;
    return (cos(x) * (t + t1)/2, sin(x) * (t1 - t)/2);
end;
(*--------------------------------------------------------*)
(*
** Function sine with complex argument
*)
function csin(s: complex): complex;
var
    x, y, t, t1: real;
begin
    (x,y) := s;
    t := exp(y);
    t1 := 1/t;
    return (sin(x) * (t + t1)/2, cos(x) * (t - t1)/2);
end;
(*--------------------------------------------------------*)
(*
** Logarithm with complex argument, principal branch
*)
function clog(s: complex): complex;
var
    x,y: real;
begin
    (x,y) := s;
    return (log(x*x + y*y)/2, arctan2(y,x));
end;
(**********************************************************************)
(*--------------------------------------------------------------------*)
function tprod(s: complex; n: integer): complex;
(* Returns  s * (s+1) * (s+2) * ... * (s+n-1) *)
var
    k: integer;
    z: complex;
begin
    z := (1.0,0.0);
    for k := 0 to n-1 do
        z := cmult(z,s);
        s[0] := s[0] + 1;
    end;
    return z;
end;
(*--------------------------------------------------------------------*)
function lntprod(s: complex; n: integer): complex;
(*
    Calculates log(s * (s+1) * (s+2) * ... * (s+n-1))

    The complex plane is cut along the negative real axis
    For s real and positive, lntprod is real

    !! It is not true in general that lntprod(s,n) = clog(tprod(s,n)) !!
*)
var
    k: integer;
    z: complex;
begin
    z := (0.0,0.0);
    for k := 0 to n-1 do
        z := z + clog(s);
        s[0] := s[0] + 1;
    end;
    return z;
end;
(*--------------------------------------------------------------------*)
function stirling(s: complex): complex;
(*
    Calculates stirling(s) = log Gamma(s) using Stirling's expansion

    log Gamma(s) = (s - 1/2)*log(x) - s + (1/2)*log(2*pi)
                + (bern(2)/(1*2)) * 1/s
                + (bern(4)/(3*4)) * 1/s**3
                + (bern(6)/(5*6)) * 1/s**5
                ...
                + (bern(2*m-2)/((2*m-3)*(2*m-2)) * 1/s**(2*m-3)
                + R(2*m)

    For Re(s) > 0 the remainder term satisfies

    |R(2*m)| <= 1/cos(phi/2)**(2*m) *|bern(2*m)|/((2*m-1)*2*m) * 1/s**(2*m-1);
                phi = arg(s); |phi| <= pi/2.

    If Re(s) >= 1 and abs(s-1) >= 11,
    we can use m=6 for a sufficiently good approximation of log Gamma(s)
*)
external
    Ber: array of real;
var
    m := 6;
    i: integer;
    b: real;
    u, s1, s2, z: complex;
begin
    u := clog(s);                           (* u = log(s) *)
    u := cmult(u,(s[0] - 0.5, s[1])) - s;   (* u = (s - 1/2)*log(s) - s *)
    z := (u[0] + log(2*pi)/2, u[1]);
        (* z = (s + 1/2)*log(s) - s + (1/2)*log(2*pi) *)

    s1 := cinv(s);                          (* s1 = 1/s *)
    s2 := cmult(s1,s1);                     (* s2 = 1/s**2 *)
    b := Ber[m-1]/((2*m-2)*(2*m-3));
    u := (b,0);
    for i := m-2 to 1 by -1 do
        u := cmult(u,s2);
        b := Ber[i]/(2*i*(2*i-1));
        u[0] := u[0] + b;
    end;
    u := cmult(u,s1);
    (*
          u = (bern(2)/(1*2))/s + (bern(4)/(3*4))/s**3
              + (bern(6)/(5*6))/s**5 + ...
              + (bern(2*m-2)/((2*m-3)*(2*m-2))/s**(2*m-3)
    *)
    return (z + u);
end;
(*--------------------------------------------------------------------*)
(*
** Auxiliary function for Gamma and logGamma
*)
function gammaprep(s: complex): integer;
(*
   Auxiliary function for logGamma();
   returns an integer t >= 0, such that  ss := s + t  satisfies:
      Re(ss) >= 1 and |ss-1| >= R := 11
*)
var
    R := 11;
    a,b: real;
    m,t: integer;
begin
    (a,b) := s;
    m := floor(a) - 1;
    if m >= R then
        t := 0;
    elsif m >= 1 then
       	if abs(b) >= R then
            t := 0;
       	else
            t := R - m;
        end;
    elsif abs(b) >= R then
        t := 1 - m;
    else
       	t := R - m;
    end;
    return t;
end;
(*--------------------------------------------------------------------*)
function logGamma(s: complex): complex;
(*
    If Re(s) < 1 or |s-1| < 11, the argument s is translated
    by an integer t > 0 and we have

        logGamma(s) = logGamma(s+t) - log(s*(s+1)*...*(s+t-1))
*)
var
    t: integer;
    u: complex;
begin
    t := gammaprep(s);
    if t = 0 then
        return stirling(s);
    else
        u := lntprod(s,t);      (* u := log(s*(s+1)*...*(s+t-1)) *)
        s[0] := s[0] + t;       (* s := s + t *)
        return (stirling(s) - u);
    end;
end;
(*--------------------------------------------------------------------*)
function Gamma(s: complex): complex;
var
    u,z: complex;
    t: integer;
begin
    t := gammaprep(s);
    u := stirling((s[0]+t,s[1]));
    z := cexp(u);
    if t /= 0 then
        u := tprod(s,t);
        z := cdiv(z,u);
    end;
    return z;
end;
(*--------------------------------------------------------------------*)
(*
** Auxiliary function for emczeta(bb,n,s)
** Coefficients for Euler-McLaurin summation formula
** calculates emc[k] for k = 1,2,...,bb
**
**      emc[k] := bern(2k)/(2k)! * s * (s+1) * (s+2) * ... * (s+2k-2);
*)
function emczetacoef(s: complex; bb: integer): array of complex;
external
    BMax: integer;
    Ber: array of real;
var
    k: integer;
    u,v: complex;
    emc: array of complex;
begin
    if bb <= 0 then 
        bb := 1; 
    elsif bb > BMax then
        bb := BMax;
    end;
    emc := alloc(array,bb+1,(0.0,0.0));
    u := s/2;
    emc[1] := u * Ber[1];
    for k := 2 to bb do
        s[0] := s[0] + 1;
        u := cmult(u,s);
        s[0] := s[0] + 1;
        u := cmult(u,s);
        u := u/((2*k-1)*2*k);
        emc[k] := u*Ber[k];
    end;
    return emc;
end;
(*--------------------------------------------------------------------*)
function emczeta(deg, n: integer; s: complex): complex;
(*
    Calculates
    z := emczeta(deg,n,s) =
         1/(s-1)*1/n**(s-1) + 1/2*1/n**s
         + emc[1]/n**(s+1) + emc[2]/n**(s+3) + ...
         + emc[deg]/n**(s+2*deg-1);
*)
external
    BMax: integer;
var
    i, n2: integer;
    s1,nms,z,z1: complex;
    emc: array of array[2];
begin
    if deg > BMax then
        writeln("emczeta(deg,n,s), deg = ",deg," too big");
        deg := 0;
    end;
    emc := emczetacoef(s,deg); 
    s1 := cinv(s-(1,0));
    z := n*s1 + (0.5,0);    (* z = n/(s-1) + 0.5 *)

    if deg <= 0 then
        z1 := 0.0;
    else
        n2 := n*n;
        z1 := emc[deg];
        for i := deg-1 to 1 by -1 do
            z1 := z1/n2 + emc[i];
        end;
        z1 := z1/n;
        (* z1 = emc[1]/n + emc[2]/n**3 + ... + emc[deg]/n**(2*deg - 1) *)
    end;
    z := z + z1;
    nms := cpow(n,-s);                 (* nms = 1/n**s *)
    return cmult(z,nms);
end;
(*--------------------------------------------------------------------*)
(*
** Partial sum of zeta series
**
**      harms(n,s) = sum( 1/k**s: 1 <= k <= n )
*)
function harms(n: integer; s: complex): complex;
var
    k: integer;
    z: complex;
begin
    z := (0,0);
    for k := 1 to n do
        z := z + cexp(-s*log(k));          (* z := z + 1/k**s *)
    end;
    return z;
end;
(*--------------------------------------------------------------------*)
function zetafuneq(s: complex): complex;
(*
    The functional equation of the zeta function is

        zeta(s) = Gamma(1-s) * (2*pi)**(s-1) * 2 * sin(pi*s/2) * zeta(1-s)

     zetafuneq(s)  calculates the factor of zeta(1-s)
     Hypothesis: Re(s) < 1
*)
var
    u,v,w: complex;
begin
    u := Gamma((1-s[0],-s[1]));             (* u = Gamma(1-s) *)
    v := cpow(2*pi,(s[0]-1, s[1]));         (* v = (2*pi)**(s-1) *)
    v := cmult(u,v);

    w := 2*csin(pi*s/2);                    (* w = 2*sin(pi*s/2) *)

    return cmult(v,w);
end;
(*--------------------------------------------------------------------*)
(*
** Calculates the zeta function of a complex argument s
*)
function zeta(s: complex; verbose := 0): complex;
var
    eps := 1.0E-10;
    bb: integer;                (* bb <= BMax *)
    u,v,ff,s1,z: complex;
    N: integer;
    funeq: boolean;
begin
    if s[0] < 0 then
        funeq := true;
        ff := zetafuneq(s);
        s := (1-s[0],-s[1]);    (* s = 1 - s *)
    else
        funeq := false;
    end;
    (bb,N) := getparams(s,eps);
    if verbose then
        writeln("calculating zeta with N = ",N," and bb = ",bb);
    end;
    u := emczeta(bb,N,s);
    v := harms(N-1,s);
    if verbose then
        writeln("emczeta(bb,N,s) = ",u);
        writeln("harms(N-1,s) = ",v);
    end;
    z := u + v;
    if funeq then
        z := cmult(z,ff);
    end;
    return z;
end;
(*--------------------------------------------------------------------*)
(*
** Error estimate for the calculation of zeta(s)
** using Euler/Maclaurin summation formula,
** if partial sum(1/n**s, 1 <= n <= N) and
** Bernoulli numbers up to B(2*r) are used.
*)
function zetaerr(s: complex; N,r: integer): real;
var
    sigma,sr,t: real;
    u,u1,u2: real;
    v1: complex;
begin
    (sigma,t) := s;
    u := 3/(2*pi)**(2*r);    (* u ~= B(2*r)/(2*r)! *)
    v1 := tprod(s,2*r);
    u1 := cnorm(v1);
    sr := sigma + 2*r - 1;
    u2 := sr * N**sr;
    return u*u1/u2;
end;
(*------------------------------------------------------------*)
function getparams(s: complex; eps: real): array[2];
external
    BMax: integer;
var
    sigma,t: real;
    r,N: integer;
    BBr,ssr,sigr,u,x: real;
begin
    (sigma,t) := s;
    t := abs(t);
    if sigma >= 10 then
        r := 1 + floor(20/sigma);
    elsif t <= 64 then
        r := 4;
    else
        r := round(t**(1/3));       (* somewhat arbitrary *)
        if r > BMax then r := BMax; end;    
    end;
    BBr := 3/(2*pi)**(2*r);    (* BBr ~= B(2*r)/(2*r)! *)
    ssr := cnorm(tprod(s,2*r));
    sigr := sigma + 2*r - 1;
    u := BBr*(ssr/sigr);
    x := (u/eps)**(1/sigr);
    N := 1 + floor(x);
    return (r,N);
end;
(*-------------------------------------------------------------zeta(*)
function Zet(t: real; verbose := 0): real;
(*
	Zet(t) := zeta(1/2 + i*t) * exp(i*zetheta(t))

  The function Zet is real valued, its zeros t correspond 
  to the zeros 1/2 + i*t of the zeta function
*)
var
    eps := 1.0E-10;
    bb: integer;                (* bb <= Bmax - 1 *)
    k,N: integer;
    s, u: complex;
    theta, x, y: real;
begin
    s := (0.5,t);
    (bb,N) := getparams(s,eps);
    if verbose then
        writeln("calculating zeta with N = ",N," and bb = ",bb);
    end;
    u := emczeta(bb,N,s);

    theta := zetheta(t);
    y := u[0]*cos(theta) - u[1]*sin(theta);

    x := 0.0;
    for k := 1 to N-1 do
        x := x + cos(theta - t*log(k))/sqrt(k);
    end;
    (* x = real part of sum(exp(i*theta)/n**(0.5+i*t): 1 <= k <= N-1) *)
    return(x + y);
end;
(*--------------------------------------------------------------------*)
function zetheta(t: real): real;
(*
    zetheta(t) = Im logGamma(1/4 + i*t/2) - (t/2)*log(pi)

  The function
       Zet(t) := zeta(1/2 + i*t) * exp(i*zetheta(t))
  is real valued

  For t -> +infinity one has

    zetheta(t) = (t/2)*log(t/(2*pi)) - t/2
                 - pi/8 + 1/(48*t) + 7/(5760*t*t*t) + ...
*)
var
    tt, theta: real;
    s,z: complex;
begin
    tt := abs(t);
    if tt <= 12 then
        s := (0.25,t/2);
        z := logGamma(s);
        theta := z[1] - (t/2)*log(pi);
    else
        theta := (t/2)*(log(tt/(2*pi)) - 1) - pi/8 + 1/(48*t);
        if tt < 160 then
            theta := theta + 7/(5760*t*t*t);
        end;
    end;
    return theta;
end;
(*--------------------------------------------------------------------*)
(*
** Calculation of a zeta zero by interval halving
** Function Zet must have opposite signs at x1 resp. x2
** otherwise return value is 0.0
** TODO: Regula falsi
*)
function zetazero(x1,x2: real; eps := 10**-8): real;
var
    z1,z2,z,x3,delta: real;
    mess := "Zet() has same sign at interval ends";
begin
    z1 := Zet(x1); z2 := Zet(x2);
    if (z1 > 0 and z2 > 0) or (z1 < 0 and z2 < 0) then
        writeln(mess);
        return(0.0);
    elsif z1 = 0.0 then
        writeln("probable zero at x1");
        return x1;
    elsif z2 = 0.0 then
        writeln("probable zero at x2");
        return x2;
    end;
    (* now z1 and z2 have opposite sign *)
    if z1 > 0 then
         (x1,x2) := (x2,x1);
    end;
    delta := x2 - x1;
    while abs(delta) > eps do
        delta := delta/2;
        x3 := x1 + delta;
        z := Zet(x3);
        if z <= 0 then
            x1 := x3;
        else
            x2 := x3;
        end;
    end;
    return x1 + delta/2;
end;
(*--------------------------------------------------------------*)
(*
** Approximate number of zeros s of the Riemann zeta function
** with 0 < Im(s) < T
*) 
function NTapprox(T: real): real;
var
    t: real;
begin
    t := T/(2*pi);
    return t*(log(t) - 1) + 7/8;
end;
(*--------------------------------------------------------------------*)
function gram(n: integer): real;
(* Calculates the n-th Gram point, which is defined by

        zetheta(gram(n)) = n*pi

    Since gram(0) > 17, the following approximation of zetheta suffices

        zetheta(t) = (t/2)*log(t/(2*pi)) - t/2 - pi/8
                      + 1/(48*t) + 7/(5760*t*t*t) + ...
*)
var
    eps := 1.0E-10;
    t0, t1: real;
    c, ln2pi: real;
begin
    if n < 0 then return 0.0; end;
    c := (2.0*n + 0.25)*pi;
    ln2pi := log(2*pi);
    t0 := 0; t1 := 4.0*n + 20.0;
    while abs(t1 - t0) > eps do
        t0 := t1;
        t1 := c + t0 - 1/(24*t0) - 7/(2880*t0*t0*t0);
        t1 := t1/(log(t0) - ln2pi);
    end;
    return t1;
end;
(*---------------------------------------------------------------*)
(*******************************************************************)
(*
** INITIALIZATION
*)
(* store Bernoulli constants with high precision *)
set_floatprec(128);

bern_ini(100);

(* default precision for calculation with zeta *)
set_floatprec(64);

(*******************************************************************)