[洛谷]P3768 简单的数学题-数论

题目地址

题意简述

给定两个整数$n,p$，求下面式子的值：

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

在模$p$意义下的值。

$n\le 10^{10},p\le 1.1\times 10^9$，$p$为质数。

---

• 分析

我们根据套路，来枚举$gcd$，原式变成：

$$\begin{gathered} \sum _{d=1}^n\sum _{i=1}^n\sum _{j=1}^{n}ijd[gcd(i,j)=d] \\ =\sum _{d=1}^{n}d\sum _{i=1}^n\sum _{j=1}^{n}ij[gcd(i,j)=d] \\ =\sum _{d=1}^{n}d\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^{\lfloor \frac{n}{d}\rfloor }id\times jd[gcd(i,j)=1] \\ =\sum _{d=1}^n{d}^3\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^{\lfloor \frac{n}{d}\rfloor }ij[gcd(i,j)=1] \end{gathered}$$

后面我们用莫比乌斯反演的套路，原式变成：

$$\begin{gathered} =\sum _{d=1}^n{d}^3\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^{\lfloor \frac{n}{d}\rfloor }ij\sum _{w|i,w|j}\mu (w) \\ =\sum _{d=1}^n{d}^3\sum _{w=1}^{\lfloor \frac{n}{d}\rfloor }\mu (w)\sum _{i=1}^{\lfloor \frac{n}{dw}\rfloor }iw\sum _{j=1}^{\lfloor \frac{n}{dw}\rfloor }jw \\ =\sum _{d=1}^n{d}^3\sum _{w=1}^{{}^{\lfloor \frac{n}{d}\rfloor }}\mu (w)w^2\sum _{i=1}^{\lfloor \frac{n}{dw}\rfloor }i\sum _{j=1}^{\lfloor \frac{n}{dw}\rfloor }j \end{gathered}$$

然后我们令$S(n)=\frac{n\times (n+1)}{2}$，原式为：

$$=\sum _{d=1}^n{d}^3\sum _{w=1}^{{}^{\lfloor \frac{n}{d}\rfloor }}\mu (w)w^{2}S(\lfloor \frac{n}{dw}\rfloor )S(\lfloor \frac{n}{dw}\rfloor )$$

但是此时$O(n\sqrt{n})$的做法显然不行，所以根据套路，我们转换枚举$dw$，令$T=dw$，原式可以变成：

$$\begin{gathered} =\sum _{T=1}^{n}S{\left(\lfloor \frac{n}{T}\rfloor \right)}^2\sum _{d|T}{d}^3\mu \left(\frac{T}{d}\right){\left(\frac{T}{d}\right)}^2 \\ =\sum _{T=1}^{n}S{\left(\lfloor \frac{n}{T}\rfloor \right)}^2\sum _{d|T}d{T}^2\mu \left(\frac{T}{d}\right) \\ =\sum _{T=1}^{n}S{\left(\lfloor \frac{n}{T}\rfloor \right)}^2{T}^2\sum _{d|T}d\mu \left(\frac{T}{d}\right) \end{gathered}$$

我们可以发现，最后那一部分就是$id⨂\mu =\phi$（$⨂$为狄利克雷卷积，其中因为$\phi ⨂1=id$,所以$\phi ⨂1⨂\mu =id⨂\mu$，而$\mu ⨂1=ϵ(单位元)$，所以$\phi ⨂ϵ=\phi =id⨂\mu$，其中$id(i)=i$）。

所以原式变成了：

$$=\sum _{T=1}^{n}S{\left(\lfloor \frac{n}{T}\rfloor \right)}^2{T}^2\phi (T)$$

接下来，前面的$S$可以数论分块$O(\sqrt{n})的求，而后面需要前缀和，而$n$很大不能线性筛出，所以我们只能用低于线性的杜教筛来求。

我们令$F(n)={\sum }_{i=1}^n{i}^2\phi (i)$，换一种写法就是：
$$F(n)=\sum _{i=1}^{n}i{d}^2\phi (i)$$

