Clarification on Calculations for $v_{1,z} $ in "A Method for Verifying the Generalized Riemann Hypothesis"

Question

Fatima Majeed 2024-12-31

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2172588-clarification-on-calculations-for-v_-1-z-in-a-method-for-verifying-the-generalized-riemann-hypo

编辑： Fatima Majeed 2025-1-6

I am reading the paper *"[A Method for Verifying the Generalized Riemann Hypothesis]( https://arxiv.org/abs/2408.00187 )"* by Ghaith A. Hiary, Summer Ireland, and Megan Kyi. On page 18,

is calculated using Corollary 9 and Lemma 10 with parameters

, $z = 10^{28} + 501675.8%$, and

. The result given is $\operatorname{Re}(v_{1,z}) \in [31.418062, 31.41806]$.

The definitions used are as follows:

- Corollary 9 states:

where ψ is the digamma function.

- Lemma 10 gives an approximation for

:

with the error term

bounded as

.

I attempted to reproduce the calculation using MATLAB. Below is my implementation and results:

---

### Attempt to Calculate $ \psi $ (Digamma Function)

matlab

% Define the complex number
z = -2 + (10^28 + 501675.8)*1i;
s = (3-z)/2;
% Euler-Mascheroni constant
gamma = 0.57721566490153286060651209;
% Number of terms in the series
N = 1e5; % Adjust as needed for precision
% Compute the series sum
series_sum = 0;
for n = 0:N
    series_sum = series_sum + (1/(n+1) - 1/(n+s));
end
% Compute the digamma value
digamma_value = -gamma + series_sum;
disp(['Digamma value at z = ', num2str(real(s)), ' + ', num2str(imag(s)), 'i is: ', num2str(digamma_value)]);

Result:

---

### Calculate $R(K, x) $

matlab

K = 10^7;
R = -(K^(-2))/2 * (2.85*(-5)/log(K) - 1);
disp(R);

Result:

---

### Approximate $ \sum_{k=1}^K \frac{\Lambda(k)}{k^{1-z}}$

matlab

K = 10^7;
% Initialize the sum
sum_value = 0;
% Loop through k from 1 to K
for k = 1:K
    % Calculate the Mangoldt function for k
    L_k = mangoldt(k);
    
    % Add the term to the sum
    sum_value = sum_value + L_k / k^(1 - (10^28 + 501675.8)*1i + 2);
end
disp(sum_value);
% Define the smallest_prime_factor function
function p = smallest_prime_factor(n)
for p = 2:sqrt(n)
    if mod(n, p) == 0
        return; % Return the first divisor
    end
end
p = n; % If no divisor is found, n is prime
end
% Define the Mangoldt function
function L = mangoldt(n)
if n == 1
    L = 0; % Von Mangoldt function is 0 for n = 1
else
    % Find the smallest prime factor of n
    p = smallest_prime_factor(n);
    
    % Check if n is a power of p
    k = log(n) / log(p);
    if abs(k - round(k)) < 1e-10 % Tolerance for floating-point errors
        L = log(p); % Logarithm of the prime
    else
        L = 0; % Not a prime power
    end
end
end

Result:

---

### Final Computation for $ v_{1,z} $

Using the components calculated above:

matlab

v = -0.5*log(pi) - 1/(-2 + (10^28 + 501675.8)*1i) + 0.0509 + 0.0094i + (9.4205e-15) + 0.5*(11.5129 - 2.00002e-23i);
disp(v);

Result:

---

Updated Attempt to Calculate ψ (Digamma Function)

z=-2 + (sym(10^28) + 501675.8)*1i
psi(sym((3-z)/2))
var = vpa(ans)

Result :

z =

vpa("63.779235423273333843086578777284") - vpa("1.5707963267948966192313216912398i")

----

#####Updated Attemp to Calculate $R(K, x) $

K = 10^7;
R = -(vpa(K^(-2)))/2 * (2.85*(-5)/vpa(log(K)) - 1);
disp(R);

Result:

vpa("0.0000000000000094204974050866702154395576247119")

------

Updated Approximate $ \sum_{k=1}^K \frac{\Lambda(k)}{k^{1-z}}$

K = 10^7;
sum_value = 0;
% Precompute primes 
limit = floor(sqrt(K));
is_prime = true(1, limit);
is_prime(1) = false; 
for i = 2:sqrt(limit)
    if is_prime(i)
        is_prime(i*i:i:limit) = false;
    end
end
primes_list = find(is_prime);
% Use parallel computing for summation
parfor k = 1:K
    L_k = mangoldt(k, primes_list);
    sum_value = sum_value + L_k / k^3;
end
disp(sum_value);
% Define the Mangoldt function with precomputed primes
function L = mangoldt(n, primes_list)
    if n == 1
        L = 0; % Von Mangoldt function is 0 for n = 1
    else
        % Find the smallest prime factor of n using precomputed primes
        for p = primes_list
            if mod(n, p) == 0
                % Check if n is a power of p
                k = log(n) / log(p);
                if abs(k - round(k)) < 1e-10 % Tolerance for floating-point errors
                    L = log(p); % Logarithm of the prime
                    return;
                else
                    L = 0; % Not a prime power
                    return;
                end
            end
        end
        L = 0; % If no prime factor found, n is not a power of any prime in the list
    end
