二维整除分块

P2260 [清华集训2012]模积和

好像是上面那道题的升级版？

题意:求$\sum\limits_{i=1}^n\sum\limits_{j=1}^m[i≠j](n\%i)(m\%j)$,对某个数取模。

默认$n<m$，下面一定要注意大小关系！

$=\sum\limits_{i=1}^n\sum\limits_{j=1}^m(n-\lfloor n/i\rfloor*i)(m-\lfloor m/j\rfloor*j)[i≠j]$

$=\sum\limits_{i=1}^n\sum\limits_{j=1}^m(n-\lfloor n/i\rfloor*i)(m-\lfloor m/j\rfloor*j)-\sum\limits_{i=1}^n(n-\lfloor n/i\rfloor*i)(m-\lfloor m/i\rfloor*i)$

前半部分:

$=\sum\limits_{i=1}^n(n-\lfloor n/i\rfloor*i)\sum\limits_{j=1}^m(m-\lfloor m/j\rfloor*j)$

设$F(n)=\sum\limits_{i=1}^n(n-\lfloor n/i\rfloor*i)$

则前半部分$=F(n)F(m)$

$F(n)=\sum\limits_{i=1}^n(n-\lfloor n/i\rfloor*i)$

$=n^2-\sum\limits_{i=1}^n\lfloor n/i\rfloor*i$

设$G(n,k)=\sum\limits_{i=1}^n\lfloor k/i\rfloor*i$,上式=$n^2-G(n,n)$

$G(n,k)$直接如上题整除分块即可。

后半部分:

$\sum\limits_{i=1}^n(n-\lfloor n/i\rfloor*i)(m-\lfloor m/i\rfloor*i)$

$=\sum\limits_{i=1}^n(nm-m\lfloor n/i\rfloor*i-n\lfloor m/i\rfloor*i+\lfloor n/i\rfloor\lfloor m/i\rfloor*i^2)$

$=n^2m-mG(n,n)-nG(n,m)+\sum\limits_{i=1}^n\lfloor n/i\rfloor\lfloor m/i\rfloor*i^2$

What?$\sum\limits_{i=1}^n\lfloor n/i\rfloor\lfloor m/i\rfloor*i^2$这个东西好像并不是很好办？

$\lfloor n/i\rfloor\lfloor m/i\rfloor$也有类似$\lfloor n/i\rfloor$的规律,也会形成块状。

令n=10,m=15,得到。

********************************************************************************(10/1)*(15/1)=150
***********************************(10/2)*(15/2)=35
***************(10/3)*(15/3)=15
******(10/ 4)*(15/ 4)=6
******(10/ 5)*(15/ 5)=6
**(10/ 6)*(15/ 6)=2
**(10/ 7)*(15/ 7)=2
*(10/ 8)*(15/ 8)=1
*(10/ 9)*(15/ 9)=1
*(10/10)*(15/10)=1

可以看到也形成块状了，但是具体怎么分块?

其实就是两个块状楼梯叠合在一起(形象理解!MC中掉落沙子)。

所以拐角的个数还是$O(\sqrt{n})$,(约4倍$\sqrt{n}$)

我们找到了一个左端点,

那么右端点=($(n/i)$楼梯的右端点,与$(m/i)$楼梯的右端点中,最近的一个)。

因为上文:$N/i'$与$N/i$相等,则$i'$的最大值为$N/(N/i)$。

所以右端点=$min(n/(n/i),m/(m/i))$;

还有一个$i^2$怎么办？，我们通过小学奥数（代码里面的S函数）可以得到$i^2$的区间和，其余同上。

综上，$ANS=F(n)F(m)-n^2m+mG(n,n)+nG(n,m)-\sum\limits_{i=1}^n\lfloor n/i\rfloor\lfloor m/i\rfloor*i^2$

复杂度$O(\sqrt{n})$

Code:

#include<algorithm>
#include<cstdio>
#define mod 19940417
#define inv6 3323403
#define inv2 9970209
using namespace std;
int n,m;
int G(int n,int lim)
{
  int l=1,r,ans=0;
  for (;l<=lim;l=r+1){
    r=min(n/(n/l),lim);
    ans=(ans+1ll*(n/l)*(r+l)%mod*(r-l+1)%mod*inv2)%mod;
  }return ans;
}
int F(int n)
{return (1ll*n*n-G(n,n))%mod;}
int S(int n)
{return 1ll*n*(n+1)%mod*(n+n+1)%mod*inv6%mod;}
int main()
{
  scanf("%d%d",&n,&m);
  if (n>m)swap(n,m);
  int ans=1ll*F(n)*F(m)%mod;
  ans=(ans+1ll*n*G(m,n)+1ll*m*G(n,n))%mod;
  ans=(ans-1ll*n*n%mod*m)%mod;
  int l=1,r;
  for (;l<=n;l=r+1){
    r=min(n/(n/l),m/(m/l));
    ans=(ans-1ll*(n/l)*(m/l)%mod*(S(r)-S(l-1)))%mod;
  }printf("%d",(ans+mod)%mod);
  return 0;
}