[BOJ 18918] 피보나치 수의 최대공약수의 합처럼 보이지만... × 25

이 문제는 두 문제를 섞은 것이다. 그들은 각각 [BOJ 13716] 피보나치 수열처럼 보이지만... 과
[BOJ 17372] 피보나치 수의 최대공약수의 합 이다. 이 문제를 도전하기 전에 위 두 문제를 해결하고 오자.
 
앞서 포스팅에서 [BOJ 17372]를 다루었다. 그때와 같은 방식으로 본 문제에서 주어진 식을 단순화 하면
\[\sum_{d=1}^{n}d^{k}F_{d}C(d)\]
이며 $C(d)=2\sum_{i=1}^{n/d}\phi(i)-1$이다. [BOJ 17372]와 근본적인 방법은 똑같다. xudyh's sieve를 이용하여 $C(d)$를 빠르게 계산하고 $C(d)$에 $n/d$항이 존재한다는 특성을 이용해 $O(\sqrt{n})$만에 문제를 해결하는 것이다. 반면, 이를 위해서는 $\sum_{d=1}^{n}d^{k}F_{d}$를 빠르게 계산할 수 있어야한다. 이것을 하기 위해 차분법을 이용한다.
 
이제 사용될 다항식들을 정의한다. 다항식 $P(X)$와 피보나치 수열 $\{F_{i}\}$그리고 $\deg(P)=k$에 대하여
\[F_{P}(X,t)=\sum_{i=0}^{t}P(X+i){t \choose i}(-1)^{i}\]
\[\Delta P(X)=P(X)-P(X+1)\]
\[h_{P}(X,c)=\sum_{i=0}^{k}F_{P}(X,i)F_{c+2i+1}\]
 
Lemma. $F_{P}(X,t) = F_{\Delta P}(X,t-1)$ for $t>0$
(pf)
\[\begin{align}F_{\Delta P}(X,t-1)&=\sum_{i=0}^{t-1}\Delta P(X+i){t-1 \choose i}(-1)^{i}=\sum_{i=0}^{t-1}(P(X+i)-P(X+i+1)) {t-1 \choose i}(-1)^{i} \\ &=\sum_{i=0}^{t-1}P(X+i) {t-1 \choose i}(-1)^{i} -\sum_{i=1}^{t}P(X+i) {t-1 \choose i-1}(-1)^{i-1}=\sum_{i=0}^{t}P(X+i) {t \choose i}(-1)^{i}=F_{P}(X,t)\ \ \   \blacksquare \end{align} \]
 
Theorem. $h_{P}(X+1,c+1)-h_{P}(X,c)=F_{c}P(X)$
(pf) 귀납법을 사용하자. 만약 $P(X)=0$ 즉, 0차식이었다면 자명하다. $P(X)$를 고정해두고 이것보다 낮은 차수에 대해서 정리가 성립한다고 가정하다. 그럼,
\[\begin{align} h_{P}(X,c)&=F_{P}(X,0)F_{c+1}+\sum_{i=1}^{k}F_{\Delta P}(X,i-1)F_{c+2i+1}\\ &=P(X)F_{c+1}+\sum_{i=0}^{k-1}F_{\Delta P}(X,i)F_{c+2i+3}=P(X)F_{c+1}+h_{\Delta P}(X,c+2)\end{align}\]
이 성립한다. 이때 아래가 성립한다. 
\[h_{\Delta P}(X+1,c+1)-h_{P}(X,c)=h_{\Delta P}(X+1,c+3)-h_{\Delta P}(X, c+2)+P(X+1)F_{c+2}-P(X)F_{c+1}\]
정의에 의해 $\deg(\Delta P)<\deg(P)$이기 때문에, 가정으로 인해 
\[h_{\Delta P}(X+1,c+3)-h_{\Delta P}(X, c+2)=F_{c+2}\Delta P(X)=F_{c+2}(P(X)-P(X+1))\]
이며 대입하면 결과를 얻는다. 
\[h_{\Delta P}(X+1,c+1)-h_{P}(X,c)=F_{c+2}(P(X)-P(X+1)) +P(X+1)F_{c+2}-P(X)F_{c+1}=F_{c}P(X) \ \ \  \blacksquare\]
 
