trying to compute Riemann's prime counting function J(x)

I am trying to compute Riemann's prime counting function J(x):

J(x) should approximate the numbers of primes <= x using this code:

function J_RiemannPrimeCount = J(x)
    if x < 2
        error("x must be >= 2");
    end
    integral_fun = @(t) (1 ./ (t.*(t.^2-1).*log(t)));
    integral_term = integral(integral_fun,x,Inf);
    
    zetaZeros = 0.5 + csvread("first 100k zeros of the Riemann zeta.txt") .* i;
    maxZero = 35;
    k = 1:1:maxZero;
    Li_term = logint(x.^zetaZeros(k)) - logint(2);
    Li_sum = sum(Li_term);
    
    J_RiemannPrimeCount = (logint(x) - logint(2)) - Li_sum - log(2) + integral_term;

end

The file "first 100k zeros of the Riemann zeta.txt" contains the Imag-values of the Riemann-zeta-function for Re = 0.5 (example: 14.134725142, 21.022039639, 25.010857580, 30.424876126, 32.935061588, ....) . I am using the first maxZero=35 of these to approximate J(x).

The periodic term "Li_sum = sum(Li_term)" is not correct - everything else should be fine. ( I am somewhat worried, that I am doing a super-stupid mistake here... )

Can anyone help nonetheless .. ??

[Thomas]

Hi , would still greatly appreciate if anyone has an idea on this here:

When calculating the explicit Riemann prime counting function J(x) for x=2,...,1000 using the program below, I get an oscillating term Li_sum that looks like in the diagram (see attachement). The local maxima are in the right place (= at the exact locations of the prime numbers): however Li_sum should not be descending in the way it is doing here, but should be distributed around zero, since it is "only" the error term in the smooth estimation. So there must be a kind of "adjustment"-factor somewhere - maybe resulting from the fact that Li_sum is the truncated version of a conditionally convergent series.

% explicit Riemann Prime Counting function J
function J_RiemannPrimeCount = J(x)

if x < 2
    error("x must be >= 2");
end
integral_fun = @(t) (1 ./ (t.*(t.^2-1).*log(t)));
integral_term = integral(integral_fun,x,Inf);
zetaZeros = 0.5 + csvread("first 100k zeros of the Riemann zeta.txt") .* i;
maxZero = 1000;
k = 1:1:maxZero;

Li_term =   logint(x .^ zetaZeros(k)) + logint(x .^ (1-zetaZeros(k)))   - 2*logint(2) ;
Li_sum = sum(Li_term);
J_RiemannPrimeCount_smooth =    (logint(x) - logint(2)) - log(2) + integral_term;   % prime estimation without oscillating term
J_RiemannPrimeCount =           (logint(x) - logint(2)) - log(2) + integral_term - Li_sum; % prime counting with oscillating term

end

Any idea, where the error could be ... ???

Answer 1

[Jyothis Gireesh]

Here are a few pointers which may help with resolving the issue:

• The second term in the Reimann’s function is computed by taking the logarithmic integral over x raised to non-trivial zeros of the zeta function.
• Since the sum is ‘conditionally convergent’, the summation should be done by taking the zeros of zeta function in a pair-wise fashion (taking ρ and 1-ρ) as follows

In the code given above, by restricting the zeros to the first 35 entries in the "first 100k zeros of the Riemann zeta.txt" file you may be violating the pair-wise summation condition causing the sum to diverge. So for ρ value in the first 35 entries you may add the corresponding 1-ρ value as well so that the sum term converges.

[Thomas]

thanks for helping me here... I already tried this. But nonetheless I do not get a converging sum. I am currently using this code (below). The more zeros I add up the more "Li_sum" grows...

function J_RiemannPrimeCount = J(x)
    if x < 2
        error("x must be >= 2");
    end
    integral_fun = @(t) (1 ./ (t.*(t.^2-1).*log(t)));
    integral_term = integral(integral_fun,x,Inf);
    
    zetaZeros = 0.5 + csvread("first 100k zeros of the Riemann zeta.txt") .* i;
    maxZero = 35;
    k = 1:1:maxZero;
    zetaZeros_maxZero = zetaZeros(1:maxZero);
    Li_term = logint(x .^ zetaZeros(k)) + logint(x .^ (1-zetaZeros(k))) - 2*logint(2);
    Li_sum = sum(Li_term) ;
    
    J_RiemannPrimeCount_test = (logint(x) - logint(2)) - 0*Li_sum - log(2) + integral_term;
    J_RiemannPrimeCount = (logint(x) - logint(2)) - Li_sum - log(2) + integral_term;
end

[Miguel Camarasa]

I have the same problem! The term Li_sum goes doesn't converge. Did you solve this? :/

[Thomas]

Unfortunately I could not resolve it. I found an alternative prime counting function, that did provide reasonable results. I published both approaches (the one which works - J2, and the one which does not work - J1) in:

two ways to compute Riemanns prime counting function - File Exchange - MATLAB Central (mathworks.com)

if you have any hints how to get J1 to work, please let me know. I am afraid it is a kind of "silly" error.

[Miguel Camarasa]

Thank you very much Thomas. I was trying to solve the previous one, but it diverges anyway, maybe because of the precision? I don't think so.

[Thomas]

I tried it with variable precision arithmetic and variables using many digits. It did not change anything. Using more Riemann zeros makes the divergence even quicker, so taking more Riemann-zeros into the summation does not help either.

[Miguel Camarasa]

I see, I don't understand why. Anyway, when I run the J2 function, it doesn't seem to get a nice approximation between the number of primes in [2,100] with 100 roots of the Zeta Function. Maybe something is missing?

[Thomas]

I think the approximation becomes better for large x in J2(x) But I did not look into this in detail...

Answer 2

[Johannes van Ek]

Very interesting. I have the same problem. Did you manage to resolve this issue?

I do the summation over the conjugate pairs and it diverges.

I found a paper in Cantor's Paradise on this and it quotes converging results with 35 roots (and 100 roots), and it presents nice graphs of J(x). I do not know how the author got these results.

I also tried the cos(alpa(ln(x)) formulation from riemann's original paper. No success either. Puzzled.

[Thomas]

I found an alternative prime counting function, that did provide reasonable results.

I published both approaches (the one which works - J2, and the one which does not work - J1) in:

two ways to compute Riemanns prime counting function - File Exchange - MATLAB Central (mathworks.com)

if you have any hints how to get J1 to work, please let me know. I am afraid it is a kind of "silly" error.

[Thomas]

btw: do you have a linkt to the paper, that you mention?

[Johannes van Ek]

https://www.cantorsparadise.com/the-riemann-hypothesis-explained-fa01c1f75d3f?source=collection_home---4------2-----------------------

The author closely follows Derbyshire's 'Prime obsession'. In the notes on chapter 21
Derbyshire mentions that het uses Ei[rho ln(x)] rather than Li(x^rho) because of details
in Riemann's original paper.
When I use the exponential integral in mpmath everything works beatifully.
I still need to look into why Li(x^rho) does not work.

[Thomas]

thanks ! I will have a look into this...

[Johannes van Ek]

It is because of branch points of li(x^rho) at 0 an 1 that care needs to ve taken. Anyway, it works now. On to the million dollars!

[Irene Chen]

Do you know how to deal with the branch points? This is so confusing......