How can I optimize and improve the accuracy of a summation involving the Mangoldt function in 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^(3);
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

[Torsten]

Could you give some explanation about what you are trying to do ?

[Walter Roberson]

It looks like another approach would be

f = factor(n);
if numel(unique(f)) == 1
    L = log(f(1));
else
    L = 0;
end

[Fatima Majeed]

@Torsten I am working on calculating the sum:

\sum_{k=1}^{10^7} \frac{\Lambda(k)}{k^{1-z}},

where \Lambda(k) is the von Mangoldt function.

This calculation is a sub-problem of

https://www.mathworks.com/matlabcentral/answers/2172588-clarification-on-calculations-for-v_-1-z-in-a-method-for-verifying-the-generalized-riemann-hypo?s_tid=answers_rc1-1_p1_MLT

[Torsten]

Estimate

lambda(k) <= k 

for all k.

Thus you get

sum_{k=N}^{k=Inf} Lambda(k)/k^3 <= sum_{k=N}^{k=Inf} 1/k^2

Now determine N such that the last sum is smaller than your desired precision.

You will see that N = 10^7 is far too large for reasonable accuracy.

[Fatima Majeed]

@Walter Roberson Thank you for your answer. However, this code gives the same accuracy as the solution I came up with. Could you please help by looking at my full question

https://www.mathworks.com/matlabcentral/answers/2172588-clarification-on-calculations-for-v_-1-z-in-a-method-for-verifying-the-generalized-riemann-hypo?s_tid=answers_rc1-1_p1_MLT

to see if there’s a way to improve my method?

I would really appreciate your help!

[Torsten]

The default for higher precision arithmetics is 32 digits, but you can enhance precision of your symbolic computations by using the "digits" function:

https://uk.mathworks.com/help/symbolic/digits.html

You didn't include your new code under

https://www.mathworks.com/matlabcentral/answers/2172588-clarification-on-calculations-for-v_-1-z-in-a-method-for-verifying-the-generalized-riemann-hypo?s_tid=answers_rc1-1_p1_MLT

so we can't tell if it makes sense what you are doing right now.

[Fatima Majeed]

@Torsten Thank you for your help! By the way, I've included my updated code.

Answer 1

[Gayathri]

Hi @Fatima Majeed,

To improve the performance of the code, you can prefer using "parfor" loops instead of a "for" loop as shown below.

parfor k = 1:K

This will considerably speed up the code.

Also, you can use a more efficient method to precompute prime number list up to a certain limit. This can speed up the determination of the smallest prime factor.

Please refer to the code below which implements the same.

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

Below is a screenshot that compares the elapsed time for the original code, the code using a "parfor" loop, and the code utilizing both a "parfor" loop and precomputation of the prime list.

For more information on "parfor" loop, please refer to the below documentation link.

https://www.mathworks.com/help/matlab/ref/parfor.html

Hope you find this information helpful!

[Fatima Majeed]

@Gayathri I would like to know how to use normalparforprime. Could you please provide a simple example of how it’s implemented?

Thank you in advance for your guidance!

[Gayathri]

@Fatima Majeed, The code I have provided in the answer is the same code which I have saved under the name "normalparforprime.m". The part of the code shown below, precomputes the prime number list to make the code run faster.

% 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);

[Walter Roberson]

Interesting, that is faster than using primes()

K = 10^7;

tic
limit = floor(sqrt(K));
primes_list1 = primes(limit);
toc

Elapsed time is 0.010151 seconds.

tic
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_list2 = find(is_prime);
toc

Elapsed time is 0.004593 seconds.

isequal(primes_list1, primes_list2)

ans = logical
   1

[Fatima Majeed]

@Walter Roberson and @Walter Roberson thanks alot for your help