Summing up the divisors raised to power k

  • MHB
  • Thread starter Saitama
  • Start date
  • Tags
    Power
In summary: Note that $F(n,k) = \sum_{d|n} d^k$.So for each $i$, we need to find $\sum_{i \le d \le j} d^k$.To do this in $O(\log n)$, we can use a segment tree.The segment tree is an array of $n$ elements, each of which is the sum of a contiguous range of numbers in the original array.That is, the $i$th element of the segment tree is the sum of $a_i, a_{i+1}, \dots, a_j$ for some pair $(i,j)$.To find the sum of a range
  • #1
Saitama
4,243
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
  • #2
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}$$
 
  • #3
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;
}
 
  • #4
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! :)
 
  • #5
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)
 
  • #6
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! :)
 
  • #7
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)
 

FAQ: Summing up the divisors raised to power k

What is the concept of "Summing up the divisors raised to power k"?

The concept involves taking each divisor of a given number and raising it to a specified power, then adding all of these values together to get a final sum.

Why is this concept important in mathematics?

This concept is important because it allows for the efficient calculation of various arithmetic functions, such as the sum of divisors, which have applications in number theory, cryptography, and other areas of mathematics.

What is the formula for summing up the divisors raised to power k?

The formula is as follows: ∑(d^k), where d represents each divisor of the given number and k represents the specified power.

Can this concept be applied to any number?

Yes, this concept can be applied to any positive integer, as long as the specified power is a positive integer as well.

What are some real-world applications of this concept?

This concept has applications in areas such as cryptography, where it can be used to generate secure keys, and in number theory, where it can be used to study perfect numbers and amicable numbers.

Similar threads

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

Hot Threads

Recent Insights

Back
Top