[BZOJ2956]模积和-根号分块

说在前面

突然觉得”说在前面”里的废话真是太多了= =

题目

BZOJ2956传送门

题目大意

给定$N$,$M$,求出$\sum\limits_{i=1}^{N}\sum\limits_{j=1}^{M}\left(N\bmod i\right)\cdot \left(M\bmod j\right)[i!=j]$ 在模19940417意义下的值。$n,m\leq 10^9$

输入输出格式

输入格式：
输入一行两个整数$N$，$M$，含义如题

输出格式：
输出一行一个整数，表示答案

解法

直接化简式子，把题目中$i\not = j$ 的拆成$\forall i,j$的答案减去$i=j$的答案。
对于$\forall i , j$ 有：
$\sum\limits_{i=1}^{N}\sum\limits_{j=1}^{M}\left(N\bmod i\right)\cdot \left(M\bmod j\right)$
$=\sum\limits_{i=1}^{N}\sum\limits_{j=1}^{M}\left(N - \left\lfloor\frac{n}{i}\right\rfloor i\right) \left(M-\left\lfloor\frac{M}{j}\right\rfloor j\right)$
$=\sum\limits_{i=1}^{N}\sum\limits_{j=1}^{M} \left(NM - N\left\lfloor\frac{M}{j}\right\rfloor j - M\left\lfloor\frac{N}{i}\right\rfloor i + ij \left\lfloor\frac{N}{i}\right\rfloor \left\lfloor\frac{M}{j}\right\rfloor \right)$
$=N^2M^2 - N^2\sum\limits_{j=1}^{M}\left\lfloor\frac{M}{j}\right\rfloor j - M^2\sum\limits_{i=1}^{n}\left\lfloor\frac{N}{i}\right\rfloor i +\sum\limits_{i=1}^{N}i*\left\lfloor\frac{N}{i}\right\rfloor\sum\limits_{j=1}^{M}j*\left\lfloor\frac{M}{j}\right\rfloor$

对于$i=j$的情况类似，可以自行推导=w=
（其实是懒得写….Latex太麻烦了，下面这张图就是上面那一串的源码，mmp）

自带大常数的代码

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std ;const int mmod = 19940417 , inv6 = 3323403 ;
long long N , M ;long long ds( long long lf , long long rg ){return ( lf + rg ) * ( rg - lf + 1 ) / 2 %mmod ;
}long long pfs( long long lf , long long rg ){lf = lf - 1 ;long long tmp1 = rg * ( rg + 1 )%mmod * ( 2 * rg + 1 ) %mmod * inv6 %mmod ,tmp2 = lf * ( lf + 1 )%mmod * ( 2 * lf + 1 ) %mmod * inv6 %mmod ;return ( tmp1 - tmp2 + mmod )%mmod ;
}long long cal( long long a , long long lim , long long div ){long long rt = 0 ;for( int i = 1 , last = 1 ; i <= lim ; i = last ){int ed = min( div / ( div / i ) , lim ) ;rt = ( rt + ( div / i ) * ds( last , ed ) ) %mmod ;last = ed + 1 ;}return rt * a %mmod ;
}long long spe(){long long rt = 0 ;for( int i = 1 , last = 1 ; i <= M ; i = last ){int ed = min( M / ( M / i ) , N / ( N / i ) ) ;rt = ( rt + ( M / i ) * ( N / i )%mmod * pfs( last , ed )%mmod )%mmod ;last = ed + 1 ;}return rt ;
}int main(){long long ans1 = 0 , ans2 = 0 ;scanf( "%lld%lld" , &N , &M ) ;if( N < M ) swap( N , M ) ;ans1 = N * N%mmod * M%mmod * M%mmod - cal( N*N%mmod , M , M ) - cal( M*M%mmod , N , N ) + cal( 1 , N , N ) * cal( 1 , M , M ) %mmod ;ans1 = ( ans1%mmod + mmod )%mmod ;ans2 = N * M%mmod * M%mmod - cal( N%mmod , M , M ) - cal( M%mmod , M , N ) + spe() ;ans2 = ( ans2%mmod + mmod )%mmod ;printf( "%lld" , ( ans1 - ans2 + mmod )%mmod ) ;
}