Help need for evaluation of latest wiki page equations

[rudolph-git-acc]

Grateful for any help to properly evaluate the sequence of equations on this Wiki page: http://michaelnielsen.org/polymath1/index.php?title=Polymath15_test_problem#Large_negative_values_of_.5Bmath.5Dt.5B.2Fmath.5D

Below is the pari/gp code I have used so far. At the end of the code is a simple rootfinder that should (approximately) find this sequence of roots:

-10.00000000, 906.72649548
-10.00000000, 931.47009390
-10.00000000, 976.80075639
-10.00000000, 990.45190417

however from eq 318 onwards it finds many more roots (also at different t, e.g. -100). Have experimented with the integral limits, the integral's fixed parameter, the real precision settings of pari and the size of the sums, but no success yet.

Have re-injected some of the non-zero multiplication factors that can be switched on/off if needed.

One thing I noticed is that after eq 307 the function doesn't seem to be real anymore (see data at the end). Not sure if this should matter though, since I use the real value of Ht for root finding.

default(realprecision, 40)

Xi(s)=return((s/2)*(s-1)*Pi^(-s/2)*gamma(s/2)*zeta(s));
M0(s)=return(s*(s-1)/2*Pi^(-s/2)*sqrt(2*Pi)*exp((s/2-1/2)*log(s/2)-s/2));

Ht302(x,t) = {
    A=intnum(v=-8,8,Xi((1+I*x)/2+I*sqrt(abs(t))*v)*exp(-v^2),1);
    return(A);
}

Ht303(x,t) = {
    A=intnum(v=-8,8,Xi((1+I*x)/2+I*sqrt(abs(t))*v-Pi*I*abs(t)/8)*exp(-(v-Pi*sqrt(abs(t))/8)^2),1);
    return(A);
}

Ht304(x,t) = {
    xt = x-Pi*abs(t)/4;
    mulfact = exp(-Pi^2*abs(t)/64);
    A=mulfact*intnum(v=-8,8,Xi((1+I*xt)/2+I*sqrt(abs(t))*v)*exp(-v^2+Pi*sqrt(abs(t))*v/4),1);
    return(A);
}

Ht307(x,t) = {
    xt = x-Pi*abs(t)/4;
    mulfact = exp(-Pi^2*abs(t)/64)*2;
    A=mulfact*intnum(v=-8,8,real(M0((1+I*xt)/2+I*sqrt(abs(t))*v))*zeta((1+I*xt)/2+I*sqrt(abs(t))*v)*exp(-v^2+Pi*sqrt(abs(t))*v/4),1);
    return(A);
}

Ht312(x,t) = {
    xt = x-Pi*abs(t)/4;
    mulfact = exp(-Pi^2*abs(t)/64)*(M0((1+I*xt)/2));
    A= mulfact*intnum(v=-8,8,exp((I*sqrt(abs(t))*v)/2*log(xt/(4*Pi))-Pi*sqrt(abs(t))*v/4+(I*abs(t)*v^2)/(2*xt))*zeta((1+I*xt)/2+I*sqrt(abs(t))*v)*exp(-v^2+Pi*sqrt(abs(t))*v/4),1);
    return(A);
}

Ht315(x,t) = {
    xt = x-Pi*abs(t)/4;
    N = sqrt(xt/(4*Pi));
    u = 4*Pi*abs(t)/xt;
    mulfact = exp(-Pi^2*abs(t)/64)*(M0((1+I*xt)/2));
    A=mulfact*intnum(v=-8,8,exp(I*sqrt(abs(t))*v*log(N)+I*u*v^2/(8*Pi))*zeta((1+I*xt)/2+I*sqrt(abs(t))*v)*exp(-v^2),1);
    return(A);
}

Ht316(x,t) = {
    xt = x-Pi*abs(t)/4;
    N = sqrt(xt/(4*Pi));
    u = 4*Pi*abs(t)/xt;
    mulfact = exp(-Pi^2*abs(t)/64)*(M0((1+I*xt)/2));
    A=mulfact*sum(n=1,1000, intnum(v=-8,8,exp(I*sqrt(abs(t))*v*log(N)+I*u*v^2/(8*Pi))*n^(-((1+I*xt)/2+I*sqrt(abs(t))*v))*exp(-v^2),1));
    return(A);
}

