Problem 487 : Sums of power sums (40%)

Description

original

Let $f_k(n)$ be the sum of the kth powers of the first n positive integers.

For example, $f_2(10) = 1^2 + 2^2 + 3^2 + 4^2 + 5^2 + 6^2 + 7^2 + 8^2 + 9^2 + 10^2 = 385$.

Let $S_k(n)$ be the sum of $f_k(i)$ for $1 \le i \le n$. For example, $S_4(100) = 35375333830$.

What is $\sum (S_{10000}(10^{12}) modp)$ over all primes $p$ between 2 ⋅ 109 and 2 ⋅ 109 + 2000?

간략해석

$f_k(n)$는 1부터 n까지 k승한 것의 합이고, $S_k(n)$은 $f_k(1)$부터 $f_k(n)$까지의 합이다. 2 ⋅ 10^9 와 2 ⋅ 10^9 + 2000사이 모든 소수 $p$에 대하여 다음을 구하시오. $\sum (S_{10000}(10^{12}) mod p)$

Idea & Algorithm

naive idea

이 문제를 풀기 위해서는 문제를 나누는 것이 중요하다.

범위 내 소수 판별
각 소수 별 함수값 계산

1의 경우에는 $O(scope \cdot \sqrt N )$이니 약 2000만에서 3000만번 연산으로 구할 수 있다. 또한 실제 소수의 개수를 세보았을 때, 범위 내 소수는 100개가 나온다.

그렇다면 이 문제를 어느정도 합리적 시간내로 풀기 위해서는 한 소수당 5초 이하의 시간을 써야 10분이내로 정답이 나올 것이다.

2의 경우에 naive하게 하면 $f_k(n)$의 시간복잡도는 $O(kn)$이고, $S_k(n)$은 시간복잡도가 $n * O(f_k(n)) \to O(kn^2)$이 되고 사실상 죽을때까지 못푼다.

그럼 식을 전개하면서 조금 더 아이디어를 생각해보자.

\[\begin{align}S_k(n) & = \sum_{i=1}^n f_k(i) \\ & = \sum_{i=1}^{n} \sum_{j=1}^{i} j^k \\ & = \sum_{i=1}^n (n+1-i) \cdot i^k \\ & = (n+1)\cdot \sum_{i=1}^{n} i^k - \sum_{i=1}^n i^{k+1} \\ & = (n+1)\cdot f_k(n) - f_{k+1}(n) \end{align}\]

다음 전개에 따라서 시간복잡도를 $O(kn)$으로 줄일 수 있지만, 여전히 부족하다. 시간복잡도에서 n연산을 log수준으로 줄이거나 아니면 n에 상관받지 않고 연산을 하느냐가 중요하다.

advanced idea

다음 식 전개에 따라서 구해야할 값은 $f_k(n)$을 빠른 시간내에 구하면 된다.

다음 식을 보자. 이 식은 이항 계수를 통한 전개이다. (Pascal’s Identity) 식에 대해 궁금하다면 아래의 파일을 참조하자.

\[(n+1)^{k+1} = _{k+1}C_{1}\cdot f_k(n) + _{k+1}C_{2}\cdot f_{k-1}(n) ... _{k+1}C_{k+1}\cdot f_0(n) + 1\] \[f_k(n) = \left\{ (n+1)^{k+1} - (_{k+1}C_{2}\cdot f_{k-1}(n) ... _{k+1}C_{k+1}\cdot f_0(n) + 1)\right\}\cdot inv(_{k+1}C_{1})\]

즉, 점화식과 이항계수 전처리에 따라 시간복잡도 $O(k^2)$으로 구할 수 있다.

이 문제보다 쉬운 버전으로 BOJ 1492 : 합도 있다.

source code

100개의 소수 * (약 2억번정도의 연산 + 많은 mod 연산 )에 의하여 3분정도 돌리면 정답이 나온다.

구현시 주의점은 MOD값이 커서, MOD^2이 LLINT MAX/2 정도 된다. 그렇기에 MOD 연산을 적절히 하여 overflow가 안나게 주의해야한다.

#include <iostream>
using namespace std;
typedef long long ll;
typedef pair<int, int> pii;

ll n,m,d[10005][10005],s[10005];
ll MOD;

void init(){
    d[0][0]=1;
    for(int i = 1 ; i <= m+2 ; i++){
        d[i][0] = d[i][i] = 1;
        for(int j = 1 ; j < i;j++){
            d[i][j] = (d[i-1][j]+d[i-1][j-1]) % MOD;
        }
    }
}

ll pw(ll x,ll y){
    ll ret = 1;
    while(y){
        if(y&1)ret=(ret*x)%MOD;
        x=(x*x)%MOD;
        y>>=1;
    }
    return ret;
}

ll inv(ll x){ return pw(x,MOD-2);}

ll ans = 0;

void process(ll mod){
    n = 1e12; m = 1e4;
    n %= mod;
    MOD = mod;
    init();
    s[0] = n;
    for(int i = 1 ; i <= m+1 ; i++){
        s[i] = (pw(n+1,i+1)-1+MOD)%MOD;
        for(int j = 2 ; j <= i+1 ; j++){
            s[i] = (s[i]+MOD-(d[i+1][j] * s[i-j+1])%MOD)%MOD;
        }
        s[i]=(s[i] * inv(d[i+1][1]))%MOD;

    }
    ans += (((n+1)%MOD) * s[m]%MOD-s[m+1]+MOD) % MOD ;
}

bool isPrime(ll n){
    for(ll i = 2 ; i*i <= n ; i++){
        if(n%i==0) return false;
    }
    return true;
}

int main(){
    ll N = 2e9;
    for(ll i = 1 ; i < 2000 ; i+=2){
        if(isPrime(N+i)){
            process(N+i);
        }
    }
    cout << ans;
}

Reference

점화식을 뭐라부르는 지 몰랐는데, 도움을 받아 자료를 찾았습니다.

Thanks to @rkm0959 : https://arxiv.org/pdf/1011.0076.pdf

Problem 487