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

题意

给定n和p(p为质数)求$$\sum _{i=1}^n\sum _{j=1}^{n}ijgcd(i,j)\%p$$
$n\le 10^{10},5\times 10^8\le p\le 1.1\times 10^9$

思路

经典的先提出gcd
$$=\sum _{d=1}^{n}d\sum _{i=1}^n\sum _{j=1}^{n}ij[gcd(i,j)==d]$$
后面的部分除以d，有
$$=\sum _{d=1}^{n}d\sum _{i=1}^{\frac{n}{d}}\sum _{j=1}^{\frac{n}{d}}didj[gcd(i,j)==1]$$
提出后面的d
$$=\sum _{d=1}^n{d}^3\sum _{i=1}^{\frac{n}{d}}\sum _{j=1}^{\frac{n}{d}}ij[gcd(i,j)==1]$$
然后进行莫比乌斯反演，将${\sum }_{d|n}\mu (d)=[n==1]$带入有
$$=\sum _{d=1}^n{d}^3\sum _{i=1}^{\frac{n}{d}}\sum _{j=1}^{\frac{n}{d}}ij\sum _{k|gcd(i,j)}\mu (k)$$
将k提出有,及枚举1到$\frac{n}{d}$内k的倍数
$$=\sum _{d=1}^n{d}^3\sum _{k=1}^{\frac{n}{d}}\mu (k)\sum _{i=1}^{\frac{n}{dk}}\sum _{j=1}^{\frac{n}{dk}}kikj$$
$$=\sum _{d=1}^n{d}^3\sum _{k=1}^{\frac{n}{d}}\mu (k)k^2\sum _{i=1}^{\frac{n}{dk}}\sum _{j=1}^{\frac{n}{dk}}ij$$
上式等价于
$$=\sum _{d=1}^n{d}^3\sum _{dk=1}^n\mu (k)k^2\sum _{i=1}^{\frac{n}{dk}}\sum _{j=1}^{\frac{n}{dk}}ij$$
令$T=dk$有

$$=\sum _{d=1}^n{d}^3\sum _{T=1}^n[T\%d==0]\mu (\frac{T}{d})(\frac{T}{d}{)}^2\sum _{i=1}^{\frac{n}{T}}\sum _{j=1}^{\frac{n}{T}}ij$$
我们将$d^3$乘入
$$=\sum _{d=1}^n\sum _{T=1}^n[T\%d==0]\mu (\frac{T}{d})d{T}^2\sum _{i=1}^{\frac{n}{T}}\sum _{j=1}^{\frac{n}{T}}ij$$
这里可以看作枚举的是1到n中每一个d的倍数，那么我们也等价于枚举1到n中每个数T的因子

$$=\sum _{T=1}^n{T}^2\sum _{d|T}d\mu (\frac{T}{d})\sum _{i=1}^{\frac{n}{T}}\sum _{j=1}^{\frac{n}{T}}ij$$
对于式子$$\sum _{i=1}^{\frac{n}{T}}\sum _{j=1}^{\frac{n}{T}}ij$$
有
$$\sum _{i=1}^{\frac{n}{T}}i\sum _{j=1}^{\frac{n}{T}}j=\sum _{i=1}^{\frac{n}{T}}i(\frac{(\lfloor \frac{n}{T}\rfloor +1)\lfloor \frac{n}{T}\rfloor }{2})=(\frac{(\lfloor \frac{n}{T}\rfloor +1)\lfloor \frac{n}{T}\rfloor }{2}{)}^2$$
对于式子${\sum }_{d|T}d\mu (\frac{T}{d})$可以发现是$id$于$\mu$的狄利克雷卷积$id\ast \mu =\phi$
所以有
$$=\sum _{T=1}^n{T}^2\phi (T)(\frac{(\lfloor \frac{n}{T}\rfloor +1)\lfloor \frac{n}{T}\rfloor }{2}{)}^2$$
然后可以对
$$\sum _{T=1}^n{T}^2\phi (T)$$用杜教筛
令$f(x)=x^2\phi (x)$，$g(x)=x^2$，$S(n)={\sum }_{i=1}^{n}f(i)$

