C program to look for large values of | zeta( 1/2 + it) | w.r.t. t

I've been doing some more work on a C program to find some large
values of Z(t) = either: | zeta(1/2 + it)| or -|zeta(1/2 + it| ;
Z(t) is sometime called the Hardy function. In my program,
only the most significant or main term of the Riemann-Siegel
formula is used; the other terms aren't included because
they are a but complicated, count mainly when one wants to
compute zeros of zeta on the critical line, and don't affect
conjectures on the asymptotic behaviour of
G(y) = max_{y in [0, t]} | Z(y) |.

Tadej Kotnik had a paper published in Math. Comp. ,
[ 73 (2004), 949-956 ] , where precise values
of increasing values of | Z(t) | were presented in a
table. He included positive and negative values of Z(t),
and his table goes up to t = 10^13.

Cf.:
http://www.ams.org/mcom/2004-73-246/S0025-5718-03-01568-0/home.html

I think a preprint is available for free from his homepage
in the "si" country domain (in postscript). He also
did computations on the growth of zeta(1 + it), with a
collaborator I believe. I'm not sure how much the range
beyond t=10^13 has been looked at for large
values of Z(t).

Based on random matrix theory and other methods,
David W. Farmer, S.M. Gonek and C.P. Hughes have conjectured that

max_{t in [0, T]} |Z(t)| = exp( (1 + o(1)) sqrt(1/2 * T log(log(T)) ) )

This is Conjecture A from the preprint
that can be found here:
http://arxiv.org/abs/math/0506218

I admit that I'm not sure what to think of what appears
on page 23 of a PDF file from a workshop with date 2004
available here:
www.aimath.org/WWN/lrmt/lrmt.pdf

(The AIM is a well-reputed organization).

<< E.1.b The maximal order of the zeta-function on the critical line.
How large is the maximum value of |zeta(1/2 + it)| for T < t < 2T? It is known that the Riemann Hypothesis implies that the maximum is at most exp(clog T/log log T) for some c > 0.

It is also known that the maximum gets as big as
exp(c_1(log T/log log T)^(1/2) ) for a sequence of T -> oo
for some c_1 > 0. It has been conjectured that the
smaller bound (the one that is known to occur) is closer
to the truth. However, the new conjectures
about moments [74] suggest that it may be the larger. >>

The above seems to say that, in 2004 at the AIM,
moments conjectures suggested that

max_{t in [0, T]} |Z(t)|
might be as large as exp(clog T/log log T) for some c > 0.
I think that would be larger than what we get from Conjecture
A in the Farmer, Gonek and Hughes preprint.

----

In any case, my program in its present form uses heuristics found in
Kotnik's paper, i.e. (unless I'm mistaken)

