bzoj2956 模积和 （分块）

bzoj2956 模积和

原题地址：http://www.lydsy.com/JudgeOnline/problem.php?id=2956

题意：
求∑∑((n mod i)*(m mod j))，其中1<=i<=n,1<=j<=m,i≠j。答案 mod 19940417

数据范围
n,m<=10^9

题解：
$∑_{i=1}^n∑_{j=1}^m((n \ mod\ i)*(m\ mod \ j)) \ \ \ (i≠j)$
$=∑_{i=1}^n∑_{j=1}^m((n \ mod\ i)*(m\ mod \ j))-∑_{i=1}^{min(n,m)}((n \ mod\ i)*(m\ mod \ i))$

前一部分
$=∑_{i=1}^n∑_{j=1}^m((n-i*⌊\frac{n}{i}⌋)*(m-j*⌊\frac{m}{j}⌋))$
$=∑_{i=1}^n(n-i*⌊\frac{n}{i}⌋)∑_{j=1}^m(m-j*⌊\frac{m}{j}⌋)$
可以分块做。

后一部分$=∑_{i=1}^{min(n,m)}((n \ mod\ i)*(m\ mod \ i))$
$=∑_{i=1}^{min(n,m)}((n-i*⌊\frac{n}{i}⌋)*(m-i*⌊\frac{m}{i}⌋))$
$=∑_{i=1}^{min(n,m)}(n*m-i*⌊\frac{n}{i}⌋*m-i*⌊\frac{m}{i}⌋*n+i^2⌊\frac{m}{i}⌋⌊\frac{n}{i}⌋)$
也可以分块做。

其中用到了 $1+2+..+n=\frac{(1+n)*n}{2}$
和$1^2+2^2+..+n^2=\frac{n*(1+n)*(2n+1)}{6}$
因为19940417和2、6互质，可以预处理逆元。

代码：

#include<cstdio> 
#include<iostream>
#include<cstring>
#include<algorithm>
using namespace std;
const int mod=19940417;
const int phi=17091780;
const int inv2=9970209;
const int inv6=3323403;
int n,m;    
int main()
{
    scanf("%d%d",&n,&m);
    int ret1=0,ret2=0;
    for(int ed,i=1;i<=m;i=ed+1)
    {
        ed=(m/(m/i)); int len=ed-i+1;
        int a=(1LL*(i+ed)*len)%mod; a=(1LL*a*inv2)%mod; a=(1LL*a*(m/i))%mod;     
        ret1=(ret1+a)%mod;
    }
    ret1=((1LL*m*m)%mod-ret1+mod)%mod;
    for(int ed,i=1;i<=n;i=ed+1)
    {
        ed=(n/(n/i)); int len=ed-i+1;
        int a=(1LL*(i+ed)*len)%mod; a=(1LL*a*inv2)%mod; a=(1LL*a*(n/i))%mod;     
        ret2=(ret2+a)%mod;
    }
    ret2=((1LL*n*n)%mod-ret2+mod)%mod;
    int ans=(1LL*ret1*ret2)%mod; 
    int ret=0;
    for(int ed,i=1;i<=min(n,m);i=ed+1)
    {
        ed=min(min((n/(n/i)),(m/(m/i))),min(m,n));  int len=ed-i+1; 
        int a=(1LL*(i+ed)*len)%mod; a=(1LL*a*inv2)%mod; 
        a=((1LL*(1LL*a*(n/i))%mod*m)%mod+(1LL*(1LL*a*(m/i))%mod*n)%mod)%mod;
        int b=(1LL*(2*ed%mod+1)%mod*((1LL*ed*(ed+1))%mod))%mod; b=(1LL*b*inv6)%mod;
        int c=(1LL*(2*(i-1)%mod+1)%mod*((1LL*i*(i-1))%mod))%mod;c=(1LL*c*inv6)%mod;
        b=(b-c+mod)%mod; b=(1LL*b*(n/i))%mod; b=(1LL*b*(m/i))%mod;
        a=(a-b+mod)%mod; ret=(ret+a)%mod;
    }
    ret=((1LL*((1LL*min(m,n)*n)%mod)*m)%mod-ret+mod)%mod;
    ans=(ans-ret+mod)%mod;
    printf("%d\n",ans);
    return 0;
}