Ht317(x,t) = {
    xt = x-Pi*abs(t)/4;
    N = sqrt(xt/(4*Pi));
    u = 4*Pi*abs(t)/xt;
    mulfact = exp(-Pi^2*abs(t)/64)*(M0((1+I*xt)/2))*N^((1+I*xt)/2)/8;
    A=mulfact*sum(n=1,1000, intnum(v=-8,8,exp(-I*sqrt(abs(t))*v*log(n/N)+I*u*v^2/(8*Pi)-((1+I*xt)/2)*log(n/N))*exp(-v^2),1));
    return(A);
}

Ht318(x,t) = {
    xt = x-Pi*abs(t)/4;
    N = sqrt(xt/(4*Pi));
    u = 4*Pi*abs(t)/xt;
    mulfact = exp(-Pi^2*abs(t)/64)*(M0((1+I*xt)/2))*N^((1+I*xt)/2)/8*sqrt(Pi);
    A=mulfact*sum(n=1,1000,exp(-abs(t)*(log(n/N))^2/(4*(1-I*u/(8*Pi)))-(1+I*xt)/2*log(n/N)));
    return(A);
}

Ht320(x,t) = {
    xt = x-Pi*abs(t)/4;
    N = sqrt(xt/(4*Pi));
    u = 4*Pi*abs(t)/xt;
    mulfact = exp(-Pi^2*abs(t)/64)*(M0((1+I*xt)/2))*N^((1+I*xt)/2)/8*sqrt(Pi);
    A = mulfact*sum(n=1,1000,exp(-(abs(t)*(n-N)^2)/(4*N^2*(1-I*u/(8*Pi)))-I*xt/2*(n-N)/N+I*xt*(n-N)^2/(4*N^2)));
    return(A);
}

Ht321(x,t) = {
    xt = x-Pi*abs(t)/4;
    N = sqrt(xt/(4*Pi));
    u = 4*Pi*abs(t)/xt;
    mulfact = exp(-Pi^2*abs(t)/64)*(M0((1+I*xt)/2))*N^((1+I*xt)/2)/8*sqrt(Pi);
    A = mulfact*sum(n=1,1000,exp(-(2*Pi*u*(n-N)^2)/(8*Pi-I*u)-2*Pi*I*N*(n-N)+Pi*I*(n-N)^2));
    return(A);
}

Ht323(x,t) = {
    xt = x-Pi*abs(t)/4;
    N = sqrt(xt/(4*Pi));
    u = 4*Pi*abs(t)/xt;
    mulfact = exp(-Pi^2*abs(t)/64)*(M0((1+I*xt)/2))*N^((1+I*xt)/2)/8*sqrt(Pi);
    A = mulfact*sum(n=1,1000,exp(-(2*Pi*u*(n-N)^2)/(8*Pi-I*u)-Pi*I*n^2+2*Pi*I*(n-N)^2));
    return(A);
}

Ht324(x,t) = {
    xt = x-Pi*abs(t)/4;
    N = sqrt(xt/(4*Pi));
    u = 4*Pi*abs(t)/xt;
    mulfact = exp(-Pi^2*abs(t)/64)*(M0((1+I*xt)/2))*N^((1+I*xt)/2)/8*sqrt(Pi);
    A = mulfact*sum(n=1,1000,exp((16*Pi^2*I*(n-N)^2)/(8*Pi-I*u))*exp(Pi*I*n));
    return(A);
}

Ht328(x,t) = {
    xt = x-Pi*abs(t)/4;
    N = sqrt(xt/(4*Pi));
    u = 4*Pi*abs(t)/xt;
    mulfact = exp(-Pi^2*abs(t)/64)*(M0((1+I*xt)/2))*N^((1+I*xt)/2)/8*sqrt(Pi);
    A = mulfact*sum(n=1,1000,exp(-Pi*I*n*(n+1)/2)*exp(2*Pi*I*(n+1/2)*N)*exp(-u*n*(n+1)/16));
    return(A);
}

print("eq.302 ",Ht302(1000, -10))
print("eq.303 ",Ht303(1000, -10))
print("eq.304 ",Ht304(1000, -10))
print("eq.307 ",Ht307(1000, -10))
print("eq.315 ",Ht315(1000, -10))
print("eq.316 ",Ht316(1000, -10))
print("eq.317 ",Ht317(1000, -10))
print("eq.318 ",Ht318(1000, -10))
print("eq.320 ",Ht320(1000, -10))
print("eq.321 ",Ht321(1000, -10))
print("eq.323 ",Ht323(1000, -10))
print("eq.324 ",Ht324(1000, -10))
print("eq.328 ",Ht324(1000, -10))

