Speed problem for summation (sum of divisors)

Question

I should implement this summation in C ++. I have tried with this code, but with very high numbers up to 10 ^ 12 it takes too long.

The summation is:

For any positive integer k, let d(k) denote the number of positive divisors of k (including 1 and k itself). For example, for the number 4: 1 has 1 divisor, 2 has two divisors, 3 has two divisors, and 4 has three divisors. So the result would be 8.

This is my code:

#include <iostream>
#include <algorithm>

using namespace std;

int findDivisors(long long n) 
{
    int c=0;
    for(int j=1;j*j<=n;j++)
    {
        if(n%j==0)
        {
            c++;
            if(j!=(n/j))
            {
                c++;
            }
        }
    }
    return c;
}

long long compute(long long n)
{
    long long sum=0;
    for(int i=1; i<=n; i++)
    {
        sum += (findDivisors(i));
    }
    return sum;
}

int main()
{
    int n, divisors;

    freopen("input.txt", "r", stdin);
    freopen("output.txt", "w", stdout);

    cin >> n;

    cout << compute(n);
}

I think it's not just a simple optimization problem, but maybe I should change the algorithm entirely. Would anyone have any ideas to speed it up? Thank you.

Answer 1

largest_prime_is_463035818's answer shows an O(N) solution, but the OP is trying to solve this problem

with very high numbers up to 10¹².

The following is an O(N^1/2) algorithm, based on some observations about the sum

n/1 + n/2 + n/3 + ... + n/n

In particular, we can count the number of terms with a specific value.

Consider all the terms n/k where k > n/2. There are n/2 of those and all are equal to 1 (integer division), so that their sum is n/2.

Similar considerations hold for the other dividends, so that we can write the following function

long long count_divisors(long long n)
{
    auto sum{ n };
    for (auto i{ 1ll }, k_old{ n }, k{ n }; i < k ; ++i, k_old = k)
    { //                                    ^^^^^ it goes up to sqrt(n)
        k = n / (i + 1);
        sum += (k_old - k) * i;
        if (i == k)
            break;
        sum += k;
    }
    
    return sum;   
}

Here it is tested against the O(N) algorithm, the only difference in the results beeing the corner cases n = 0 and n = 1.

Edit

Thanks again to largest_prime_is_463035818, who linked the Wikipedia page about the divisor summatory function, where both an O(N) and an O(sqrt(N)) algorithm are mentioned.

An implementation of the latter may look like this

auto divisor_summatory(long long n)
{
    auto sum{ 0ll };
    auto k{ 1ll };
    for ( ; k <= n / k; ++k )
    {
        sum += n / k;
    }
    --k;
    return 2 * sum - k * k;
}

They also add this statement:

Finding a closed form for this summed expression seems to be beyond the techniques available, but it is possible to give approximations. The leading behavior of the series is given by

D(x) = xlogx + x(2γ - 1) + Δ(x)

where γ is the Euler–Mascheroni constant, and the error term is Δ(x) = O(sqrt(x)).

Answer 2

Let's start off with some math and reduce the O(n * sq(n)) factorization to O(n * log(log(n))) and for counting the sum of divisors the overall complexity is O(n * log(log(n)) + n * n^(1/3)).

For instance:

In Codeforces himanshujaju explains how we can optimize the solution of finding divisors of a number. I am simplifying it a little bit.

Let, n as the product of three numbers p, q, and r.

so assume p * q * r = n, where p <= q <= r. 
The maximum value of p = n^(1/3).

Now we can loop over all prime numbers in a range [2, n^(1/3)]
and try to reduce the time complexity of prime factorization.

We will split our number n into two numbers x and y => x * y = n. 

And x contains prime factors up to n^(1/3) and y deals with higher prime factors greater than n^(1/3).

Thus gcd(x, y) = 1.
Now define F(n) as the number of prime factors of n. 

From multiplicative rules, we can say that 

F(x * y) = F(x) * F(y), if gcd(x, y) = 1.

For finding F(n) => F(x * y) = F(x) * F(y)

So first find F(x) then F(y) will F(n/x)

And there will 3 cases to cover for y:
1. y is a prime number: F(y) = 2.
2. y is the square of a prime number: F(y) = 3.
3. y is a product of two distinct prime numbers: F(y) = 4.

So once we are done with finding F(x) and F(y), we are also done with finding F(x * y) or F(n).

