P4240 毒瘤之神的考验

Lanthanmum

于 2024-10-04 22:50:40 发布

阅读量142

点赞数 4

分类专栏：数论文章标签：算法数论根号分治

本文链接：https://blog.csdn.net/Lanthamum/article/details/142706666

版权

数论专栏收录该内容

1 篇文章 0 订阅

订阅专栏

毒瘤之神的考验 - 洛谷

$\sum_{i=1}^{n}\sum_{j=1}^{m}\varphi (ij) \qquad 1-1$

定义 $\varphi (i)=i*\prod_{p|i}^{}\frac{p-1}{p}$

猜想 $\varphi (ij)$ 与 $\varphi (i)\varphi (j)$ 有关

$\varphi (ij)=\frac{i*\prod_{p|i}^{}\frac{p-1}{p}*j*\prod_{p|j}^{}\frac{p-1}{p}}{\prod_{p|gcd(i,j)}^{}*\frac{p-1}{p}} \qquad 1-2$

$\varphi (i)*\varphi(j)=i*\prod_{p|i}^{}*\frac{p-1}{p}*j*\prod_{p|j}^{}*\frac{p-1}{p} \qquad 1-3$

发现上式1-1 上下两边乘gcd(i,j)有

$\varphi(ij)=\frac{\varphi(i)*\varphi(j)*gcd(i,j)}{\varphi(gcd(i,j))}\qquad 1-4$

带入1-1有

$\sum_{i=1}^{n}\sum_{j=1}^{m}\frac{\varphi(i)*\varphi(j)*gcd(i,j)}{\varphi(gcd(i,j))}\qquad 2-1$

化简 n<m

$\sum_{k=1}^{n}\sum_{i=1}^{n}\sum_{j=1}^{m}\frac{\varphi(i)*\varphi(j)*k*[gcd(i,j)=k]}{\varphi(k)}\qquad 2-2$

$\sum_{k=1}^{n}\sum_{ik=1}^{n}\sum_{jk=1}^{m}\frac{\varphi(ik)*\varphi(jk)*k*[gcd(i,j)=1]}{\varphi(k)}\qquad 2-3$

$\sum_{k=1}^{n}\frac{k}{\varphi(k)}\sum_{i=1}^{n/k}\sum_{j=1}^{m/k}\varphi(ik)*\varphi(jk)*[gcd(i,j)=1]\qquad 2-4$

$\sum_{k=1}^{n}\frac{k}{\varphi(k)}\sum_{e=1}^{n/k}\sum_{ie=1}^{n/k}\sum_{je=1}^{m/k}\varphi(ike)*\varphi(jke)*\mu (e)\qquad 2-5$

$\sum_{k=1}^{n}\frac{k}{\varphi(k)}\sum_{e=1}^{n/k}\sum_{i=1}^{n/ke}\sum_{j=1}^{m/ke}\varphi(ike)*\varphi(jke)*\mu (e)\qquad 2-6$

经典代换T=ke e=T/k

$\sum_{T=1}^{n}\sum_{k|T}^{}\frac{k}{\varphi(k)}*\mu (\frac{T}{k})\sum_{i=1}^{n/T}\sum_{j=1}^{m/T}\varphi(iT)*\varphi(jT)\qquad 2-7$

然后化简不了了

这个时候我们可以把一部分看出一个整体

$f(x)=\sum_{k|x}^{}\frac{k*\mu(\frac{x}{k})}{\varphi(k)}\qquad 2-8$

$g(x,y)=\sum_{i=1}^{y}\varphi(ix) \qquad2-9-1$

分析这两个函数，发现f(x) 可以在 $O(nlnn)$ 下预处理出来

g(x,y)有以下递推式

$g(x,y)=g(x-1,y)+\varphi(xy)\qquad 2-9-2$

因此也可以在 $O(nlnn)$ 下处理出来

于是莫队 $O(nlnn+qsqrt(n)lnn)$

#include <bits/stdc++.h>