$$\sum _{i=1}^{n}g\ast f=\sum _{i=1}^n\sum _{d|i}g(d)f(\frac{i}{d})=\sum _{i=1}^n\sum _{d|i}{d}^2(\frac{i}{d}{)}^2\phi (\frac{i}{d})$$
$$=\sum _{i=1}^n{i}^2\sum _{d|i}\phi (\frac{i}{d})$$
有${\sum }_{d|n}\phi (d)=n$，故
$$=\sum _{i=1}^n{i}^3$$
再看${\sum }_{i=1}^n{\sum }_{d|i}g(d)f(\frac{i}{d})$
考虑枚举每个因子的倍数
$$\sum _{i=1}^n\sum _{d|i}g(d)f(\frac{i}{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 )$$
$$g(1)S(n)+\sum _{d=2}^{n}g(d)S(\lfloor \frac{n}{d}\rfloor )=\sum _{i=1}^n{i}^3$$
$g(x)=x^2$带入移项得
$$S(n)=\sum _{i=1}^n{i}^3-\sum _{d=2}^n{d}^{2}S(\lfloor \frac{n}{d}\rfloor )$$
有
$$\sum _{i=1}^n{i}^3=(\frac{n(n+1)}{2}{)}^2,\sum _{i=1}^n{i}^2=\frac{n(n+1)(2n+1)}{6}$$
用杜教筛就可以做了
注意n大于int，预处理时可能会爆，注意取模

#include <bits/stdc++.h>
using namespace std;
const int N=1e7+5;
int prime[N];
long long phi[N];
int vis[N];
int cnt;
long long mod;
long long inv2;
long long inv4;
long long inv6;
long long quickmod(long long a,long long b)
{
    long long ans=1;
    while(b)
    {
        if(b%2==1)
            ans=ans*a%mod;
        b=b/2;
        a=a*a%mod;
    }
    return ans;
}
void init()
{
    cnt=0;
    phi[1]=1;
    for(int i=2;i<N;i++)
    {
        if(!vis[i])
        {
            prime[cnt++]=i;
            phi[i]=i-1;
        }
        for(int j=0;j<cnt&&i*prime[j]<N;j++)
        {
            vis[i*prime[j]]=1;
            if(i%prime[j]==0)
            {
                phi[i*prime[j]]=phi[i]*prime[j];
                break;
            }
            else
                phi[i*prime[j]]=phi[i]*(prime[j]-1);
        }
    }
    for(int i=1;i<N;i++)
        phi[i]=(phi[i-1]+1ll*i*i%mod*phi[i]%mod)%mod;
}
map<long long,long long>P;
long long Sum2(long long n)
{
    return n%mod*((n+1)%mod)%mod*(2*n%mod+1)%mod*inv6%mod;
}
long long Sum(long long n)
{
    if(n<N) return phi[n];
    if(P[n]) return P[n];
    long long sum=n%mod*((n+1)%mod)%mod*(n%mod)%mod*((n+1)%mod)%mod*inv4%mod;
    for(long long i=2,last;i<=n;i=last+1)
    {
        last=n/(n/i);
        sum=(sum-(Sum2(last)-Sum2(i-1)+mod)%mod*Sum(n/i)%mod+mod)%mod;
    }
    return P[n]=sum;
}

int main()
{
    long long n;
    scanf("%lld%lld",&mod,&n);
    init();
    inv2=quickmod(2,mod-2);
    inv4=quickmod(4,mod-2);
    inv6=quickmod(6,mod-2);
    long long sum=0;
    for(long long i=1,last;i<=n;i=last+1)
    {
        last=n/(n/i);
        sum=(sum+(Sum(last)-Sum(i-1)+mod)%mod*((n/i+1)%mod)%mod*((n/i)%mod)%mod*((n/i+1)%mod)%mod*((n/i)%mod)%mod*inv4%mod)%mod;
    }
    printf("%lld\n",sum);
    return 0;
}