이항정리를 이용한 거듭제곱의 합 1^k+2^k+3^k+...+n^k 을 구하는 방법

1. 이항정리(binomial theorem)

 

0이상의 정수 n과 음이 아닌 정수 x,y에 대하여, $$(x+y)^{n} = \sum_{k = 0}^{n}\binom{n}{k}x^{k}y^{n-k}$$

 

https://en.wikipedia.org/wiki/Binomial_theorem

Binomial theorem - Wikipedia

 

y = 1로 두면 다음과 같은 식을 얻는다.

 

$$(x+1)^{n} = \sum_{k=0}^{n} \binom{n}{k} x^{k}$$

 

 

2. 점화식 유도

 

이항정리로부터, $$(m+1)^{k+1} = \sum_{x = 0}^{k+1} \binom{k+1}{x} m^{x}$$

 

우변에 x = k+1인 경우를 제외하면,

 

$$(m+1)^{k+1} - m^{k+1} = \sum_{x=0}^{k} \binom{k+1}{x}m^{x}$$

 

이제 양변에 m = 1부터 n까지 sigma를 취하면,

 

$$\sum_{m = 1}^{n} ((m+1)^{k+1} - m^{k+1}) = \sum_{m=1}^{n} \sum_{x=0}^{k} \binom{k+1}{x}m^{x}$$

 

좌변은 전개하면 망원급수형태인데,

 

$$2^{k+1} - 1^{k+1} + 3^{k+1} - 2^{k+1} + ... +n^{k+1} - (n-1)^{k+1} + (n+1)^{k+1} - n^{k} = (n+1)^{k+1} - 1$$

 

우변에 대해 정리를 해보자면, "푸비니의 정리"는 다중적분의 순서를 바꿀 수 있다는 정리이다.

 

이에 따라 우변의 두 이산합의 순서를 바꿔보면...

 

https://en.wikipedia.org/wiki/Fubini%27s_theorem

Fubini's theorem - Wikipedia

 

$$\sum_{m=1}^{n} \sum_{x=0}^{k} \binom{k+1}{x}m^{x} = \sum_{x=0}^{k} \sum_{m=1}^{n} \binom{k+1}{x}m^{x}$$

 

이항계수 부분은 m과 무관하므로,

 

$$\sum_{x = 0}^{k} \binom{k+1}{x} (\sum_{m = 1}^{n} m^{x})$$

 

$S_{x}(n) = 1^{x} + 2^{x} + ... + n^{x} = \sum_{m = 1}^{n} m^{x}$라고 정의하자. 

 

모든 결과를 종합하면 다음 식을 얻는다.

 

$$(n+1)^{k+1} - 1 = \sum_{x = 0}^{k} \binom{k+1}{x} S_{x}(n)$$

 

 

3. 연습문제

 

25974번: 거듭제곱의 합 1 (acmicpc.net)

25974번: 거듭제곱의 합 1

 

 

4. 풀이

 

n과 k가 주어질때, $S_{k}(n) = 1^{k} + 2^{k} + ... + n^{k}$를 $10^{9} + 7$로 나눈 나머지를 구하는 문제.

 

위에서 유도한 $$(n+1)^{k+1} - 1 = \sum_{x = 0}^{k} \binom{k+1}{x} S_{x}(n)$$을 이용하면 문제를 풀 수 있을 것 같다.

 

n은 일정하니까, x = 0,1,2,...,k에 대하여 $S_{x}(n)$을 다이나믹 프로그래밍을 활용하여 차례대로 구한다.

 

i번째 $S_{i}(n)$을 구해놓는다면, 점화식을 이용해서 $S_{i+1}(n)$을 구할 수 있을 것 같다.

 

결국에 문제는 이항계수를 구하는게 문제인데, k!까지 모든 팩토리얼을 계산해놓고, $(k!)^{-1}$까지 구해놓는다.

 

def mod_factorial(p,mod):
    
    factorial = [0]*(p+2)
    factorial[0] = 1

    for i in range(1,p+2):
    
        factorial[i] = factorial[i-1]*i
        factorial[i] %= mod
    
    return factorial

 

 

팩토리얼 역원을 구하는 방법이 정리가 안되어있네..? (나중에 정리할 것 같지만.. 여기선 간단하게)

 

매우 큰 소수 $p = 10^{9} + 7$에 대하여, n이 최대 $10^{9}$까지니까, 모든 n!은 p와 서로소이다.

 

그래서 p에 대한 모듈로 곱셈의 역원이 존재하는데 페르마의 소정리로부터,

 

$$(n!)^{p-1} ≡ 1 (mod p)$$

 

n!을 빼보면, $n! * (n!)^{p-2} ≡ 1 (mod p)$

 

그러면, $n!^{-1} ≡ (n!)^{p-2} (mod p)$

 

