BZOJ 2956 模积和 (数学推导+数论分块)
2018-02-23 16:47
337 查看
题目链接: http://www.lydsy.com/JudgeOnline/problem.php?id=2956
题目大意: 求∑i=1n∑j=1,j≠im(nmodi)(mmodj)∑i=1n∑j=1,j≠im(nmodi)(mmodj)对19940417取模的值。
思路分析:
从heheda神犇的博客“heheda的数论专题练习”上翻到这样一道比较有趣的题。
一开始想歪了,一直在考虑n mod i与n mod (i+1)之间的递推关系,后来发现行不通。
既然是取模,那就可以采用“化模为除”的方法,即利用公式nmodi=n−[n/i]∗inmodi=n−[n/i]∗i化简原式,暂且不考虑i≠ji≠j的条件可得ans=∑i=1n∑j=1m(n−i[n/i])(m−j[m/j])=∑i=1n∑j=1m(nm−nj[m/j]−mi[n/i]+ij[n/i][m/j])=n2m2−n2∑j=1mj[m/j]−m2∑i=1ni[m/i]+∑i=1ni[n/i]∗∑j=1mj[m/j]ans=∑i=1n∑j=1m(n−i[n/i])(m−j[m/j])=∑i=1n∑j=1m(nm−nj[m/j]−mi[n/i]+ij[n/i][m/j])=n2m2−n2∑j=1mj[m/j]−m2∑i=1ni[m/i]+∑i=1ni[n/i]∗∑j=1mj[m/j]
做到这一步,既然只有i和[n/i]出现在式子里(j同理),显然可以用数论分块加速这个过程了。(至于数论分块的实现可以参见我的另一篇博客http://blog.csdn.net/suncongbo/article/details/78819470) 分成(2n−−√+2m−−√)(2n+2m)个块,每个块内的[n/i][n/i]与[m/i][m/i]值均不变,块内求和直接用求和公式getsum(n)=n(n+1)2getsum(n)=n(n+1)2, 块内求和等于块右端点的getsum返回值减去块左端点的getsum返回值。
然后去掉i=ji=j的情况:∑i=1min(n,m)(n−i[n/i])(m−i[m/i])=∑i=1min(n,m)(nm−ni[m/i]−mi[n/i]+i2[n/i][m/i])=nmmin(n,m)−∑i=1min(n,m)i(n[m/i]+m[n/i])+∑i=1min(n,m)i2[n/i][m/i]∑i=1min(n,m)(n−i[n/i])(m−i[m/i])=∑i=1min(n,m)(nm−ni[m/i]−mi[n/i]+i2[n/i][m/i])=nmmin(n,m)−∑i=1min(n,m)i(n[m/i]+m[n/i])+∑i=1min(n,m)i2[n/i][m/i], 前两项仍然可以用前面的方法进行计算。棘手的问题来了:
对于最后一项,牵扯到求连续自然数的平方和。我们知道sqrsum(i)=∑i=1ni2=n(n+1)(2n+1)6sqrsum(i)=∑i=1ni2=n(n+1)(2n+1)6, 但是如果把n(n+1)(2n+1)n(n+1)(2n+1)算出来将会很大,long long无法支持。因此只能中途对1994041719940417取模。但是又要除以6,牵扯到求逆元,而19940417又不是质数,无法用费马小定理求逆元。此时做法一是通过其他方法(如exgcd)求出6对于1994041719940417的逆元,二是直接算出6的逆元并作为常量使用。其实,这个数的值是33234033323403.
而我使用的是第三种做法——分类讨论。(其实是因为懒得写逆元)
观察到当n≡0n≡0或n≡2(mod3)n≡2(mod3)时,n(n+1)n(n+1)的值可以被6整除,因此先算出n(n+1)6n(n+1)6即可。当n≡1(mod3)n≡1(mod3)时,2n+12n+1被3整除,n(n+1)n(n+1)被2整除,故计算2n+132n+13与n(n+1)2n(n+1)2相乘即可。(当然如果模数是24而不是6我甘愿老老实实地求逆元)
注意运算符优先级问题,不可以有任何连着三个19940417或1e9级别的数不加取模地连乘,否则必爆。
代码实现
题目大意: 求∑i=1n∑j=1,j≠im(nmodi)(mmodj)∑i=1n∑j=1,j≠im(nmodi)(mmodj)对19940417取模的值。
思路分析:
从heheda神犇的博客“heheda的数论专题练习”上翻到这样一道比较有趣的题。
一开始想歪了,一直在考虑n mod i与n mod (i+1)之间的递推关系,后来发现行不通。
既然是取模,那就可以采用“化模为除”的方法,即利用公式nmodi=n−[n/i]∗inmodi=n−[n/i]∗i化简原式,暂且不考虑i≠ji≠j的条件可得ans=∑i=1n∑j=1m(n−i[n/i])(m−j[m/j])=∑i=1n∑j=1m(nm−nj[m/j]−mi[n/i]+ij[n/i][m/j])=n2m2−n2∑j=1mj[m/j]−m2∑i=1ni[m/i]+∑i=1ni[n/i]∗∑j=1mj[m/j]ans=∑i=1n∑j=1m(n−i[n/i])(m−j[m/j])=∑i=1n∑j=1m(nm−nj[m/j]−mi[n/i]+ij[n/i][m/j])=n2m2−n2∑j=1mj[m/j]−m2∑i=1ni[m/i]+∑i=1ni[n/i]∗∑j=1mj[m/j]
做到这一步,既然只有i和[n/i]出现在式子里(j同理),显然可以用数论分块加速这个过程了。(至于数论分块的实现可以参见我的另一篇博客http://blog.csdn.net/suncongbo/article/details/78819470) 分成(2n−−√+2m−−√)(2n+2m)个块,每个块内的[n/i][n/i]与[m/i][m/i]值均不变,块内求和直接用求和公式getsum(n)=n(n+1)2getsum(n)=n(n+1)2, 块内求和等于块右端点的getsum返回值减去块左端点的getsum返回值。
然后去掉i=ji=j的情况:∑i=1min(n,m)(n−i[n/i])(m−i[m/i])=∑i=1min(n,m)(nm−ni[m/i]−mi[n/i]+i2[n/i][m/i])=nmmin(n,m)−∑i=1min(n,m)i(n[m/i]+m[n/i])+∑i=1min(n,m)i2[n/i][m/i]∑i=1min(n,m)(n−i[n/i])(m−i[m/i])=∑i=1min(n,m)(nm−ni[m/i]−mi[n/i]+i2[n/i][m/i])=nmmin(n,m)−∑i=1min(n,m)i(n[m/i]+m[n/i])+∑i=1min(n,m)i2[n/i][m/i], 前两项仍然可以用前面的方法进行计算。棘手的问题来了:
对于最后一项,牵扯到求连续自然数的平方和。我们知道sqrsum(i)=∑i=1ni2=n(n+1)(2n+1)6sqrsum(i)=∑i=1ni2=n(n+1)(2n+1)6, 但是如果把n(n+1)(2n+1)n(n+1)(2n+1)算出来将会很大,long long无法支持。因此只能中途对1994041719940417取模。但是又要除以6,牵扯到求逆元,而19940417又不是质数,无法用费马小定理求逆元。此时做法一是通过其他方法(如exgcd)求出6对于1994041719940417的逆元,二是直接算出6的逆元并作为常量使用。其实,这个数的值是33234033323403.
而我使用的是第三种做法——分类讨论。(其实是因为懒得写逆元)
观察到当n≡0n≡0或n≡2(mod3)n≡2(mod3)时,n(n+1)n(n+1)的值可以被6整除,因此先算出n(n+1)6n(n+1)6即可。当n≡1(mod3)n≡1(mod3)时,2n+12n+1被3整除,n(n+1)n(n+1)被2整除,故计算2n+132n+13与n(n+1)2n(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; }
相关文章推荐
- bzoj 2956: 模积和 分块+数学
- BZOJ3505 & 洛谷P3166 [Cqoi2014]数三角形 【数学、数论】
- 【BZOJ1257】余数之和(数论分块,暴力)
- BZOJ 4403: 序列统计 (组合数 Lucas 数论推导)
- Bzoj3834:[Poi2014]Solar Panels:数论,分块
- BZOJ 2142 礼物 组合数学+数论
- 【bzoj 3209】花神的数论题(组合数学)
- 数学(数论)BZOJ 3309:DZY Loves Math
- 【BZOJ】4804 欧拉心算 莫比乌斯函数+欧拉函数+数论分块
- BZOJ 1045 HAOI 2008 糖果传递 数学推导
- Bzoj 4403: 序列统计 Lucas定理,组合数学,数论
- BZOJ 4173 数学 数论
- [BZOJ4173]数学(数论)
- BZOJ 2956 模积和 (分块加速)
- 洛谷 P3317 [SDOI2014]重建(矩阵树定理+数学推导) [bzoj3534]
- bzoj2154: Crash的数字表格/2693: jzptab [莫比乌斯反演、数论推导]
- [BZOJ1432][ZJOI2009]Funtion(数学推导)
- bzoj 2956 数学展开,分段处理
- 【bzoj3994】[SDOI2015]约数个数和 线性筛法+莫比乌斯反演+数论分块
- BZOJ_2956_模积和_数学