Description

洛谷 P6156 簡單題

Solution

題意非常清晰明瞭。

首先，我們不難發現，\(f(x) = \mu^2(x)\)。

然後就是一波基礎而不失難度的推式子。

\[\begin{aligned} & \sum\limits_{i = 1}^n\sum_{j = 1}^n(i + j)^k\mu^2(\gcd(i, j))\gcd(i, j) \\ =& \sum_{i = 1}^n\sum_{j = 1}^n(i + j)^k\sum_{d = 1}^n\mu^2(d)d[\gcd(i, j) = d] \\ =& \sum_{d = 1}^nd\mu^2(d)\sum_{i = 1}^{\lfloor \frac{n}{d} \rfloor}\sum_{j = 1}^{\lfloor \frac{n}{d} \rfloor}(id + jd)^k[\gcd(i, j) = 1] \\ =& \sum_{d = 1}^nd^{k + 1}\mu^2(d)\sum_{i = 1}^{\lfloor \frac{n}{d} \rfloor}\sum_{j = 1}^{\lfloor \frac{n}{d} \rfloor}(i + j)^k[\gcd(i, j) = 1] \\ =& \sum_{d = 1}^nd^{k + 1}\mu^2(d)\sum_{i = 1}^{\lfloor \frac{n}{d} \rfloor}\sum_{j = 1}^{\lfloor \frac{n}{d} \rfloor}\sum_{p \mid (i, j)}\mu(p)(i + j)^k \\ =& \sum_{d = 1}^nd^{k + 1}\mu^2(d)\sum_{p = 1}^{\lfloor \frac nd \rfloor}\mu(p) \sum_{i = 1}^{\lfloor \frac{n}{dp} \rfloor}\sum_{j = 1}^{\lfloor \frac{n}{dp} \rfloor}\sum_{p \mid (i, j)}(ip + jp)^k \\ =& \sum_{d = 1}^nd^{k + 1}\mu^2(d)\sum_{p = 1}^{\lfloor \frac nd \rfloor}\mu(p)p^k \sum_{i = 1}^{\lfloor \frac{n}{dp} \rfloor}\sum_{j = 1}^{\lfloor \frac{n}{dp} \rfloor}\sum_{p \mid (i, j)}(i + j)^k \\ \end{aligned} \]

經典令 \(T = dp\)

，則：

\[\begin{aligned} & \sum_{T = 1}^n\sum_{i = 1}^{\lfloor \frac{n}{T} \rfloor}\sum_{j = 1}^{\lfloor \frac{n}{T} \rfloor}(i + j)^k\sum_{d | T}d^{k + 1}\mu^2(d)\mu(\frac Td)(\frac Td) ^k \\ =& \sum_{T = 1}^nT^k\sum_{i = 1}^{\lfloor \frac{n}{T} \rfloor}\sum_{j = 1}^{\lfloor \frac{n}{T} \rfloor}(i + j)^k\sum_{d \mid T}d\mu^2(d)\mu(\frac Td) \\ \end{aligned} \]

我們再令

\[S(n) = \sum\limits_{i = 1}^n\sum\limits_{j = 1}^n(i + j)^k \\ f(n) = \sum\limits_{d \mid n}d\mu^2(d)\mu(\frac nd) \]

那麼原式為：

\[\sum_{T = 1}^nT^kS(\frac nT)f(T) \]

我們只要能夠快速求出 \(S(n)\) 和 \(f(n)\) 即可。

---

先來看 \(S(n)\)。

設 \(F(n) = \sum\limits_{i = 1}^ni^k\)，\(G(n) = \sum\limits_{i = 1}^nF(i)\)。

結論：

\[S(n) = G(2n) - 2G(n) \]

證明：

考慮使用數學歸納法。

\(S(n) = G(2n) - 2G(n)\)，那麼有 \(S(n + 1) = G(2n + 2) - 2G(n + 1)\)。

\[\begin{aligned} S(n + 1) = & \sum_{i = 1}^{n + 1}\sum_{j = 1}^{n + 1}(i + j)^k\\ =& S(n) + 2\sum_{i = 1}^n(i + n + 1)^k + (2n + 2)^k \\ =& S(n) + 2F(2n + 1) - 2F(n + 1) + F(2n + 2) - F(2n + 1) \\ =& G(2n) - 2G(n) + F(2n + 1) + F(2n + 2) - 2F(n + 1) \\ =& \sum_{i = 1}^{2n}F(i) - 2\sum_{i = 1}^nF(i) + F(2n + 1) + F(2n + 2) - 2F(n + 1) \\ =& G(2n + 2) - 2G(n + 1) \end{aligned} \]

證畢.