그런데, 식을 살짝만 더 바꿔보면 다음과 같은 점화식을 얻는다.

 

$$(n-1)! * (n * (n!)^{p-2}) ≡ 1 (mod p)$$

 

즉, (n-1)!의 역원은 n!의 역원에 n을 곱하면 된다.

 

그러므로, inverse[n-1] = inverse[n] * n이라는 점화식을 얻는다.

 

이 점화식을 이용하면, 가장 큰 k에 대하여 inverse[k]는 $(k!)^{p-2} = (factorial[k])^{p-2}$로 구해놓고,

 

다이나믹 프로그래밍을 이용해서 거꾸로 순회해서 구하면 되겠다.

 

#(n)!*(n!)^(p-2) == 1
#(n-1)!*n*(n!)^(p-2) == 1
def mod_inv_factorial(factorial,p,mod):
    
    inverse = [0]*(p+2)

    inverse[p+1] = pow(factorial[p+1],mod-2,mod)

    for i in range(p+1,0,-1):
        
        inverse[i-1] = inverse[i]*i
        inverse[i-1] %= mod
    
    return inverse

 

factorial과 inverse를 모두 구해놓는다면, 이항계수를 소수 p로 나눈 나머지는..

 

$$\binom{n}{r} ≡ factorial[n] * inverse[r] * inverse[n-r] (mod p)$$로 구할 수 있을 것이다.

 

S[x]를 $S_{x}(n)$을 계산해놓은 값을 저장한다고 하자.

 

S[0]는 $1^{0} + 2^{0} + ... + n^{0} = n$이다.

 

이제 i는 1부터 p까지 순회하여, 배열을 채워넣는다.

 

점화식 $(n+1)^{i+1} - 1 = \sum_{x=0}^{i} \binom{i+1}{x} S_{x}(n)$을 이용하자.

 

매 반복문마다 구해야하는 값은 $S_{i}(n)$

 

그러면, $$(n+1)^{i+1} - 1 - \sum_{x=0}^{i-1} \binom{i+1}{x} S_{x}(n) = \binom{i+1}{i} S_{i}(n)$$

 

0~i-1까지 S배열에 $S_{x}(n)$의 값이 이미 계산되어 있고,

 

factorial, inverse 배열을 이용하면, 모든 $\binom{i+1}{x}$를 계산할 수 있다.

 

또한 $(n+1)^{i+1}$은 분할정복을 이용한 거듭제곱으로 빠르게 계산 가능하다.

 

좌변은 모두 계산가능하고, 다음 우변에 이항계수 $\binom{i+1}{i}$가 문제인데, .. 이제 보니까.. 

 

$\binom{i+1}{i} = i+1$이네...??

 

그런데, $i+1 < p \geq 1000$이므로, $10^{9} + 7$과 서로소이다. 

 

그러므로 모듈로 곱셈 역원이 존재하고, 페르마의 소정리에 의해 $(i+1)^{10^{9}+7-2}$로 구할 수 있다.

 

따라서, 좌변을 모두 계산한 다음, $(i+1)^{10^{9}+7-2}$를 곱하면, $S_{i}(n)$을 얻게된다.

 

i = 1부터 p까지 반복하면, 구하고자 하는 값 $S_{p}(n) = 1^{p} + 2^{p} + ... + n^{p}$은 S[p]에 들어있다.

 

#(n+1)^(p+1) - 1 = sigma(t=0)^(p) (p+1Ct)St(n)
#nCr = n!/(r!)*(n-r)!
def sum_of_power(factorial,inverse,n,p):
    
    S = [0]*(p+1)

    S[0] = n % mod

    for i in range(1,p+1):
        
        sigma = 0

        #sum(j = 0)^(i-1) (i+1)Cj*Sj(n)
        for j in range(i+1):
            
            sigma += (S[j] * factorial[i+1]*inverse[j]*inverse[i+1-j])
            sigma %= mod

        #Si(n) = ((n+1)^(i+1)-1 - sigma)*(i+1)Ci^(-1)

        S[i] = pow(n+1,i+1,mod) -1 - sigma
        
        S[i] %= mod

        S[i] *= pow(i+1,mod-2,mod)

        S[i] %= mod
    
    return S[p]

n,p = map(int,input().split())

mod = 10**9+7

factorial = mod_factorial(p,mod)
inverse = mod_inv_factorial(factorial,p,mod)

print(sum_of_power(factorial,inverse,n,p))

 

놀랍게도, 최대 $10^{9}$인 n과는 전혀 무관하게 최대 1000인 p에 대하여 $O(p^{2})$으로 문제를 해결할 수 있다.

 

당연히 합동식 성질을 이용해 mod로 계속 나눠줘서 수를 압축해주는건 필수다.