How can I optimize and improve the accuracy of a summation involving the Mangoldt function in MATLAB?

2 次查看（过去 30 天）
Fatima Majeed
Fatima Majeed 2025-1-1
评论： Fatima Majeed 2025-1-6
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
  7 个评论
显示 5更早的评论

请先登录，再进行评论。

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

采纳的回答

Gayathri
Gayathri 2025-1-2
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!
  4 个评论
显示 2更早的评论
Walter Roberson
Walter Roberson 2025-1-3
Ran in:
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

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Creating and Concatenating Matrices 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by