Join the Stack Overflow Community
Stack Overflow is a community of 6.6 million programmers, just like you, helping each other.
Join them; it only takes a minute:
Sign up

Is there any algorithm to calculate (1^x + 2^x + 3^x + ... + n^x) mod 1000000007?
Note: a^b is the b-th power of a.

The constraints are 1 <= n <= 10^16, 1 <= x <= 1000. So the value of N is very large.

I can only solve for O(m log m) if m = 1000000007. It is very slow because the time limit is 2 secs.

Do you have any efficient algorithm?

There was a comment that it might be duplicate of this question, but it is definitely different.

share|improve this question
    
@TobySpeight No, it is definitely different. 1^x, 2^x..., n^x is mod exponential, and I want to find FAST ALGORITHM because n <= 10^16 and mod = 10^9+7. – square1001 3 hours ago
    
n^x mod m is the same as ((n^(x-1) mod m) * (n mod m)) mod m; (a + b) mod m is the same as ((a mod m) + (b mod m)) mod m. – Vatine 3 hours ago
    
@TobySpeight In my slow algorithm, I have to take mod about 10^10 times. But it is very slow because the time limit is 2sec. I only want to find fast algorithm. – square1001 3 hours ago
1  
@CiaPan Please note that 10^9+7 is a prime. – square1001 3 hours ago
2  
Give a look to en.wikipedia.org/wiki/Faulhaber's_formula. If your problem comes from a challenge site like HackerRank, I bet the inner product always simplifies to an integer and that there is a way to get rid of the global factor as well. – Fabian Pijcke 3 hours ago
up vote 9 down vote accepted

You can sum up the serie

1**X + 2**X + ... + N**X

with a help of Faulhaber's formula and you'll get a polynome with X + 1 power to compute for arbitrary N.

If you don't want to compute Bernoulli Numbers, you can find the the polynom by solving X + 2 linear equations (for N = 1, N = 2, N = 3, ..., N = X + 2) which is a slower method but easier to implement.

Let's have an example for X = 2. In this case we have X + 1 = 3 order polynom

    A*N**3 + B*N**2 + C**N + D

The linear equations are

      A +    B +   C + D = 1              =  1 
    A*8 +  B*4 + C*2 + D = 1 + 4          =  5
   A*27 +  B*9 + C*3 + D = 1 + 4 + 9      = 14
   A*64 + B*16 + C*4 + D = 1 + 4 + 9 + 16 = 30 

having solved the equations we'll get

  A = 1/3
  B = 1/2
  C = 1/6
  D =   0 

The final formula is

  1**2 + 2**2 + ... + N**2 == N**3 / 3 + N**2 / 2 + N / 6

Now, all you have to do it to put an arbitrary large N into the formula. So far the algorithm has O(X**2) complexity (doesn't depend on N)

share|improve this answer
    
But unfortunately x might be 1000... Solving linear equation can be O(X^3), but I think it is slow. I am checking about Faulhaber's Formula right now. – square1001 3 hours ago
    
@square1001: if you have to be very fast, I suggest precomputing Bernoulli Numbers (up to 1000); it's not an easy task (like linear equations solving) but in this case you'll get the polynom almost ready to use. – Dmitry Bychenko 2 hours ago
    
I checked about Bernoulli number, and finally understood that B[1], B[2],..., B[x] for O(x^2). So this problem can solve for O(N^2) ! – square1001 2 hours ago
    
@square1001: In O(X**2), to be true, it doesn't depend N (number of items to sum), but of power X they should be raised – Dmitry Bychenko 2 hours ago
1  
Implementated. I got Accepted. – square1001 2 hours ago

There are a few ways of speeding up modular exponentiation. From here on, I will use ** to denote "exponentiate" and % to denote "modulus".

First a few observations. It is always the case that (a * b) % m is ((a % m) * (b % m)) % m. It is also always the case that a ** n is the same as (a ** floor(n / 2)) * (a ** (n - floor(n/2)). This means that for an exponent <= 1000, we can always complete the exponentiation in at most 20 multiplications (and 21 mods).

We can also skip quite a few calculations, since (a ** b) % m is the same as ((a % m) ** b) % m and if m is significantly lower than n, we simply have multiple repeating sums, with a "tail" of a partial repeat.

share|improve this answer
    
But your algorithm has to calculate 1^x, 2^x, ..., 1000000006^x. It has over 10^9 terms, so you have to iterate over 10^10 times. So, I think it can't pass because the time limit is 2 secs. – square1001 3 hours ago
1  
@square1001 Just tested a NOP loop on my laptop, takes ~0.1 seconds to iterate 10 ** 10 times and that is likely dominated by "fork" and "exec". – Vatine 3 hours ago
    
This becomes much easier if x is odd, because, for example, 1000000006 is congruent to -1, so 1^x and 1000000006^x cancel out, and so do many other terms. I don’t see a shortcut if x is even, though. – Tom Zych 3 hours ago
    
And if x is known in advance, you could build a closed formula for the (non-mod) sum of exponentials. But, yes, you can probably collapse this to "calculate 500000003 ** x % m; multiply by floor(n / m); loop for i from 1 to n % m, exponentiate and sum" – Vatine 3 hours ago

I think Vatine’s answer is probably the way to go, but I already typed this up and it may be useful, for this or for someone else’s similar problem.

I don’t have time this morning for a detailed answer, but consider this. 1^2 + 2^2 + 3^2 + ... + n^2 would take O(n) steps to compute directly. However, it’s equivalent to (n(n+1)(2n+1)/6), which can be computed in O(1) time. A similar equivalence exists for any higher integral power x.

There may be a general solution to such problems; I don’t know of one, and Wolfram Alpha doesn’t seem to know of one either. However, in general the equivalent expression is of degree x+1, and can be worked out by computing some sample values and solving a set of linear equations for the coefficients.

However, this would be difficult for large x, such as 1000 as in your problem, and probably could not be done within 2 seconds.

Perhaps someone who knows more math than I do can turn this into a workable solution?

Edit: Whoops, I see Fabian Pijcke had already posted a useful link about Faulhaber's formula before I posted.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.