BZOJ 2956 模积和 (数学推导+数论分块)

题目链接: http://www.lydsy.com/JudgeOnline/problem.php?id=2956

题目大意: 求

$$\sum^{n}_{i=1} \sum^{m}_{j=1, j\ne i} (n \mod i)(m \mod j)$$ 对19940417取模的值。

思路分析:
从heheda神犇的博客“heheda的数论专题练习”上翻到这样一道比较有趣的题。
一开始想歪了，一直在考虑n mod i与n mod (i+1)之间的递推关系，后来发现行不通。
既然是取模，那就可以采用“化模为除”的方法，即利用公式

$$n \mod i = n-[n/i]*i$$ 化简原式，暂且不考虑 $i\ne j$的条件可得 $$ans=\sum^{n}_{i=1} \sum^{m}_{j=1} (n-i[n/i])(m-j[m/j]) = \sum^{n}_{i=1} \sum^{m}_{j=1} (nm-nj[m/j]-mi[n/i]+ij[n/i][m/j]) = n^2m^2-n^2\sum^{m}_{j=1}j[m/j]-m^2\sum^{n}_{i=1}i[m/i]+\sum^{n}_{i=1}i[n/i]*\sum^{m}_{j=1}j[m/j]$$
做到这一步，既然只有i和[n/i]出现在式子里(j同理)，显然可以用数论分块加速这个过程了。(至于数论分块的实现可以参见我的另一篇博客 http://blog.csdn.net/suncongbo/article/details/78819470) 分成 $(2\sqrt n+2\sqrt m)$个块，每个块内的 $[n/i]$与 $[m/i]$值均不变，块内求和直接用求和公式 $$getsum(n)=\frac{n(n+1)}{2}$$ , 块内求和等于块右端点的getsum返回值减去块左端点的getsum返回值。

然后去掉$i=j$的情况:

$$\sum^{\min(n,m)}_{i=1} (n-i[n/i])(m-i[m/i]) =\sum^{\min(n,m)}_{i=1}(nm-ni[m/i]-mi[n/i]+i^2[n/i][m/i]) =nm\min(n,m)-\sum^{min(n,m)}_{i=1}i(n[m/i]+m[n/i])+\sum^{min(n,m)}_{i=1}i^2[n/i][m/i]$$ , 前两项仍然可以用前面的方法进行计算。棘手的问题来了：
对于最后一项，牵扯到求连续自然数的平方和。我们知道 $$sqrsum(i)=\sum^{n}_{i=1}i^2=\frac{n(n+1)(2n+1)}{6}$$ , 但是如果把 $n(n+1)(2n+1)$算出来将会很大，long long无法支持。因此只能中途对 $19940417$取模。但是又要除以6，牵扯到求逆元，而19940417又不是质数，无法用费马小定理求逆元。此时做法一是通过其他方法(如exgcd)求出6对于 $19940417$的逆元，二是直接算出6的逆元并作为常量使用。其实，这个数的值是 $3323403$.
而我使用的是第三种做法——分类讨论。(其实是因为懒得写逆元)
观察到当 $n\equiv 0$或 $n\equiv 2 (\mod 3)$时， $n(n+1)$的值可以被6整除，因此先算出 $\frac{n(n+1)}{6}$即可。当 $n\equiv 1 (\mod 3)$时， $2n+1$被3整除， $n(n+1)$被2整除，故计算 $\frac{2n+1}{3}$与 $\frac{n(n+1)}{2}$相乘即可。(当然如果模数是24而不是6我甘愿老老实实地求逆元)
注意运算符优先级问题， 不可以有任何连着三个19940417或1e9级别的数不加取模地连乘，否则必爆。

代码实现

#include<cstdio>
#include<algorithm>
#include<cmath>
using namespace std;

const int NM = 31623;
const long long P = 19940417ll;
long long bl[(NM<<2)+4];
long long tmp[(NM<<2)+4];
long long n,m;
int nm,nn,mm;

void modinc(long long &a)
{
    if(a<0) a+=P;
    if(a>=P) a-=P;
}

void merge(int lb1,int lb2,int rb2)
{
    int i = lb1,j = lb2,k = lb1;
    while(i<=lb2-1 && j<=rb2)
    {
        if(bl[i]<bl[j]) {tmp[k] = bl[i]; i++;}
        else {tmp[k] = bl[j]; j++;}
        k++;
    }
    while(i<=lb2-1) {tmp[k] = bl[i]; i++; k++;}
    while(j<=rb2) {tmp[k] = bl[j]; j++; k++;}
    for(int i=lb1; i<=rb2; i++) bl[i] = tmp[i];
}

long long sqrsum(long long rb)
{
    if(rb%3==1) return ((((rb+rb+1)/3)%P)*(rb*(rb+1)/2)%P)%P;
    return ((rb*(rb+1)/6)%P*(rb+rb+1))%P;
}

long long getsum(long long rb)
{
    return (rb*(rb+1)/2)%P;
}

long long min_ll(long long x,long long y)
{
    return x<y ? x : y;
}

int main()
{
    scanf("%lld%lld",&n,&m); nm = 0;
    nn = sqrt(n); mm = sqrt(m);
    for(int i=1; i<=nn; i++) bl[++nm] = i;
    for(int i=nn; i>=1; i--) bl[++nm] = n/i;
    for(int i=1; i<=mm; i++) bl[++nm] = i;
    for(int i=mm; i>=1; i--) bl[++nm] = m/i;
    merge(1,1+nn+nn,nm);
    long long a = 0ll,b = 0ll,c = 0ll,d = 0ll;
    for(int i=1; i<=nm && bl[i]<=n; i++) {a += (getsum(bl[i])-getsum(bl[i-1])+P)*(n/bl[i])%P; modinc(a);}
    for(int i=1; i<=nm && bl[i]<=m; i++) {b += (getsum(bl[i])-getsum(bl[i-1])+P)*(m/bl[i])%P; modinc(b);}
    c = a*b%P; d = (min_ll(n,m)*n%P)*m%P;
    for(int i=1; i<=nm && bl[i]<=min_ll(n,m); i++)
    {
        long long tmp = (((sqrsum(bl[i])-sqrsum(bl[i-1])+P)*(n/bl[i])%P)*(m/bl[i]))%P;
        d += tmp; modinc(d);
    }
    for(int i=1; i<=nm && bl[i]<=min_ll(n,m); i++)
    {
        long long tmp = ((getsum(bl[i])-getsum(bl[i-1])+P)*((m*(n/bl[i])+n*(m/bl[i]))%P))%P;
        d -= tmp; modinc(d);
    }
    long long ans = ((n*m)%P)*((n*m)%P)%P;
    ans -= (((n*n)%P)*b)%P; modinc(ans);
    ans -= (((m*m)%P)*a)%P; modinc(ans);
    ans += c; modinc(ans);
    ans -= d; modinc(ans);
    printf("%lld\n",ans);
    return 0;
}