【数论分块】bzoj2956: 模积和

数论分块并不精通……第一次调了一个多小时才搞到60pts；因为不会处理i==j的情况，只能枚举了……

Description

$sum_{i=1}^{n}sum_{j=1 land i ot = j}^{m}(n mod i)(m mod j)$

Input

第一行两个数n,m。

Output

一个整数表示答案mod 19940417的值

Sample Input

3 4

Sample Output

样例说明

答案为(3 mod 1)*(4 mod 2)+(3 mod 1) * (4 mod 3)+(3 mod 1) * (4 mod 4) + (3 mod 2) * (4 mod 1) + (3 mod 2) * (4 mod 3) + (3 mod 2) * (4 mod 4) + (3 mod 3) * (4 mod 1) + (3 mod 3) * (4 mod 2) + (3 mod 3) * (4 mod 4) = 1

数据规模和约定

30%: n,m <= 1000

60%: n,m <= 10^6

100% n,m <= 10^9

题目分析

我们有

$原式=sum_{i=1}^{n}sum_{j=1}^{m}(n-{left lfloor frac{n}{i} ight floor}i)(m-{left lfloor frac{m}{j} ight floor}j)-sum_{i=1}^{min(n,m)}(n-{left lfloor frac{n}{i} ight floor}i)(m-{left lfloor frac{m}{i} ight floor}i)$

$=sum_{i=1}^{n}(n-{left lfloor frac{n}{i} ight floor}i)sum_{j=1}^{m}(m-{left lfloor frac{m}{j} ight floor}j)-sum_{i=1}^{min(n,m)}(nm+{left lfloor frac{n}{i} ight floor}{left lfloor frac{m}{i} ight floor}i^2-(m{left lfloor frac{n}{i} ight floor}+n{left lfloor frac{m}{i} ight floor})i)$

化出来的后一项$sum_{i=1}^{min(n,m)}(nm+{left lfloor frac{n}{i} ight floor}{left lfloor frac{m}{i} ight floor}i^2-(m{left lfloor frac{n}{i} ight floor}+n{left lfloor frac{m}{i} ight floor})i)$不是很常规。但注意到$left lfloor frac{n}{i} ight floor$和$left lfloor frac{m}{i} ight floor$都是单调的，那么就可以从小到大枚举的时候顺带取一个min来做。这样的复杂度就是$O(sqrt n+sqrt m)$的了。
大概是这样的：

早上被这最后一步卡住了……

然后就是一些细节上注意取模

 1 #include<bits/stdc++.h>
 2 typedef long long ll;
 3 const int MO = 19940417;
 4 const int inv6 = 3323403;
 5 
 6 ll n,m,ans,del;
 7 
 8 inline void Add(ll &x, ll y){x = ((x+y)%MO+MO)%MO;}
 9 ll sum(ll x){return x*(x+1)%MO*(2*x+1)%MO*inv6%MO;}
10 ll calc(ll x)
11 {
12     ll ret = 0;
13     for (ll i=1, j=0; i<=x; i=j+1)
14     {
15         j = x/(x/i);
16         Add(ret, 1ll*(x/i)*(i+j)*(j-i+1)/2%MO);
17     }
18     return ((x%MO*x%MO-ret)+MO)%MO;
19 }
20 int main()
21 {
22     scanf("%lld%lld",&n,&m);
23     if (n > m) std::swap(n, m);
24     ans = calc(n)*calc(m)%MO;
25     del = n*m%MO*n%MO;
26     for (ll i=1, j=0; i<=n; i=j+1)
27     {
28         j = std::min(n/(n/i), m/(m/i));
29         ll s1 = (sum(j)-sum(i-1))*(n/i)%MO*(m/i)%MO;
30         ll s2 = (n*(m/i)%MO+m*(n/i)%MO)%MO*((i+j)*(j-i+1)/2%MO);
31         Add(del, (s1-s2)%MO);
32     }
33     Add(ans, -del);
34     printf("%lld
",ans);
35     return 0;
36 }

END