牛客小白月赛 18 I：Forsaken遇到了毒瘤（数论，整除分块）

---

$n\%d+m\%d\ge d$ 可以写成 $n-\lfloor \frac{n}{d}\rfloor \ast d+m-\lfloor \frac{m}{d}\rfloor \ast d\ge d$

移项一下可以得出：$\lfloor \frac{n}{d}\rfloor +\lfloor \frac{m}{d}\rfloor +1\ge \lfloor \frac{n+m}{d}\rfloor$

枚举 $d$ ，$check$ 这个条件，这个条件可以整除分块。注意到 $d$ 的 上限是 $n+m$，设 $n\le m$，将答案分成三个部分进行分块统计，即：$1\le d\le n，n+1\le d\le m，m+1\le d\le n+m$

最后一个块整个块都是答案。
现在只需要能快速处理块间数字的因数和即可。
$\sum _{x=1}^n\sigma (x)=\sum _{x=1}^n\sum _{d|x}d=\sum _{d=1}^{n}d\sum _{x=1}^{\lfloor \frac{n}{d}\rfloor }1=\sum _{d=1}^{n}d\ast \lfloor \frac{n}{d}\rfloor$

这个式子也可以分块，但这两个分块并不是嵌套的形式，直接暴力做复杂度会达到$10^9$，而不是$n^{\frac{3}{4}}$
用线性筛先预处理一部分（大概 $10^7$），把复杂度降下来，然后就可以过了

---

代码：

#include<bits/stdc++.h>
using namespace std;
const int mod = 1e9 + 7;
typedef long long ll;
ll n,m;

const int N=1e7 + 5;
bool mark[N];
int prim[N];
ll sd[N],sp[N];
int cnt;
void initial()
{
    cnt=0;
    sd[1]=1;
    for (int i=2 ; i<N ; ++i)
    {
        if (!mark[i])
        {
            prim[cnt++]=i;
            sd[i]=i+1;
            sp[i]=i+1;
        }
        for (int j=0 ; j<cnt && i*prim[j]<N ; ++j)
        {
            mark[i*prim[j]]=1;
            if (!(i%prim[j]))
            {
                sp[i*prim[j]]=sp[i]*prim[j]+1;
                sd[i*prim[j]]=sd[i]/sp[i]*sp[i*prim[j]];
                break;
            }
            sd[i*prim[j]]=sd[i]*sd[prim[j]] % mod;
            sp[i*prim[j]]=1+prim[j];
        }
    }
    for(int i = 1; i < N; i++)
    	sd[i] = (sd[i] + sd[i - 1]) % mod;
}
ll fpow(ll a,ll b) {
	ll r = 1;
	while(b) {
		if(b & 1) r = r * a % mod;
		b >>= 1;
		a = a * a % mod;
	}
	return r;
}
ll inv2;
ll cal(ll x) {
	return 1ll * x * (x + 1) % mod * inv2 % mod;
}
ll getsum(ll x) {
	if(x < N) return sd[x];
	int l,r;
	ll ans = 0;
	for(l = 1; l <= x; l = r + 1) {
		r = x / (x / l);
		ans = ans +  (cal(r) - cal(l - 1) + mod) % mod * (x / l) % mod;
		ans %= mod;
	}
	return ans;
}
int main() {
	initial();
	inv2 = fpow(2,mod - 2);
	scanf("%lld%lld",&n,&m);
	if(n > m) swap(n,m);
	ll l,r;
	ll res = 0;
	for(l = 1; l <= n; l = r + 1) {
		r = min(n / (n / l),m / (m / l));
		r = min(r,(n + m) / ((n + m) / l));
		if((n / l + m / l + 1) <= (n + m) / l) {
			res = (res + getsum(r) - getsum(l - 1) + mod) % mod;
		}
	}
	for(l = n + 1; l <= m; l = r + 1) {
		r = min(m / (m / l),(n + m) / ((n + m) / l));
		if((m / l + 1) <= (m + n) / l)
			res = (res + getsum(r) - getsum(l - 1) + mod) % mod;
	}
	res = (res + getsum(n + m) - getsum(m) + mod) % mod;
	printf("%lld\n",res);
	return 0;
}