특히 $P(X)=X^{k}$그리고 $c = X$라고 두면, $h_{P}(X+1,X+1)-h_{P}(X,X)=F_{X}P(X)=F_{X}X^{k}$이다. $h(X):=h_{P}(X+1,X+1)$이라 정의하면, $h(X+1)-h(1)=\sum_{i=1}^{X}F_{i}i^{k}$가 된다. 이제 $h(X)$를 빠르게 계산해보자.
$h(X)$에 $F_{P}$를 넣고 전개하면 다음과 같다. 
\[\begin{align}h(X)&=\sum_{i=0}^{k}\sum_{j=0}^{i}(X+j)^{k}{i \choose j}(-1)^{j} F_{X+2i+1}=\sum_{j=0}^{k}(X+j)^{k}(-1)^{j}\sum_{i=j}^{k}{i\choose j}F_{X+2i+1} \\ &=\sum_{j=0}^{k}(X+j)^{k}(F_{X-1}\left[(-1)^{j}\sum_{i=j}^{k}F_{2i+1}{i\choose j}\right]+F_{X} \left[(-1)^{j}\sum_{i=j}^{k}F_{2i+2}{i\choose j}\right])  \end{align}\]
유도 과정에서 항등식 $F_{a+b}=F_{a}F_{b+1}+F_{a-1}F_{b}$를 이용하였다. 
$F_{X-1},F_{X}$의 계수를 각각 $C_{j},D_{j}$라고 두면 피보나치 수열의 값과 조합의 값을 미리 전처리 해두었다는 가정 하에 $O(k^{2})$로 $C_{j},D_{j}$들을 전처리 할 수 있다. 그렇게 된다면 우리는 $O(k\log k)$의 시간으로 $h(X)$를 계산할 수 있다. 
 
이를 바탕으로 구현한 코드가 아래에 있다. 
 

# 코드 보기

#include <bits/stdc++.h>
#define MAX 500003
using namespace std;
typedef long long ll;
const int mod = 1e9 + 7;
ll fib[1000001];
ll combi[4001][4001];
ll C[4001], D[4001];
ll k;
ll phi[MAX];
unordered_map<ll, ll> mp_phi;
 
ll fpow(ll n, ll r){
    ll s = 1;
    for(;r;r >>= 1){
        if(r & 1) s = s * n % mod;
        n = n * n % mod;
    }
    return s;
}
 
ll s_phi(int x){
    int i,j;
    if(x < MAX) return phi[x];
    if(mp_phi.find(x) != mp_phi.end()) return mp_phi[x];
    ll ret = x * (x + 1LL) / 2 %mod;
    for(i = 2;i <= x;i = j + 1){
        j = x / (x / i);
        ret = (ret - (j - i + 1) * s_phi(x / i)) % mod;
    }
    return mp_phi[x] = ret;
}
 
void mat_mul(ll A[2][2], ll B[2][2]){
    ll C[2][2];
    C[0][0] = ((A[0][0] * B[0][0]) % mod + (A[0][1] * B[1][0]) % mod) % mod;
    C[0][1] = ((A[0][0] * B[0][1]) % mod + (A[0][1] * B[1][1]) % mod) % mod;
    C[1][0] = ((A[1][0] * B[0][0]) % mod + (A[1][1] * B[1][0]) % mod) % mod;
    C[1][1] = ((A[1][0] * B[0][1]) % mod + (A[1][1] * B[1][1]) % mod) % mod;
    A[0][0] = C[0][0]; A[0][1] = C[0][1]; A[1][0] = C[1][0]; A[1][1] = C[1][1];
}
 
