【数论】bzoj4174tty的求助

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

计算

$$\sum_{n=1}^N\sum_{m=1}^M\sum_{k=0}^{m-1}\lfloor\frac{nk+x}m\rfloor$$
其中 $N，M<=500000，0<x<=100000$， $x$精确到小数点后$8$
位小数。

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

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

$$\sum_{k=0}^{m-1}\lfloor\frac{nk+x}m\rfloor$$
$$=\sum_{k=0}^{m-1}(\lfloor \frac{nk(mod)m+x}m\rfloor+\frac{nk}{m}-\frac{nk(mod)m}m)$$
上式是考虑 $nk$中 $m$的倍数部分。

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

$$\sum_{k=0}^{m-1}\lfloor \frac{nk(mod)m+x}m\rfloor$$
这相当于在枚举 $m$的剩余系，然后很明显剩余系重复了$d$次，所以只需要枚举一个剩余系即可。
$$=d*\sum_{k=0}^{\frac{m}d-1}\lfloor \frac{kd+x}m\rfloor$$
和之前同样的考虑方式，考虑 $x$中$m$的倍数部分。
$$=d*(\frac{m}d*\frac{x-x(mod)m}m+\sum_{k=0}^{\frac{m}d-1}\lfloor\frac{kd+x(mod)m}m\rfloor)$$
很明显， $kd+x(mod)m$是不会比 $2*m$大的…
$$=d*(\frac{x-x(mod)m}d+\sum_{k=0}^{\frac{m}d-1}[kd+x(mod)m≥m])$$
移项
$$=d*(\frac{x-x(mod)m}d+\sum_{k=0}^{\frac{m}d-1}[k≥\lfloor\frac{m-x(mod)m}d\rfloor])$$
$$=d*(\frac{x-x(mod)m}d+\lfloor\frac{m}d\rfloor-\lfloor\frac{m-x(mod)m}d\rfloor)$$
$$=d*(\frac{x-x(mod)m}d+\lfloor\frac{x(mod)m}d\rfloor)$$
$$=d*\lfloor\frac{x}d\rfloor$$

对于第二个和式：

$$\sum_{k=0}^{m-1}\frac{nk}{m}$$
根据和式的性质，可以将常数提取出来。
$$=n*\sum_{k=0}^{m-1}\frac{k}m=n*\frac{m*(m-1)}{2*m}$$
$$=\frac{nm-n}2$$

对于第三个和式，化简的过程与第一个和式的部分相似。

$$\sum_{k=0}^{m-1}\frac{nk(mod)m}m$$
这相当于在枚举 $m$的剩余系，然后很明显剩余系重复了$d$次，所以只需要枚举一个剩余系即可。
$$=d*\sum_{k=0}^{\frac{m}d-1}\frac{kd}m=\frac{d^2}m*\frac{\frac{m}d*(\frac{m}d-1)}2$$
$$=\frac{m-d}2$$

接着，我们来整理一下这些柿子：

$$\sum_{n=1}^N\sum_{m=1}^M(d*\lfloor\frac{x}d\rfloor+\frac{nm-n}2-\frac{m-d}2)$$
$$=\frac{1}2\sum_{n=1}^N\sum_{m=1}^M(2*d*\lfloor\frac{x}d\rfloor+nm-n-m+d)$$
把和式展开，其中 $Sum(i)$表示从 $1$到$i$的公差为 $1$的等差数列的前$n$项和，即 $Sum(i)=\frac{i*(i+1)}2$。
$$=\frac{1}2(Sum(N)*Sum(M)-N*Sum(M)-M*Sum(N)+\sum_{n=1}^N\sum_{m=1}^M(2*d*\lfloor\frac{x}d\rfloor+d))$$
对于最后的和式中的 $d=gcd(n,m)$，我们考虑每一个 $d$出现的次数，这样就可以计算每一个相同的$d$对答案的贡献。
$$=\frac{1}2(Sum(N)*Sum(M)-N*Sum(M)-M*Sum(N)+$$
$$\sum_{d=1}^{min(N,M)}(2*d*\lfloor\frac{x}d\rfloor+d)*\sum_{g=1}^{min(\frac{n}d,\frac{m}d)}*\mu(g)*\lfloor\frac{N}{dg}\rfloor*\lfloor\frac{M}{dg}\rfloor))$$

然后我们就枚举$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 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;
}