题目链接
(ps:这是一道权限题…)
计算
∑n=1N∑m=1M∑k=0m−1⌊nk+xm⌋
其中 N,M<=500000,0<x<=100000 , x 精确到小数点后
位小数。
根本不会的说…
然后去看po姐博客就一脸懵逼辣…
以下摘自popoqqq博客,并加上蒟蒻的辣鸡理解…
单独考虑最后一个和式
∑k=0m−1⌊nk+xm⌋
=∑k=0m−1(⌊nk(mod)m+xm⌋+nkm−nk(mod)mm)
上式是考虑 nk 中 m 的倍数部分。
然后分开考虑每一个独立的和式。
令
∑k=0m−1⌊nk(mod)m+xm⌋
这相当于在枚举 m 的剩余系,然后很明显剩余系重复了
=d∗∑k=0md−1⌊kd+xm⌋
和之前同样的考虑方式,考虑 x 中
=d∗(md∗x−x(mod)mm+∑k=0md−1⌊kd+x(mod)mm⌋)
很明显, kd+x(mod)m 是不会比 2∗m 大的…
=d∗(x−x(mod)md+∑k=0md−1[kd+x(mod)m≥m])
移项
=d∗(x−x(mod)md+∑k=0md−1[k≥⌊m−x(mod)md⌋])
=d∗(x−x(mod)md+⌊md⌋−⌊m−x(mod)md⌋)
=d∗(x−x(mod)md+⌊x(mod)md⌋)
=d∗⌊xd⌋
对于第二个和式:
∑k=0m−1nkm
根据和式的性质,可以将常数提取出来。
=n∗∑k=0m−1km=n∗m∗(m−1)2∗m
=nm−n2
对于第三个和式,化简的过程与第一个和式的部分相似。
∑k=0m−1nk(mod)mm
这相当于在枚举 m 的剩余系,然后很明显剩余系重复了
=d∗∑k=0md−1kdm=d2m∗md∗(md−1)2
=m−d2
接着,我们来整理一下这些柿子:
∑n=1N∑m=1M(d∗⌊xd⌋+nm−n2−m−d2)
=12∑n=1N∑m=1M(2∗d∗⌊xd⌋+nm−n−m+d)
把和式展开,其中 Sum(i) 表示从 1 到
=12(Sum(N)∗Sum(M)−N∗Sum(M)−M∗Sum(N)+∑n=1N∑m=1M(2∗d∗⌊xd⌋+d))
对于最后的和式中的 d=gcd(n,m) ,我们考虑每一个 d 出现的次数,这样就可以计算每一个相同的
=12(Sum(N)∗Sum(M)−N∗Sum(M)−M∗Sum(N)+
∑d=1min(N,M)(2∗d∗⌊xd⌋+d)∗∑g=1min(nd,md)∗μ(g)∗⌊Ndg⌋∗⌊Mdg⌋))
然后我们就枚举
d
,然后用经典的莫比乌斯反演求一下
#include <iostream>
#include <cstdio>
#define LL long long int
#define mod 998244353
#define MAXN 500005
using namespace std;
int n, m, prime[MAXN], tot;
LL u[MAXN];
bool flag[MAXN];
double x;
inline void init()
{
u[1]=1;
for(int i=2, j, k;i<MAXN;++i)
{
if(!flag[i])prime[++tot]=i, u[i]=-1;
for(j=1;j<=tot&&(k=i*prime[j])<MAXN;++j)
{
flag[k]=1;
if(i%prime[j]==0){u[k]=0;break;}
u[k]=-u[i];
}
}
for(int i=1;i<MAXN;++i)
u[i]=(u[i]+u[i-1]+mod+mod)%mod;
}
LL sum(LL n){return n*(n+1)/2%mod;}
int main()
{
init();
scanf("%d%d%lf",&n,&m,&x);
LL ans=((sum(n)*sum(m)-n*sum(m)-m*sum(n))%mod+mod)%mod, k;
if(n>m)swap(n,m);
for(int i=1, nn, mm, last;i<=n;++i)
{
k=i+int(x/i)*i*2;
nn=n/i, mm=m/i;
for(int j=1;j<=nn;j=last+1)
{
last=min(nn/(nn/j),mm/(mm/j));
ans=(ans+k*(u[last]-u[j-1])%mod*(nn/j)%mod*(mm/j))%mod;
}
}
cout<<ans*499122177%mod<<'\n';
return 0;
}