pair<ll, ll> fibo(ll n){
    ll A[2][2] = {1, 1, 1, 0};
    ll B[2][2] = {1, 0, 0, 1};
    n--;
    for(;n;n >>= 1){
        if(n & 1) mat_mul(B, A);
        mat_mul(A, A);
    }
    return {B[0][0], B[1][0]};
}
 
ll h(ll X){
    ll ret = 0;
    pair<ll, ll> fib_XX = fibo(X);
    ll fib_X = fib_XX.first, fib_X_1 = fib_XX.second;
    for(ll i = 0;i <= k;i++){
        ll t = (C[i] * fib_X_1 % mod + D[i] * fib_X % mod) % mod;
        ret = (ret + t * fpow(X + i, k) % mod) % mod;
    }
    return ret;
}
 
void init(){
    ll a = 1, b = 0;
    ll i, j;
    for(i = 1;i <= 1e6;i++){
        fib[i] = (a + b) % mod;
        a = b; b = fib[i];
    }
    
    combi[0][0] = 1;
    for(i = 1; i <= 4000;i++){
        for(j = 0;j <= i;j++){
            if(j == 0){
                combi[i][j] = 1;
                continue;
            }
            combi[i][j] = (combi[i - 1][j - 1] + combi[i - 1][j]) % mod; 
        }
    }
    
    for(i = 0;i <= k;i++){
        for(j = i;j <= k;j++){
            D[i] = (D[i] + fib[2 * j + 2] * combi[j][i] % mod) % mod;
            C[i] = (C[i] + fib[2 * j + 1] * combi[j][i] % mod) % mod;
        }
        if(i % 2) C[i] = mod - C[i];
        if(i % 2) D[i] = mod - D[i];
    }
    for(i = 1;i < MAX;i++){
        phi[i] = (phi[i] + i) % mod;
        for(j = 2;j * i < MAX;j++) phi[j * i] = (phi[j * i] - phi[i]) % mod;
        phi[i] = (phi[i] + phi[i - 1]) % mod;
    }
}
 
int main(){
    ll n, i, j;
    ll pi, fi, ans = 0;
    scanf("%lld %lld", &n, &k);
    init();
    for(i = 1;i <= min((ll)1e6, n);i++){
        pi = 2LL * s_phi(n / i) % mod - 1LL;
        ans = (ans + fib[i] * fpow(i, k) % mod * pi % mod) % mod;
    }
    for(;i <= n;i = j + 1){
        j = n / (n / i);
        fi = ((h(j + 1) - h(i)) % mod + mod) % mod;
        pi = 2LL * s_phi(n / i) % mod - 1LL;
        ans = (ans + fi * pi % mod) % mod;
    }
    printf("%lld", (ans + mod) % mod);
}

# 코드 닫기

 
이거의 업그레이드 버전으로 $\times 25$문제가 있다. 이거는 FFT를 이용해서 $C_{j}, D_{j}$를 빠르게 구해주면 되는데 그 방법을 모르겠다ㅎ
 
#2024-01-24 edited.
이제 안다. 먼저 $C_{j}$항을 생각해보자. 부호를 제외하면 $ \sum_{i=j}^{k}{i\choose j}F_{X+2i+1} $이다. 이 식을 좀 변형시켜보자. 
\[ \sum_{i=j}^{k}{i\choose j}F_{X+2i+1} = \frac{1}{j!}\sum_{i=j}^{k}F_{X+2i+1}\frac{i!}{(i-j)!}=\frac{1}{j!}\sum_{d=0}^{k-j}F_{2d+2+2j}\frac{(d+j)!}{d!}\]
이때 수열 $\{a_{n}\}$, $\{b_{n}\}$, $\{c_{n}\}=\{a_{n}\}*\{b_{n}\}$(convolution) 를 정의하자. 
\[a_{n}=F_{2k+2-2n}(k-n)!\]
\[b_{n}=\frac{1}{n!}\]
그럼,
\[c_{k-j}=\sum_{d=0}^{k-j}a_{k-j-d}b_{d}=\sum_{d=0}^{k-j}F_{2+2j+2d}(j+d)!\frac{1}{d!}\]
이고 이것이 정확히 $C_{j}$ 항의 일부임을 알 수 있다. $\frac{1}{j!}$는 전처리로 쉽게 계산할 수 있다. 수열 $c$는 쉽게 전처리 가능한 두 수열의 convolution이므로 FFT를 이용해 구할 수 있다. 즉, 아래의 다항식의 계수가 $c$가 된다.
\[(\frac{1}{0!}+ \frac{1}{1!}x+ \frac{1}{2!}x^{2}+\cdots)(F_{2k+2}k!+F_{2k}(k-1)!x+F_{2k-2}(k-2)!x^{2}+\cdots)\]
 