xs = 900; xe = 1010;
ts = -10; te = -1;
f=Ht318;
forstep(t=ts, te, 10, {
    xc = xs;
    while(xc <= xe,
        if(real(f(xc,t))*real(f(xc+1,t)) < 0, root=solve(a=xc, xc+1, real(f(a,t))); printf("%3.8f, %3.8f\n",t, root));
        xc=xc+1);
})

eq.302 1.201485926838305713580758240853488704456 E-165 + 3.514303326279329051 E-221*I
eq.303 1.201485926838305713582303927841938470920 E-165 + 3.522103826014908131 E-221*I
eq.304 1.201485926838305713582303927841938470920 E-165 + 3.522103826014922598 E-221*I
eq.307 1.201477416607903787215861199656638682224 E-165 + 3.992587328505261451855248837227037420869 E-169*I
eq.312 1.200242494404236645648308584282451391036 E-165 + 4.014384445003979486997306873566995218158 E-169*I
eq.315 1.200242494404236645648308584282451391036 E-165 + 4.014384445003979486997306873566995218158 E-169*I
eq.316 1.200242494404236645648308630967087228353 E-165 + 4.014384445003979487000178658198142256762 E-169*I
eq.317 1.210630457223022714918883693478442584040 E-165 - 5.581422172204223838770738013241198104262 E-166*I
eq.318 1.209227903617935860382758702133562335051 E-165 - 5.611945144292097211229689587084475461023 E-166*I
eq.320 2.791674127351571571455964690537038072397 E-167 - 2.126109533152847673586778474944806978351 E-167*I
eq.321 2.791674127351571571455964690537038072397 E-167 - 2.126109533152847673586778474944806978351 E-167*I
eq.323 -3.076893561828845392623560432368639943370 E-167 + 1.687160985046708488717820653125387727510 E-167*I
eq.324 -3.076893561828845392623560432368639943370 E-167 + 1.687160985046708488717820653125387727510 E-167*I
eq.328 -3.076893561828845392623560432368639943370 E-167 + 1.687160985046708488717820653125387727510 E-167*I

[teorth]

Thanks for this - it's easy for me to step through it myself now. The multiplicative factor for 317 should have N^(-(1+I*xt)/2) rather than N^((1+I*xt)/2)/8 (I guess the 1/8 factor was just dropped at the very start, as it wasn't important), which seems to fix that step at least.

Similarly, the multiplicative factor for 318 should be N^(-(1+I*xt)/2)*sqrt(Pi/(1-I*u/(8*Pi))) rather than N^((1+I*xt)/2)/8*sqrt(Pi) .

Now checking 3.20...

[teorth]

OK, it looks like even after fixing the multiplicative factor, the 3.20 approximation is not very accurate unless x and t are extremely large (of size 10^5 or so). However they do become reasonably good at that scale. What this could mean is that the wiki asymptotics are for a regime which is not the one seen in the numerics.

On the other hand, now that 3.18 remains reasonably accurate and significantly faster to calculate than the preceding expressions due to the lack of integration, perhaps it could be used to explore the zeroes at a larger range than previously done? e.g. for x,t of size around 10^5 or 10^6. I don't know right now whether the patterns we are seeing are transient; it could be that a different dynamics takes over at larger scales.

[rudolph-git-acc]

Many thanks for your help! I now indeed get the correct zeros for eq 3.18, however the equations beyond that remain unstable.

Do you know where the imaginary parts are coming from beyond 3.07?

Great idea to use the fast 3.18 to explore very large t, x. Will try this out.

P.S.:
print("eq.328 ",Ht324(1000, -10)) should have been print("eq.328 ",Ht328(1000, -10)) of course.

[teorth]

I think the small imaginary parts are coming from the lower order terms in the Stirling approximation that we are omitting (they are of order about O(1/12s), which is consistent with the fact that the imaginary part is four orders of magnitude smaller than the real part from 3.07 onwards). It can be fixed by adding additional terms to definition of M_0, but given that we have many more serious accuracy issues, it probably isn't worth expending the effort to fix it at this stage.

[teorth]

