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
题解:
∑ni=1∑mj=1((n mod i)∗(m mod j)) (i≠j)
=∑ni=1∑mj=1((n mod i)∗(m mod j))−∑min(n,m)i=1((n mod i)∗(m mod i))
前一部分
=∑ni=1∑mj=1((n−i∗⌊ni⌋)∗(m−j∗⌊mj⌋))
=∑ni=1(n−i∗⌊ni⌋)∑mj=1(m−j∗⌊mj⌋)
可以分块做。
后一部分
=∑min(n,m)i=1((n mod i)∗(m mod i))
=∑min(n,m)i=1((n−i∗⌊ni⌋)∗(m−i∗⌊mi⌋))
=∑min(n,m)i=1(n∗m−i∗⌊ni⌋∗m−i∗⌊mi⌋∗n+i2⌊mi⌋⌊ni⌋)
也可以分块做。
其中用到了
1+2+..+n=(1+n)∗n2
和
12+22+..+n2=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;
}