1.P3704-[SDOI2017]数字表格[莫比乌斯反演]
Problem
有一个\(n\times m\)的表格，坐标\((i,j)\)处的数字是\(f_{gcd(i,j)}\)，其中\(f\)是斐波那契数列，要求计算\(\prod_{i=1}^n\prod_{j=1}^mf_{gcd(i,j)}\)，答案对\(10^9+7\)取模
Solve
括号内的可以预处理，幂次记得用费马小定理计算
Code
#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int N=1e6+10;
const int mod=1e9+7;
int primes[N],st[N],cnt;
int mu[N];
ll f[N],invf[N],g[N],invg[N];
void get_primes()
{
mu[1]=1;
for(int i=2;i<=N-10;i++)
{
if(!st[i]) primes[++cnt]=i,mu[i]=-1;
for(int j=1;primes[j]*i<=N-10;j++)
{
st[primes[j]*i]=1;
if(i%primes[j]==0) break;
mu[i*primes[j]]=-mu[i];
}
}
}
ll power(ll x,int y)
{
ll res=1;
while(y)
{
if(y&1) res=res*x%mod;
x=x*x%mod;
y>>=1;
}
return res;
}
int main()
{
ios::sync_with_stdio(false);
cin.tie(nullptr);
get_primes();
f[1]=1;
invf[0]=invf[1]=1;
for(int i=2;i<=1000000;i++) f[i]=(f[i-1]+f[i-2])%mod,invf[i]=power(f[i],mod-2);
for(int i=0;i<=1000000;i++) g[i]=1;
for(int i=1;i<=1000000;i++)
{
if(!mu[i]) continue;
for(int j=i;j<=1000000;j+=i)
{
if(mu[i]<0) g[j]=g[j]*invf[j/i]%mod;
else g[j]=g[j]*f[j/i]%mod;
}
}
invg[1]=invg[0]=1;
for(int i=2;i<=1000000;i++)
g[i]=g[i]*g[i-1]%mod,invg[i]=power(g[i],mod-2);
int T;
cin>>T;
while(T--)
{
int n,m;
cin>>n>>m;
if(n>m) swap(n,m);
ll ans=1;
for(int l=1,r;l<=n;l=r+1)
{
r=min(n/(n/l),m/(m/l));
ll t=g[r]*invg[l-1]%mod;
ans=ans*power(t,1LL*(n/l)*(m/l)%(mod-1))%mod;//欧拉定理，求幂
}
cout<<(ans+mod)%mod<<'\n';
}
return 0;
}
2.[SDOI2014]数表[莫比乌斯反演,树状数组]
Problem
\(Q\)个询问，每次给定\(n\)，\(m\)，\(a\)，计算\(\sum_{i=1}^n\sum_{j=1}^m\sigma(gcd(i,j))[\sigma(gcd(i,j))\le a]\)(原问题是求\(ij\)的自然数因子的和)
\(1\le n,m\le 10^5\) \(|a|\le 10^9\)
Solve
先不考虑\(a\)的限制
在不考虑\(a\)的情况下，后面部分直接\(O(nlogn)\)预处理即可。现在考虑小于等于\(a\)的情况，考虑把所有询问离线按照\(a\)的值从小到大排序，对于\(\sum_{d|T}\sigma(d)\mu(\frac{T}{d})\)，发现对于一个固定的\(a\)，计算的是\(\sigma(d)\le a\)的贡献，那么考虑\(\sigma(d)\)对于\(T\)的贡献是多少，显然就是\(d\)的每个倍数\(T\)会对\(a\)产生贡献，然后分块查询的时候，是一个求前缀的过程，在确定\(d\)的情况下，统计\(T\)对\(a\)的贡献是单点修改，于是考虑用树状数组来维护。
Code
#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int N=1e5+10;
const ll mod=1LL<<31;
int primes[N],st[N],cnt;
int mu[N];
ll d[N],g[N];
void get_primes()
{
mu[1]=1;
for(int i=2;i<=N-10;i++)
{
if(!st[i]) primes[++cnt]=i,mu[i]=-1;
for(int j=1;primes[j]*i<=N-10;j++)
{
st[primes[j]*i]=1;
if(i%primes[j]==0) break;
mu[i*primes[j]]=-mu[i];
}
}
}
struct Query{
int n,m,a,id;
bool operator < (const Query &t)const{
return a<t.a;
}
};
ll tr[N];
void add(int p,ll x)
{
for(;p<=100000;p+=p&-p) tr[p]=(tr[p]+x)%mod;
}
ll ask(int p)
{
ll res=0;
for(;p;p-=p&-p) res=(res+tr[p])%mod;
return res;
}
int main()
{
ios::sync_with_stdio(false);
cin.tie(nullptr);
get_primes();
for(int i=1;i<=100000;i++)
for(int j=i;j<=100000;j+=i)
d[j]+=i;
vector<pair<int,int>>p(100000);
for(int i=0;i<100000;i++) p[i].first=d[i+1],p[i].second=i+1;
sort(p.begin(), p.end());
int t;
cin>>t;
vector<Query>q;
for(int i=0;i<t;i++)
{
int n,m,a;
cin>>n>>m>>a;
if(n>m) swap(n,m);
q.push_back({n,m,a,i});
}
vector<ll>ans(t);
sort(q.begin(), q.end());
auto solve=[&](int n,int m,int id)->void{
ll res=0;
for(int l=1,r;l<=n;l=r+1)
{
r=min(n/(n/l),m/(m/l));
res=(res+(n/l)*(m/l)%mod*(ask(r)-ask(l-1)+mod)%mod)%mod;
}
ans[id]=res;
};
for(int i=0,j=0;i<t;i++)
{
auto [n,m,a,id]=q[i];
while(p[j].first<=a&&p[j].second<=100000)
{
int x=p[j].second;
ll delta=p[j].first;
for(int k=1;k*x<=100000;k++) add(k*x,1LL*mu[k]*d[x]);
j++;
}
solve(n,m,id);
}
for(int i=0;i<t;i++) cout<<ans[i]<<"\n";
return 0;
}
3.P6156 简单题
Problem
给定\(n\)，\(k\)，计算\(\sum_{i=1}^n\sum_{j=1}^n(i+j)^kf(gcd(i,j))gcd(i,j)\)
其中\(f(k)\)的定义如下：如果\(k\)有平方因子，则\(f(k)=0\)，否则\(f(k)=1\)
答案对\(998244353\)取模
\(1\le n\le 5\times 10^6\) \(k\le 10^{18}\)
Solve
可以发现\(f(k)=\mu(k)^2\)
考虑计算\(S(n)=\sum_{i=1}^n\sum_{j=1}^n=(i+j)^k\)
记\(P(n)=\sum_{i=1}^ni^k\)，则
考虑计算\(f(T)=\sum_{d|T}d\mu(d)^2\mu(\frac{T}{d})=\)，这是一个卷积的形式，由于只包含积性函数，所以\(f\)也是积性函数，因此可以用线性筛计算积性函数
素数\(p\)的\(f(p)=\mu(1)^2\mu(p)+p\mu(p)^2\mu(1)=p-1\)
考虑再线性筛的过程计算\(ip\)的\(f(ip)\)
若\(i\%p\neq 0\)，说明\(ip\)不会包含\(p\)的二重因子，所以\(gcd(i,p)=1\)，所以\(f(ip)=f(i)f(p)=f(i)(p-1)\)
若\(i\%p=0\)，说明\(i\)和\(p\)不互素：
- 若\(i/p\%p\neq 0\)，说明\(ip\)最多包含\(p\)的二重因子，\(gcd(i/p,p^2)=1\)，\(f(ip)=f(i/p)f(p^2)=-pf(i/p)\)
- 若\(i/p\%p = 0\)，假设\(ip\)包含\(p\)的\(k\)重因子(\(k\ge 3\))，那么\(gcd(i/p^k,p^k)=1\)，\(f(ip)=f(i/p^k)f(p^k)=0\)
Code
#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int N=5e6+10,NN=5e6;
const int mod=998244353;
int primes[2*N],st[2*N],cnt;
ll s[2*N],f[2*N],pw[2*N],k;
ll power(ll x,ll y)
{
ll res=1;
while(y)
{
if(y&1) res=res*x%mod;
x=x*x%mod;
y>>=1;
}
return res;
}
void get_primes()
{
int MM=NN*2;
f[1]=1,pw[1]=1,s[1]=power(2,k);
for(int i=2;i<=MM;i++)
{
if(!st[i]) primes[++cnt]=i,f[i]=i-1,pw[i]=power(i,k);
for(int j=1;primes[j]*i<=N-10;j++)
{
st[primes[j]*i]=1;
pw[primes[j]*i]=pw[i]*pw[primes[j]]%mod;
if(i%primes[j]==0)
{
//f(T)=mu(d) x mu(T/d)[d|T]
if(i/primes[j]%primes[j]!=0)
f[i*primes[j]]=f[i/primes[j]]*(mod-primes[j])%mod; //a_i=2的情况,-p_i
else f[i*primes[j]]=0;
break;
}
f[i*primes[j]]=f[i]*(primes[j]-1)%mod; //a_i=1的情况,p_i-1
}
}
for(int i=1;i<=NN;i++)f[i]=f[i]*pw[i]%mod;
for(int i=1;i<=NN;i++) f[i]=(f[i]+f[i-1])%mod;
for(int i=1;i<=MM;i++) pw[i]=(pw[i]+pw[i-1])%mod;
for(int i=1;i<NN;i++) s[i+1]=(s[i]+pw[(i<<1)+2]+pw[(i<<1)+1]-2LL*pw[i+1]%mod+mod)%mod;
}
int main()
{
ios::sync_with_stdio(false);
cin.tie(nullptr);
int n;
cin>>n>>k;
k=k%(mod-1);
get_primes();
ll ans=0;
for(int l=1,r;l<=n;l=r+1)
{
r=n/(n/l);
ans=(ans+s[n/l]*(f[r]-f[l-1]+mod)%mod)%mod;
}
cout<<ans<<'\n';
return 0;
}
4.P5221 Product
Problem
计算\(\prod_{i=1}^n\prod_{j=1}^n\frac{lcm(i,j)}{gcd(i,j)}(mod\ 104857601)\)
\(1\le n\le 1000000\)
Solve
分看开，前面\(\prod_{i=1}^n\prod_{j=1}^nij=\prod_{i=1}^ni^n\prod_{j=1}^nj=\prod_{i=1}^ni^nn!=(n!)^n(\prod_{i=1}^i)^n=(n!)^n(n!)^n=(n!)^{2n}\)
再看后面\(\sum_{i=1}^n\sum_{j=1}^ngcd(i,j)^{-2}\)，先不考虑幂次，考虑如何计算\(\sum_{i=1}^n\sum_{j=1}^ngcd(i,j)\)
考虑用欧拉定理减小幂次\(a^k\equiv a^{k\%\phi(p)}(mod\ p)\)，\(p\in primes\)
于是\((\prod_{d=1}^nd^{\sum_{i=1}^{\lfloor \frac{n}{d} \rfloor}\sum_{j=1}^{\lfloor \frac{n}{d} \rfloor}[gcd(i,j)==1]})^{-2}=(\prod_{d=1}^nd^{(\sum_{i=1}^{\lfloor \frac{n}{d} \rfloor}\sum_{j=1}^{\lfloor \frac{n}{d} \rfloor}[gcd(i,j)==1])\%(mod-1)})^{-2}\)
\(\sum_{i=1}^{\lfloor \frac{n}{d} \rfloor}\sum_{j=1}^{\lfloor \frac{n}{d} \rfloor}[gcd(i,j)==1]\)这个用欧拉函数求前缀和即可，即\(s({\lfloor \frac{n}{d} \rfloor})=\sum_{i=1}^{\lfloor \frac{n}{d} \rfloor}\sum_{j=1}^{\lfloor \frac{n}{d} \rfloor}[gcd(i,j)==1]=2\times \sum_{i=1}^{\lfloor \frac{n}{d} \rfloor}\sum_{j=i}^{\lfloor \frac{n}{d} \rfloor}\phi(i)=2\phi({\lfloor \frac{n}{d} \rfloor})+s({\lfloor \frac{n}{d} \rfloor}-1)\)
负数幂次用除法逆元，并且本题卡空间，不能开long long的数组，st数组用bool类型，primes数组事先计算\(N\)范围的素数个数来开
Code
#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int N=1000005,NN=1000000;
const int mod=104857601;
int primes[80000],phi[N],cnt;
bool st[N];
int power(int x,int y)
{
int res=1;
while(y)
{
if(y&1) res=1LL*res*x%mod;
x=1LL*x*x%mod;
y>>=1;
}
return res;
}
void get_primes()
{
phi[1]=1;
for(int i=2;i<=NN;i++)
{
if(!st[i]) primes[++cnt]=i,phi[i]=i-1;
for(int j=1;primes[j]*i<=N-10;j++)
{
st[primes[j]*i]=1;
if(i%primes[j]==0){
phi[i*primes[j]]=phi[i]*primes[j];
break;
}
phi[i*primes[j]]=phi[i]*phi[primes[j]];
}
}
for(int i=1;i<=NN;i++) phi[i]=2*phi[i]+phi[i-1]%(mod-1);
}
int main()
{
ios::sync_with_stdio(false);
cin.tie(nullptr);
get_primes();
int n;
cin>>n;
int ans1=1,ans2=1,ans=0;
for(int i=1;i<=n;i++) ans1=1LL*ans1*i%mod;
ans1=power(ans1,2*n);
for(int i=2;i<=n;i++) ans2=(1LL*ans2*power(i,phi[n/i]-1)%mod);
ans=(1LL*ans1*power(1LL*ans2*ans2%mod,mod-2)%mod);
cout<<ans<<'\n';
return 0;
}
5.2020 ICPC JiNan F. Gcd Product
Problem
给定两个长度\(N\)的序列\(A\)和\(B\)，对于每个\(k\)(\(1\le k\le N\))，\(C_k=\sum_{i=1}^kA_{gcd(i,k)}B_{gcd(k+1-i,k)}\)，由于答案可能很大，对\(998244353\)取模，最后输出答案\(C_1\oplus C_2\oplus...\oplus C_N\)
\(1\le N\le 5\times 10^5\)
\(0\le A_i,B_i\le 998244353\)
Solve
考虑合并\(d_1,t_1\)和\(d_2,t_2\)，也就是把后面的限制条件\([d_1|i],[t_1|\frac{i}{d_1}]\)合并
令\(f(T_1)=\sum_{d_1|T_1}A_{d_1}\mu(\frac{T_1}{d_1}),g(T_2)=\sum_{d_2|T_2}B_{d_2}\mu(\frac{T_2}{d_2})\)
注意到\(d\frac{T}{d}=T\)，很明显的一个卷积形式，所以这部分可以用狄利克雷卷积，时间复杂度是\(O(nlogn)\)
考虑计算后面\(\sum_{i=1}^k[T_1|i][T_2|k+1-i]\)
下面考虑如何快速求\((k_1,k_2)\)正整数解的个数
证明：
假设\(k_2\)在模\(T_1\)下不止有\(1\)个解，假设有\(2\)个解分别设为\(x_1,x_2\)，那么有\(r_1T_1=x_1T_2-1,r_2T_1=x_2T_2-1\)，两式作差得\((r_1-r_2)T_1=(x_1-x_2)T_2\)，由于\(gcd(T_1,T_2)=1\)，所以有\((r_1-r_2)=kT_2,(x_1-x_2)=kT_1\)，不妨设\(k=1\)，那么\((x_1-x_2)=T_1\)，在模\(T_1\)意义下，这样迭代下去，一定可以找到一个小于\(T_1\)的唯一满足条件的\(x_1\)。
设\(k_1,k_2\)的解为\(y_1,y_2\)，那么\(y_1T_1<T_2T_1,y_2T_2<T_1T_2\)，设\(k_1=r_1T_2+y_1,k_2=r_2T_1+y_2\),带入\(k_1T_1+k_2T_2=k+1\)得，\(y_1T_1+y_2T_2=k+1-(r_1+r_2)T_1T_2\)，两边同时对\(T_1T_2\)取模，得到\(y_1T_1+y_2T_2\equiv1(mod\ T_1T_2)\)，进一步得到\(y_1T_1+y_2T_2=T_1T_2+1\)(因为\(1\le y_1T_1+y_2T_2\lt T_2T_1+T_1T_2=2T_2T_2\))，之后可以认为将剩余的\(k-T_1T_2\)分配到\(y_1T_1,y_2T_2\)上，所以解的个数就是\(\frac{k}{T_1T_2}\)。
另一种理解形式：将\(y_1T_1+y_2T_2=T_1T_2+1\)带入\(y_1T_1+y_2T_2=k+1-(r_1+r_2)T_1T_2\)，得到\(r_1+r_2=\frac{k}{T_1T_2}-1\)，由于\(x+y=n\)共有\(n+1\)对非负整数解，所以\(r_1+r_2\)共有非负\(\frac{k}{T_1T_2}\)对整数解，由于\(y_1,y_2\)唯一，而\(k_1=r_1T_1+y_1,k_2=r_2T_2+y_2\)，所以\((k_1,k_2)\)的解数和\((r_1,r_2)\)相同。
对于\(T_1T_2|k\)的证明：\(r_1T_1=k,r_2T_2=k\)，因此\(r_1T_1=r_2T_2\)，所以\(T_1|r_2T_2\)，由于\(gcd(T_1,T_2)=1\)，所以\(T_1|r_2\)，所以\(r_2=r_3T_1\)，所以\(r_3T_1T_2=k\)
总⽽⾔之, 答案就是\(id*(a*\mu)\times (b*\mu)\) , 其中\(\times\)是互质狄利克雷卷积,\(*\) 是普通狄利克雷卷积， 时间复杂度: \(O(nlogn)\)
#include <bits/stdc++.h>
using namespace std;
const int P=998244353;
const int N=5e5+10;
int f[N],g[N],w[N],A[N],B[N];
int primes[N],st[N],phi[N],mu[N],cnt;
void init(int n)
{
mu[1]=phi[1]=1;
for(int i=2;i<=n;i++)
{
if(!st[i])
{
primes[++cnt]=i;
mu[i]=-1;
phi[i]=i-1;
}
for(int j=1;primes[j]*i<=n;j++)
{
st[i*primes[j]]=1;
if(i%primes[j]==0)
{
phi[i*primes[j]]=phi[i]*primes[j];
break;
}
mu[i*primes[j]]=-mu[i];
phi[i*primes[j]]=phi[i]*(primes[j]-1);
}
}
}
int C[N];
int main()
{
ios::sync_with_stdio(false);
cin.tie(nullptr);
int n;
cin>>n;
init(n);
for(int i=1;i<=n;i++) cin>>A[i];
for(int i=1;i<=n;i++) cin>>B[i];
for(int i=1;i<=n;i++) //O(nlogn)
for(int j=i;j<=n;j+=i)
if(mu[j/i]!=0)
f[j]=((f[j]+mu[j/i]*A[i])%P+P)%P,
g[j]=((g[j]+mu[j/i]*B[i])%P+P)%P;
for(int i=1;i<=n;i++) //O(nlogn)
for(int j=i;j<=n;j+=i)
if(phi[i]*phi[j/i]==phi[j]) w[j]=(w[j]+1LL*f[i]*g[j/i]%P)%P;
for(int i=1;i<=n;i++)//O(nlogh)
for(int j=i;j<=n;j+=i) C[j]=(C[j]+1LL*w[i]*(j/i)%P)%P;
int ans=0;
for(int i=1;i<=n;i++) ans^=C[i];
cout<<ans<<'\n';
return 0;
}