In Cp-Algorithm there is also a nice explanation of how to count the number of divisors on a number. And also in GeeksForGeeks a nice coding example of how to count the number of divisors of a number in an efficient way. One can check the articles and can generate a nice solution to this problem.

C++ implementation

#include <bits/stdc++.h>
using namespace std;
const int maxn = 1e6 + 11;
bool prime[maxn];
bool primesquare[maxn];
int table[maxn]; // for storing primes

void SieveOfEratosthenes()
{

    for(int i = 2; i < maxn; i++){
        prime[i] = true;
    }

    for(int i = 0; i < maxn; i++){
        primesquare[i] = false;
    }

    // 1 is not a prime number
    prime[1] = false;

    for(int p = 2; p * p < maxn; p++){
        // If prime[p] is not changed, then
        // it is a prime
        if(prime[p] == true){
            // Update all multiples of p
            for(int i = p * 2; i < maxn; i += p){
                prime[i] = false;
            }
        }
    }

    int j = 0;
    for(int p = 2; p < maxn; p++) {
        if (prime[p]) {
            // Storing primes in an array
            table[j] = p;

            // Update value in primesquare[p * p],
            // if p is prime.
            if(p < maxn / p) primesquare[p * p] = true;
            j++;
        }
    }
}

// Function to count divisors
int countDivisors(int n)
{
    // If number is 1, then it will have only 1
    // as a factor. So, total factors will be 1.
    if (n == 1)
        return 1;

    // ans will contain total number of distinct
    // divisors
    int ans = 1;

    // Loop for counting factors of n
    for(int i = 0;; i++){
        // table[i] is not less than cube root n
        if(table[i] * table[i] * table[i] > n)
            break;

        // Calculating power of table[i] in n.
        int cnt = 1; // cnt is power of prime table[i] in n.
        while (n % table[i] == 0){ // if table[i] is a factor of n
            n = n / table[i];
            cnt = cnt + 1; // incrementing power
        }

        // Calculating the number of divisors
        // If n = a^p * b^q then total divisors of n
        // are (p+1)*(q+1)
        ans = ans * cnt;
    }

    // if table[i] is greater than cube root of n

    // First case
    if (prime[n])
        ans = ans * 2;

    // Second case
    else if (primesquare[n])
        ans = ans * 3;

    // Third case
    else if (n != 1)
        ans = ans * 4;

    return ans; // Total divisors
}

int main()
{
    SieveOfEratosthenes();
    int sum = 0;
    int n = 5;
    for(int i = 1; i <= n; i++){
        sum += countDivisors(i);
    }
    cout << sum << endl;
    return 0;
}

Output

n = 4 => 8
n = 5 => 10

Complexity

Time complexity: O(n * log(log(n)) + n * n^(1/3))

Space complexity: O(n)

Thanks, @largest_prime_is_463035818 for pointing out my mistake.

Answer 3

I used your brute force approach as reference to have test cases. The ones I used are

compute(12) == 35
cpmpute(100) == 482

Don't get confused by computing factorizations. There are some tricks one can play when factorizing numbers, but you actually don't need any of that. The solution is a plain simple O(N) loop:

#include <iostream>
#include <limits>

long long compute(long long n){
    long long sum = n+1;
    for (long long i=2; i < n ; ++i){
        sum += n/i;
    }
    return sum;
}

int main()
{
    std::cout << compute(12) << "\n";
    std::cout << compute(100) << "\n";
}

Output:

35
482

Why does this work?

The key is in Marc Glisse's comment:

As often with this kind of problem, this sum actually counts pairs x, y where x divides y, and the sum is arranged to count first all x corresponding to a fixed y, but nothing says you have to keep it that way.

I could stop here, because the comment already explains it all. Though, if it didn't click yet...

The trick is to realize that it is much simpler to count divisors of all numbers up to n rather than n-times counting divisors of individual numbers and take the sum.

You don't need to care about factorizations of eg 123123123 or 52323423 to count all divisors up to 10000000000. All you need is a change of perspective. Instead of trying to factorize numbers, consider the divisors. How often does the divisor 1 appear up to n? Simple: n-times. How often does the divisor 2 appear? Still simple: n/2 times, because every second number is divisible by 2. Divisor 3? Every 3rd number is divisible by 3. I hope you can see the pattern already.

You could even reduce the loop to only loop till n/2, because bigger numbers obviously appear only once as divisor. Though I didn't bother to go further, because the biggest change is from your O(N * sqrt(N)) to O(N).