#define int long long
#define INF (1ll<<60)
#define eps 1e-6
using namespace std;
typedef long long ll;
typedef long long LL;
typedef pair<int, int> pii;
typedef vector<int> vi;
typedef unsigned long long ull;
typedef vector<vector<int> > vii;
typedef vector<vector<vector<int> > > viii;
typedef vector<ll> vl;
typedef vector<vector<ll> > vll;
typedef vector<double> vd;
typedef vector<vector<double> > vdd;
#define time mt19937_64 rnd(chrono::steady_clock::now().time_since_epoch().count());//稳定随机卡牛魔 ull
const int N = 1e5 + 9;
const int mod = 998244353;
int phi[N], mu[N], vis[N], p[N], cnt;
ll t, ans[N], res;
ll f[N];
int inv[N];
vi g[N], d[N];
int L[N],R[N],pos[N];

struct Q {
    int n, m, id;

    friend bool operator<(Q x, Q y) {
        return x.n / t == y.n / t ? ((x.n / t) & 1) ? y.m < x.m : x.m < y.m : x.n < y.n;
    }
} que[N];

void init(int n) {
    //mu,phi
    mu[1] = 1, phi[1] = 1;
    for (int i = 2; i <= n; i++) {
        if (!vis[i]) p[++cnt] = i, mu[i] = -1, phi[i] = i - 1;
        for (int j = 1; j <= cnt && i * p[j] <= n; j++) {
            int t = i * p[j];
            vis[t] = 1;
            if (i % p[j]) mu[t] = -mu[i], phi[t] = phi[i] * (p[j] - 1);
            else {
                phi[t] = phi[i] * p[j];
                break;
            }
        }
    }
    //inv
    inv[1] = 1;
    for (ll i = 2; i <= n; i++) {
        inv[i] = inv[mod % i] * (mod - mod / i) % mod;
    }
    //g(x,y)
    for (ll i = 1; i <= n; i++) {
        g[i].resize(n / i + 2);
        for (ll j = 1; j <= n / i; j++) {
            //前缀和推出
            g[i][j] = (g[i][j - 1] + phi[i * j]) % mod;
        }
    }
    //f(x)
    for (ll i = 1; i <= n; i++) {
        for (ll j = i; j <= n; j += i) {
            f[j] = f[j] % mod + i * inv[phi[i]] % mod * mu[j / i] % mod;
            d[j].push_back(i);//j的约数
        }
    }
}

void update(ll a, ll b, ll w) {
    for (ll i = 0; i < d[a].size(); ++i) {
        res = (res + 1ll * w % mod * phi[a] % mod * g[d[a][i]][b / d[a][i]] % mod * f[d[a][i]] % mod) % mod;
    }
    res = (res + mod) % mod;
}

signed main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr), cout.tie(nullptr);
    init(N - 9);
    int q;
    cin >> q;
    t = sqrt(q);
    for (int i = 1; i <= q; i++) {
        int n, m;
        cin >> n >> m;
        if (n > m) {
            swap(n, m);
        }
        que[i] = {n, m, i};
    }
    sort(que + 1, que + 1 + q);
    ll l = 0, r = 0;
    for (int i = 1; i <= q; i++) {
        while (l < que[i].n) update(++l, r, 1);
        while (r < que[i].m) update(++r, l, 1);
        while (l > que[i].n) update(l--, r, -1);
        while (r > que[i].m) update(r--, l, -1);
        ans[que[i].id] = res % mod;
    }
    for (int i = 1; i <= q; i++) {
        cout << ans[i] << '\n';
    }
    return 0;
}

接着推

把两个函数带入2-7可以得到

$\sum_{T=1}^{n}f(T)*g(T,n/T)*g(T,m/T) \qquad 3-1$

设整体一个函数有

$h(n,m,T)=\sum_{k=1}^{T}f(k)*g(k,n/k)*g(k,m/k) \qquad 3-1-1$

整除分块得到部分答案

$h(n,m,r)-h(n,m,l-1)=ansi\qquad 3-1-2$

$O(n^{2}lnn)$ 前缀和处理出答案