- vtheta(t) close to zero in Riemann-Siegel,
- t/(2pi/log(2)) close to 0 (mod 1) , either as 0.0129 ...
or 0.9928....
- t/(2pi/(log(3)) also close to 0 "" """ ""
- idem to a lesser degree for primes 5, 7.

There is no attempt to find all increasing local maxima of Z(t).
Some factors of merit (quality figures) are computed, and
some maxima will no doubt be missed.

- A feature I haven't seen elsewhere is a Monte Carlo estimate of
the sum of the terms passed #20,000 in the RS formula. If
the Monte Carlo estimate, added to the "exact" sum of
terms 1 to 20,000 is large enough, then an "exact"
sum from first to last term in R.S. is done.

- Double precision with 64-bits is used. That way, I can't
use the program for 10^20 . I've tried converting to
a library with arbitrary precision arithmetic, but
the expressions for functions are much longer.
I've managed to reproduce my 64-bit results, but
everything is slow, and I have bugs in some further
versions. The library is MPFR. No Monte Carlo is
used here.

I'm now passed 10^13 in the 64-bit version based only on
math.h , but this was from 1.3 x 10^13 .

t Z(t) MC Exact
14762037972713.855469 447.85 14.56 14.77
............
14772754195582.750000 404.66 13.27 13.34

MC is the Monte Carlo approximation of the RS sum,
divided by log(t), and is used to decide whether to
do a full-term "exact" R.S. sum . "Exact" refers
to the full term R.S. sum, divided by log(t).

The source code in full (including Mersenne Twister)
appears below, in case someone is interested.

David Bernier


P.S.: some lines may be wrapped below. I can
email it as an attached file if needed.
========================================================

#include <stdio.h>
#include <math.h>



#define Nmersenne 624
#define Mmersenne 397
#define MATRIX_A 0x9908b0dfUL
#define UPPER_MASK 0x80000000UL
#define LOWER_MASK 0x7fffffffUL

#define G8 genrand_int8()

static unsigned long mt[Nmersenne];
static int mti=Nmersenne+1;


void init_genrand(unsigned long s)
{
mt[0]= s & 0xffffffffUL;
for (mti=1; mti<Nmersenne; mti++) {
mt[mti] =
(1812433253UL * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti);
mt[mti] &= 0xffffffffUL;

}
}



void init_by_array(unsigned long[], int);
void init_by_array(unsigned long init_key[], int key_length)
{
int i, j, k;
init_genrand(19650218UL);
i=1; j=0;
k = (Nmersenne>key_length ? Nmersenne : key_length);
for (; k; k--) {
mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525UL))
+ init_key[j] + j;
mt[i] &= 0xffffffffUL;
i++; j++;
if (i>=Nmersenne) { mt[0] = mt[Nmersenne-1]; i=1; }
if (j>=key_length) j=0;
}
for (k=Nmersenne-1; k; k--) {
mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941UL)) - i;
mt[i] &= 0xffffffffUL;
i++;
if (i>=Nmersenne) { mt[0] = mt[Nmersenne-1]; i=1; }
}

mt[0] = 0x80000000UL;
}


unsigned long genrand_int8(void)
{
unsigned long y;
static unsigned long mag01[2]={0x0UL, MATRIX_A};


if (mti >= Nmersenne) {
int kk;

if (mti == Nmersenne+1)
init_genrand(5489UL);

for (kk=0;kk<Nmersenne-Mmersenne;kk++) {
y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
mt[kk] = mt[kk+Mmersenne] ^ (y >> 1) ^ mag01[y & 0x1UL];
}
for (;kk<Nmersenne-1;kk++) {
y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
mt[kk] = mt[kk+(Mmersenne-Nmersenne)] ^ (y >> 1) ^ mag01[y & 0x1UL];
}
y = (mt[Nmersenne-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
mt[Nmersenne-1] = mt[Mmersenne-1] ^ (y >> 1) ^ mag01[y & 0x1UL];

mti = 0;
}

y = mt[mti++];


y ^= (y >> 11);
y ^= (y << 7) & 0x9d2c5680UL;
y ^= (y << 15) & 0xefc60000UL;
y ^= (y >> 18);

return y;
}



static double small[10200];
static double smallprlog[8];
static double sroot[10200];



int main(void)
{
double pi=3.14159265358979323846;
double dt;
double t, t0;
double vtheta;
double r7;
unsigned long multiplier;
unsigned long Nbill;
double tau;
double maxz;
double Z;
double Z20;
double Z100;
double Z300;
double Z1000;
double Z3000;
double Z10000;
double Z20000;

double maxt;
double d1, d2, d8, d1000, d3000, d10000;
double dBill;
double dhalf;
double C;
double m, msafe;
double tmptmax;
double sigma;
double mlog2, mlog3, mlog5, sum;
double m2, m3, m5;

unsigned long initmersenne[4]={0x3251, 0x6504, 0xa948, 0xbedc};
int lengthmersenne=4;
unsigned long random;
double div, unif;

double t0sin, t0cos;
double sv;
double cv;
double sigrnd;



long i;
long counter=0;
long indexer;
int excess;
int numcounter;
int newmax;
int dofullsum;
int prompt;
C = d2*pi/log((double)2);



small[0] = (double)0;
sroot[0] = (double)0;
for(i=1;i<10200;i++)
{
small[i] = log( (double) i);
sroot[i] = sqrt( (double) i);
}
smallprlog[0] = small[2];
smallprlog[1] = small[3];
smallprlog[2] = small[5];
smallprlog[3] = small[7];
smallprlog[4] = small[11];
smallprlog[5] = small[13];
smallprlog[6] = small[17];
smallprlog[7] = small[19];

maxz = (double) 0;
d1 = (double) 1;
d2 = (double) 2;
d8 = (double) 8;
d1000 = (double) 1000;
d3000 = (double) 3000;
d10000 = (double) 10000;
dhalf = ((double)1)/d2;
dBill = ((double) 1000000000) ;
newmax = 0;
C = d2*pi/small[2];

multiplier = 92767630;

Nbill = 524;

Nbill = Nbill + 100;

Nbill = Nbill + 137;
Nbill = Nbill + 100;


Nbill = Nbill + 600;


init_by_array(initmersenne, lengthmersenne);
div = (double)4294967295UL;




while(1 == 1)
{


multiplier++;

if((multiplier%100000000)==0)
{
printf(".");
fflush(stdout);
}


if(multiplier == 1000000000)
{
multiplier=0;
Nbill++;

}



t0 = ((double)multiplier);
t0 = ((double)Nbill)*dBill + t0;
t0 = t0*C;
t = t0;




vtheta = (t/d2)*log(t/(d2*pi)) - t/d2 - pi/d8;

dt=vtheta/(d2*pi);
dt = dt-floor(dt);

if(dt>dhalf)
{
t = t + d2*pi*(d1-dt)/( dhalf*log(t/(d2*pi)) );
}
else
{
t = t - d2*pi*dt/( dhalf*log(t/(d2*pi)) );
}
vtheta = (t/d2)*log(t/(d2*pi)) - t/d2 - pi/d8;


dt=vtheta/(d2*pi);
dt = dt-floor(dt);

if(dt>dhalf)
{
t = t + d2*pi*(d1-dt)/( dhalf*log(t/(d2*pi)) );
}
else
{
t = t - d2*pi*dt/( dhalf*log(t/(d2*pi)) );
}

vtheta = (t/d2)*log(t/(d2*pi)) - t/d2 - pi/d8;


/*** test for diophantine log2, log3 ***/

excess = 0;

for(i=0;i<3;i++)
{
mlog2 = t*smallprlog[i]/(d2*pi);
mlog2 = mlog2 - floor(mlog2);
if(mlog2>0.5)
{
mlog2 = ((double)1) - mlog2;
}

if(mlog2 > 0.12)
{
excess = 1;
}
if( i == 0)
{
m2 = mlog2;
if(m2 > 0.025)
{
excess = 1;
}

}
if( i == 1)
{
m3 = mlog2;
if(m3 > 0.035)
{
excess = 1;
}

}

if( i == 2)
{
m5 = mlog2;


}



}


if( excess == 0)
{


vtheta = (t/d2)*log(t/(d2*pi)) - t/d2 - pi/d8;
cv = cos(vtheta);
sv = sin(vtheta);
if(cv>0)
{
if( fabs(sv)<0.25 )
{


Z = (double) 0;
i = (long) 1;
while( ((double)i) < 20.5 )
{
Z = Z + cos(t*small[i] - vtheta)/sroot[i];
i++;
}
Z20 = Z + Z;


if(Z20>14.0 )
{// do test ....

Z = (double) 0;
i = (long) 1;
while( ((double)i) < 100.5 )
{
Z = Z + cos(t*small[i] - vtheta)/sroot[i];
i++;
}
Z100 = Z + Z;






if(Z100> ((double)33) )
{



Z = (double) 0;

i = (long) 1;
while( ((double)i) < 300.5 )
{
Z = Z + cos(t*small[i] - vtheta)/sroot[i];
i++;
}
Z300 = Z + Z;




if(Z300>52.0 )
{
Z = (double) 0;

i = (long) 1;
while( ((double)i) < 1000.5 )
{
Z = Z + cos(t*small[i] - vtheta)/sroot[i];
i++;
}
Z1000 = Z + Z;


if(Z1000>80.0 )
{

Z = (double) 0;
i = (long) 1;

while( ((double)i) < 3000.5 )
{
Z = Z + cos(t*small[i] - vtheta)/sroot[i];
i++;
}

Z3000 = Z + Z;

if(Z3000 > 113.0)
{
Z = (double) 0;
i = (long) 1;

while( ((double)i) < 10000.5 )
{
Z = Z + cos(t*small[i] - vtheta)/sroot[i];
i++;
}

Z10000 = Z + Z;

if(Z10000>154.0)
{
Z = (double) 0;
i = (long) 1;

while( ((double)i) < 20000.5 )
{
Z = Z + cos(t*log( (double) i) - vtheta)/sqrt( (double) i);
i++;
}

Z20000 = Z + Z;

if(Z20000 > 183.0)
{
Z = (double) 0;
i = (long) 1;


tau = sqrt(t/(d2*pi));
m= floor(tau);
Z = (double) 0;


for(indexer=0;indexer<10000;indexer++)
{
random = (G8);
unif = ((double)random)/div;
i = floor(unif*(m-d10000*d2)) + d1 + d2*d10000;
Z = Z + cos(t*log( i) - vtheta)/sqrt(i);

}

Z = (m-d2*d10000)/((double)10000)*Z;

Z = Z + Z+ Z20000;
sigma = Z/log(t);
if(sigma>2.8) // 11.8
{
sigrnd = sigma;

tau = sqrt(t/(d2*pi));
m= floor(tau);



msafe = m + dhalf;
Z = (double) 0;

i = (long) 1;



while( ((double)i) < msafe )
{
Z = Z + cos(t*log( (double) i) - vtheta)/sqrt( (double) i);
i++;
}

Z = Z + Z;
sigma = Z/log(t);

if(sigma>2.0)
{
printf("\n%16lf %.2lf %.2lf %.2lf %.4lf %.4lf\n",t, Z, sigrnd, sigma, fabs(sv), m5);
}


} // if sigma

} // if Z20000

} // if Z10000

} // if Z3000

} // if Z1000

} // if Z300

} // if Z100

} // if Z20

} // if sv

} // if cv

} // if excess


} // while ..



printf("\n\n");
return 0;
}
=====================================================


.

Relevant Pages

  • Re: Performance of hand-optimised assembly
    ... Why not use unsigned for sum in the test harness? ... int is better yet because the output will then be same across machines. ... gcc -O3 seems ... but taking an average by eye my best asm version was *slower* than ...
    (comp.lang.c)
  • Re: [PATCH] input: Introduce light-weight contact tracking
    ... seems the patch was damaged somehow. ... static int illegal(int nslot, int npos, unsigned x) ... int i, j, sum; ...
    (Linux-Kernel)
  • Re: Algorithm
    ... Wont sum of all positive numbers will be the largest sub-array? ... int getint ... struct sofar *next; ... struct sofar *discard(struct sofar *trail) ...
    (comp.lang.c)
  • Re: Arrays vs Buffers
    ... and the next array access is dependent on the value ... and then return the overall sum to the test harness ... nextPointer(int[] array, int point) ... pointer = nextPointer; ...
    (comp.lang.java.programmer)
  • Re: sum(2,4,3,5,-1,...)
    ... The only example I have to imitate is the one on p.289 ... of an array related to the first argument of the function sum I'm ... int sum ... no args, and get 0, like you can do in Scheme or C++. ...
    (comp.lang.c.moderated)