所以線性篩時篩出 \(i^k\) 求兩遍字首和得到 \(G(n)\)，就可以 \(O(1)\) 求出 \(S(n)\) 了。

---

接下來看 \(f(n)\) 怎麼求，\(f(n) = \sum\limits_{d \mid n}d\mu^2(d)\mu(\frac nd)\)。

我們發現 \(f(n)\) 是由好幾個積性函式乘起來的，所以它也是積性函式。

對於質數 \(p\) ：\(f(p) = \mu^2(1)\mu(p) \times p\mu^2(p)\mu(1) = p - 1\)

回憶一下線性篩的過程，我們對一個數 \(i\) 配上一個質數 \(p\)。

下面我們進行分類討論：

• \(p \nmid i\) ：\(f(ip) = f(i) \times f(p)\)

• \(p \mid i\) ：

  • \(p^2 \mid i\) 且 \(p^3 \nmid i\) ：

    \(f(p^2) = \mu^2(1)\mu(p^2) \times p\mu^2(p)\mu(p) \times p^2\mu^2(p^2)\mu(1) = -p\)

    那麼 \(f(ip) = f(p^2) \times f(i) = (-p)f(i)\)。

  • \(p^k \mid i \ \ (k \geq 3)\) ：

    此時，每一項中的 \(d\) 和 \(\frac nd\) 中必有一個有二次項，使得 \(\mu^2(d)\) 和 \(\mu(\frac nd)\) 中必有一個為 0，所以 \(f(p) = 0\)。

附上線性篩的程式碼：

\(Code\)

inline void euler(){
    f[1] = 1;
    for(int i = 2; i < N; ++i){
        if(!vis[i]) p[++tot] = i, f[i] = i - 1;
        for(int j = 1; j <= tot && i * p[j] < N; ++j){
            vis[i * p[j]] = 1;
            if(i % p[j]) f[i * p[j]] = f[i] * f[p[j]] % mod;
            else{
                if((i / p[j]) % p[j]) f[i * p[j]] = f[i / p[j]] * (mod - p[j]) % mod;
                break;
            }
        }
    }
}

---

至此，這道題就可以愉快的解決啦，計算答案的時候整除分塊走一個即可。

我在預處理時把 \(T^k\) 和 \(f(T)\) 乘到了一起，具體為線性篩後面的第一個字首和。

Code

#include <bits/stdc++.h>
#define ll long long

using namespace std;

namespace IO{
    inline ll read(){
        ll x = 0;
        char ch = getchar();
        while(!isdigit(ch)) ch = getchar();
        while(isdigit(ch)) x = (x << 3) + (x << 1) + ch - '0', ch = getchar();
        return x;
    }

    template <typename T> inline void write(T x){
        if(x > 9) write(x / 10);
        putchar(x % 10 + '0');
    }
}
using namespace IO;

const ll N = 1e7 + 10;
const ll mod = 998244353;
ll n, k, tot;
ll p[N], f[N];
ll F[N], G[N];
bool vis[N];

inline ll qpow(ll a, ll b){
    ll res = 1;
    while(b){
        if(b & 1) res = res * a % mod;
        a = a * a % mod, b >>= 1;
    }
    return res;
}

inline void euler(){
    f[1] = F[1] = 1;
    for(int i = 2; i < N; ++i){
        if(!vis[i]) p[++tot] = i, f[i] = i - 1, F[i] = qpow(i, k);
        for(int j = 1; j <= tot && i * p[j] < N; ++j){
            vis[i * p[j]] = 1;
            F[i * p[j]] = F[i] * F[p[j]] % mod;
            if(i % p[j]) f[i * p[j]] = f[i] * f[p[j]] % mod;
            else{
                if((i / p[j]) % p[j]) f[i * p[j]] = f[i / p[j]] * (mod - p[j]) % mod;
                break;
            }
        }
    }
    for(int i = 1; i < N; ++i) f[i] = (f[i - 1] + f[i] * F[i] % mod) % mod, F[i] = (F[i] + F[i - 1]) % mod;
    for(int i = 1; i < N; ++i) F[i] = (F[i] + F[i - 1]) % mod;
}

inline ll S(ll n){
    return (F[n << 1] - (F[n] << 1) % mod + mod) % mod;
}

inline ll solve(ll n, ll k){
    ll res = 0;
    for(int l = 1, r; l <= n; l = r + 1){
        r = n / (n / l);
        res = (res + S(n / l) * (f[r] - f[l - 1] + mod) % mod) % mod;
    }
    return res;
}

signed main(){
    n = read(), k = read();
    euler();
    write(solve(n, k)), puts("");
    return 0;
}

\[\_EOF\_ \]

本文來自部落格園，作者：xixike，轉載請註明原文連結：https://www.cnblogs.com/xixike/p/15699746.html