尤拉函式&莫比烏斯反演
阿新 • • 發佈：2019-01-01
近幾天做了幾道有關反演的問題，在此集合一下吧。
1、[BZOJ 2301]HAOI2011 Problem b
2、[BZOJ 2440]中山市選2011 完全平方數
3、gcd
4、[BZOJ 2186]SDOI2008 莎拉公主的困惑
5、[BZOJ 3529]SDOI2014 數表
（蒟蒻自認為反演一類的的題目重要的就是記住兩個重要的公式：1、sigma(mu[i] , i|n ) = [n==1] 2、sigma(phi[i] , i|n ) = n）
題解：
1、[BZOJ 2301]HAOI2011 Problem b
這基本上算是入門題了吧。。。。
gcd(a,b)==k等價於gcd(a/k,b/k)==1，這樣我們將範圍除k，等價到求區間內互質的二元組數。
假設query(n,m)得到gcd(a,b)==1的數量，那麼答案很明顯就是gcd(b,d) - gcd(a-1,d) - gcd(b,c-1) + gcd(a-1,c-1)
簡單了，下面簡述一下推導過程吧
為了不被多組資料卡掉，我們觀察到n/k的答案只有不超過根號n個不同答案，對mu[]字首和優化，那麼每次詢問時O(sqrt(n))的複雜度
#include <cmath> #include <cstdio> #include <cstring> #include <algorithm> using namespace std; int l,r,ans,mid,INF,K,Case,tot; int prime[41005],mu[41005],i,j; bool check[41005]; void init(){ mu[1]=1; for (i=2;i<41000;i++){ if (!check[i]) mu[i]=-1, prime[++tot]=i; for (j=1;j<=tot;j++){ if (i*prime[j]>41000) break; check[i*prime[j]]=1; if (i%prime[j]==0) {mu[i*prime[j]]=0;break;} else mu[i*prime[j]] = -mu[i]; } } } bool Judge(){ int ret=0; for (int i=1;i*i<=mid;i++) ret+=(mid/i/i)*mu[i]; return (ret>=K); } int main(){ //freopen("2440.in","r",stdin); //freopen("2440.out","w",stdout); init(); INF=41000; INF=INF*INF; scanf("%d",&Case); while (Case--){ scanf("%d",&K); l=1; r=INF; while (l<=r){ mid=((long long)l+r)>>1; if (Judge()) ans=mid, r=mid-1; else l=mid+1; } printf("%d\n",ans); } return 0; }
2、[BZOJ 2440]中山市選2011 完全平方數
好題啊~~~~~~
首先必須想到一個利器——二分！
我們二分出一個數怎樣檢驗有多少個完全平方數比他小呢？
從莫比烏斯函式的定義來看
接下來利用容斥從n中減去為平方數倍數的數。
如果一個數有平方因子，莫比烏斯函式值為0，不必考慮。
如果不為0，實質上就是容斥，減去mu[i]*n/(i*i)即可
#include <cmath> #include <cstdio> #include <cstring> #include <algorithm> using namespace std; int l,r,ans,mid,INF,K,Case,tot; int prime[41005],mu[41005],i,j; bool check[41005]; void init(){ mu[1]=1; for (i=2;i<41000;i++){ if (!check[i]) mu[i]=-1, prime[++tot]=i; for (j=1;j<=tot;j++){ if (i*prime[j]>41000) break; check[i*prime[j]]=1; if (i%prime[j]==0) {mu[i*prime[j]]=0;break;} else mu[i*prime[j]] = -mu[i]; } } } bool Judge(){ int ret=0; for (int i=1;i*i<=mid;i++) ret+=(mid/i/i)*mu[i]; return (ret>=K); } int main(){ //freopen("2440.in","r",stdin); //freopen("2440.out","w",stdout); init(); INF=41000; INF=INF*INF; scanf("%d",&Case); while (Case--){ scanf("%d",&K); l=1; r=INF; while (l<=r){ mid=((long long)l+r)>>1; if (Judge()) ans=mid, r=mid-1; else l=mid+1; } printf("%d\n",ans); } return 0; }
3、gcd
（由於網上沒找到原題，就簡述一下題意吧）
給你n個正整數，以及一個正整數k。然後有q個詢問，每次詢問一個x。要求回答從這n個數中選出正好k個數並且使得他們的最大公約數是x的方案數是多少。輸出答案mod 10^9+7值。n,q,x<=10^6,k<=n。時限：5s
還是容斥。。。。。
sum[]記錄數字i的約數有多少個，這樣方便計算。
預處理計算對於每一個x的答案ans[x]，用到組合之類的東西。
複雜度大約是nlogn。
詢問直接輸了，具體還是見程式碼吧。
#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
typedef long long LL;
const int Maxn=1000001, Mod=(1e9)+7;;
#define C(n,m) ((LL)jc[n]*ny[m]%Mod*ny[(n)-(m)]%Mod)
int check[Maxn],jc[Maxn],ny[Maxn],mu[Maxn],prime[Maxn];
int sum[Maxn],i,j,tot,x,n,k,Case;
LL ans[Maxn];
int qck(int a,int b){
int ret=1;
for (;b>0;b>>=1){
if (b&1) ret=(LL)ret*a%Mod;
a=(LL)a*a%Mod;
}
return ret;
}
void init(){
mu[1]=1;
for (i=2;i<Maxn;i++){
if (check[i]==0) { mu[i]=-1;prime[++tot]=i; }
for (j=1;j<=tot;j++){
if (prime[j]*i>=Maxn) break;
check[prime[j]*i]=prime[j];
if (i%prime[j]) mu[prime[j]*i]=-mu[i];
else { mu[prime[j]*i]=0;break; }
}
if (check[i]==0) ny[i]=qck(i,Mod-2);
else ny[i]=(LL)ny[i/check[i]]*ny[check[i]]%Mod;
}
jc[0]=ny[0]=1;
jc[1]=ny[1]=1;
for (i=2;i<Maxn;i++){
jc[i]=(LL)jc[i-1]*i%Mod;
ny[i]=(LL)ny[i]*ny[i-1]%Mod;
}
}
int main(){
freopen("gcd.in","r",stdin);
freopen("gcd.out","w",stdout);
init();
scanf("%d%d",&n,&k);
for (i=1;i<=n;i++){
scanf("%d",&x);
sum[x]++;
}
for (i=1;i<Maxn;i++)
for (j=i+i;j<Maxn;j+=i)
sum[i]=sum[j]+sum[i];
for (i=1;i<Maxn;i++)
for (j=1;i*j<Maxn;j++){
if (sum[i*j]<k) continue;
ans[i]=(ans[i]+(LL)mu[j]*C(sum[i*j],k)+Mod)%Mod;
}
scanf("%d",&Case);
while (Case--){
scanf("%d",&x);
printf("%I64d\n",ans[x]);
}
return 0;
}
鳴謝zy給我的啟發！
如果0<x<n，x與n互質，那麼x+n也一定與n互質。
否命題也成立。
那麼我們只有求出在m！中與m！互質的數，也就是phi(m!)。
答案就是n!/m!*phi(m!)。
關於計算答案我用了些什麼線性逆元亂七八糟的東西，理論上覆雜度還可以，實際上被卡常數了，跑得很慢的。
trick：雖然題目上沒說，但是資料上是這麼給的：R是一個比較大的質數（大於所有n，m）
#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
typedef long long LL;
const int Maxn=10000001;
int Case,Mod,check[Maxn],prime[Maxn],tot,i,j;
int phi[Maxn],jc[Maxn],ny[Maxn],n,m,l,r,mid,t,ans;
int qck(int a,int b){
int ret=1;
for (;b>0;b>>=1){
if (b&1) ret=(LL)ret*a%Mod;
a=(LL)a*a%Mod;
}
return ret;
}
void init(){
jc[0]=1; ny[0]=1;
jc[1]=1; ny[1]=1;
for (i=2;i<Maxn;i++){
if (check[i]==0) prime[++tot]=i;
for (j=1;j<=tot;j++){
if (i*prime[j]>=Maxn) break;
check[i*prime[j]]=prime[j];
if (i%prime[j]==0) break;
}
jc[i]=(LL)jc[i-1]*i%Mod;
if (check[i]==0) ny[i]=qck(i,Mod-2);
else ny[i]=(LL)ny[check[i]]*ny[i/check[i]]%Mod;
}
for (i=1,phi[0]=1;i<=tot;i++)
phi[i]=(LL)phi[i-1]*(prime[i]-1)%Mod*ny[prime[i]]%Mod;
for (i=1;i<Maxn;i++)
ny[i]=(LL)ny[i-1]*ny[i]%Mod;
}
int main(){
//freopen("2186.in","r",stdin);
//freopen("2186.out","w",stdout);
scanf("%d%d",&Case,&Mod);
init();
while (Case--){
scanf("%d%d",&n,&m);
if (m>n) {printf("0\n");continue;}
l=1; r=tot; t=0;
while (l<=r){
mid=(l+r)>>1;
if (prime[mid]<=m) t=mid, l=mid+1;
else r=mid-1;
}
ans=(LL)jc[m]*phi[t]%Mod*jc[n]%Mod*ny[m]%Mod;
printf("%d\n",ans);
}
return 0;
}5、[BZOJ 3529]SDOI2014 數表
（終於到最後一題了
）
其實題目裡把a給很大是無用的，算一下表格內最大的數字不會超過50W
以sum[]為關鍵字排序，逐步更新G[]，用Problem b的方法分段優化時間，用樹狀陣列維護G[]，複雜度大約就是O( q*logn*sqrt(n) )
貼上醜醜的程式碼：
#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
#define lowbit(x) ((x)&(-(x)))
int ans[20005],tr[100005],prime[100005];
int i,j,k,kk,N,tot,mu[100005];
bool check[100005];
struct query
{
int n,m,k,num;
bool operator <(const query &a)const
{ return k<a.k; }
} q[20005];
struct divsum
{
int x,y;
bool operator <(const divsum &a)const
{ return y<a.y; }
}sum[100005];
void init(){
for (i=1;i<=100000;i++){
for (j=1;j*j<=i;j++)
if (i%j==0){
sum[i].y+=j;
if (i/j!=j) sum[i].y+=i/j;
}
sum[i].x=i;
}
sort(sum+1,sum+100001);
mu[1]=1;
for (i=2;i<=100000;i++){
if (check[i]==0) {prime[++tot]=i;mu[i]=-1;}
for (j=1;j<=tot;j++){
if (prime[j]*i>100000) break;
check[prime[j]*i]=1;
if (i%prime[j]==0) {mu[i*prime[j]]=0;break;}
else mu[i*prime[j]]=-mu[i];
}
}
}
void ins(int x,int y){
for (int i=x;i<=100000;i+=lowbit(i))
tr[i]+=y;
}
int get(int x){
int ret=0;
for (int i=x;i>0;i-=lowbit(i))
ret+=tr[i];
return ret;
}
int get(int l,int r)
{ return get(r)-get(l-1);}
int main(){
//freopen("table.in","r",stdin);
//freopen("table.out","w",stdout);
init();
scanf("%d",&N);
for (i=1;i<=N;i++){
scanf("%d%d%d",&q[i].n,&q[i].m,&q[i].k);
q[i].num=i;
}
sort(q+1,q+N+1);
for (i=1,j=1;i<=N;i++){
for (;j<=100000 && sum[j].y<=q[i].k;j++)
for (k=1;k*sum[j].x<=100000;k++)
ins(k*sum[j].x,sum[j].y*mu[k]);
for (k=1;k<=q[i].n&&k<=q[i].m;k=kk+1){
kk=min( q[i].n/(q[i].n/k), q[i].m/(q[i].m/k) );
ans[q[i].num]+=get(k,kk)*(q[i].n/k)*(q[i].m/k);
}
if (ans[q[i].num]<0){
ans[q[i].num]+=2147483647;
ans[q[i].num]++;
}
}
for (i=1;i<=N;i++)
printf("%d\n",ans[i]);
return 0;
}——————————————————————————————————————————
最後戳兩個題目，以後做吧：
[bzoj 2154]Crash的數字表格
[bzoj 3434]wc2014時空穿梭