【做题记录】数论（马思博）
LG P3312 [SDOI2014] 数表
首先不考虑 \(a\) 的限制，有：
\[\begin{aligned}
&\sum_{i=1}^{n}\sum_{j=1}^{m}\operatorname{s}(\gcd(i,j))\\
=&\sum_{i=1}^{n}\sum_{j=1}^{m}\sum_{d=1}^{n}\operatorname{s}(d)[\gcd(i,j)=d]\\
=&\sum_{d=1}^{n}\operatorname{s}(d)\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}[\gcd(i,j=1]\\
=&\sum_{d=1}^{n}\operatorname{s}(d)\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}\sum_{k\mid\gcd(i,j)}\mu(k)\\
=&\sum_{d=1}^{n}\operatorname{s}(d)\sum_{k=1}^{n}\mu(k)\lfloor\frac{n}{dk}\rfloor\lfloor\frac{m}{dk}\rfloor
\end{aligned}
\]
改变一下枚举顺序，先枚举 \(dk\)，于是有：
\[\sum_{k=1}^{n}\lfloor\frac{n}{k}\rfloor\lfloor\frac{m}{k}\rfloor\sum_{d\mid k}\operatorname{s}(d)\mu(\frac{k}{d})
\]
考虑 \(a\) 的限制，显然 \(d\) 能产生贡献当且仅当 \(\operatorname{s}(d)\le a\)。考虑离线，按照 \(a\) 排序，每次对于能贡献的 \(d\)，将 \(k\,\operatorname{s.t.}d\mid k\) 的贡献增加 \(\operatorname{s}(d)\mu(\frac{k}{d})\)。于是对于询问 \((n,m)\) 可以对 \(k\) 数论分块。需要单点加区间查，树状数组即可。时间复杂度 \(O(n\ln n\log n+q\sqrt{n}\log n)\)。
Code
#include<bits/stdc++.h>
#define int long long
#define il inline
#define pii pair<int,int>
#define fir first
#define sec second
#define mp make_pair
using namespace std;
namespace asbt{
const int n=1e5,maxn=1e5+5,mod=1<<31;
il int pls(int x,int y){
return x+y<mod?x+y:x+y-mod;
}
il void add(int &x,int y){
x=pls(x,y);
}
il int sub(int x,int y){
return x<y?x-y+mod:x-y;
}
il void mns(int &x,int y){
x=sub(x,y);
}
int q,prn,prm[maxn],mu[maxn],g[maxn],ans[maxn];
bool npr[maxn];
pii f[maxn];
struct node{
int n,m,a,id;
il bool operator<(const node &x)const{
return a<x.a;
}
}a[maxn];
il void euler(){
f[1]=mp(1,1),mu[1]=1;
for(int i=2;i<=n;i++){
if(!npr[i]){
prm[++prn]=i;
f[i].fir=g[i]=i+1;
f[i].sec=i;
mu[i]=-1;
}
for(int j=1;j<=prn&&i*prm[j]<=n;j++){
npr[i*prm[j]]=1;
if(i%prm[j]==0){
// cout<<i*prm[j]<<'\n';
mu[i*prm[j]]=0;
g[i*prm[j]]=g[i]*prm[j]+1;
f[i*prm[j]]=mp(f[i].fir/g[i]*g[i*prm[j]],i*prm[j]);
// if(i*prm[j]==4){
// cout<<4<<' '<<g[i]<<' '<<prm[j]<<' '<<g[i*prm[j]]<<'\n';
// }
break;
}
mu[i*prm[j]]=-mu[i];
g[i*prm[j]]=prm[j]+1;
f[i*prm[j]]=mp(f[i].fir*(prm[j]+1),i*prm[j]);
}
}
// for(int i=1;i<=10;i++){
// cout<<g[i]<<' ';
// }
// cout<<'\n';
// for(int i=1;i<=10;i++){
// cout<<f[i].fir<<' ';
// }
// cout<<'\n';
sort(f+1,f+n+1);
}
struct{
#define lowbit(x) (x&-x)
int tr[maxn];
il void ad(int p,int v){
for(;p<=n;p+=lowbit(p)){
add(tr[p],v);
}
}
il int query(int p){
int res=0;
for(;p;p-=lowbit(p)){
add(res,tr[p]);
}
return res;
}
#undef lowbit
}F;
il int solve(int n,int m){
int res=0;
for(int i=1,j;i<=n;i=j+1){
j=min(n/(n/i),m/(m/i));
res=((F.query(j)-F.query(i-1)+mod)*(n/i)%mod*(m/i)+res)%mod;
}
return res;
}
int main(){
ios::sync_with_stdio(0),cin.tie(0);
cin>>q;
for(int i=1;i<=q;i++){
cin>>a[i].n>>a[i].m>>a[i].a;
if(a[i].n>a[i].m){
swap(a[i].n,a[i].m);
}
a[i].id=i;
}
euler();
sort(a+1,a+q+1);
for(int i=1,j=1;i<=q;i++){
for(;j<=n&&f[j].fir<=a[i].a;j++){
for(int k=f[j].sec;k<=n;k+=f[j].sec){
F.ad(k,((f[j].fir*mu[k/f[j].sec])%mod+mod)%mod);
}
}
ans[a[i].id]=solve(a[i].n,a[i].m);
}
for(int i=1;i<=q;i++){
cout<<ans[i]<<'\n';
}
return 0;
}
}
signed main(){return asbt::main();}
D. [国家集训队] Crash的数字表格 / JZPTAB
\[\begin{aligned}
&\sum_{i=1}^{n}\sum_{j=1}^{m}\operatorname{lcm}(i,j)\\
=&\sum_{i=1}^{n}\sum_{j=1}^{m}\frac{i\times j}{\gcd(i,j)}\\
=&\sum_{d=1}^{\min(n,m)}\frac{1}{d}\sum_{i=1}^{n}\sum_{j=1}^{m}i\times j\times[\gcd(i,j)=d]\\
\end{aligned}
\]
不妨令 \(n\le m\)，上式即
\[\sum_{d=1}^{n}\frac{1}{d}\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}(i\times d)\times(j\times d)\times[\gcd(i,j)=1]
\]
莫比乌斯反演，得
\[\begin{aligned}
&\sum_{d=1}^{n}d\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}i\times j\sum_{k\mid\gcd(i,j)}\mu(k)\\
=&\sum_{d=1}^{n}d\sum_{k=1}^{\lfloor\frac{n}{d}\rfloor}\mu(k)\sum_{i=1}^{n}i\times[k\mid i]\sum_{j=1}^{m}j\times[k\mid j]\\
=&\sum_{d=1}^{n}d\sum_{k=1}^{\lfloor\frac{n}{d}\rfloor}\mu(k)\times(k\times\frac{\lfloor\frac{n}{d\times k}\rfloor\times(\lfloor\frac{n}{d\times k}\rfloor+1)}{2})\times(k\times\frac{\lfloor\frac{m}{d\times k}\rfloor\times(\lfloor\frac{m}{d\times k}\rfloor+1)}{2})
\end{aligned}
\]
于是可以调和级数求。预处理 \(i^2\) 和 \(\frac{i\times(i+1)}{2}\) 减少取模即可。
Code
#include<bits/stdc++.h>
#define ll long long
#define il inline
using namespace std;
namespace asbt{
const int maxn=1e7+5,mod=20101009;
int n,m,prn,prm[maxn],mu[maxn],f[maxn],g[maxn];
bool npr[maxn];
il void euler(int n=1e7){
mu[1]=npr[1]=1;
for(int i=2;i<=n;i++){
if(!npr[i]){
prm[++prn]=i,mu[i]=-1;
}
for(int j=1;j<=prn&&i*1ll*prm[j]<=n;j++){
npr[i*prm[j]]=1;
if(i%prm[j]==0){
mu[i*prm[j]]=0;
break;
}
mu[i*prm[j]]=-mu[i];
}
}
}
il void init(int n=1e7){
for(int i=1;i<=n;i++){
f[i]=i*1ll*i%mod;
g[i]=i*1ll*(i+1)/2%mod;
}
}
int main(){
ios::sync_with_stdio(0),cin.tie(0);
euler(),init();
cin>>n>>m;
if(n>m){
swap(n,m);
}
int ans=0;
for(int i=1;i<=n;i++){
int res=0;
for(int j=1;j<=n/i;j++){
int tmp=g[n/i/j]*1ll*g[m/i/j]%mod*f[j]%mod;
if(mu[j]==1){
(res+=tmp)%=mod;
}
else if(mu[j]==-1){
(res+=mod-tmp)%=mod;
}
}
ans=(res*1ll*i+ans)%mod;
}
cout<<ans;
return 0;
}
}
int main(){return asbt::main();}
然而这还不够优秀。观察这个式子：
\[\begin{aligned}
&\sum_{d=1}^{n}d\sum_{k=1}^{\lfloor\frac{n}{d}\rfloor}\mu(k)\times(k\times\frac{\lfloor\frac{n}{d\times k}\rfloor\times(\lfloor\frac{n}{d\times k}\rfloor+1)}{2})\times(k\times\frac{\lfloor\frac{m}{d\times k}\rfloor\times(\lfloor\frac{m}{d\times k}\rfloor+1)}{2})\\
=&\sum_{d=1}^{n}d\sum_{k=1}^{\lfloor\frac{n}{d}\rfloor}\mu(k)\times k^2\times\frac{\lfloor\frac{n}{d\times k}\rfloor\times(\lfloor\frac{n}{d\times k}\rfloor+1)}{2}\times\frac{\lfloor\frac{m}{d\times k}\rfloor\times(\lfloor\frac{m}{d\times k}\rfloor+1)}{2}
\end{aligned}
\]
考虑设 \(\operatorname{sum}(n,m)=\sum_{k=1}^{n}\mu(k)\times k^2\times\frac{\lfloor\frac{n}{k}\rfloor\times(\lfloor\frac{n}{k}\rfloor+1)}{2}\times\frac{\lfloor\frac{m}{k}\rfloor\times(\lfloor\frac{m}{k}\rfloor+1)}{2}\)，给 \(\mu(k)\times k^2\) 做前缀和，可以数论分块做。于是原式变成：
\[\sum_{d=1}^{n}d\times\operatorname{sum}(\lfloor\frac{n}{d}\rfloor,\lfloor\frac{m}{d}\rfloor)
\]
这个东西还可以数论分块，于是时间复杂度线性。
Code
#include<bits/stdc++.h>
#define ll long long
#define il inline
using namespace std;
namespace asbt{
const int maxn=1e7+5,mod=20101009;
int n,m,prn,prm[maxn],mu[maxn],f[maxn],g[maxn],h[maxn];
bool npr[maxn];
il void euler(int n=1e7){
mu[1]=npr[1]=1;
for(int i=2;i<=n;i++){
if(!npr[i]){
prm[++prn]=i,mu[i]=-1;
}
for(int j=1;j<=prn&&i*1ll*prm[j]<=n;j++){
npr[i*prm[j]]=1;
if(i%prm[j]==0){
mu[i*prm[j]]=0;
break;
}
mu[i*prm[j]]=-mu[i];
}
}
}
il void init(int n=1e7){
for(int i=1;i<=n;i++){
f[i]=(f[i-1]+i)%mod;
g[i]=((g[i-1]+mu[i]*1ll*i*i)%mod+mod)%mod;
h[i]=i*1ll*(i+1)/2%mod;
}
}
il int calc(int n,int m){
int res=0;
for(int i=1;i<=n;i++){
int j=min(n/(n/i),m/(m/i));
res=((g[j]-g[i-1]+mod)*1ll*h[n/i]%mod*h[m/i]+res)%mod;
i=j;
}
return res;
}
int main(){
ios::sync_with_stdio(0),cin.tie(0);
euler(),init();
cin>>n>>m;
if(n>m){
swap(n,m);
}
int ans=0;
for(int i=1;i<=n;i++){
int j=min(n/(n/i),m/(m/i));
ans=((f[j]-f[i-1]+mod)*1ll*calc(n/i,m/i)+ans)%mod;
i=j;
}
cout<<ans;
return 0;
}
}
int main(){return asbt::main();}
LG P3327 [SDOI2015] 约数个数和
考虑将 \(i\) 和 \(j\) 分解质因数，分别记为 \(\prod p_k^{\alpha_k}\) 和 \(\prod p_k^{\beta_k}\)。于是 \(\operatorname{d}(ij)=\prod(\alpha_k+\beta_k+1)\)。考虑怎么转化这个式子，考虑 \(x\mid i\) 和 \(y\mid j\) 且 \(x\perp y\)，对于每个 \(p_k\)，显然 \(\alpha_{x,k}=0\lor\alpha_{y,k}=0\)，于是这样的 \((x,y)\) 共有 \(\prod(\alpha_k+\beta_k+1)\) 种取值。于是 \(\operatorname{d}(ij)=\sum_{x\mid i}\sum_{y\mid j}[\gcd(x,y)=1]\)。于是：
\[\begin{aligned}
\sum_{i=1}^{n}\sum_{j=1}^{m}\operatorname{d}(ij)=&\sum_{i=1}^{n}\sum_{j=1}^{m}\sum_{x\mid i}\sum_{y\mid j}[\gcd(x,y)=1]\\
=&\sum_{i=1}^{n}\sum_{j=1}^{m}\sum_{x\mid i}\sum_{y\mid j}\sum_{d\mid\gcd(x,y)}\mu(d)\\
=&\sum_{i=1}^{n}\sum_{j=1}^{m}\sum_{d\mid\gcd(i,j)}\mu(d)\sum_{x\mid i}[d\mid x]\sum_{y\mid j}[d\mid y]\\
=&\sum_{i=1}^{n}\sum_{j=1}^{m}\sum_{d\mid\gcd(i,j)}\mu(d)\operatorname{d}(\frac{i}{d})\operatorname{d}(\frac{j}{d})\\
=&\sum_{d=1}^{n}\mu(d)\sum_{i=1}^{n}[d\mid i]\operatorname{d}(\frac{i}{d})\sum_{j=1}^{m}[d\mid j]\operatorname{d}(\frac{j}{d})\\
=&\sum_{d=1}^{n}\mu(d)\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\operatorname{d}(i)\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}\operatorname{d}(j)
\end{aligned}
\]
给 \(\operatorname{d}\) 和 \(\mu\) 做前缀和，数论分块即可。
Code
#include<bits/stdc++.h>
#define int long long
#define il inline
using namespace std;
namespace asbt{
const int maxn=5e4+5;
int T,n,m,prn,prm[maxn],g[maxn],f[maxn];
bool npr[maxn];
il void init(int n=5e4){
for(int i=1;i<=n;i++){
for(int j=1;j<=i/j;j++){
if(i%j==0){
f[i]++;
if(j!=i/j){
f[i]++;
}
}
}
}
g[1]=npr[1]=1;
for(int i=2;i<=n;i++){
if(!npr[i]){
prm[++prn]=i;
g[i]=-1;
}
for(int j=1;j<=prn&&i*1ll*prm[j]<=n;j++){
npr[i*prm[j]]=1;
if(i%prm[j]==0){
g[i*prm[j]]=0;
break;
}
g[i*prm[j]]=-g[i];
}
}
for(int i=1;i<=n;i++){
f[i]+=f[i-1];
g[i]+=g[i-1];
}
}
int main(){
ios::sync_with_stdio(0),cin.tie(0);
init();
cin>>T;
while(T--){
cin>>n>>m;
if(n>m){
swap(n,m);
}
int ans=0;
for(int i=1;i<=n;i++){
int j=min(n/(n/i),m/(m/i));
ans+=(g[j]-g[i-1])*f[n/i]*f[m/i];
i=j;
}
cout<<ans<<"\n";
}
return 0;
}
}
signed main(){return asbt::main();}
浙公网安备 33010602011771号