The equation (3.18), by the way, is almost already entirely in (N,u) coordinates (replace xt by 4PiN^2) after discarding the multiplicative factor, so we should for instance now be able to fill in the missing portions of https://drive.google.com/file/d/101wtC7VyM9aHTst2F-X1Xf62FQT0m6e9/view
relatively easily.

[rudolph-git-acc]

Have played quite a bit with 3.18. It is really fast and also appears to be sufficiently accurate as long as x / |t| > 10 (-ish). For other values it quickly drifts away from the correct value of Ht.

I see what you mean with the (N,u) coordinates, however once the multiplicative factor is removed, the zeros of 3.18 immediately become very 'noisy' again. I have seen a similar effect when I first tested the (unbounded) Ht-integral for larger negative t and noticed that a 'normaliser' (like e.g. multiplying by 10^166) yielded more stable results.

The Ht-zeros used for the graph you refer to were in the range t=-200...0, x=20..600 and we have to go higher to test the new equation. Therefore just kicked off a real root finding run for t=100..0, x=10000...10600 and will overlay this with the Ht data (did some test runs and expect this to all match).

[rudolph-git-acc]

There is a perfect fit between the zeros of Ht-actual and eq.3.18 (red + green turns into brownish):

Will now expand the data in the t-direction to see if the pattern is re-emerging, however expect divergence to start when t decreases < -800.

[teorth]

The perfect fit is very promising! It suggests that much of the numerical deviation between H_t and 3.18 is coming from the multiplier rather than the sum.

By the way, I think I know why your root finding algorithm creates spurious zeroes when the multiplier is removed - it's because the multiplier and sum are both complex numbers, with approximately opposite phase; only the product is (approximately) real. So any attempt to find zeroes based on detecting sign changes in the sum is going to create artificial zeroes when the phase is close to \pm \pi/2. On the other hand this means that one can remove any purely real components from the multiplier (in particular M0 can be replaced by its phase M0/abs(M0)), which will remove the exponentially small magnitude in these sums and perhaps lead to slightly more accurate numerics.

I'm curious to see exactly where the breakdown in accuracy from (3.18) to (3.20) is coming from. There are three substitutions being made here:

• replacing 1 + I*xt/2 with I*xt/2;
• replacing the first log(n/N) with (n-N)/N;
• replacing the second log(n/N) with (n-N)/N - (n-N)^2/2N^2.

From my own primitive testing it seemed that the first change was fairly accurate but the second one was not for smallish values of t,x, but it should get more accurate for larger values of t,x. (Amusingly, replacing the first log(n/N) with the seemingly more accurate (n-N)/N - (n-N)^2/2N^2 gives far worse results due to the fact that the latter approximation has an additional zero leading to some new very large terms in the sum.)

[rudolph-git-acc]

The phase M0/abs(M0)) works great and no longer a need for any other multiplicative factor than that. Have also tried Newton-Raphson instead of sign-changed based root finding, however this worked only successfully for 3.18 with the new factor (still drift in 3.20).

The graph has been extended to t=-300 and still has a perfect fit with the Ht-actual zeros. At t=-1000 the drift becomes ~0.3 at each zero (i.e. Ht-actual is 0.3 higher than 3.18). This increases to ~0.7 at t=-1800 (here only a single 'straight line' is left to detect). Will generate the roots for x=10000...10600, t=-1800...-10 in steps of 10, so we have that ready when we break through the 3.20 Barrier :)

[km-git-acc]

Rudolph, just to confirm, is this the eqn 3.18 function you are using post the changes?

Ht318_new(x,t) = {
    xt = x-Pi*abs(t)/4;
    N = sqrt(xt/(4*Pi));
    u = 4*Pi*abs(t)/xt;
    m0val = M0((1+I*xt)/2);
    mulfact = (m0val/abs(m0val))*N^(-(1+I*xt)/2)*sqrt(Pi/(1-I*u/(8*Pi)));
    A=mulfact*sum(n=1,1000,exp(-abs(t)*(log(n/N))^2/(4*(1-I*u/(8*Pi)))-(1+I*xt)/2*log(n/N)));
    return(A);
}

Ht318_new(1000,-10) -> 1.55003 + 0.00051*I

Using the root finder at the end, the output is close to the one in the first post.

-10.00000000, 906.76541539
-10.00000000, 931.50799535
-10.00000000, 976.83695739
-10.00000000, 990.48752071

[rudolph-git-acc]

