Summing up the divisors raised to power k

  • Context: MHB 
  • Thread starter Thread starter Saitama
  • Start date Start date
  • Tags Tags
    Power
Click For Summary
SUMMARY

The discussion focuses on optimizing the solution for the SPOJ problem IITD4, which involves summing divisors raised to a power k. The original implementation using a brute-force approach resulted in Time Limit Exceeded (TLE) errors. Suggestions include limiting iterations to prime factors and utilizing the divisor function σ_k(n) to transform the summation into a product of geometric series. Additionally, implementing a sieve for prime generation and optimizing the calculation of modular powers are recommended to enhance performance.

PREREQUISITES
  • Understanding of modular arithmetic and modular exponentiation
  • Familiarity with the divisor function σ_k(n) and its properties
  • Experience with the Sieve of Eratosthenes for prime number generation
  • Knowledge of C++ programming, particularly with STL and performance optimization techniques
NEXT STEPS
  • Implement the Sieve of Eratosthenes to generate primes up to 100,000 efficiently
  • Research the divisor function σ_k(n) and its applications in number theory
  • Learn about lazy evaluation techniques to optimize repeated calculations
  • Explore advanced modular arithmetic techniques, including the Euclidean algorithm for modular inverses
USEFUL FOR

Competitive programmers, algorithm enthusiasts, and anyone looking to optimize mathematical computations in C++ for performance-critical applications.

Saitama
Messages
4,244
Reaction score
93
Here is the problem I am trying to solve: SPOJ.com - Problem IITD4

Following is what I tried:

Code: 
#include <iostream>
#include <cstdio>
#include <cmath>
#define MOD 1000000007
using namespace std;

typedef long long ll;

ll modularPower(ll base, ll exponent)
{
	ll res=1;
	while(exponent)
	{
		if(exponent%2)
			res=(res*base)%MOD;
		exponent/=2;
		base=(base*base)%MOD;
	}
	return res;
}

ll f(ll i, ll k)
{
	ll j, u, sum=0;
	u=sqrt(i);
	for(j=1; j<=u; ++j)
	{
		if(i%j==0)
		{
			sum=(sum+modularPower(j, k));
			sum=(sum+modularPower(i/j, k));
		}
	}
	if(u*u==i)
	{
		sum=(sum-modularPower(u, k)+MOD)%MOD;
	}
	return sum;
}

int main()
{
	int T, a, b, k, i;
	scanf("%d", &T);
	while(T--)
	{
		ll ans=0;
		scanf(" %d %d %d", &a, &b, &k);
		for(i=a; i<=b; ++i)
			ans=(ans+f(i,k))%MOD;
		cout<<ans<<endl;
	}
}

..but this gives TLE. What optimisations should I try? :confused:

Thanks!
 
Technology news on Phys.org
Pranav said:
Here is the problem I am trying to solve: SPOJ.com - Problem IITD4

Following is what I tried:

Code: 
#include <iostream>
#include <cstdio>
#include <cmath>
#define MOD 1000000007
using namespace std;

typedef long long ll;

ll modularPower(ll base, ll exponent)
{
	ll res=1;
	while(exponent)
	{
		if(exponent%2)
			res=(res*base)%MOD;
		exponent/=2;
		base=(base*base)%MOD;
	}
	return res;
}

ll f(ll i, ll k)
{
	ll j, u, sum=0;
	u=sqrt(i);
	for(j=1; j<=u; ++j)
	{
		if(i%j==0)
		{
			sum=(sum+modularPower(j, k));
			sum=(sum+modularPower(i/j, k));
		}
	}
	if(u*u==i)
	{
		sum=(sum-modularPower(u, k)+MOD)%MOD;
	}
	return sum;
}

int main()
{
	int T, a, b, k, i;
	scanf("%d", &T);
	while(T--)
	{
		ll ans=0;
		scanf(" %d %d %d", &a, &b, &k);
		for(i=a; i<=b; ++i)
			ans=(ans+f(i,k))%MOD;
		cout<<ans<<endl;
	}
}

..but this gives TLE. What optimisations should I try? :confused:

Thanks!

Hi Pranav,

The number of iterations in the loop in your function [m]f[/m] can be reduced.
Instead of iterating over all divisors, we can limit ourselves to the primes.

Note that $F(n,k)$ is the divisor function $\sigma_k(n)$.
In particular, if we break up $n$ into its primes, the summation turns into a summation of a set of geometric series.
So suppose we have the primes and their powers in $n$:
$$n=\prod_{i=1}^r p_i^{a_i}$$
Then:
$$\sigma_k(n) = F(n,k) = \prod_{i=1}^r \frac {p_i^{(a_i+1)k} - 1 }{p_i^k - 1}$$
 
I like Serena said:
Hi Pranav,

The number of iterations in the loop in your function [m]f[/m] can be reduced.
Instead of iterating over all divisors, we can limit ourselves to the primes.

