题目
题目描述
由于出题人懒得写背景了，题目还是简单一点好。
输入一个整数n和一个整数p，你需要求出$(\sum_{i=1}^n\sum_{j=1}^n ij\gcd(i,j))\mod p$，其中gcd(a,b)表示a与b的最大公约数。
输入输出格式
输入格式：
一行两个整数p、n。
输出格式：
一行一个整数$(\sum_{i=1}^n\sum_{j=1}^n ij\gcd(i,j))\mod p$。
输入输出样例
输入样例#1：
998244353 2000
输出样例#1：
883968974
说明
对于20%的数据，$n \leq 1000$。
对于30%的数据，$n \leq 5000$。
对于60%的数据，$n \leq 10^6$，时限1s。
对于另外20%的数据，$n \leq 10^9$，时限3s。
对于最后20%的数据，$n \leq 10^{10}$，时限6s。
对于100%的数据，$5 \times 10^8 \leq p \leq 1.1 \times 10^9$且p为质数。
题解
不对这种一不小心就会爆long long的题才是最该死的
题目要求：$\sum_{i=1}^n\sum_{j=1}^n ij\gcd(i,j)$
先化简：$\sum_{d=1}^{n}\sum_{i=1}^{n}\sum_{j=1}^{n}[\gcd(i,j)==d]ijd$
可得：$\sum_{d=1}^{n}\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{n}{d}}[\gcd(i,j)==1]ijd^3$
$\sum_{d=1}^{n}d^3\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{n}{d}}[\gcd(i,j)==1]ij$
后面的部分很莫比乌斯$\sum_{d=1}^nd^3\sum_{i=1}^{\frac{n}{d}}\mu(i)i^2(1+2+…\left\lfloor\frac{n}{id}\right\rfloor)^2$
套路啦，设$T=id$，则有$\sum_{T=1}^{n}(1+2+…\left\lfloor\frac{n}{T}\right\rfloor)\sum_{d\vert T}\mu(\frac{T}{d})(\frac{T}{d})^2d^3$
化简：$\sum_{T=1}^{n}(1+2+…\left\lfloor\frac{n}{T}\right\rfloor)T^2\sum_{d\vert T}\mu(\frac{T}{d})d$
因为$\sum_{d\vert T}\mu(\frac{T}{d})d=\phi(T)$
所以式子为$\sum_{T=1}^{n}(1+2+…\left\lfloor\frac{n}{T}\right\rfloor)T^2\phi(T)$
这样可以$O(\sqrt n)$的整除分块做了，但是数据范围觉得你太naive了
尝试用杜教筛搞出$T^2\phi(T)$，设个函数$t(i)=i^2\phi(i)$，显然这是个积性函数
设$S(n)=\sum_{i=1}^{n}t(i)$
则$f(1)S(n)=\sum_{i=1}^{n}(t*f)(i)-\sum_{i=2}^{n}f(i)S(\left\lfloor\frac{n}{i}\right\rfloor)$
$(t*f)(i)=\sum_{d\vert i}t(d)f(\frac{i}{d})=\sum_{d\vert i}\phi(d)d^2f(\frac{i}{d})$
因为$\sum_{d\vert i}\phi(d)=i$，可以设$f(i)=i^2$
所以$\sum_{d\vert i}\phi(d)d^2f(\frac{i}{d})=\sum_{d\vert i}\phi(d)d^2(\frac{i}{d})^2=i^2\sum_{d\vert i}\phi(d)=i^3$
所以$S(n)=\sum_{i=1}^{n}i^3-\sum_{i=1}^{n}i^2S(\left\lfloor\frac{n}{i}\right\rfloor)$
这个就可以杜教筛处理了
然后这题最该死的地方就是到处爆long long，mdzz我查了一晚上TAT
代码
# include<iostream>
# include<cstring>
# include<cstdio>
# include<algorithm>
# define LL long long
using namespace std;
const int MAX=5e6+1;
LL n,mod,tot,inv6,inv2;
int pm[MAX];
LL phi[MAX],ph[MAX];
bool use[MAX];
LL pow(LL b,LL a=6)
{
LL ans=1;
for(;b;b>>=1,a=a*a%mod)
if(b%2) ans=ans*a%mod;
return ans;
}
LL SS(LL x)
{
x%=mod;
return x*(x+1)%mod*(x+x+1)%mod*inv6%mod;
}
LL S(LL x)
{
x%=mod;
return x*(x+1)%mod*inv2%mod;
}
LL dfs(LL x)
{
if(x<MAX) return phi[x];
LL tt=n/x,ttt;
if(ph[tt]) return ph[tt];
ph[tt]=S(x),ph[tt]=ph[tt]*ph[tt]%mod;
for(LL l=2,r;l<=x;l=r+1)
{
r=x/(x/l),ttt=(SS(r)-SS(l-1))%mod;
ph[tt]=(ph[tt]-dfs(x/l)*ttt%mod)%mod;
}
return ph[tt]=(ph[tt]+mod)%mod;
}
int main()
{
scanf("%lld%lld",&mod,&n);
phi[1]=1,inv6=pow(mod-2),inv2=pow(mod-2,2);
for(int i=2;i<MAX;++i)
{
if(!use[i]) pm[++tot]=i,phi[i]=i-1;
for(int j=1;j<=tot&&i*pm[j]<MAX;++j)
{
use[i*pm[j]]=1;
if(i%pm[j]==0)
{
phi[i*pm[j]]=phi[i]*pm[j]%mod;
break;
}
phi[i*pm[j]]=phi[i]*phi[pm[j]]%mod;
}
}
for(LL i=1;i<MAX;++i)
phi[i]=(phi[i]*i%mod*i%mod+phi[i-1])%mod;
if(n>=MAX) dfs(n);
LL ans=0,L,R,U;
for(LL l=1,r;l<=n;l=r+1)
{
r=n/(n/l),U=S(n/l),L=l-1<MAX?phi[l-1]:ph[n/(l-1)],R=((r<MAX?phi[r]:ph[n/r])-L)%mod;
U=U*U%mod,ans=(ans+R*U%mod)%mod;
}
return printf("%lld",(ans+mod)%mod),0;
}
