数论分块是一种非常重要的思想。就是对于一些表达式,它的值只有sqrt(n)种,那么我们就对于这sqrt(n)种数值进行分块,然后暴力算即可。我们通过下面一道例题来了解数论分块:BZOJ1257
Description
给出正整数n和k,计算j(n, k)=k mod 1 + k mod 2 + k mod 3 + … + k mod n的值
其中k mod i表示k除以i的余数。
例如j(5, 3)=3 mod 1 + 3 mod 2 + 3 mod 3 + 3 mod 4 + 3 mod 5=0+1+0+3+3=7
Input
输入仅一行,包含两个整数n, k。
1<=n ,k<=10^9
Output
输出仅一行,即j(n, k)。
Sample Input
5 3
Sample Output
7
题目要求的就是∑(n mod i)(i∈[1,k]),我们可以先进行分类讨论(k大于n的情况),对于i∈[n+1,k]答案就是
(k-(n+1)+1)*n,而对于区间[1,n],我们就可以进行数论分块。(n mod i)
可以表示为(n-(n/i)*i)
,而(n/i)
总共只有2*sqrt(n)种取值,而且这些取值是连续的,从i到n/(n/l)
这段区间内(n/i)
都是相同的。(大家可以思考一下为什么),所以这道题我们就可以统计[i,n/(n/i)],将O(n)的复杂度变为O(sqrt(n))。
#include<bits/stdc++.h>
#define ll long long
using namespace std;
ll read(){
char c;ll x=0,y=1;while(c=getchar(),(c<'0'||c>'9')&&c!='-');
if(c=='-') y=-1;else x=c-'0';while(c=getchar(),c>='0'&&c<='9')
x=x*10+c-'0';return x*y;
}
ll n,k,l,ans,tot,r;
int main()
{
n=read();k=read();l=1;
register int i;
while(l<=n){
if(k<l) break;
r=k/(k/l);r=min(r,n);
ll res=0;
tot+=k*(r-l+1);ans+=(k/l)*(l+r)*(r-l+1)/2; //(k/l)*(l+r)*(r-l+1)/2相当于(k/l)*∑i (i∈[l,r])
l=r+1;
}
if(l<=n) tot+=(n-l+1)*k;
printf("%lld",tot-ans);
return 0;
}
然后我们继续扩展,现在来到了二维:BZOJ2956
Description
求∑∑((n mod i)*(m mod j))其中1<=i<=n,1<=j<=m,i≠j。
Input
第一行两个数n,m。
Output
一个整数表示答案mod 19940417的值
Sample Input
3 4
Sample Output
1
样例说明
答案为(3 mod 1)(4 mod 2)+(3 mod 1) (4 mod 3)+(3 mod 1) * (4 mod 4) + (3 mod 2) * (4 mod 1) + (3 mod 2) * (4 mod 3) + (3 mod 2) * (4 mod 4) + (3 mod 3) * (4 mod 1) + (3 mod 3) * (4 mod 2) + (3 mod 3) * (4 mod 4) = 1
数据规模和约定
对于100%的数据n,m<=10^9。
这题要求的是 ∑ni=1∑mj=1(n mod i)(m mod i)−∑min(n,m)i=1(n mod i)(m mod i) ∑ i = 1 n ∑ j = 1 m ( n m o d i ) ( m m o d i ) − ∑ i = 1 m i n ( n , m ) ( n m o d i ) ( m m o d i )
那我们前面部分可以直接用数论分块做,但是关于 ∑min(n,m)i=1(n mod i)(m mod i) ∑ i = 1 m i n ( n , m ) ( n m o d i ) ( m m o d i ) 我们需要推一下式子
∑min(n,m)i=1(n mod i)(m mod i)=∑min(n,m)i=1(n−(n/i)∗i)(m−(m/i)∗i) ∑ i = 1 m i n ( n , m ) ( n m o d i ) ( m m o d i ) = ∑ i = 1 m i n ( n , m ) ( n − ( n / i ) ∗ i ) ( m − ( m / i ) ∗ i )
=∑min(n,m)i=1(nm−ni(m/i)−mi(n/i)+i2(n/i)(m/i)) = ∑ i = 1 m i n ( n , m ) ( n m − n i ( m / i ) − m i ( n / i ) + i 2 ( n / i ) ( m / i ) )
化成这样之后我们就可以进行数论分块了,我们依旧分区间l,r,只不过这次算r不一样r=min(n/(n/l),m/(m/l));
然后我们在数论分块种求出这个表达式的值即可。
#include<bits/stdc++.h>
#define MD 19940417
#define ll long long
#define rev 3323403 //6关于MD的逆元
using namespace std;
ll read(){
char c;ll x;while(c=getchar(),c<'0'||c>'9');x=c-'0';
while(c=getchar(),c>='0'&&c<='9') x=x*10+c-'0';return x;
}
ll n,m,l,r,p,totn,totm,cut;
ll pows(ll a,ll b){
ll base=1;
while(b){
if(b&1) base=base*a%MD;
a=a*a%MD;b/=2;
}
return base;
}
ll calc(ll x){
return x*(x+1)%MD*(2*x+1)%MD*rev%MD;
}
int main()
{
n=read();m=read();l=1;p=min(n,m);
while(l<=n){
r=n/(n/l);
totn=(totn+(r-l+1)*n%MD-(n/l)*(l+r)*(r-l+1)/2+MD)%MD;
l=r+1;
}
l=1;
while(l<=m){
r=m/(m/l);
totm=(totm+(r-l+1)*m%MD-(m/l)*(l+r)*(r-l+1)/2+MD)%MD;
l=r+1;
}
l=1;
while(l<=p){
r=min(n/(n/l),m/(m/l));
cut=(cut+n*m%MD*(r-l+1)%MD-((n/l)*m+(m/l)*n)%MD*((r-l+1)*(l+r)/2%MD)%MD+(n/l)*(m/l)%MD*(calc(r)-calc(l-1)+MD)%MD+MD*5)%MD; //核心代码
l=r+1;
}
printf("%lld",(totn*totm%MD-cut+MD)%MD);
return 0;
}