Note that $F(n,k)$ is the divisor function $\sigma_k(n)$.
In particular, if we break up $n$ into its primes, the summation turns into a summation of a set of geometric series.
So suppose we have the primes and their powers in $n$:
$$n=\prod_{i=1}^r p_i^{a_i}$$
Then:
$$\sigma_k(n) = F(n,k) = \prod_{i=1}^r \frac {p_i^{(a_i+1)k} - 1 }{p_i^k - 1}$$

So I guess I have to implement a sieve and make an array of primes <=10^5. Umm...I wonder if it will pass all the cases but still, I will implement this tomorrow morning. Meanwhile, I found the following solution which does not require implementing sieve. I tried to understand it but I am stuck at the ndiv function used. Any idea what it is used for? :confused:

Code: 
#include <iostream>
#include <math.h>
using namespace std;
#define ULL unsigned long long
# define MOD 1000000007;
int ndiv(int n, int m, int p)
{
    return m/p - (n-1)/p;

}

unsigned long long fsum(unsigned long long a, unsigned long long b)
{
    return (a+b)%1000000007;
}
unsigned long long fprod(unsigned long long a,unsigned long long b)
{
    return ((a)*(b))%1000000007;
}
ULL Pow(ULL a, ULL b)
{
    if(a==1)
    {
        return 1;
    }

    ULL x=1,y=a;
    while(b != 0)
    {
        if(b&1)
        {
            x=(x*y);
            if(x>1000000007) x%=1000000007;
        }
        y = (y*y);
        if(y>1000000007) y%=1000000007;
        b /= 2;
    }
    return x;
}
ULL G(ULL a, ULL b, ULL k)
{
    long long z = 0;
    for(int i = 1; i<=b; i++)
    {
        z += fprod(ndiv(a,b,i),Pow(i,k));
    }
    z%=MOD;

    return z;
}
int main()
{
    int T;

    cin>>T;
    ULL a,b,k;
    for(int q = 1; q<=T; q++)
    {
        ((cin>>a)>>b)>>k;
        cout<<G(a,b,k);
        if(q!=T)
        {
            cout<<endl;
        }
    }
    return 0;
}
 
It is as I feared, I am still getting TLE. Following is the code I have written:

Code: 
#include <iostream>
#include <vector>
#define MOD 1000000007
using namespace std;

typedef long long ll;

bool isPrime[100010];
vector<int> p;

void sieve() {
    
    for(int i = 0; i <= 100000;++i) {
        isPrime[i] = true;
    }
    isPrime[0] = false;
    isPrime[1] = false;
    for(int i = 2; i * i <= 100000; ++i) {
         if(isPrime[i] == true) {
             // Mark all the multiples of i as composite numbers
             for(int j = i * i; j <= 100000 ;j += i)
                 isPrime[j] = false;
        }
    }
    for(int i=0; i<=100000; ++i)
    	if(isPrime[i])
    		p.push_back(i);
}

ll modularPower(ll base, ll exponent)
{
	ll res=1;
	while(exponent)
	{
		if(exponent%2)
			res=(res*base)%MOD;
		exponent/=2;
		base=(base*base)%MOD;
	}
	return res;
}

ll f(int i, int k)
{
	int j;
	ll product=1, count=0, aux1, aux2;
	for(j=0; j<p.size(); ++j)
	{
		if(p[j]>i)
			break;
		if(!(i%p[j]))
		{
			count=0;
			while(i%p[j]==0)
			{
				i/=p[j];
				count++;
			}
			aux1=modularPower(p[j],(count+1)*k)-1;
			aux2=modularPower(p[j],k)-1;
			aux2=modularPower(aux2, MOD-2);
			product=(((product*aux1)%MOD)*aux2)%MOD;
		}
	}
	return product;
}

int main() {
	sieve();
	int T, a, b, k;
	cin>>T;
	while(T--)
	{
		cin>>a>>b>>k;
		ll ans=0;
		for(int i=a; i<=b; ++i)
			ans=(ans+f(i,k))%MOD;
		cout<<(ans+MOD)%MOD<<endl;
	}
	return 0;
}

Have I implemented your idea incorrectly? Please let me know.

Also, any clue about the ndiv function in the code I posted in my previous post? Thanks! :)
 
Pranav said:
It is as I feared, I am still getting TLE. Following is the code I have written:

Have I implemented your idea incorrectly? Please let me know.

Your code looks correct except for one thing.

The numerator is guaranteed to be divisible by the denominator.
However, I'm afraid that's no longer guaranteed once we take them modulo $m$.
To evaluate $\frac a b \bmod m$, I'm afraid we will have to find $(b \bmod m)^{-1}$ with the Euclidean algorithm. (Sweating)And as an observation [m]p.push_back(i)[/m] is pretty expensive, since it triggers a reallocation for every iteration.
I've just noticed that the time limit is 0.53 seconds and just these push_back's might already cause an timeout.
Also, any clue about the ndiv function in the code I posted in my previous post? Thanks! :)

With a timeout of 0.53 seconds, we'll have to be smarter.
One of the loops has to be collapsed.
So let's see what happens if we reverse them.

Suppose we want to find G(1,6,2).

