题目描述
由于出题人懒得写背景了，题目还是简单一点好。
输入一个整数n和一个整数p，你需要求出((sum_{i=1}^nsum_{j=1}^n ijgcd(i,j))~mod~p)，其中gcd(a,b)表示a与b的最大公约数。
输入输出格式
输入格式：
一行两个整数p、n。
输出格式：
一行一个整数((sum_{i=1}^nsum_{j=1}^n ijgcd(i,j))~mod~p)。
输入输出样例
输入样例#1：
998244353 2000
输出样例#1：
883968974
说明
对于20%的数据，(n leq 1000)。
对于30%的数据，(n leq 5000)。
对于60%的数据，(nleq 10^6)，时限1s。
对于另外20%的数据，(nleq 10^9)，时限3s。
对于最后20%的数据，(n leq 10^{10})，时限6s。
对于100%的数据，(5 imes 10^8 leq p leq 1.1 imes 10^9)且p为质数。
题解
同样的莫比乌斯反演，加上杜教筛
[ans=sum_{i=1}^nsum_{j=1}^nijcdot gcd(i,j)
]
[=sum_{d=1}^nsum_{i=1}^{lfloor frac{n}{i}
floor}sum_{j=1}^{lfloor frac{n}{i}
floor}dicdot djcdot d[gcd(i,j)=1]
]
[=sum_{d=1}^nd^3sum_{i=1}^{lfloor frac{n}{i}
floor}sum_{j=1}^{lfloor frac{n}{i}
floor}ij[gcd(i,j)=1]
]
[=sum_{d=1}^nd^3sum_{i=1}^{lfloor frac{n}{d}
floor}mu(i)cdot i^2cdot s(lfloor frac{n}{id}
floor) (s(n)=(frac{n*(n+1)}{2})^2)
]
[=sum_{T=1}^ns(lfloor frac{n}{T}
floor)sum_{d|T}d^3cdot mu(frac{T}{d})cdot (frac{T}{d})^2
]
[=sum_{T=1}^ns(lfloor frac{n}{T}
floor)T^2sum_{d|T}dcdot mu(frac{T}{d})
]
[=sum_{T=1}^ns(lfloor frac{n}{T}
floor)T^2varphi(T)
]
最后一步与欧拉函数有关，也与卷积有关
对于前面(sum_{T=1}^ns(lfloor frac{n}{T} floor))可以整除分块，我们考虑后面部分的前缀和
设(S(n)=sum_{i=1}^ni^2varphi(i))
然后上杜教筛
[S(n)=sum_{i=1}^ni^2varphi(i)=sum_{i=1}^ni^2(sum_{d|i}varphi(d)-sum_{d|i,d
e i}varphi(d))
]
[=sum_{i=1}^ni^2sum_{d|i}varphi(d)-sum_{i=1}^ni^2sum_{d|i,d
e i}varphi(d)
]
[=sum_{i=1}^ni^3-sum_{i=1}^nsum_{d|i,d
e i}i^2varphi(d)
]
[=sum_{i=1}^ni^3-sum_{frac{i}{d}=2}^nsum_{frac{i}{d}|i}^ni^2varphi(d)
]
[=sum_{i=1}^ni^3-sum_{x=2}^nsum_{x|xd}^n(xd)^2varphi(d)
]
[=sum_{i=1}^ni^3-sum_{x=2}^nsum_{d=1}^{lfloor frac{n}{x}
floor}x^2d^2varphi(d)
]
[=sum_{i=1}^ni^3-sum_{x=2}^nx^2sum_{d=1}^{lfloor frac{n}{x}
floor}d^2varphi(d)
]
[=sum_{i=1}^ni^3-sum_{x=2}^nx^2S(lfloor frac{n}{x}
floor)
]
成功杜教筛，复杂度(O(n^frac{2}{3}))？
先预处理前(1e7)项，后面的用杜教筛的式子求
[ans=sum_{T=1}^ns(lfloor frac{n}{T}
floor )S(T)
]
整除分块了，杜教筛了，然后就过了
这一题涉及了欧拉函数和卷积，但这些还没学通，式子中有些东西是强背的
以后还会系统地学
在代码实现的过程中用到了平方与立方数列求和公式，不然复杂度不对，详细公式可以见这里（以前从未听说过这么强的公式。。。）
#include<bits/stdc++.h>
#define ll long long
#define db double
#define ld long double
const int MAXN=1000000+10;
int Mod,cnt,prime[MAXN],vis[MAXN];
ll phi[MAXN],f[MAXN],six,two;
std::map<ll,ll> M;
template<typename T> inline void read(T &x)
{
T data=0,w=1;
char ch=0;
while(ch!='-'&&(ch<'0'||ch>'9'))ch=getchar();
if(ch=='-')w=-1,ch=getchar();
while(ch>='0'&&ch<='9')data=((T)data<<3)+((T)data<<1)+(ch^'0'),ch=getchar();
x=data*w;
}
template<typename T> inline void write(T x,char c=' ')
{
if(x<0)putchar('-'),x=-x;
if(x>9)write(x/10);
putchar(x%10+'0');
if(c!=' ')putchar(c);
}
template<typename T> inline void chkmin(T &x,T y){x=(y<x?y:x);}
template<typename T> inline void chkmax(T &x,T y){x=(y>x?y:x);}
template<typename T> inline T min(T x,T y){return x<y?x:y;}
template<typename T> inline T max(T x,T y){return x>y?x:y;}
inline ll qexp(ll a,ll b)
{
ll res=1;
while(b)
{
if(b&1)res=res*a%Mod;
a=a*a%Mod;
b>>=1;
}
return res;
}
inline void init()
{
two=qexp(2,Mod-2);
six=qexp(6,Mod-2);
memset(vis,1,sizeof(vis));
vis[0]=vis[1]=0;
phi[1]=1;
for(register int i=2;i<MAXN;++i)
{
if(vis[i])
{
prime[++cnt]=i;
phi[i]=i-1;
}
for(register int j=1;j<=cnt&&i*prime[j]<MAXN;++j)
{
vis[i*prime[j]]=0;
if(i%prime[j])phi[i*prime[j]]=phi[i]*phi[prime[j]]%Mod;
else
{
phi[i*prime[j]]=phi[i]*(ll)prime[j]%Mod;
break;
}
}
}
for(register ll i=1;i<MAXN;++i)f[i]=(f[i-1]+i*i%Mod*phi[i]%Mod)%Mod;
}
inline ll s3(ll x)
{
x%=Mod;
ll res=x*(x+1)%Mod*two%Mod;
return res*res%Mod;
}
inline ll s2(ll x)
{
x%=Mod;
return x*(x+1)%Mod*(x+x+1)%Mod*six%Mod;
}
inline ll Phis(ll x)
{
if(x<MAXN)return f[x];
if(M[x])return M[x];
ll res=s3(x);
for(register ll i=2;;)
{
if(i>x)break;
ll j=x/(x/i);
(res-=Phis(x/i)*(s2(j)-s2(i-1))%Mod)%=Mod;
i=j+1;
}
return M[x]=(res+Mod)%Mod;
}
inline ll solve(ll n)
{
ll res=0;
for(register ll i=1;;)
{
if(i>n)break;
ll j=n/(n/i);
(res+=s3(n/i)*(Phis(j)-Phis(i-1))%Mod)%=Mod;
i=j+1;
}
return (res+Mod)%Mod;
}
int main()
{
ll n;
read(Mod);read(n);
init();
write(solve(n),'
');
return 0;
}