由于$i{d}^k⨂i{d}^k(n)=\sum _{d|n}i{d}^k(d)i{d}^k(\frac{n}{d})=n^k{\sum }_{d|n}1$

所以我们令$f(n)=n^2\phi (n),g(n)=n^2$，那么:
$$\begin{gathered} f⨂g(n)=\sum _{d|n}\phi (d)d^2\frac{n^2}{d^2} \\ =n^2\sum _{d|n}\phi (d) \end{gathered}$$

由于${\sum }_{d|n}\phi (d)=n$，所以$f⨂g=i{d}^3$

套入杜教筛公式，我们可以得到：

$$F(n)=\sum _{i=1}^n{i}^3-\sum _{i=2}^n{i}^{2}F(\frac{n}{i})$$

根据公式${\sum }_{i=1}^n{i}^3={\left(\frac{n\times (n+1)}{2}\right)}^2$

${\sum }_{i=1}^n{i}^2=\frac{n\times (n+1)\times (2n+1)}{6}$

所以每次用杜教筛算前缀和即可。

复杂度大概为$O(\sqrt{n}+2n^{\frac{2}{3}})$

下面上代码：

#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
#define ll long long
using namespace std;
const int M=12000001,N=1e4+1;
ll n,P;
ll A[M],B[N],F[M];
ll prime[N],phi[M],cnt,inv_2,inv_6;
bool vis[M];
ll Sqr(ll a){return a*a%P;}
ll fpow(ll a,ll b){
	ll ans=1;
	for(;b;b>>=1,a=(a*a)%P){
		if(b&1)ans=(ans*a)%P;
	}
	return ans;
}
void init(ll w){
    phi[1]=1;
    inv_2=fpow(2,P-2);
    inv_6=fpow(6,P-2);
    for(ll i=2;i<=w;i++){
        if(!vis[i])prime[++cnt]=i,phi[i]=i-1;
        for(ll j=1,v;j<=cnt&&i*prime[j]<=w;j++){
            v=i*prime[j];
            vis[v]=1;
            if(!(i%prime[j])){
                phi[v]=phi[i]*prime[j];
                break;
            }
            phi[v]=phi[i]*phi[prime[j]];
        }
    }
    A[0]=0;
    for(int i=1;i<=w;i++)A[i]=(A[i-1]+Sqr(i)*phi[i]%P)%P;
}
ll &Rec(ll a){
    if(a<M) return A[a];
    else return B[n/a];//不用map，储存空间的小技巧
}
ll S(ll x){if(x>=P)x%=P;return (x*(x+1)%P)*inv_2%P;}
ll S2(ll x){if(x>=P)x%=P;return ((x*(x+1))%P)*((x+x+1)%P)%P*inv_6%P;}
ll SS(ll L,ll R){return (S(R)-S(L-1))%P;}
ll SS2(ll L,ll R){return (S2(R)-S2(L-1))%P;}
ll S3(ll x){return Sqr(S(x));}
ll calc(ll x){
    if(!x) return 0;
    ll &ans=Rec(x);
    if(ans!=-1) return ans;
    ans=S3(x);
    for(ll i=2,j;i<=x;i=j+1){
        j=x/(x/i);
        (ans-=SS2(i,j)*calc(x/i))%=P;
    }
    return (ans+P)%P;
}
ll area(ll L,ll R){
    return ((calc(R)-calc(L-1))%P+P)%P;
}
ll solve(){
    ll ans=0;
    for(ll i=1,j;i<=n;i=j+1){
        j=n/(n/i);
        (ans+=Sqr(S(n/i))*area(i,j)%P)%=P;
    }
    return (ans+P)%P;
}
int main(){
    memset(A,-1,sizeof(A));
    memset(B,-1,sizeof(B));
    scanf("%lld%lld",&P,&n);
    init(min(n,M-1ll));
    printf("%lld\n",solve());
    return 0;
}