P6825 「EZEC-4」求和

求：

\[\sum_{i=1}^{n}\sum_{j=1}^{n}\gcd(i,j)^{i+j} \]

先化簡原式，有：

\[\begin{aligned} &\sum_{i=1}^{n}\sum_{j=1}^{n}\gcd(i,j)^{i+j}\\ &=\sum_{d=1}^{n}\sum_{i=1}^{n}\sum_{j=1}^{n}d^{i+j}[\gcd(i,j)=d]\\ &=\sum_{d=1}^{n}\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}d^{d\cdot(i+j)}[\gcd(i,j)=1]\\ &=\sum_{d=1}^{n}\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}d^{d\cdot(i+j)}\sum_{p|\gcd(i,j)}\mu(p)\\ &=\sum_{d=1}^{n}\sum_{p=1}^{\lfloor\frac{n}{d}\rfloor}\mu(p)\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}d^{d\cdot(i+j)}[p|i][p|j]\\ &=\sum_{d=1}^{n}\sum_{p=1}^{\lfloor\frac{n}{d}\rfloor}\mu(p)\sum_{i=1}^{\lfloor\frac{n}{pd}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{pd}\rfloor}d^{dp\cdot(i+j)}\\ &=\sum_{d=1}^{n}\sum_{p=1}^{\lfloor\frac{n}{d}\rfloor}\mu(p)\sum_{i=1}^{\lfloor\frac{n}{T}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{T}\rfloor}d^{T\cdot(i+j)} & T=dp\\ &=\sum_{d=1}^{n}\sum_{p=1}^{\lfloor\frac{n}{d}\rfloor}\mu(p)\sum_{i=1}^{\lfloor\frac{n}{T}\rfloor}(d^T)^i\sum_{j=1}^{\lfloor\frac{n}{T}\rfloor}(d^T)^j\\ &=\sum_{d=1}^{n}\sum_{p=1}^{\lfloor\frac{n}{d}\rfloor}\mu(p)\left(\sum_{i=1}^{\lfloor\frac{n}{T}\rfloor}(d^T)^i\right)^2 \end{aligned} \]

再看後面那部分，令 \(x=\frac{T}{p}^T\)

。

設 \(S_n\) 為等比數列前 \(n\) 項和，則有：

\[S_n=S_{\lfloor\frac{n}{2}\rfloor}(1+a^{\lfloor\frac{n}{2}\rfloor})+a^n[2 \nmid n] \]

即可 \(\mathcal O(\log n)\) 求解。

時間複雜度 \(\mathcal O(n \log n)\)。

#include <bits/stdc++.h>

using namespace std;

#define ll long long

#define pll pair<ll, ll>

#define mpr make_pair

#define fi first

#define se second

const int _ = 1.5e6 + 10;

int mod;

int read()
{
    int x = 0, f = 1;
    char c = getchar();
    while (c < '0' || c > '9')
    {
        if (c == '-')
            f = -1;
        c = getchar();
    }
    while (c >= '0' && c <= '9')
    {
        x = x * 10 + c - '0';
        c = getchar();
    }
    return x * f;
}

int qpow(int x, int y)
{
    int res = 1;
    while (y)
    {
        if (y & 1)
            res = 1ll * res * x % mod;
        x = 1ll * x * x % mod;
        y >>= 1;
    }
    return res;
}

pll qsum(int a, int n)
{
    if (n == 0)
        return mpr(1, 1);
    if (n == 1)
        return mpr(a, a);
    pll tmp = qsum(a, n / 2);
    if (n & 1)
        return mpr((tmp.fi * (tmp.se + 1) % mod + tmp.se * tmp.se % mod * a % mod) % mod, tmp.se * tmp.se % mod * a % mod);
    else
        return mpr(tmp.fi * (tmp.se + 1) % mod, tmp.se * tmp.se % mod);
}

int T, n, m, k;

bitset<_> vis;

int cnt, pri[_], mu[_];

void init(int n)
{
    mu[1] = 1;
    for (int i = 2; i <= n; i++)
    {
        if (!vis[i])
            pri[++cnt] = i, mu[i] = -1;
        for (int j = 1; j <= cnt && i * pri[j] <= n; j++)
        {
            vis[i * pri[j]] = 1;
            if (i % pri[j] == 0)
                break;
            mu[i * pri[j]] = -mu[i];
        }
    }
}

int solve(int n)
{
    int ans = 0;
    for (int d = 1; d <= n; d++)
    {
        int a = qpow(d, d);
        int b = a;
        for (int p = 1; p <= n / d; p++)
        {
            pll tmp = qsum(b, n / (d * p));
            ans = (ans + tmp.fi * tmp.fi % mod * (mu[p] + mod) % mod) % mod;
            b = 1ll * b * a % mod;
        }
    }
    return ans;
}

signed main()
{
    T = read();
    init(_ - 10);
    while (T--)
    {
        n = read(), mod = read();
        printf("%d\n", solve(n));
    }
    return 0;
}
 

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