end

Result:

Starting parallel pool (parpool) using the 'Processes' profile ...

Connected to parallel pool with 2 workers.

0.1648

____

Updated Final Computation for $ v_{1,z} $

v =vpa( -0.5*log(pi) - 1/(-2 + (sym(10)^28 + 501675.8)*1i) - R -0.0509-0.0094i + 0.5*sym(var));
disp(v);

Result:

-----------

My result does not match the range provided in the paper ($ \operatorname{Re}(v_{1,z}) \in [31.418062, 31.41806] $). Could someone help identify where my calculations or implementation might be incorrect? Additionally, are there specific precision considerations or assumptions in the paper that I may have overlooked?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

John D'Errico 2024-12-31

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2172588-clarification-on-calculations-for-v_-1-z-in-a-method-for-verifying-the-generalized-riemann-hypo#answer_1556731

编辑：John D'Errico 2024-12-31

在 MATLAB Online 中打开

I think, from almost line 1 of what you did, you don't understand double precision arithmetic, and that will cause everything you do in this enterprise to fail.

z = -2 + (10^28 + 501675.8)*1i;

OOPS. You did it right there. Mistake 1. And everything thereafter will fail, or at least be suspect, because that number is not what you think it is.

A complex number, thus what you have in z, is just two doubles, the real and imaginary parts. Now look at what you wrote in the imaginary part:

10^28 + 501675.8

so you are adding 501675.8 to a number on the order of 1e28. What is eps(1e28)?

eps(1e28)
ans = 2.1990e+12

eps gives you the distance from 1e28, and the nest larger number you can store in a double. You can think of it as effectively the smallest number you can add to 1e28, and get a different result.

Trying to add a number smaller than 2.19e12 to 1e28, this will get you nowhere. In fact, try this:

(1e28 +  501675.8) == 1e28
ans = logical
   1

Do you see that anything you now do with z is going to be potentially mathematical gibberish, since z is not what you think it is?

What else? Lets see.

gamma = 0.57721566490153286060651209;

Most of the digits in that number are useless, since MATLAB will only be able to store the first 16 digits or so as a double. (And even worse, you just overloaded the highly useful function gamma, so the gamma function will not work for you, in case you need it at some point. Don't name variables in a way that overloads existing function names. You will be sorry soon, and then you will be asking on Answers how to make the gamma function work.)

Until you understand that a double has limits in what it can do, you are wasting your time. That is not to say that what you want to do is a waste of imte. But you need to learn to use computational tools properly, else nothing you do will succeed.

So, COULD you do anything with this computation, even if using high precision tools? Possibly. I won't even try to get into that argument. But you MIGHT conceivably then have a chance, using syms and vpa. Even at that, you will need to learn how to work with high precision tools, as it is very easy to still have a problem, if you mix numbers that are entered into MATLAB as doubles, and then turn them into syms.

There are other things you do, that are just a bit silly, not making use of the tools in MATLAB.

function p = smallest_prime_factor(n)
for p = 2:sqrt(n)
    if mod(n, p) == 0
        return; % Return the first divisor
    end
end
p = n; % If no divisor is found, n is prime
end

So, a loop to find the smallest prime factor? What does this do?

smallestPrimeFactor = @(n) min(factor(n));

Learn to use the tools already in MATLAB. I suppose that might even extend to learning to use the function psi. But you might not want to go that far, apparently trying to do at least that piece on your own.

lookfor digamma
psi                            - Digamma and polygamma functions
sym.psi                        - Digamma function

4 个评论
显示 2更早的评论隐藏 2更早的评论

Fatima Majeed 2024-12-31

编辑：Fatima Majeed 2024-12-31

Thank you so much for your detailed explanation and for introducing me to syms and vpa! I now understand the limitations of double precision in MATLAB and how symbolic computation helps achieve the accuracy needed for such calculations.

I calculated Re(v1,z) using vpa and your suggestion regarding the digamma function and found the value:

Re(v1,z)=31.26635276871197626010195638099

This is much closer to the expected value, but it’s still an approximation and not exactly the same as the value in the paper. Could this discrepancy be due to numerical methods or the symbolic computation itself? Or is there a particular adjustment I might be missing in my approach?

Walter Roberson 2025-1-2

By the way, the output of factor() is always sorted in practice (though that is not documented.)

请先登录，再进行评论。

Clarification on Calculations for $v_{1,z} $ in "A Method for Verifying the Generalized Riemann Hypothesis"

Updated Approximate $ \sum_{k=1}^K \frac{\Lambda(k)}{k^{1-z}}$

Updated Final Computation for $ v_{1,z} $

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

回答（1 个）

4 个评论
显示 2更早的评论隐藏 2更早的评论

另请参阅

类别

标签

Community Treasure Hunt

Clarification on Calculations for $v_{1,z} $ in "A Method for Verifying the Generalized Riemann Hypothesis"

Updated Approximate $ \sum_{k=1}^K \frac{\Lambda(k)}{k^{1-z}}$

Updated Final Computation for $ v_{1,z} $

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

回答（1 个）

4 个评论 显示 2更早的评论隐藏 2更早的评论

另请参阅

类别

标签

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

4 个评论
显示 2更早的评论隐藏 2更早的评论