이것을 그대로 구현하고 나면 시간이 촉박할 것이다. 적어도 필자는 그랬다. 최대한 최적화를 해주자. 그 방법은 먼저 Barret reduction을 모든 곱셈에 대해서 하는 것이다. 또 $(a+b)%mod$의 사용을 지양하고 $a+b$를 먼저 계산한 뒤 $a+b$가 $mod$보다 크면 $a+b-mod$를 계산해주는 식으로 하자. 이렇게까지 최적화를 했음에도 불구하고 필자는 시간이 부족했다. 
 
다음 방법으로는 직접 계산하는 항의 개수를 늘리는 것이다. 앞선 코드에서 $min(1e6, n)$을 이용해서 $10^6$까지의 수들에 대해서는 직접 그 합을 계산했다. 이때, 이 10^6을 적당히 큰 수로 바꾸는 것이다. 필자는 $14084507$이라는 수를 선택했는데, 그 과정은 이러하다. 일단 가장 큰 테스트케이스는 $n=10^9$, $k=10^5$일 것이다. 우리의 알고리즘에서는 적절한 구간 $(i,j)$를 잡고 이 구간에서 $\lfloor n/i\rfloor$의 값이 일정함으로 한번에 $d^{k}F_{i}$의 부분합을 계산했다. 이때 $h$ 함수를 2번 호출하였고 $h$ 내부에서는 총 $2k$번 연산이 일어난다. 따라서, 구간을 잡아서 한번에 부분합을 계산하는 위와 같은 방법은 구간의 길이가 적어도 $2k$여야 효율적이라고 생각했기에 구간의 길이가 적어도 $2k$가 되는 최소의 $i$를 찾았고 이것이 아까 찾은 수, $14084507$ 이다. 
 
하지만 이것으로는 부족하다. FFT를 돌리는 과정에서 펙토리얼 계수의 다항식은 2번 연산된다. 그 이유는 $C_{i},D_{i}$각각을 계산하기 때문이다. 그럼 같은 다항식이 2번 FFT 되는 것이기 때문에, 이 과정을 1번으로 압축하였다. 
 
IDE에서 코드를 테스트해보고 싶다면 $n=1000000000$,  $k = 98303$을 테스트해보는 것이 효율적이다. $98303$은 10만 아레에서 비트 수가 가장 많은 수로, 이것은 거듭제곱 함수에 치명적이다. 게다가, 거의 모든 거듭제곱은 지수가 $k$이기 때문에 항상 16번(98303의 비트 개수) 이상은 곱셈을 하고 있는 것이다.
 
이것을 하면서 알게 된 사실은; ideone 을 사용하는데 돌릴 때마다 시간이 약간 차이나게 나오는데 거기서 나온 시간의 최솟값이 백준에서 돌렸을 때의 결과가 되는 것 같다. 
 
코드는 아래에 있다. 
 

# 코드 보기

