【数论】bzoj4174tty的求助

题目链接
（ps：这是一道权限题…）

计算

∑n=1N∑m=1M∑k=0m−1⌊nk+xm⌋

其中N，M<=500000，0<x<=100000，x精确到小数点后8
位小数。

根本不会的说…
然后去看po姐博客就一脸懵逼辣…

以下摘自popoqqq博客，并加上蒟蒻的辣鸡理解…
单独考虑最后一个和式

∑k=0m−1⌊nk+xm⌋

=∑k=0m−1(⌊nk(mod)m+xm⌋+nkm−nk(mod)mm)

上式是考虑nk中m的倍数部分。

然后分开考虑每一个独立的和式。
令d=gcd(n,m)，这个变量将贯穿整篇博客。

∑k=0m−1⌊nk(mod)m+xm⌋

这相当于在枚举m的剩余系，然后很明显剩余系重复了d次，所以只需要枚举一个剩余系即可。

=d∗∑k=0md−1⌊kd+xm⌋

和之前同样的考虑方式，考虑x中m的倍数部分。

=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次，所以只需要枚举一个剩余系即可。

=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到i的公差为1的等差数列的前n项和，即Sum(i)=i∗(i+1)2。

=12(Sum(N)∗Sum(M)−N∗Sum(M)−M∗Sum(N)+∑n=1N∑m=1M(2∗d∗⌊xd⌋+d))

对于最后的和式中的d=gcd(n,m)，我们考虑每一个d出现的次数，这样就可以计算每一个相同的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，然后用经典的莫比乌斯反演求一下gcd(i,j)==d(1<=i<=N,1<=j<=M)的个数

#include <iostream>#include <cstdio>#define LL long long int#define mod 998244353#define MAXN 500005using 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;}