题意
题意:求∑i=1n∑j=1m[i=j](nmodi)⋅(mmodj)。
数据范围:1⩽n,m⩽109。
分析
这是一道简单的推式子题,但是实现比较恶心。
首先不妨设n⩽m(如果n>m交换一下就好了)
然后可以用容斥将题目拆成两个部分:
i=1∑nj=1∑m(nmodi)⋅(mmodj)−i=1∑n(nmodi)⋅(mmodi)
将两个求和拆开:
i=1∑n(nmodi)⋅j=1∑m(mmodj)−i=1∑n(nmodi)⋅(mmodi)
因为取模很难搞,所以可以用一个性质amodb=a−⌊ba⌋⋅b(取模的定义),将上式转换为:
i=1∑n(n−⌊in⌋⋅i)⋅j=1∑m(m−⌊jm⌋⋅j)−i=1∑n(n−⌊in⌋⋅i)⋅(m−⌊im⌋⋅i)
将括号拆开,可以得到:
(n2−i=1∑ni⋅⌊in⌋)⋅(m2−i=1∑mi⋅⌊im⌋)−i=1∑n(nm−mi⋅⌊in⌋−ni⋅⌊im⌋+i2⋅⌊in⌋⋅⌊im⌋)
此时我们的复杂度已经是O(n)了,然而这个复杂度仍然不足以通过本题。
要做这道题需要一个简单的技巧——整除分块。
我们很容易发现很多⌊in⌋的值是一样的,且所有的⌊in⌋为一个不下降子序列,呈块状分布
通过简单的计算可以得到,对于每个起点为i的块,它的值为⌊in⌋,终点为⌊⌊in⌋n⌋,然后我们就可以用O(n
)的算法计算∑i=1ni⋅⌊in⌋了:
主要步骤,对于每一个块[l,r],它乘号前面前缀和差分得到,即(1+2+⋯+r)−(1+2+⋯+(l−1)),乘号后面的就是⌊ln⌋。
代码:
inline long long sum1(long long x){
return x*(x+1)%mod*inv2%mod;
}
l=1,sum=n*n%mod;
while(l<=n){
r=n/(n/l);
sum=(sum-(sum1(r)-sum1(l-1)+mod)%mod*(n/l)%mod+mod)%mod;
l=r+1;
}
同样,∑i=1mi⋅⌊im⌋的值也可以用相同的方法求得。
考虑求∑i=1n(nm−mi⋅⌊in⌋−ni⋅⌊im⌋+i2⋅⌊in⌋⋅⌊im⌋),发现用整除分块完成它需要使区间[l,r]中所有的i都满足⌊in⌋与⌊im⌋相同。
可以用类似的方法,将上面代码中的r=n/(n/l);改为r=min(n/(n/l),m/(m/l));,然后就可以直接求值了!
且慢,虽然前三项nm,mi⋅⌊in⌋,ni⋅⌊im⌋都很容易求得,但是i2⋅⌊in⌋⋅⌊im⌋不是很容易求,因为(12+22+⋯+i2)不好处理。
这里有一个简单的结论:12+22+⋯+n2=6n(n+1)(2n+1),会在文后证明,现在先给出这一部分的代码:
inline long long sum2(long long x){
return x*(x+1)%mod*(2*x+1)%mod*inv6%mod;
}
l=1,tmp3=0;
while(l<=n){
long long a,b,c;
r=min(n/(n/l),m/(m/l));
a=(r-l+1)*n%mod*m%mod;
b=(sum1(r)-sum1(l-1)+mod)%mod*((n/l)*m%mod+(m/l)*n%mod)%mod;
c=(sum2(r)-sum2(l-1)+mod)%mod*(n/l)%mod*(m/l)%mod;
tmp3=(tmp3+a-b+c+mod)%mod;
l=r+1;
}
代码
记得开long long取模很恶心,如果错了多半是这种原因因为三个整数相乘会爆long long,因此除法用逆元实现(逆元提前求得)
#include<stdio.h>
const long long mod=19940417,inv2=9970209,inv6=3323403;
long long i,j,k,m,n,l,r,ans,tmp1,tmp2,tmp3;
inline long long min(long long a,long long b){
return a<b? a:b;
}
inline void swp(long long &a,long long &b){
a+=b,b=a-b,a-=b;
}
inline long long sum1(long long x){
return x*(x+1)%mod*inv2%mod;
}
inline long long sum2(long long x){
return x*(x+1)%mod*(2*x+1)%mod*inv6%mod;
}
int main(){
scanf("%lld%lld",&n,&m);
if(n>m)
swp(n,m);
l=1,tmp1=n*n%mod;
while(l<=n){
r=n/(n/l);
tmp1=(tmp1-(sum1(r)-sum1(l-1)+mod)%mod*(n/l)%mod+mod)%mod;
l=r+1;
}
l=1,tmp2=m*m%mod;
while(l<=m){
r=m/(m/l);
tmp2=(tmp2-(sum1(r)-sum1(l-1)+mod)%mod*(m/l)%mod+mod)%mod;
l=r+1;
}
l=1,tmp3=0;
while(l<=n){
long long a,b,c;
r=min(n/(n/l),m/(m/l));
a=(r-l+1)*n%mod*m%mod;
b=(sum1(r)-sum1(l-1)+mod)%mod*((n/l)*m%mod+(m/l)*n%mod)%mod;
c=(sum2(r)-sum2(l-1)+mod)%mod*(n/l)%mod*(m/l)%mod;
tmp3=(tmp3+a-b+c+mod)%mod;
l=r+1;
}
ans=(tmp1*tmp2%mod-tmp3+mod)%mod;
printf("%lld\n",ans);
return 0;
}
简单的证明
证明:
12+22+⋯+n2=6n(n+1)(2n+1)
(这个式子可以用简单构造法与数学归纳法证明,由于大部分题解都是用的数学归纳法,且用数学归纳法读者自证不难,因此这里使用简单构造法)
由于有这样一个式子:i2=i2−i+i(i−1)⋅i+i,因此我们可以将这个拆开:
12+22+⋯+n2=i=1∑ni2=i=1∑n(i−1)⋅i+i=1∑ni
然后乘上31:
i=1∑ni2=31i=1∑n(i−1)⋅i⋅3+i=1∑ni
构造一下:
i=1∑ni2=31i=1∑n(i−1)⋅i⋅((i+1)−(i−2))=31i=1∑n(−(i−2)⋅(i−1)⋅i+(i−1)⋅i⋅(i+1))+i=1∑ni
把它们全部展开:
i=1∑ni2=31(−(−1)⋅0⋅1+0⋅1⋅2−0⋅1⋅2+1⋅2⋅3−1⋅2⋅3+⋯(n−1)⋅n⋅(n+1))+2n⋅(n+1)
发现括号里的很多项都可以抵消:
i=1∑ni2=3(n−1)⋅n⋅(n+1)+2n⋅(n+1)=6n⋅(n+1)⋅(2(n−1)+3)=6n⋅(n+1)⋅(2n+1)
证毕。