Yes, that's the one. You could simplify the mulfact by only using M0/abs(M0)). Eq. 3.18 does deliver very well as long as x / |t| > 10; it is eqs. 3.20 onwards that are still unstable when establishing the zeros.

[km-git-acc]

It seems mulfact still requires the N^() term. On removing it,

Ht318_new(x,t) = {
    xt = x-Pi*abs(t)/4;
    N = sqrt(xt/(4*Pi));
    u = 4*Pi*abs(t)/xt;
    m0val = M0((1+I*xt)/2);
    mulfact = (m0val/abs(m0val));
    A=mulfact*sum(n=1,1000,exp(-abs(t)*(log(n/N))^2/(4*(1-I*u/(8*Pi)))-(1+I*xt)/2*log(n/N)));
    return(A);
}

Ht318_new(1000,-10) -> -2.54492 + 0.56468*I
the imaginary term becomes significant unlike in the previous version

The root finder also gives different and a lot more zeros,

-10.00000000, 902.09670707
-10.00000000, 904.48290536
-10.00000000, 909.25164595
-10.00000000, 911.63424081
-10.00000000, 914.01563999
..
..

Also while dealing with large x, it may be useful to change the number of terms in the sum dynamically (for eg sqrt(x) instead of 1000),

Ht318_new2(x,t) = {
    xt = x-Pi*abs(t)/4;
    N = sqrt(xt/(4*Pi));
    u = 4*Pi*abs(t)/xt;
    m0val = M0((1+I*xt)/2);
    mulfact = (m0val/abs(m0val))*N^(-(1+I*xt)/2)*sqrt(Pi/(1-I*u/(8*Pi)));
    A=mulfact*sum(n=1,floor(sqrt(x)),exp(-abs(t)*(log(n/N))^2/(4*(1-I*u/(8*Pi)))-(1+I*xt)/2*log(n/N)));
    return(A);
}

for eg. check the behavior of the new vs new2 functions at x=10^7 or higher.

[rudolph-git-acc]

That is strange. I am currently running x=10000...11000, t=-1800...-10 in steps of 10, with the code below (at realprecision=40) and don't have issues with the root finding. Note that I already use sqrt(x) and count to 100 rather than 1000 (could that explain the difference?):

EDIT 1:
Oops. You are correct. I use log(n) in the last part and therefore have implicitly done the N^() multiplication !

Ht318(x,t) = {
    xt = x-Pi*abs(t)/4;
    N = sqrt(xt/(4*Pi));
    u = 4*Pi*abs(t)/xt;
    mulfact = exp(-Pi^2*abs(t)/64)*(M0((1+I*xt)/2))*N^(-0*(1+I*xt)/2)*sqrt(Pi);
    A=M0((1+I*xt)/2)/abs(M0((1+I*xt)/2))*sum(n=1,100,exp(-abs(t)*(log(n/N))^2/(4*(1-I*u/(8*Pi)))-(1+I*xt)/2*log(n)));
    return(A);
}

EDIT 2:
If I have done the math correctly, moving from 3.16 to 3.17 just gives log(n) as the last log in 3.17. The non-zero multiplicative N^() is then brought in for the later steps. However, when I skip this multiplication and just continue the last log(n) in 3.20, it does produce the correct zeros (this works all the way to 3.23):

Ht320(x,t) = {
    xt = x-Pi*abs(t)/4;
    N = sqrt(xt/(4*Pi));
    u = 4*Pi*abs(t)/xt;
    A = M0((1+I*xt)/2)/abs(M0((1+I*xt)/2))*sum(n=1,100,exp(-(abs(t)*(n-N)^2)/(4*N^2*(1-I*u/(8*Pi)))-(1+I*xt)/2*log(n)));
    return(A);
}

[km-git-acc]

Great. I was able to match the output with the above Ht318 function after removing the sqrt(Pi/(1-I*u/(8*Pi))) factor. Had kept it since it is complex, but the imaginary part is small.

The Ht320 function is giving the correct zeros in the test range, although the imaginary part here is not insignificant which may require caution.

Ht320(1000,-10) -> 0.69424 - 0.20181*I

[rudolph-git-acc]

Here is a plot of both eq 3.18 and Ht-actual in the range: x=10000..15000, t=-3000...-10 (step 10). The brownish color indicates the surprisingly good fit (< 1 difference) even at t=-3000 and using n=1...100.

Normalized version: