T i p s Tips Tips:本文默认 n ≤ m , a ≤ b n\leq m,\: a\leq b n≤m,a≤b
[POI2007] ZAP-Queries
(题目传送门)
题意
给出 a , b , d a,b,d a,b,d,求满足 1 ≤ x ≤ a 1 \leq x \leq a 1≤x≤a, 1 ≤ y ≤ b 1 \leq y \leq b 1≤y≤b,且 gcd ( x , y ) = d \gcd(x,y)=d gcd(x,y)=d 的二元组 ( x , y ) (x,y) (x,y) 的数量。
1 ≤ n ≤ 5 × 1 0 4 1 \leq n \leq 5 \times 10^4 1≤n≤5×104, 1 ≤ d ≤ a , b ≤ 5 × 1 0 4 1 \leq d \leq a,b \leq 5 \times 10^4 1≤d≤a,b≤5×104。
解题思路
-
原题要求的是
∑ i = 1 a ∑ j = 1 b [ gcd ( i , j ) = d ] \sum\limits_{i=1}^a\sum\limits_{j=1}^b{[\gcd(i,j)=d]} i=1∑aj=1∑b[gcd(i,j)=d] -
将 i , j i,j i,j 除以 d d d 得
= ∑ i = 1 ⌊ a d ⌋ ∑ j = 1 ⌊ b d ⌋ [ gcd ( i , j ) = 1 ] =\sum\limits_{i=1}^{\left\lfloor\frac{a}{d}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{b}{d}\right\rfloor}{[\gcd(i,j)=1]} =i=1∑⌊da⌋j=1∑⌊db⌋[gcd(i,j)=1] -
由莫比乌斯反演得
= ∑ i = 1 ⌊ a d ⌋ ∑ j = 1 ⌊ b d ⌋ ∑ k ∣ gcd ( i , j ) μ ( k ) =\sum\limits_{i=1}^{\left\lfloor\frac{a}{d}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{b}{d}\right\rfloor}\sum\limits_{k|\gcd(i,j)}{\mu(k)} =i=1∑⌊da⌋j=1∑⌊db⌋k∣gcd(i,j)∑μ(k) -
再将 i , j i,j i,j 除以 k k k,这样就不关心 k k k 是否为 i , j i,j i,j 的约数,可以直接枚举 k k k :
= ∑ k = 1 ⌊ a d ⌋ μ ( k ) ⌊ a k d ⌋ ⌊ b k d ⌋ =\sum\limits_{k=1}^{\left\lfloor\frac{a}{d}\right\rfloor}\mu(k){\left\lfloor\frac{a}{kd}\right\rfloor}{\left\lfloor\frac{b}{kd}\right\rfloor} =k=1∑⌊da⌋μ(k)⌊kda⌋⌊kdb⌋
-
这个式子的后半部分可以整除分块来优化,前半部分预处理前缀和即可。在代码实现时,可以预先让 a , b a,b a,b 同除以 d d d,方便实现
代码
#include<bits/stdc++.h>
#define LL long long
using namespace std;
const int N=50010;
int T,a,b,d,n,m;
int prime[N],v[N],tot,mu[N],sum[N];
LL ans;
void primes()
{
mu[1]=1;
for(int i=2; i<=50000; i++)
{
if(!v[i])
{
v[i]=i;
prime[++tot]=i;
mu[i]=-1;
}
for(int j=1; j<=tot; j++)
{
if(prime[j]>v[i] || prime[j]>50000/i)
break;
v[i*prime[j]]=prime[j];
if(i%prime[j]==0)
break;
mu[prime[j]*i]=-mu[i];
}
}
for(int i=1; i<=50000; i++)
sum[i]=sum[i-1]+mu[i];
}
int main()
{
primes();
scanf("%d",&T);
while(T--)
{
scanf("%d%d%d",&a,&b,&d);
n=a/d; m=b/d;
if(n>m)
swap(n,m);
ans=0;
for(int l=1,r; l<=n; l=r+1)
{
r=min(n/(n/l),m/(m/l));
ans+=1LL*(sum[r]-sum[l-1])*1LL*(n/l)*1LL*(m/l);
}
printf("%lld\n",ans);
}
return 0;
}
[HAOI2011] Problem b
(题目传送门)
题意
对于给出的 n n n 个询问,每次求有多少个数对 ( x , y ) (x,y) (x,y),满足 a ≤ x ≤ b a \le x \le b a≤x≤b, c ≤ y ≤ d c \le y \le d c≤y≤d,且 gcd ( x , y ) = k \gcd(x,y) = k gcd(x,y)=k, gcd ( x , y ) \gcd(x,y) gcd(x,y) 函数为 x x x 和 y y y 的最大公约数。
1 ≤ n , k ≤ 5 × 1 0 4 1 \le n,k \le 5 \times 10^4 1≤n,k≤5×104, 1 ≤ a ≤ b ≤ 5 × 1 0 4 1 \le a \le b \le 5 \times 10^4 1≤a≤b≤5×104, 1 ≤ c ≤ d ≤ 5 × 1 0 4 1 \le c \le d \le 5 \times 10^4 1≤c≤d≤5×104。
解题思路
-
类比上面一题
ZAP
,设 f ( n , m ) f(n,m) f(n,m) 表示满足 1 ≤ x ≤ n 1 \leq x \leq n 1≤x≤n, 1 ≤ y ≤ m 1 \leq y \leq m 1≤y≤m,且 gcd ( x , y ) = k \gcd(x,y)=k gcd(x,y)=k 的二元组 ( x , y ) (x,y) (x,y) 的数量 -
类似二维前缀和,容易得出
a n s = f ( b , d ) − f ( a − 1 , d ) − f ( b , c − 1 ) + f ( a − 1 , c − 1 ) ans=f(b,d)-f(a-1,d)-f(b,c-1)+f(a-1,c-1) ans=f(b,d)−f(a−1,d)−f(b,c−1)+f(a−1,c−1) -
同样,可以一开始先令 a , b , c , d a,b,c,d a,b,c,d 除以 k k k,方便代码实现
代码
#include<bits/stdc++.h>
#define LL long long
using namespace std;
const int N=50010;
int T,a,b,c,d,k;
int prime[N],v[N],tot,mu[N],sum[N];
LL ans;
void primes()
{
mu[1]=1;
for(int i=2; i<=50000; i++)
{
if(!v[i])
{
v[i]=i;
prime[++tot]=i;
mu[i]=-1;
}
for(int j=1; j<=tot; j++)
{
if(prime[j]>v[i] || prime[j]>50000/i)
break;
v[i*prime[j]]=prime[j];
if(i%prime[j]==0)
break;
mu[prime[j]*i]=-mu[i];
}
}
for(int i=1; i<=50000; i++)
sum[i]=sum[i-1]+mu[i];
}
LL f(int n,int m)
{
if(n>m)
swap(n,m);
LL res=0;
for(int l=1,r; l<=n; l=r+1)
{
r=min(n/(n/l),m/(m/l));
res+=1LL*(sum[r]-sum[l-1])*1LL*(n/l)*1LL*(m/l);
}
return res;
}
int main()
{
primes();
scanf("%d",&T);
while(T--)
{
scanf("%d%d%d%d%d",&a,&b,&c,&d,&k);
ans=0;
ans=f(b/k,d/k)-f((a-1)/k,d/k)-f(b/k,(c-1)/k)+f((a-1)/k,(c-1)/k);
printf("%lld\n",ans);
}
return 0;
}
YY的GCD
(题目传送门)
题意
T T T 组数据,每组数据给定 N , M N, M N,M,求 1 ≤ x ≤ N 1 \leq x \leq N 1≤x≤N, 1 ≤ y ≤ M 1 \leq y \leq M 1≤y≤M 且 gcd ( x , y ) \gcd(x, y) gcd(x,y) 为质数的 ( x , y ) (x, y) (x,y) 有多少对。
T = 1 0 4 T = 10^4 T=104, N , M ≤ 1 0 7 N, M \leq 10^7 N,M≤107。
解题思路
-
设 p r i m e s primes primes 表示质数集合
-
原题要求的是
a n s = ∑ k ∈ p r i m e s ∑ i = 1 n ∑ j = 1 m [ gcd ( i , j ) = k ] ans=\sum\limits_{k\in primes}\sum\limits_{i=1}^{n}\sum\limits_{j=1}^m{[\gcd(i,j)=k]} ans=k∈primes∑i=1∑nj=1∑m[gcd(i,j)=k] -
同除以 k k k 得:
= ∑ k ∈ p r i m e s ∑ i = 1 ⌊ n k ⌋ ∑ j = 1 ⌊ m k ⌋ [ gcd ( i , j ) = k ] =\sum\limits_{k\in primes}\sum\limits_{i=1}^{\left\lfloor\frac{n}{k}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{m}{k}\right\rfloor}{[\gcd(i,j)=k]} =k∈primes∑i=1∑⌊kn⌋j=1∑⌊km⌋[gcd(i,j)=k] -
由莫比乌斯反演得
= ∑ k ∈ p r i m e s ∑ i = 1 ⌊ n k ⌋ ∑ j = 1 ⌊ m k ⌋ ∑ d ∣ gcd ( i , j ) μ ( d ) =\sum\limits_{k\in primes}\sum\limits_{i=1}^{\left\lfloor\frac{n}{k}\right\rfloor}\sum\limits_{j=1}^{\left\lfloor\frac{m}{k}\right\rfloor}\sum\limits_{d|\gcd(i,j)}{\mu(d)} =k∈primes∑i=1∑⌊kn⌋j=1∑⌊km⌋d∣gcd(i,j)∑μ(d) -
我们枚举 d d d 得:
= ∑ k ∈ p r i m e s ∑ d = 1 ⌊ n k ⌋ μ ( d ) ⌊ n k d ⌋ ⌊ m k d ⌋ =\sum\limits_{k\in primes}\sum\limits_{d=1}^{\left\lfloor\frac{n}{k}\right\rfloor}{\mu(d)}\left\lfloor\frac{n}{kd}\right\rfloor\left\lfloor\frac{m}{kd}\right\rfloor =k∈primes∑d=1∑⌊kn⌋μ(d)⌊kdn⌋⌊kdm⌋ -
设 T = k d T=kd T=kd,把 T T T 提到前面去
= ∑ T = 1 n ⌊ n T ⌋ ⌊ m T ⌋ ∑ k ∣ T , k ∈ p r i m e s μ ( T k ) =\sum\limits_{T=1}^n\left\lfloor\frac{n}{T}\right\rfloor\left\lfloor\frac{m}{T}\right\rfloor\sum\limits_{k|T,k\in primes}{\mu(\frac{T}{k})} =T=1∑n⌊Tn⌋⌊Tm⌋k∣T,k∈primes∑μ(kT) -
对于后半部分式子,我们可以在线性筛后预处理出来,再做一次整除分块即可求出答案
代码
#include<bits/stdc++.h>
#define LL long long
using namespace std;
const int N=10000010;
int T,n,m;
int prime[N],v[N],tot,mu[N];
LL sum[N],g[N];
void primes()
{
mu[1]=1;
for(int i=2; i<=1e7; i++)
{
if(!v[i])
{
v[i]=i;
prime[++tot]=i;
mu[i]=-1;
}
for(int j=1; j<=tot; j++)
{
if(prime[j]>v[i] || prime[j]>1e7/i)
break;
v[i*prime[j]]=prime[j];
if(i%prime[j]==0)
break;
mu[prime[j]*i]=-mu[i];
}
}
for(int i=1; i<=tot; i++)
{
for(int j=1; j*prime[i]<=1e7; j++)
g[prime[i]*j]+=mu[j];
}
for(int i=1; i<=1e7; i++)
sum[i]=sum[i-1]+g[i];
}
LL cacl(int n,int m)
{
if(n>m)
swap(n,m);
LL res=0;
for(int l=1,r; l<=n; l=r+1)
{
r=min(n/(n/l),m/(m/l));
res+=1LL*(sum[r]-sum[l-1])*(n/l)*(m/l);
}
return res;
}
int main()
{
primes();
scanf("%d",&T);
while(T--)
{
scanf("%d%d",&n,&m);
printf("%lld\n",cacl(n,m));
}
return 0;
}
于神之怒加强版
(题目传送门)
题意
给定 n , m , k n,m,k n,m,k,计算
∑ i = 1 n ∑ j = 1 m gcd ( i , j ) k \sum_{i=1}^n \sum_{j=1}^m \gcd(i,j)^k i=1∑nj=1∑mgcd(i,j)k
对 1 0 9 + 7 10^9 + 7 109+7 取模的结果。
对于全部的测试点,保证 1 ≤ T ≤ 2 × 1 0 3 1 \leq T \leq 2 \times 10^3 1≤T≤2×103, 1 ≤ n , m , k ≤ 5 × 1 0 6 1 \leq n, m, k \leq 5 \times 10^6 1≤n,m,k≤5×106。
解题思路
-
直接开始推式子
a n s = ∑ d = 1 n ∑ i = 1 n ∑ j = 1 m d k [ gcd ( i , j ) = d ] ans=\sum_{d=1}^n\sum_{i=1}^n\sum_{j=1}^m{d^k[\gcd(i,j)=d]} ans=d=1∑ni=1∑nj=1∑mdk[gcd(i,j)=d] -
除以 d d d 得:
= ∑ d = 1 n ∑ i = 1 ⌊ n d ⌋ ∑ j = 1 ⌊ m d ⌋ d k [ gcd ( i , j ) = 1 ] =\sum_{d=1}^n\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor}{d^k[\gcd(i,j)=1]} =d=1∑ni=1∑⌊dn⌋j=1∑⌊dm⌋dk[gcd(i,j)=1]
= ∑ d = 1 n d k ∑ i = 1 ⌊ n d ⌋ ∑ j = 1 ⌊ m d ⌋ ∑ x ∣ gcd ( i , j ) μ ( x ) =\sum_{d=1}^n{d^k}\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor}\sum_{x|\gcd(i,j)}{\mu(x)} =d=1∑ndki=1∑⌊dn⌋j=1∑⌊dm⌋x∣gcd(i,j)∑μ(x) -
枚举 x x x 得:
∑ d = 1 n d k ∑ x = 1 ⌊ n d ⌋ μ ( x ) ⌊ n d x ⌋ ⌊ m d x ⌋ \sum_{d=1}^n{d^k}\sum_{x=1}^{\left\lfloor\frac{n}{d}\right\rfloor}{\mu(x)\left\lfloor\frac{n}{dx}\right\rfloor\left\lfloor\frac{m}{dx}\right\rfloor} d=1∑ndkx=1∑⌊dn⌋μ(x)⌊dxn⌋⌊dxm⌋ -
类似上一题的优化,设 T = d x T=dx T=dx 得:
∑ T = 1 n ⌊ n T ⌋ ⌊ m T ⌋ ∑ d ∣ T d k μ ( T d ) \sum_{T=1}^n{\left\lfloor\frac{n}{T}\right\rfloor\left\lfloor\frac{m}{T}\right\rfloor}\sum_{d|T}{d^k\mu(\frac{T}{d})} T=1∑n⌊Tn⌋⌊Tm⌋d∣T∑dkμ(dT) -
观察到, ∑ d ∣ T d k μ ( T d ) \sum\limits_{d|T}{d^k\mu(\frac{T}{d})} d∣T∑dkμ(dT) 是狄利克雷卷积的形式,而 f ( d ) = d k f(d)=d^k f(d)=dk 和 μ ( T d ) \mu(\frac{T}{d}) μ(dT) 又是 积性函数,所以 g ( T ) = ∑ d ∣ T d k μ ( T d ) g(T)=\sum\limits_{d|T}{d^k\mu(\frac{T}{d})} g(T)=d∣T∑dkμ(dT) 也是积性函数。因此,我们可以考虑线性筛预处理出 g g g 的值
-
设 T T T 质因数分解为 T = ∏ i = 1 q p i c i T=\prod\limits_{i=1}^q{p_i^{c_i}} T=i=1∏qpici,则
g ( T ) = ∏ i = 1 q g ( p i c i ) g(T)=\prod\limits_{i=1}^q{g(p_i^{c_i})} g(T)=i=1∏qg(pici)
= ∏ i = 1 q ∑ d ∣ p i c i d k μ ( p i c i d ) =\prod\limits_{i=1}^q\sum_{d|p_i^{c_i}}{d^k\mu(\frac{p_i^{c_i}}{d})} =i=1∏qd∣pici∑dkμ(dpici) -
观察到,当 d = 1 , p i , p i 2 … p i c i − 2 d=1,p_i,pi^2\dots p_i^{c_i-2} d=1,pi,pi2…pici−2 时, μ ( p i c i d ) = 0 \mu(\frac{p_i^{c_i}}{d})=0 μ(dpici)=0,所以
g ( T ) = ∏ i = 1 q p i ( c i − 1 ) k × μ ( p i ) + p i c i k × μ ( 1 ) g(T)=\prod\limits_{i=1}^q{p_i^{(c_i-1)k}\times\mu(p_i)}+p_i^{c_ik}\times \mu(1) g(T)=i=1∏qpi(ci−1)k×μ(pi)+picik×μ(1)
= ∏ i = 1 q p i ( c i − 1 ) k × ( p i k − 1 ) =\prod\limits_{i=1}^q{p_i^{(c_i-1)k}\times(p_i^k-1)} =i=1∏qpi(ci−1)k×(pik−1) -
这个式子是可以预处理出来的
-
所以
a n s = ∑ T = 1 n ⌊ n T ⌋ ⌊ m T ⌋ g ( T ) ans=\sum_{T=1}^n{\left\lfloor\frac{n}{T}\right\rfloor\left\lfloor\frac{m}{T}\right\rfloor}g(T) ans=T=1∑n⌊Tn⌋⌊Tm⌋g(T) -
利用整除分块即可求解
代码
#include<bits/stdc++.h>
#define LL long long
using namespace std;
const int N=5000010;
const int MOD=1e9+7;
int T,k,n,m;
int prime[N],v[N],tot;
LL g[N],mi[N],sum[N];
LL ksm(int a,int b)
{
if(b==0)
return 1;
if(b==1)
return a;
LL tmp=ksm(a,b/2);
if(b%2==0)
return (tmp*tmp)%MOD;
else
return ((tmp*tmp)%MOD*(LL)a)%MOD;
}
void primes()
{
g[1]=1;
for(int i=2; i<=5e6; i++)
{
if(!v[i])
{
v[i]=i;
prime[++tot]=i;
mi[i]=ksm(i,k);
g[i]=(mi[i]-1+MOD)%MOD; //i为质数时g(i)=i^k-1
}
for(int j=1; j<=tot; j++)
{
if(prime[j]>v[i] || prime[j]>5e6/i)
break;
v[i*prime[j]]=prime[j];
if(i%prime[j]==0)
g[i*prime[j]]=(g[i]*mi[prime[j]])%MOD; //当i和prime[j]不互质时,读者可自行证明该结论正确性
else
g[i*prime[j]]=(g[i]*g[prime[j]])%MOD;
}
}
for(int i=1; i<=5e6; i++)
sum[i]=sum[i-1]+g[i];
}
LL cacl(int n,int m)
{
if(n>m)
swap(n,m);
LL res=0;
for(int l=1,r; l<=n; l=r+1)
{
r=min(n/(n/l),m/(m/l));
res=(res+((sum[r]-sum[l-1]+MOD)%MOD)*(n/l)%MOD*(m/l)%MOD)%MOD;
}
return res;
}
int main()
{
scanf("%d%d",&T,&k);
primes();
while(T--)
{
scanf("%d%d",&n,&m);
printf("%lld\n",cacl(n,m));
}
return 0;
}
[国家集训队] Crash的数字表格 / JZPTAB
(题目传送门)
题意
给定
n
,
m
n,m
n,m,计算
∑
i
=
1
n
∑
j
=
1
m
lcm
(
i
,
j
)
\sum_{i=1}^n\sum_{j=1}^m{\operatorname{lcm}(i,j)}
i=1∑nj=1∑mlcm(i,j)
1
≤
n
,
m
≤
1
0
7
1\leq n,m \leq 10^7
1≤n,m≤107
解题思路
-
对 lcm ( i , j ) \operatorname{lcm}(i,j) lcm(i,j) 变换一下得:
a n s = ∑ i = 1 n ∑ j = 1 m i j gcd ( i , j ) ans=\sum_{i=1}^n\sum_{j=1}^m{\frac{ij}{\gcd(i,j)}} ans=i=1∑nj=1∑mgcd(i,j)ij -
枚举 gcd ( i , j ) \gcd(i,j) gcd(i,j) 得:
= ∑ i = 1 n ∑ j = 1 m ∑ d = 1 n i j d [ gcd ( i , j ) = d ] =\sum_{i=1}^n\sum_{j=1}^m\sum_{d=1}^n{\frac{ij}{d}[\gcd(i,j)=d]} =i=1∑nj=1∑md=1∑ndij[gcd(i,j)=d] -
除以 d d d 得:
= ∑ d = 1 n ∑ i = 1 ⌊ n d ⌋ ∑ j = 1 ⌊ m d ⌋ i j × d [ gcd ( i , j ) = 1 ] =\sum_{d=1}^n\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor}{ij\times d[\gcd(i,j)=1]} =d=1∑ni=1∑⌊dn⌋j=1∑⌊dm⌋ij×d[gcd(i,j)=1]
= ∑ d = 1 n d ∑ i = 1 ⌊ n d ⌋ ∑ j = 1 ⌊ m d ⌋ i j [ gcd ( i , j ) = 1 ] =\sum_{d=1}^n{d}\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor}{ij[\gcd(i,j)=1]} =d=1∑ndi=1∑⌊dn⌋j=1∑⌊dm⌋ij[gcd(i,j)=1] -
把后半部分式子单独提出来,设:
f ( a , b ) = ∑ i = 1 a ∑ j = 1 b i j [ gcd ( i , j ) = 1 ] f(a,b)=\sum\limits_{i=1}^{a}\sum\limits_{j=1}^{b}{ij[\gcd(i,j)=1]} f(a,b)=i=1∑aj=1∑bij[gcd(i,j)=1] -
化简得:
f ( a , b ) = ∑ i = 1 a ∑ j = 1 b i j ∑ k ∣ gcd ( i , j ) μ ( k ) f(a,b)=\sum\limits_{i=1}^{a}\sum\limits_{j=1}^{b}{ij\sum_{k|\gcd(i,j)}{\mu(k)}} f(a,b)=i=1∑aj=1∑bijk∣gcd(i,j)∑μ(k)
= ∑ k = 1 a μ ( k ) ∑ i = 1 ⌊ a k ⌋ ∑ j = 1 ⌊ b k ⌋ i j × k 2 =\sum_{k=1}^{a}{\mu(k)}\sum_{i=1}^{\left\lfloor\frac{a}{k}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{b}{k}\right\rfloor}{ij\times k^2} =k=1∑aμ(k)i=1∑⌊ka⌋j=1∑⌊kb⌋ij×k2
= ∑ k = 1 a μ ( k ) × k 2 ∑ i = 1 ⌊ a k ⌋ ∑ j = 1 ⌊ b k ⌋ i j =\sum_{k=1}^{a}{\mu(k)\times k^2}\sum_{i=1}^{\left\lfloor\frac{a}{k}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{b}{k}\right\rfloor}{ij} =k=1∑aμ(k)×k2i=1∑⌊ka⌋j=1∑⌊kb⌋ij -
设
g ( a , b ) = ∑ i = 1 a ∑ j = 1 b i j g(a,b)=\sum_{i=1}^{a}\sum_{j=1}^{b}{ij} g(a,b)=i=1∑aj=1∑bij -
化简得:
= ( 1 + a ) a 2 × ( b + 1 ) b 2 =\frac{(1+a)a}{2}\times \frac{(b+1)b}{2} =2(1+a)a×2(b+1)b -
所以,求出 g ( a , b ) g(a,b) g(a,b) 的时间复杂度为 O ( 1 ) O(1) O(1)
-
将 g ( a , b ) g(a,b) g(a,b) 代回 f ( a , b ) f(a,b) f(a,b) 中得:
f ( a , b ) = ∑ k = 1 a μ ( k ) × k 2 × g ( ⌊ a k ⌋ , ⌊ b k ⌋ ) f(a,b)=\sum_{k=1}^{a}{\mu(k)\times k^2\times g(\left\lfloor\frac{a}{k}\right\rfloor,\left\lfloor\frac{b}{k}\right\rfloor)} f(a,b)=k=1∑aμ(k)×k2×g(⌊ka⌋,⌊kb⌋) -
我们利用整除分块即可求出 f ( a , b ) f(a,b) f(a,b)
-
将 f ( a , b ) f(a,b) f(a,b) 代回 a n s ans ans 得:
a n s = ∑ d = 1 n d × f ( ⌊ n d ⌋ , ⌊ m d ⌋ ) ans=\sum_{d=1}^n{d}\times f(\left\lfloor\frac{n}{d}\right\rfloor,\left\lfloor\frac{m}{d}\right\rfloor) ans=d=1∑nd×f(⌊dn⌋,⌊dm⌋) -
再做一次整除分块即可求解
代码
#include<bits/stdc++.h>
#define LL long long
using namespace std;
const int N=10000010;
const int MOD=20101009;
int n,m;
int prime[N],v[N],tot,mu[N];
LL sum[N];
void primes()
{
mu[1]=1;
for(int i=2; i<=1e7; i++)
{
if(!v[i])
{
v[i]=i;
prime[++tot]=i;
mu[i]=-1;
}
for(int j=1; j<=tot; j++)
{
if(prime[j]>v[i] || prime[j]>1e7/i)
break;
v[i*prime[j]]=prime[j];
if(i%prime[j]==0)
break;
mu[i*prime[j]]=-mu[i];
}
}
for(int i=1; i<=1e7; i++)
sum[i]=(sum[i-1]+1LL*(mu[i]+MOD)%MOD*i%MOD*i%MOD)%MOD;
}
LL g(int a,int b)
{
return (1LL*a*(a+1)/2%MOD)*(1LL*b*(b+1)/2%MOD)%MOD;
}
LL f(int a,int b)
{
if(a>b)
swap(a,b);
LL res=0;
for(int l=1,r; l<=a; l=r+1)
{
r=min(a/(a/l),b/(b/l));
res=(res+(sum[r]-sum[l-1]+MOD)%MOD*g(a/l,b/l)%MOD)%MOD;
}
return res;
}
LL calc(int a,int b)
{
if(a>b)
swap(a,b);
LL res=0;
for(int l=1,r; l<=a; l=r+1)
{
r=min(a/(a/l),b/(b/l));
res=(res+1LL*(r-l+1)*(l+r)/2%MOD*f(a/l,b/l)%MOD)%MOD;
}
return res;
}
int main()
{
primes();
scanf("%d%d",&n,&m);
printf("%lld",calc(n,m));
return 0;
}
[SDOI2014] 数表
(题目传送门)
题意
有一张 n × m n\times m n×m 的数表,其第 i i i 行第 j j j 列( 1 ≤ i ≤ n 1\le i\le n 1≤i≤n, 1 ≤ j ≤ m 1\le j\le m 1≤j≤m)的数值为能同时整除 i i i 和 j j j 的所有自然数之和。给定 a a a,计算数表中不大于 a a a 的数之和。
1 ≤ n , m ≤ 1 0 5 1\le n,m\le 10^5 1≤n,m≤105, 1 ≤ Q ≤ 2 × 1 0 4 1\le Q\le 2\times 10^4 1≤Q≤2×104。
解题思路
-
设 f ( x ) f(x) f(x) 表示 x x x 的因数之和
-
则
a n s = ∑ i = 1 n ∑ j = 1 m f ( gcd ( i , j ) ) × [ f ( gcd ( i , j ) ) ≤ a ] ans=\sum_{i=1}^n\sum_{j=1}^m{f(\gcd(i,j))\times[f(\gcd(i,j))\leq a]} ans=i=1∑nj=1∑mf(gcd(i,j))×[f(gcd(i,j))≤a] -
枚举 gcd ( i , j ) \gcd(i,j) gcd(i,j) 得:
= ∑ d = 1 n ∑ i = 1 n ∑ j = 1 m f ( d ) × [ gcd ( i , j ) = d ] × [ f ( d ) ≤ a ] =\sum_{d=1}^n\sum_{i=1}^n\sum_{j=1}^m{f(d)\times [\gcd(i,j)=d]\times[f(d)\leq a]} =d=1∑ni=1∑nj=1∑mf(d)×[gcd(i,j)=d]×[f(d)≤a]
= ∑ d = 1 n f ( d ) × [ f ( d ) ≤ a ] ∑ i = 1 ⌊ n d ⌋ ∑ j = 1 ⌊ m d ⌋ [ gcd ( i , j ) = 1 ] =\sum_{d=1}^nf(d)\times[f(d)\leq a]\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor}{[\gcd(i,j)=1]} =d=1∑nf(d)×[f(d)≤a]i=1∑⌊dn⌋j=1∑⌊dm⌋[gcd(i,j)=1] -
由莫比乌斯反演得:
= ∑ d = 1 n f ( d ) × [ f ( d ) ≤ a ] ∑ i = 1 ⌊ n d ⌋ ∑ j = 1 ⌊ m d ⌋ ∑ k ∣ gcd ( i , j ) μ ( k ) =\sum_{d=1}^nf(d)\times[f(d)\leq a]\sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{j=1}^{\left\lfloor\frac{m}{d}\right\rfloor}\sum_{k|\gcd(i,j)}{\mu(k)} =d=1∑nf(d)×[f(d)≤a]i=1∑⌊dn⌋j=1∑⌊dm⌋k∣gcd(i,j)∑μ(k)
= ∑ d = 1 n f ( d ) × [ f ( d ) ≤ a ] ∑ k = 1 ⌊ n d ⌋ μ ( k ) ⌊ n k d ⌋ ⌊ m k d ⌋ =\sum_{d=1}^nf(d)\times[f(d)\leq a]\sum_{k=1}^{\left\lfloor\frac{n}{d}\right\rfloor}{\mu(k)}\left\lfloor\frac{n}{kd}\right\rfloor\left\lfloor\frac{m}{kd}\right\rfloor =d=1∑nf(d)×[f(d)≤a]k=1∑⌊dn⌋μ(k)⌊kdn⌋⌊kdm⌋ -
套路设 T = k d T=kd T=kd 得:
= ∑ T = 1 n ⌊ n T ⌋ ⌊ m T ⌋ ∑ x ∣ T μ ( T x ) × f ( x ) × [ f ( x ) ≤ a ] =\sum_{T=1}^n{\left\lfloor\frac{n}{T}\right\rfloor\left\lfloor\frac{m}{T}\right\rfloor}\sum_{x|T}{\mu(\frac{T}{x})\times f(x)\times [f(x)\leq a]} =T=1∑n⌊Tn⌋⌊Tm⌋x∣T∑μ(xT)×f(x)×[f(x)≤a] -
f ( x ) f(x) f(x) 可以在 O ( n ln n ) O(n\ln n) O(nlnn) 的时间预处理出来
-
因为有 a a a 的限制,当 a a a 变大时,后半部分式子也在改变。所以考虑将询问离线,按 a a a 的值从小到大排序,每次 a a a 变大时都类似插入操作,考虑用树状数组维护
-
设 g ( T ) = ∑ x ∣ T μ ( T x ) × f ( x ) g(T)=\sum\limits_{x|T}{\mu(\frac{T}{x})\times f(x)} g(T)=x∣T∑μ(xT)×f(x),对于所有 f ( x ) ≤ a f(x)\leq a f(x)≤a 来说,它会对 g ( x ) , g ( 2 x ) … g(x),g(2x)\dots g(x),g(2x)… 产生贡献,所以给每个 g g g 都加上 μ ( T x ) × f ( x ) \mu(\frac{T}{x})\times f(x) μ(xT)×f(x) 即可
代码
#include<bits/stdc++.h>
#define LL long long
using namespace std;
const int N=100010,Q=20010;
const LL MOD=2147483648;
struct node
{
int dat;
int id;
}f[N];
struct que
{
int n,m,s,id;
}a[Q];
int T;
int prime[N],v[N],tot,mu[N];
LL c[N],ans[N];
bool cmp(que xx,que yy)
{
return xx.s<yy.s;
}
bool cmp2(node xx,node yy)
{
return xx.dat<yy.dat;
}
void add(int x,LL y)
{
for(; x<=1e5; x+=(x&-x))
c[x]=(c[x]+y)%MOD;
}
LL ask(int x)
{
LL s=0;
for(; x; x-=(x&-x))
s=(c[x]+s)%MOD;
return s;
}
void update(int x,LL y)
{
for(int i=1; i*x<=1e5; i++)
add(i*x,1LL*mu[i]*y);
}
void primes()
{
mu[1]=1;
for(int i=2; i<=1e5; i++)
{
if(!v[i])
{
v[i]=i;
prime[++tot]=i;
mu[i]=-1;
}
for(int j=1; j<=tot; j++)
{
if(prime[j]>v[i] || prime[j]>1e5/i)
break;
v[i*prime[j]]=prime[j];
if(i%prime[j]==0)
break;
mu[i*prime[j]]=-mu[i];
}
}
for(int i=1; i<=1e5; i++)
for(int j=1; i*j<=1e5; j++)
f[i*j].dat=(f[i*j].dat+i)%MOD;
for(int i=1; i<=1e5; i++)
f[i].id=i;
sort(f+1,f+1+100000,cmp2);
}
LL calc(int a,int b)
{
if(a>b)
swap(a,b);
LL res=0;
for(int l=1,r; l<=a; l=r+1)
{
r=min(a/(a/l),b/(b/l));
res=(res+(ask(r)-ask(l-1)+MOD)%MOD*1LL*(a/l)%MOD*(b/l)%MOD)%MOD;
}
return res;
}
int main()
{
primes();
scanf("%d",&T);
for(int i=1; i<=T; i++)
scanf("%d%d%d",&a[i].n,&a[i].m,&a[i].s),a[i].id=i;
sort(a+1,a+1+T,cmp);
int p=1;
for(int i=1; i<=T; i++)
{
while(f[p].dat<=a[i].s && p<=1e5)
update(f[p].id,f[p].dat),p++;
ans[a[i].id]=calc(a[i].n,a[i].m);
}
for(int i=1; i<=T; i++)
printf("%lld\n",ans[i]);
return 0;
}