题目描述
由于出题人懒得写背景了,题目还是简单一点好。
输入一个整数n和一个整数p,你需要求出sum(i*j*gcd(i,j)) {1<=i,j<=n}。gcd(a,b)表示a与b的最大公约数。
刚才题面打错了,已修改
输入输出格式
输入格式:
一行两个整数p、n。
输出格式:
一行一个整数
输入输出样例
输入样例#1:
998244353 2000
输出样例#1:
883968974
说明
对于20%的数据, n≤1000 。
对于30%的数据, n≤5000 。
对于60%的数据, n≤10^6
,时限1s。
对于另外20%的数据, n≤10^9
,时限3s。
对于最后20%的数据, n≤10^10
,时限6s。
对于100%的数据, 5×10^8 ≤p≤1.1×10^9,且p为质数。
分析:https://blog.csdn.net/qq_33229466/article/details/80093762
容易反演得到ans=sum(g(n/i)*i^2*phi(i)) 其中n/i为整除
后面这个可以杜教筛解决了。
注意:map是平衡树,自带log,不可以把预处理的压入map中(跑了13s);另外10^10时计算时要先取模。
代码:
// luogu-judger-enable-o2
#include <iostream>
#include <cstdio>
#include <ctime>
#include <cmath>
#include <map>
#define LL long long
const int maxn=1e7+7;
using namespace std;
LL n,p,ans,inv6,inv4;
LL cnt;
LL f[maxn];
int phi[maxn];
int prime[maxn];
bool not_prime[maxn];
map <LL,LL> h;
void getf(LL n)
{
phi[1]=1;
for (LL i=2;i<=n;i++)
{
if (!not_prime[i])
{
prime[++cnt]=i;
phi[i]=i-1;
}
for (LL j=1;j<=cnt;j++)
{
if (i*prime[j]>n) break;
not_prime[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);
}
}
LL sum=0;
for (LL i=1;i<=n;i++)
{
f[i]=(f[i-1]+i*i%p*phi[i])%p;
}
}
LL power(LL x,LL y)
{
if (y==1) return x;
LL c=power(x,y/2);
c=(c*c)%p;
if (y%2) c=(c*x)%p;
return c;
}
LL sum_power2(LL n)
{
return n%p*(n+1)%p*(2*n%p+1)%p*inv6%p;
}
LL sum_power3(LL n)
{
return n%p*n%p*(n+1)%p*(n+1)%p*inv4%p;
}
LL getsum(LL n)
{
if (n<=10000000) return f[n];
LL c=h[n];
if (c) return c;
LL sum=0;
for (LL i=2,last;i<=n;i=last+1)
{
last=n/(n/i);
sum=(sum+getsum(n/i)*(sum_power2(last)+p-sum_power2(i-1)))%p;
}
c=sum_power3(n)-sum;
h[n]=c;
return c;
}
int main()
{
scanf("%lld%lld",&p,&n);
inv6=power(6,p-2);
inv4=power(4,p-2);
getf(10000000);
LL x=0,y=0;
for (LL i=1,last;i<=n;i=last+1)
{
last=n/(n/i);
y=getsum(last);
ans=(ans+sum_power3(n/i)*(y+p-x)%p)%p;
x=y;
}
printf("%lld",ans);
}