Luogu3768 简单的数学题

https://www.luogu.com.cn/problem/P3768

欧拉反演/杜教筛

欧拉反演公式：

\[\sum_{d|n} \varphi(d)=n \]

本题：

\[\sum_{i=1}^n \sum_{j=1}^n ij\gcd(i,j) \\ 对\gcd(i,j)进行欧拉反演 \\ \sum_{i=1}^n \sum_{j=1}^n ij\sum_{k|\gcd(i,j)}\varphi(k)\\ =\sum_{i=1}^n \sum_{j=1}^n ij\sum_{k|i,k|j}\varphi(k)\\ =\sum_{k=1}^n \varphi(k)k^2 (\sum_{i=1}^{\lfloor \frac{n}{k}\rfloor}i)^2\\ 进行杜教筛\\ 设f(n)=\varphi(i)i^2 ,g(n)=n^2\\ (f*g)(n)=\sum_{d|n}f(d)g(\frac{n}{d})\\ =\sum_{d|n}\varphi(d)d^2(\frac{n}{d})^2\\ =\sum_{d|n}\varphi(d)n^2\\ =n^3\\ \sum_{i=1}^n (f*g)(i)\\ =\sum_{i=1}^n \sum_{d|i} f(d)*g(\frac{i}{d})\\ =\sum_{i=1}^n \sum_{d|i} f(\frac{i}{d})*g(d)\\ =\sum_{d=1}^n g(d) \sum_{i=1}^{\lfloor \frac{n}{d} \rfloor}f(i)\\ =\sum_{d=1}^n g(d)S(\lfloor \frac{n}{d}\rfloor)\\ 根据公式\sum_{i=1}^n i^3=\frac{n^2(n+1)^2}{4}\\ \sum_{i=1}^n (f*g)(i)\\ =\sum_{i=1}^n i^3\\ =\frac{n^2(n+1)^2}{4}\\ 列出杜教筛核心公式,当d=1时\\ g(1)S(n)=\sum_{i=1}^n (f*g)(i)-\sum_{i=2}^n g(i)*S(\lfloor \frac{n}{i}\rfloor)\\ 代入g,(f*g)\\ S(n)=\frac{n^2(n+1)^2}{4}-\sum_{i=2}^n i^2*S(\lfloor \frac{n}{i}\rfloor)\\ 再观察原式\\ \sum_{k=1}^n \varphi(k)k^2 (\sum_{i=1}^{\lfloor \frac{n}{k}\rfloor}i)^2\\ 在数论分块的同时，设分块区间为[l,r]，那么计算(S(r)-S(l-1))*(\sum_{i=1}^{\lfloor \frac{n}{\lfloor \frac{n}{l} \rfloor }\rfloor}i)^2即可 \]

\(map\)可以跑过去，\(unordered\_map\)略快

\(C++ Code:\)

#include<iostream>
#include<cstdio>
#include<cmath>
#include<map>
#define ll long long
#define N 5000000
using namespace std;
ll inv6,cnt,prime[N+10],p,n,phi[N+10],s[N+10];
bool pri[N+10];
map<ll,ll>a;
ll exgcd(ll a,ll b,ll &x,ll &y)
{
    if (!b)
    {
        x=1;
        y=0;
        return a;
    }
    ll z=exgcd(b,a%b,x,y);
    ll t=x;
    x=y;
    y=t-a/b*y;
    return z;
}
ll inv(ll w)
{
    ll x,y;
    ll c=exgcd(w,p,x,y);
    ll v=p/c;
    x=(x%v+v)%v;
    return x;
}
void Pre()
{
    phi[1]=1;
    for (ll i=2;i<=N;i++)
    {
        if (!pri[i])
        {
            prime[++cnt]=i;
            phi[i]=i-1;
        }
        for (ll j=1;j<=cnt;j++)
        {
            ll g=i*prime[j];
            if (g>N)
                break;
            pri[g]=true;
            if (i%prime[j]==0)
            {
                phi[g]=phi[i]*prime[j];
                break;
            }
            phi[g]=phi[i]*phi[prime[j]];
        }
    }
    for (ll i=1;i<=N;i++)
        phi[i]%=p;
    for (ll i=1;i<=N;i++)
        s[i]=(s[i-1]+phi[i]*i%p*i%p)%p;
}
ll _s(ll x)
{
    x%=p;
    return x*(x+1)%p*(x+x+1)%p*inv6%p;
}
ll S(ll n)
{
    if (n<=N)
        return s[n];
    if (a[n])
        return a[n];
    ll z=n%p;
    ll ans=z*(z+1)/2;
    ans%=p;
    ans=ans*ans%p;
    ll l,r;
    for (ll l=2;l<=n;l=r+1)
    {
        r=n/(n/l);
        ll o=S(n/l);
        ll e=_s(r)-_s(l-1);
        e=(e%p+p)%p;
        ans=(ans-o*e%p)%p;
        ans=(ans%p+p)%p;
    }
    ans=(ans%p+p)%p;
    return a[n]=ans;
}
int main()
{
    scanf("%lld%lld",&p,&n);
    Pre();
    inv6=inv(6);
    ll ans=0;
    ll l,r;
    for (l=1;l<=n;l=r+1)
    {
        r=n/(n/l);
        ll w=n/l;
        w%=p;
        w=w*(w+1)/2;
        w%=p;
        w=w*w%p;
        ll e=S(r),o=S(l-1);
        e=(e-o)%p;
        e=(e%p+p)%p;
        ans=(ans+e*w%p)%p;
    }
    ans=(ans%p+p)%p;
    cout << ans << endl;
    return 0;
}