P3768 简单的数学题
题目描述
由于出题人懒得写背景了，题目还是简单一点好。
输入一个整数(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)。
说明
对于(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 imes 10^8 leq p leq 1.1 imes 10^9)且(p)为质数。
从各种方向推推式子，你会差不多发现有
[sum_{i=1}^nsum_{j=1}^n ijgcd(i,j)
]
[=sum_{T=1}^nF(lfloorfrac{n}{T}
floor)^2T^2varphi(T)
]
其中(F(n)=sumlimits_{i=1}^ni)
然后上杜教筛设(mathbf f(n)=n^2varphi(n))，则有
[mathbf {Id}^3=mathbf f *mathbf {Id}^2
]
带进去杜教筛得到
[mathbf s(n)=sum_{i=1}^n i^3-sum_{i=2}^n i^2 mathbf s(lfloorfrac{n}{i}
floor)
]
然后小学奥数一波算前缀和就行了
小心爆(long long)
Code:
#include <cstdio>
#include <unordered_map>
#define ll long long
const int N=5e6;
ll n,mod,phi[N+10],inv2,inv6;
int pri[N+10],ispri[N+10],cnt;
ll qp(ll d,ll k){ll re=1;while(k){if(k&1)re=re*d%mod;d=d*d%mod,k>>=1;}return re;}
ll f(ll x){x%=mod;return x*(x+1)%mod*inv2%mod;}
ll g(ll x){x%=mod;return x*(x+1)%mod*(2*x%mod+1)%mod*inv6%mod;}
void init()
{
phi[1]=1;
for(int i=2;i<=N;i++)
{
if(!ispri[i])
{
phi[i]=i-1;
pri[++cnt]=i;
}
for(int j=1;j<=cnt&&i*pri[j]<=N;j++)
{
ispri[i*pri[j]]=1;
if(i%pri[j]==0){phi[i*pri[j]]=phi[i]*pri[j]%mod;break;}
else phi[i*pri[j]]=phi[i]*(pri[j]-1)%mod;
}
}
for(int i=1;i<=N;i++)
phi[i]=(phi[i]*i%mod*i%mod+phi[i-1])%mod;
}
std::unordered_map <ll,ll> Phi;
ll calphi(ll n)
{
if(n<=N) return phi[n];
if(Phi.find(n)!=Phi.end()) return Phi[n];
ll ret=f(n)*f(n)%mod;
for(ll l=2,r;l<=n;l=r+1)
{
r=n/(n/l);
(ret-=(calphi(n/l)*(g(r)-g(l-1))%mod))%=mod;
}
ret=(ret%mod+mod)%mod;
return Phi[n]=ret;
}
int main()
{
scanf("%lld%lld",&mod,&n);
init();
ll ans=0;inv6=qp(6,mod-2);inv2=qp(2,mod-2);
for(ll l=1,r;l<=n;l=r+1)
{
r=n/(n/l);
(ans+=f(n/l)*f(n/l)%mod*(calphi(r)-calphi(l-1))%mod)%=mod;
}
ans=(ans%mod+mod)%mod;
printf("%lld
",ans);
return 0;
}
2018.11.26