$h(n,m,x)=h(n,m,x-1)+f(x)*g(x,n/x)*g(x,m/x) \,3-2$

现在我们有两种方案

1. $O(nlnn)$ 预处理， $O(n)$ 回答 TLE

2. $O(n^{2}lnn)$ 预处理， $O(sqrt(n))$ 回答 MLE+TLE

因此用根号分治来平衡时间

设定一个阈值 k

1.<=t预处理暴力 $O(nlnn+\frac{tn}{k})$

2.>t 处理过的整块算整除分块 $O(nklnk+tsqrt(n))$

最优的k $(klnk+\frac{t}{k})min$ t=1e4 k=47

k=47-500 都能过还#define int long long 了！

#include <bits/stdc++.h>
#define int long long
#define INF (1ll<<60)
#define eps 1e-6
using namespace std;
typedef long long ll;
typedef long long LL;
typedef pair<int, int> pii;
typedef vector<int> vi;
typedef unsigned long long ull;
typedef vector<vector<int> > vii;
typedef vector<vector<vector<int> > > viii;
typedef vector<ll> vl;
typedef vector<vector<ll> > vll;
typedef vector<double> vd;
typedef vector<vector<double> > vdd;
#define time mt19937_64 rnd(chrono::steady_clock::now().time_since_epoch().count());//稳定随机卡牛魔 ull
const int N = 1e5 + 9;
const int mod = 998244353;
const int B = 47;
int phi[N], mu[N], vis[N], p[N], cnt;
int f[N], inv[N];
vi g[N], h[B + 5][B + 5];

void add(int &x, int y) { x = x + y >= mod ? x + y - mod : x + y; }

void init(int n) {
    //mu,phi
    mu[1] = 1, phi[1] = 1;
    for (int i = 2; i <= n; i++) {
        if (!vis[i]) p[++cnt] = i, mu[i] = -1, phi[i] = i - 1;
        for (int j = 1; j <= cnt && i * p[j] <= n; j++) {
            int t = i * p[j];
            vis[t] = 1;
            if (i % p[j]) mu[t] = -mu[i], phi[t] = phi[i] * (p[j] - 1);
            else {
                phi[t] = phi[i] * p[j];
                break;
            }
        }
    }
    //inv
    inv[1] = 1;
    for (ll i = 2; i <= n; i++) {
        inv[i] = inv[mod % i] * (mod - mod / i) % mod;
    }
    //g(x,y)
    for (ll i = 1; i <= n; i++) {
        g[i].resize(n / i + 5);
        for (ll j = 1; j <= n / i; j++) {
            //前缀和推出
            g[i][j] = (g[i][j - 1] + phi[i * j]) % mod;
        }
    }
    //f(x)
    for (ll i = 1; i <= n; i++) {
        for (ll j = i; j <= n; j += i) {
            f[j] = f[j] % mod + i * inv[phi[i]] % mod * mu[j / i] % mod;
        }
    }
    //h(n,m,t)
    for (int i = 1; i <= B; i++) {
        for (int j = i; j <= B; j++) {
            h[i][j].resize(n / j + 5);
            for (int t = 1; t <= n / j; t++) {
                h[i][j][t] = (h[i][j][t - 1] + 1ll * f[t] * g[t][i] % mod * g[t][j] % mod) % mod;
            }
        }
    }

}


signed main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr), cout.tie(nullptr);
    init(N - 9);
    auto solve = [&](int n, int m) {
        if (n > m) {
            swap(n, m);
        }
        int res = 0;
        //根号分治
        for (int i = 1; i <= m / B; i++) {
            add(res, 1ll * f[i] * g[i][n / i] % mod * g[i][m / i] % mod);
        }
        for (int l = m / B + 1, r; l <= n; l = r + 1) {
            r = min(n / (n / l), m / (m / l));
            add(res, (h[n / l][m / l][r] - h[n / l][m / l][l - 1] + mod) % mod);
        }
        return res;
    };
    int q;
    cin >> q;
    while (q--) {
        int n, m;
        cin >> n >> m;
        cout << solve(n, m) << '\n';
    }
    return 0;
}