#include <bits/stdc++.h>
#define MAX 500003
#define LIM 100011
#define TOP 14084507
using namespace std;
typedef long long ll;
using base = complex<double>;
const int mod = 1e9 + 7;
ll fib[TOP + 1];
ll C[LIM + 1], D[LIM + 1];
ll fact[LIM + 1];
ll inv_fact[LIM + 1];
ll k;
ll phi[MAX];
unordered_map<ll, ll> mp_phi;
ll power[LIM]; //n^k
ll inv_power[LIM]; //n^(-k)
 
typedef complex<double> base;
void fft(vector<base> &a, bool inv){
    int n = a.size(), j = 0;
    vector<base> roots(n/2);
    for(int i=1; i<n; i++){
        int bit = (n >> 1);
        while(j >= bit){
            j -= bit;
            bit >>= 1;
        }
        j += bit;
        if(i < j) swap(a[i], a[j]);
    }
    double ang = 2 * acos(-1) / n * (inv ? -1 : 1);
    for(int i=0; i<n/2; i++){
        roots[i] = base(cos(ang * i), sin(ang * i));
    }
    for(int i=2; i<=n; i<<=1){
        int step = n / i;
        for(int j=0; j<n; j+=i){
            for(int k=0; k<i/2; k++){
                base u = a[j+k], v = a[j+k+i/2] * roots[step * k];
                a[j+k] = u+v;
                a[j+k+i/2] = u-v;
            }
        }
    }
    if(inv) for(int i=0; i<n; i++) a[i] /= n; // skip for OR convolution.
}
 
void multiply(vector<ll> &v, vector<ll> &w, vector<ll> &s){
    int n = 2; while(n < v.size() + w.size()) n <<= 1;
    vector<base> v1(n), v2(n), v3(n), r1(n), r2(n), r3(n), r4(n);
    for(int i=0; i<v.size(); i++){
        v1[i] = base(v[i] >> 15, v[i] & 32767);
    }
    for(int i=0; i<w.size(); i++){
        v2[i] = base(w[i] >> 15, w[i] & 32767);
    }
    for(int i=0; i<w.size(); i++){
        v3[i] = base(s[i] >> 15, s[i] & 32767);
    }
    fft(v1, 0);
    fft(v2, 0);
    fft(v3, 0);
    for(int i=0; i<n; i++){
        int j = (i ? (n - i) : i);
        base ans1 = (v1[i] + conj(v1[j])) * base(0.5, 0);
        base ans2 = (v1[i] - conj(v1[j])) * base(0, -0.5);
        base ans3 = (v2[i] + conj(v2[j])) * base(0.5, 0);
        base ans4 = (v2[i] - conj(v2[j])) * base(0, -0.5);
        base ans5 = (v3[i] + conj(v3[j])) * base(0.5, 0);
        base ans6 = (v3[i] - conj(v3[j])) * base(0, -0.5);
        r1[i] = (ans1 * ans3) + (ans1 * ans4) * base(0, 1);
        r2[i] = (ans2 * ans3) + (ans2 * ans4) * base(0, 1);
        r3[i] = (ans1 * ans5) + (ans1 * ans6) * base(0, 1);
        r4[i] = (ans2 * ans5) + (ans2 * ans6) * base(0, 1);
    }
    fft(r1, 1);
    fft(r2, 1);
    fft(r3, 1);
    fft(r4, 1);
    s.resize(n); w.resize(n);
    for(int i=0; i<n; i++){
        ll av = (ll)round(r1[i].real());
        ll bv = (ll)round(r1[i].imag()) + (ll)round(r2[i].real());
        ll cv = (ll)round(r2[i].imag());
        av %= mod, bv %= mod, cv %= mod;
        w[i] = (av << 30) + (bv << 15) + cv;
        w[i] %= mod;
        w[i] += mod;
        w[i] %= mod;
    }
    for(int i=0; i<n; i++){
        ll av = (ll)round(r3[i].real());
        ll bv = (ll)round(r3[i].imag()) + (ll)round(r4[i].real());
        ll cv = (ll)round(r4[i].imag());
        av %= mod, bv %= mod, cv %= mod;
        s[i] = (av << 30) + (bv << 15) + cv;
        s[i] %= mod;
        s[i] += mod;
        s[i] %= mod;
    }
}
 
