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

10 Ansichten (letzte 30 Tage)
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 Kommentare
Torsten
Torsten am 2 Jan. 2025
Bearbeitet: Torsten am 2 Jan. 2025
The default for higher precision arithmetics is 32 digits, but you can enhance precision of your symbolic computations by using the "digits" function:
You didn't include your new code under
so we can't tell if it makes sense what you are doing right now.

Melden Sie sich an, um zu kommentieren.

Akzeptierte Antwort

Gayathri
Gayathri am 2 Jan. 2025
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.
Hope you find this information helpful!
  4 Kommentare
Walter Roberson
Walter Roberson am 3 Jan. 2025
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

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (0)

Kategorien

Mehr zu Matrix Indexing finden Sie in Help Center und File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by