$$\DeclareMathOperator{\ndiv}{ndiv}
G(1,6,2) = \sum_{n=1}^6 \sigma_2(n) = \sum_{n=1}^6 \sum_{d|n} d^2 \\
= 1^2 \\
+ 1^2 + 2^2 \\
+ 1^2 \phantom{+ 2^2} + 3^2 \\
+ 1^2 + 2^2 \phantom{+ 3^2} + 4^2 \\
+ 1^2 \phantom{ + 2^2} \phantom{ + 3^2} \phantom{ + 4^2} + 5^2 \\
+ 1^2 + 2^2 + 3^2 \phantom{ + 4^3} \phantom{ + 5^2}+ 6^2 \\
= \ndiv(1, 6,1) \cdot 1^2 + \ndiv(1, 6,2)\cdot2^2 + ... + \ndiv(1,6,6)\cdot 6^2
$$
where $\ndiv(a,b,i)$ is the number of times that a number between $a$ and $b$ is divisible by $i$.

More generally:
$$G(a,b,k) = \sum_{n=1}^b \ndiv(a,b,n) \cdot n^k$$

Now how many numbers between $1$ and $b$ are divisible by, say, $2$? (Wondering)
 
I like Serena said:
And as an observation [m]p.push_back(i)[/m] is pretty expensive, since it triggers a reallocation for every iteration.
I've just noticed that the time limit is 0.53 seconds and just these push_back's might already cause an timeout.
I tried replacing the vector with an array but I still get TLE.

With a timeout of 0.53 seconds, we'll have to be smarter.
One of the loops has to be collapsed.
So let's see what happens if we reverse them.

Sorry, I am a little confused here, which loop are you talking about? :confused:

Suppose we want to find G(1,6,2).

$$\DeclareMathOperator{\ndiv}{ndiv}
G(1,6,2) = \sum_{n=1}^6 \sigma_2(n) = \sum_{n=1}^6 \sum_{d|n} d^2 \\
= 1^2 \\
+ 1^2 + 2^2 \\
+ 1^2 \phantom{+ 2^2} + 3^2 \\
+ 1^2 + 2^2 \phantom{+ 3^2} + 4^2 \\
+ 1^2 \phantom{ + 2^2} \phantom{ + 3^2} \phantom{ + 4^2} + 5^2 \\
+ 1^2 + 2^2 + 3^2 \phantom{ + 4^3} \phantom{ + 5^2}+ 6^2 \\
= \ndiv(1, 6,1) \cdot 1^2 + \ndiv(1, 6,2)\cdot2^2 + ... + \ndiv(1,6,6)\cdot 6^2
$$
where $\ndiv(a,b,i)$ is the number of times that a number between $a$ and $b$ is divisible by $i$.

More generally:
$$G(a,b,k) = \sum_{n=1}^b \ndiv(a,b,n) \cdot n^k$$

Now how many numbers between $1$ and $b$ are divisible by, say, $2$? (Wondering)

Thanks a lot, I finally understand it! :)
 
Pranav said:
I tried replacing the vector with an array but I still get TLE.

For the record, you could make another optimization by storing the result in a table whenever you calculate $p_i^k$ and use the table entry if you already have it (so called lazy evaluation).

But I'm pretty sure it won't be enough to get below 0.53 seconds.
Pranav said:
Sorry, I am a little confused here, which loop are you talking about? :confused:

The outer loop is the summation from $a$ to $b$.
Within that loop we have a summation of the k-powers of the divisors (or primes in your latest version).
The most inner loop is to calculate the powers.

I'm talking about the order of the 2 summations - the outer loops.
They live in different functions, but I still consider them nested loops.

Effectively we've reversed the order of those loops and then collapsed the loop that runs from $a$ to $b$.

Thanks a lot, I finally understand it! :)

Good! (Mmm)
 

Similar threads

C/C++ How can I optimize my C++ code for a sum problem on CodeChef?
  • · Replies 6 ·
Replies
6
Views
3K
Debugging a Backprop Neural Network: Issues & Solutions
  • · Replies 3 ·
Replies
3
Views
2K
MHB Writing a number as sum of squares
  • · Replies 2 ·
Replies
2
Views
2K
MHB How to Efficiently Compute Euler's Totient Function for Large Ranges?
  • · Replies 2 ·
Replies
2
Views
2K
Python The Wimol-Banlue attractor via the Taylor Series Method
  • · Replies 6 ·
Replies
6
Views
2K
C/C++ Optimizing K^(2^N) Mod 10^9+7 in C++
  • · Replies 19 ·
Replies
19
Views
5K
Java Solve Confused Java HW: Generate 225 Samples Size 625 Random Nos U(10,22)
  • · Replies 9 ·
Replies
9
Views
2K
Review C Code: Find & Resolve Issues in Visual Studio & repl.it
  • · Replies 47 ·
2
Replies
47
Views
5K
How to obtain a function from gsl library integration?
  • · Replies 1 ·
Replies
1
Views
2K
Trying to compile a source code in C
  • · Replies 5 ·
Replies
5
Views
2K