ll fpow(ll n, ll r){
    ll s = 1, x;
    for(;r;r >>= 1){
        if(r & 1){
            x = s * n;
            s = x - (uint64_t)(((__uint128_t)(-1ULL / mod) * x) >> 64) * mod;
        }
        x = n * n;
        n = x - (uint64_t)(((__uint128_t)(-1ULL / mod) * x) >> 64) * mod;
    }
    return s;
}
 
ll s_phi(int x){
    int i,j;
    if(x < MAX) return phi[x];
    if(mp_phi.find(x) != mp_phi.end()) return mp_phi[x];
    ll ret = x * (x + 1LL) / 2;
    ret = ret - (uint64_t)(((__uint128_t)(-1ULL / mod) * ret) >> 64) * mod;
    for(i = 2;i <= x;i = j + 1){
        j = x / (x / i);
        ret = (ret - (j - i + 1) * s_phi(x / i));
        ret = ret - (uint64_t)(((__uint128_t)(-1ULL / mod) * ret) >> 64) * mod;
    }
    return mp_phi[x] = ret;
}
 
void mat_mul(ll A[2][2], ll B[2][2]){
    ll C[2][2];
    ll x;
    x = A[0][0] * B[0][0];
    C[0][0] = x - (uint64_t)(((__uint128_t)(-1ULL / mod) * x) >> 64) * mod;
    x = A[0][1] * B[1][0];
    C[0][0] += x - (uint64_t)(((__uint128_t)(-1ULL / mod) * x) >> 64) * mod;
    if(C[0][0] >= mod) C[0][0] -= mod;
    
    x = A[0][0] * B[0][1];
    C[0][1] = x - (uint64_t)(((__uint128_t)(-1ULL / mod) * x) >> 64) * mod;
    x = A[0][1] * B[1][1];
    C[0][1] += x - (uint64_t)(((__uint128_t)(-1ULL / mod) * x) >> 64) * mod;
    if(C[0][1] >= mod) C[0][1] -= mod;
    
    x = A[1][0] * B[0][0];
    C[1][0] = x - (uint64_t)(((__uint128_t)(-1ULL / mod) * x) >> 64) * mod;
    x = A[1][1] * B[1][0];
    C[1][0] += x - (uint64_t)(((__uint128_t)(-1ULL / mod) * x) >> 64) * mod;
    if(C[1][0] >= mod) C[1][0] -= mod;
    
    x = A[1][0] * B[0][1];
    C[1][1] = x - (uint64_t)(((__uint128_t)(-1ULL / mod) * x) >> 64) * mod;
    x = A[1][1] * B[1][1];
    C[1][1] += x - (uint64_t)(((__uint128_t)(-1ULL / mod) * x) >> 64) * mod;
    if(C[1][1] >= mod) C[1][1] -= mod;
    
    
    A[0][0] = C[0][0]; A[0][1] = C[0][1]; A[1][0] = C[1][0]; A[1][1] = C[1][1];
}
 
pair<ll, ll> fibo(ll n){
    ll A[2][2] = {1, 1, 1, 0};
    ll B[2][2] = {1, 0, 0, 1};
    n--;
    for(;n;n >>= 1){
        if(n & 1) mat_mul(B, A);
        mat_mul(A, A);
    }
    return {B[0][0], B[1][0]};
}
 
ll h(ll X){
    ll ret = 0, x;
    pair<ll, ll> fib_XX = fibo(X);
    ll fib_X = fib_XX.first, fib_X_1 = fib_XX.second;
    for(ll i = 0;i <= k;i++){
        x = C[i] * fib_X_1;
        ll t = x - (uint64_t)(((__uint128_t)(-1ULL / mod) * x) >> 64) * mod;
        x = D[i] * fib_X;
        t += x - (uint64_t)(((__uint128_t)(-1ULL / mod) * x) >> 64) * mod;
        x = t * fpow(X + i, k);
        ret += x - (uint64_t)(((__uint128_t)(-1ULL / mod) * x) >> 64) * mod;
        if(ret >= mod) ret-= mod;
    }
    return ret;
}
 
