P3768 简单的数学题 (莫比乌斯反演+杜教筛)

传送门

题意

求解\((\sum\limits^n_{i=1}\sum\limits^n_{j=1}ijgcd(i,j))mod~p\)

解题过程

\(\sum\limits^n_{i=1}\sum\limits^n_{j=1}ijgcd(i,j)\)

枚举gcd的值

\(=\sum\limits^n_{d=1}d\sum\limits^n_{i=1}\sum\limits^n_{j=1}ij[gcd(i,j)=d]\)

\(=\sum\limits^n_{d=1}d^3\sum\limits^{\lfloor \frac{n}{d}\rfloor}_{i=1}\sum\limits^{\lfloor \frac{n}{d}\rfloor}_{j=1}ij[gcd(i,j)=1]\)

先把后面的部分莫比乌斯反演一下

\(=\sum\limits^n_{d=1}d^3\sum\limits^{\lfloor \frac{n}{d}\rfloor}_{i=1}\sum\limits^{\lfloor \frac{n}{d}\rfloor}_{j=1}ij\sum\limits_{t|gcd(i,j)}\mu(t)\)

把\(t\)提前

\(=\sum\limits^n_{d=1}d^3\sum\limits^{\lfloor \frac{n}{d}\rfloor}_{t=1}\mu(t)*t^2\sum\limits^{\lfloor \frac{n}{td}\rfloor}_{i=1}\sum\limits^{\lfloor \frac{n}{td}\rfloor}_{j=1}ij\)

令\(1+2+3+\cdots+n=S(n)\)

那么\(=\sum\limits^n_{d=1}d^3\sum\limits^{\lfloor \frac{n}{d}\rfloor}_{t=1}\mu(t)*t^2\sum\limits^{\lfloor \frac{n}{td}\rfloor}_{i=1}\sum\limits^{\lfloor \frac{n}{td}\rfloor}_{j=1}ij\)相当于

\(=\sum\limits^n_{d=1}d\sum\limits^{\lfloor \frac{n}{d}\rfloor}_{t=1}\mu(t)*t^2d^2S(\lfloor \frac{n}{td}\rfloor)^2\)

令\(T=td\)

\(=\sum\limits^n_{d=1}d\sum\limits^{\lfloor \frac{n}{d}\rfloor}_{t=1}\mu(t)*T^2S(\lfloor \frac{n}{T}\rfloor)^2\)

把T提到最外面

\(=\sum\limits_{T=1}^nT^2S(\lfloor \frac{n}{T}\rfloor)^2\sum\limits_{d|T}\frac{\mu(d)T}{d}\)

由\(\sum\limits_{d|n}\frac{\mu(d)}{d}=\frac{\varphi(n)}{n}\)

得\(=\sum\limits_{T=1}^nT^2S(\lfloor \frac{n}{T}\rfloor)^2\sum\limits_{d|T}\frac{\mu(d)T}{d}=\sum\limits_{T=1}^nT^2S(\lfloor \frac{n}{T}\rfloor)^2\varphi(T)\)

好像还是不太行

我们令\(f(i)=i^2\varphi(i)\)

注意到\(\sum\limits_{d|n}\varphi(d)=n\)

那么设一个\(g(i)=i^2\)

那么再看一下杜教筛的套路式子

\(g(1)S(n)=\sum\limits_{i=1}^{n}(f*g)(i)-\sum\limits_{d=2}^{n}g(d)\cdot S(\lfloor\frac{n}{d}\rfloor)=\sum\limits_{i=1}^{n}i^3-\sum\limits_{d=2}^{n}g(d)\cdot S(\lfloor\frac{n}{d}\rfloor)\)

这里的\(S\)是\(f\)的前缀和,为了方别区分将原来的\(1+2+3+\cdots+n\)称为\(Sum(n)\)

\(S(n)=\sum\limits_{i=1}^{n}i^3-\sum\limits_{d=2}^{n}d^2\cdot S(\lfloor\frac{n}{d}\rfloor)\)

