洛谷 P3768 简单的数学题（莫比乌斯反演+杜教筛）

传送门

---

题目大意

求${\sum }_{i=1}^n{\sum }_{j=1}^{n}i\ast j\ast gcd(i,j)$

莫比乌斯反演

在莫比乌斯反演那里学到了这个式子：$n={\sum }_{d|n}\phi (d)$

那么将其代入上式可以得到：${\sum }_{i=1}^n{\sum }_{j=1}^{n}i\ast j{\sum }_{d|gcd(i,j)}\phi (d)$

那么我们可以枚举$d\in [1,n]$，然后设$i=x\ast d,j=y\ast d$，枚举$x,y$并替换$i$，我们可以得到：

${\sum }_{d=1}^n\phi (d)\ast d^2{\sum }_{x=1}^{\lfloor \frac{n}{d}\rfloor }x{\sum }_{y=1}^{\lfloor \frac{n}{d}\rfloor }y$

令$g(x)=\frac{x(x+1)}{2}$，那么上式等价于：

${\sum }_{d=1}^n\phi (d)\ast d^2[g(\lfloor \frac{n}{d}\rfloor ){]}^2$

如果我们能预处理出${\sum }_{d=1}^n\phi (d)\ast d^2$，但是$n\le 1e^{10}$，显然无法线性筛，此时就需要杜教筛了

杜教筛

对于${\sum }_{d=1}^n\phi (d)\ast d^2$，是比较简单的杜教筛，令$f=\phi \ast i{d}_2,g=i{d}_2$，由狄利克雷卷积得$h=f\ast g={\sum }_{d|n}(\phi (d)\ast d^2)\ast (\frac{n}{d}{)}^2$

展开得到${\sum }_{d|n}\phi (d)\ast n^2$，根据莫比乌斯反演的结论$n={\sum }_{d|n}\phi (d)$，代入得$h=n^3$

设$f(n)$的前缀和为$S(n)$，那么根据杜教筛的常用套路：

$S(n)={\sum }_{i=1}^n{i}^3-{\sum }_{i=2}^n\phi (i)\ast i^2$

预处理$n^{\frac{2}{3}}$以内的$f$的前缀和，然后套用这个公式，至于${\sum }_{i=1}^n{i}^3$，显然不难得到：$(\frac{n(n+1)}{2}{)}^2$

代码

一开始线性递推求${\sum }_{i=1}^n{i}^3$，但是不知道为什么样例都不对，调了好久也不知道为什么

//
// Created by Happig on 2020/9/27
//
#include <bits/stdc++.h>
#include <unordered_map>
#include <unordered_set>

using namespace std;
#define fi first
#define se second
#define pb push_back
#define ins insert
#define Vector Point
#define ENDL "\n"
#define lowbit(x) (x&(-x))
#define mkp(x, y) make_pair(x,y)
#define mem(a, x) memset(a,x,sizeof a);
typedef long long ll;
typedef long double ld;
typedef unsigned long long ull;
typedef pair<int, int> pii;
typedef pair<ll, ll> pll;
typedef pair<double, double> pdd;
const double eps = 1e-8;
const double pi = acos(-1.0);
const int inf = 0x3f3f3f3f;
const double dinf = 1e300;
const ll INF = 1e18;
ll Mod = 1e9 + 7;
const int maxn = 5e6 + 10;

vector<int> prime;
bitset<maxn> vis;
int N;
ll sum[maxn], phi[maxn];
unordered_map<ll, ll> mp;
ll inv;

ll qkp(ll x, ll n, ll p) {
    x %= p;
    ll ans = 1;
    while (n) {
        if (n & 1) ans = ans * x % p;
        x = x * x % p;
        n >>= 1;
    }
    return ans;
}

void init(ll n) {
    sum[1] = phi[1] = 1;
    vis.reset(), prime.clear();
    //N = pow(n, 2 / 3);
    inv = qkp(6, Mod - 2, Mod);
    for (int i = 2; i < maxn; i++) {
        if (!vis[i]) {
            prime.push_back(i);
            phi[i] = i - 1;
        }
        for (int j = 0; j < prime.size() && 1LL * i * prime[j] < maxn; j++) {
            vis[i * prime[j]] = 1;
            if (i % prime[j]) {
                phi[i * prime[j]] = phi[i] * phi[prime[j]] % Mod;
            } else {
                phi[i * prime[j]] = phi[i] * prime[j] % Mod;
                break;
            }
        }
        ll x = 1LL * i * i % Mod;
        sum[i] = (sum[i - 1] + phi[i] * x % Mod) % Mod;
        //pre[i] = (pre[i - 1] + x * i % Mod) % Mod;
    }
}

ll g(ll x) {
    x %= Mod;  //注意
    return x * (x + 1) / 2 % Mod;
}

ll f(ll x) {
    x %= Mod;  //注意
    return x * (x + 1) % Mod * (x + x + 1) % Mod * inv % Mod;
}

ll getSum(ll n) {
    if (n < maxn) return sum[n];
    if (mp[n]) return mp[n];
    ll x = g(n);
    ll ans = x * x % Mod;
    for (ll l = 2, r; l <= n; l = r + 1) {
        r = n / (n / l);
        if (r > n) r = n;
        ans = (ans - (f(r) - f(l - 1) + Mod) % Mod * getSum(n / l) % Mod + Mod) % Mod;
    }
    return mp[n] = ans;
}

int main() {
    //freopen("in.txt","r",stdin);
    //freopen("out.txt","w",stdout);
    ios_base::sync_with_stdio(0), cin.tie(0), cout.tie(0);
    ll n;
    cin >> Mod >> n;
    init(n);
    ll ans = 0;
    for (ll l = 1, r; l <= n; l = r + 1) {
        r = n / (n / l);
        if (r > n) r = n;
        ll res = g(n / l) * g(n / l) % Mod;
        ans += (getSum(r) - getSum(l - 1) + Mod) % Mod * res % Mod;
        ans %= Mod;
    }
    cout << ans << ENDL;
    return 0;
}