void init(){
    ll a = 1, b = 0;
    ll i, j;
    for(i = 1;i <= TOP;i++){
        fib[i] = (a + b);
        if(fib[i] >= mod) fib[i] -= mod;
        a = b; b = fib[i];
    }
    
    ll t = 1;
    fact[0] = 1;
    for(i = 1;i <= LIM;i++){
        t = t * i;
        t = t - (uint64_t)(((__uint128_t)(-1ULL / mod) * t) >> 64) * mod;
        fact[i] = t;
    }
    t = fpow(t, mod - 2);
    inv_fact[LIM] = 1;
    for(i = LIM;i >= 1;i--){
        t = t * i;
        t = t - (uint64_t)(((__uint128_t)(-1ULL / mod) * t) >> 64) * mod;
        inv_fact[i - 1] = t;
    }
    
    vector<ll> r(LIM + 1), d(LIM + 1), c(LIM + 1);
    for(i = 0;i <= k;i++){
        r[i] = inv_fact[i];
        d[i] = fib[2 * k + 2 - 2 * i] * fact[k - i];
        d[i] = d[i] - (uint64_t)(((__uint128_t)(-1ULL / mod) * d[i]) >> 64) * mod;
        c[i] = fib[2 * k + 1 - 2 * i] * fact[k - i];
        c[i] = c[i] - (uint64_t)(((__uint128_t)(-1ULL / mod) * c[i]) >> 64) * mod;
    }
    
    multiply(r, d, c);
    
    for(i = 0;i <= k;i++){
        C[i] = inv_fact[i] * c[k - i];
        C[i] = C[i] - (uint64_t)(((__uint128_t)(-1ULL / mod) * C[i]) >> 64) * mod;
        D[i] = inv_fact[i] * d[k - i];
        D[i] = D[i] - (uint64_t)(((__uint128_t)(-1ULL / mod) * D[i]) >> 64) * mod;
        if(i % 2){
            C[i] = mod - C[i];
            D[i] = mod - D[i];
        }
    }
    
    for(i = 1;i < MAX;i++){
        phi[i] = (phi[i] + i);
        if(phi[i] >= mod) phi[i] -= mod;
        for(j = 2;j * i < MAX;j++) phi[j * i] = (phi[j * i] - phi[i]) % mod;
        phi[i] = (phi[i] + phi[i - 1]);
        if(phi[i] >= mod) phi[i] -= mod;
    }
}
 
int main(){
    ll n, i, j;
    ll pi, fi, ans = 0, x;
    scanf("%lld %lld", &n, &k);
    init();
    for(i = 1;i <= TOP && i <= n;i++){
        pi = 2LL * s_phi(n / i) - 1LL;
        x = fib[i] * fpow(i, k);
        x = x - (uint64_t)(((__uint128_t)(-1ULL / mod) * x) >> 64) * mod;
        x = x * pi;
        ans += x - (uint64_t)(((__uint128_t)(-1ULL / mod) * x) >> 64) * mod;
        if(ans >= mod) ans -= mod;
    }
    for(;i <= n;i = j + 1){
        j = n / (n / i);
        fi = mod - h(i) + h(j + 1);
        if(fi >= mod) fi -= mod;
        pi = 2LL * s_phi(n / i) - 1LL;
        x = fi * pi;
        ans += x - (uint64_t)(((__uint128_t)(-1ULL / mod) * x) >> 64) * mod;
        if(ans >= mod) ans -= mod;
    }
    printf("%lld", (ans + mod) % mod);
}

# 코드 닫기