其实\(\sum\limits_{i=1}^{n}i^3=Sum(i)^2\)

最后搞出来的是\(\sum\limits_{T=1}^nSum(\lfloor \frac{n}{T}\rfloor)^2f(T)\)

前半部分能数论分块，后半部分能杜教筛，那么应该差不多了

代码

#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/hash_policy.hpp>
#include <ext/pb_ds/tree_policy.hpp>
#include <ext/pb_ds/trie_policy.hpp>
using namespace __gnu_pbds;
using namespace std;
// freopen("k.in", "r", stdin);
// freopen("k.out", "w", stdout);
// clock_t c1 = clock();
// std::cerr << "Time:" << clock() - c1 <<"ms" << std::endl;
//#pragma comment(linker, "/STACK:1024000000,1024000000")
mt19937 rnd(time(NULL));
#define de(a) cout << #a << " = " << a << endl
#define rep(i, a, n) for (int i = a; i <= n; i++)
#define per(i, a, n) for (int i = n; i >= a; i--)
#define ls ((x) << 1)
#define rs ((x) << 1 | 1)
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int, int> PII;
typedef pair<double, double> PDD;
typedef pair<char, char> PCC;
typedef pair<ll, ll> PLL;
typedef vector<int> VI;
#define inf 0x3f3f3f3f
const ll INF = 0x3f3f3f3f3f3f3f3f;
const ll MAXN = 1e7 + 7;
const ll MAXM = 4e5 + 7;
const int MOD = 1e9 + 7;
const double eps = 1e-7;
ll p, n;
gp_hash_table<ll, ll> mp;
int phi[MAXN], pri[MAXN], tot = 0;
int vis[MAXN];
ll sum[MAXN];
void init()
{
    phi[1] = 1;
    for (int i = 2; i < MAXN; i++)
    {
        if (!vis[i])
            pri[++tot] = i, phi[i] = i - 1;
        for (int j = 1; j <= tot && pri[j] * i < MAXN; j++)
        {
            vis[i * pri[j]] = 1;
            if (i % pri[j] == 0)
            {
                phi[i * pri[j]] = phi[i] * pri[j];
                break;
            }
            phi[i * pri[j]] = phi[i] * (pri[j] - 1);
        }
    }
    sum[0] = 0;
    for (int i = 1; i < MAXN; i++)
        sum[i] = (sum[i - 1] + 1LL * phi[i] * i % p * i % p) % p;
}
void add(ll &x, ll &y) { x = x + y >= p ? x + y - p : x + y; }
void del(ll &x, ll &y) { x = x - y < 0 ? x - y + p : x - y; }
ll Inv(ll a) { return a == 1 ? 1 : (p - p / a) * Inv(p % a) % p; }
ll inv6;
ll calSum(ll x)
{
    x %= p;
    ll temp = x * (x + 1) / 2 % p;
    temp %= p;
    return temp * temp % p;
}
ll cal_sumg(ll x)
{
    x %= p;
    return x * (x << 1 | 1) % p * (x + 1) % p * inv6 % p;
}
ll cal_S(ll x)
{
    if (x < MAXN)
        return sum[x];
    if (mp.find(x) != mp.end())
        return mp[x];
    ll ans = calSum(x);
    for (ll l = 2, r = 0; l <= x; l = r + 1)
    {
        r = x / (x / l);
        ans -= (cal_sumg(r) - cal_sumg(l - 1)) % p * cal_S(x / l) % p;
        ans = ((ans % p) + p) % p;
    }
    return mp[x] = ans;
}
int main()
{
    mp.clear();
    scanf("%lld%lld", &p, &n);
    init();
    inv6 = Inv(6);
    ll ans = 0;
    for (ll l = 1, r = 0; l <= n; l = r + 1)
    {
        r = n / (n / l);
        ll s = calSum(n / l);
        ans += (cal_S(r) - cal_S(l - 1)) % p * s % p;
        ans = ((ans % p) + p) % p;
    }
    printf("%lld\n", ans);
    return 0;
}