【做題記錄】數論（馬思博）
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}w0obha2h00\rfloor}\sum_{j=1}^{\lfloor\frac{m}w0obha2h00\rfloor}[\gcd(i,j=1]\\
=&\sum_{d=1}^{n}\operatorname{s}(d)\sum_{i=1}^{\lfloor\frac{n}w0obha2h00\rfloor}\sum_{j=1}^{\lfloor\frac{m}w0obha2h00\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}w0obha2h00)
\]
考慮 \(a\) 的限制，顯然 \(d\) 能產生貢獻當且僅當 \(\operatorname{s}(d)\le a\)?？紤]離線，按照 \(a\) 排序，每次對于能貢獻的 \(d\)，將 \(k\,\operatorname{s.t.}d\mid k\) 的貢獻增加 \(\operatorname{s}(d)\mu(\frac{k}w0obha2h00)\)。于是對于詢問 \((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}w0obha2h00\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}w0obha2h00\sum_{i=1}^{\lfloor\frac{n}w0obha2h00\rfloor}\sum_{j=1}^{\lfloor\frac{m}w0obha2h00\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}w0obha2h00\rfloor}\sum_{j=1}^{\lfloor\frac{m}w0obha2h00\rfloor}i\times j\sum_{k\mid\gcd(i,j)}\mu(k)\\
=&\sum_{d=1}^{n}d\sum_{k=1}^{\lfloor\frac{n}w0obha2h00\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}w0obha2h00\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}w0obha2h00\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}w0obha2h00\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}w0obha2h00\rfloor,\lfloor\frac{m}w0obha2h00\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}\)。于是 \(\operatornamew0obha2h00(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)\) 種取值。于是 \(\operatornamew0obha2h00(ij)=\sum_{x\mid i}\sum_{y\mid j}[\gcd(x,y)=1]\)。于是：
\[\begin{aligned}
\sum_{i=1}^{n}\sum_{j=1}^{m}\operatornamew0obha2h00(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)\operatornamew0obha2h00(\frac{i}w0obha2h00)\operatornamew0obha2h00(\frac{j}w0obha2h00)\\
=&\sum_{d=1}^{n}\mu(d)\sum_{i=1}^{n}[d\mid i]\operatornamew0obha2h00(\frac{i}w0obha2h00)\sum_{j=1}^{m}[d\mid j]\operatornamew0obha2h00(\frac{j}w0obha2h00)\\
=&\sum_{d=1}^{n}\mu(d)\sum_{i=1}^{\lfloor\frac{n}w0obha2h00\rfloor}\operatornamew0obha2h00(i)\sum_{j=1}^{\lfloor\frac{m}w0obha2h00\rfloor}\operatornamew0obha2h00(j)
\end{aligned}
\]
給 \(\operatornamew0obha2h00\) 和 \(\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號