zoukankan      html  css  js  c++  java
  • 莫比乌斯反演

    重新学习了一遍莫比乌斯反演,整理一下。

    莫比乌斯函数

    莫比乌斯函数(mu)是一个积性函数。

    [mu(x)=egin{cases}1 &(x=1)\ (-1)^k & x=p_1p_2...p_k\ 0 & elseend{cases}]

    即对于一个数(x)的莫比乌斯函数分三种情况讨论
    (1.)(x=1)时,(mu(x)=1)
    (2.)(x)为互异素数的乘积时(mu(x)=(-1)^k)(k)(x)分解后的素数个数。
    (3.)其他情况(mu(x)=0)
    莫比乌斯函数是个积性函数,所以可以线性筛。

    性质

    (1.) $$sum_{d|n}mu(d)=[n=1]$$
    (2.) $$sum_{d|n}frac{mu(d)}{d}=frac{phi(n)}{n}$$

    线性筛莫比乌斯函数代码

    	for(int i = 2;i <= END;++i) {
    		if(!vis[i]) {
    			mu[i] = -1;
    			dis[++js] = i;
    		}
    		for(int j = 1;j <= js && dis[j] * i <= END;++j) {
    			vis[i * dis[j]] = 1;
    			if(i % dis[j]) mu[i * dis[j]] = -mu[i];
    			else {
    				mu[i * dis[j]] = 0;
    				break;
    			}
    		}
    	}
    

    莫比乌斯反演

    其实莫比乌斯反演就是两个公式
    (F(n))(f(n))是定义在非负整数上的两个函数
    莫比乌斯反演仅仅靠这两个公式是很难做题的,做完后面两道例题应该就算入门了

    公式一

    (F(n))(f(n))满足

    [F(n)=sumlimits_{d|n}f(d) ]

    则有$$f(n)=sumlimits_{d|n}mu(d)F(lfloorfrac{n}{d} floor)$$

    公式二

    (F(n))(f(n))满足

    [F(n)=sumlimits_{n|d}f(d) ]

    则有

    [f(n)=sumlimits_{n|d}mu(frac{d}{n})F(d) ]

    证明

    留坑待填吧(233)

    题目

    例题1

    求出有多少(x,y)满足(1le xle Bquad 1le y le D)(gcd(x,y) = k)

    评测戳这里
    这是非常基础的一道题目。为了便于讨论,我们下面都假定(Ble D)
    我们设(f(n))表示(gcd(x,y)=n)的数量,(F(n))表示(gcd(x,y)\% n= 0)的数量
    显然有下面的式子

    [F(n)=sumlimits_{n|d}^Bf(d) ]

    然后根据公式二,可以反演出下面的式子

    [f(n)=sumlimits_{n|d}^Bmu(frac{d}{n})F(d) ]

    现在(f(k))就是答案,如果我们让(B)(D)都先除以(k),那么(f(1))就是答案。
    所以下面我们假定(B=frac{B}{k},D=frac{D}{k})
    现在就是求$$f(1)=sumlimits_{d=1}^Bmu(d)F(d)$$
    然后现在只要能够快速的计算出(F(d))就能够(O(B))的计算了。
    显然有$$F(i)=lfloorfrac{B}{i} floorlfloorfrac{D}{i} floor$$
    所以$$f(1)=sumlimits_{d=1}^Bmu(d)lfloorfrac{B}{d} floorlfloorfrac{D}{d} floor$$
    然后计算就行了。因为题目中说要去重,发现被计算两遍的那些数字是(x,y) 都小于(B)的情况,所以最后在减去就行了。
    代码如下

    /*
    * @Author: wxyww
    * @Date:   2019-02-22 11:00:34
    * @Last Modified time: 2019-02-22 11:24:13
    */
    #include<cstdio>
    #include<iostream>
    #include<cstdlib>
    #include<cstring>
    #include<algorithm>
    #include<queue>
    #include<vector>
    #include<ctime>
    using namespace std;
    typedef long long ll;
    #define int ll
    const int N = 100000 + 100;
    ll read() {
    	ll x=0,f=1;char c=getchar();
    	while(c<'0'||c>'9') {
    		if(c=='-') f=-1;
    		c=getchar();
    	}
    	while(c>='0'&&c<='9') {
    		x=x*10+c-'0';
    		c=getchar();
    	}
    	return x*f;
    }
    int mu[N],dis[N],vis[N];
    void pre() {
    	mu[1] = 1;
    	int js = 0;
    	int END = 100000;
    	for(int i = 2;i <= END;++i) {
    		if(!vis[i]) {
    			mu[i] = -1;
    			dis[++js] = i;
    		}
    		for(int j = 1;j <= js && dis[j] * i <= END;++j) {
    			vis[i * dis[j]] = 1;
    			if(i % dis[j]) mu[i * dis[j]] = -mu[i];
    			else {
    				mu[i * dis[j]] = 0;
    				break;
    			}
    		}
    	}
    }
    int solve(int x,int y) {
    	int ans = 0;
    	for(int i = 1;i <= min(x,y);++i)
    		ans += mu[i] * (x / i) * (y / i);
    	return ans;
    }
    signed main() {
    	int T = read();
    	pre();
    	for(int i = 1;i <= T;++i) {
    		read();int n = read();read();int m = read();int K = read();
    		int z = min(n,m);	
    		if(K != 0) 
    		printf("Case %lld: %lld
    ",i,solve(n / K,m / K) - solve(z / K,z / K) / 2);
    		else printf("Case %lld: 0
    ",i);
    	}
    
    	return 0;
    }
    

    例题2

    直接看原题面吧

    首先一个明显的转化就是,转化成对于一个数字(k),有多少(1le x le n,1le yle m)满足(gcd(x,y)=k)。然后对于每个(k),计算答案就行了。
    发现上面这个式子不就是上面那道题么。同样设(n le m)
    我们设(f(i))表示(gcd(x,y)=i)时的答案。(fb_i)斐波那契数列第(i)项那么这道题最终的答案就是。$$prodlimits_{i=1}^nfb_i^{f(i)}$$
    根据上一题推出的式子,我们可以推出来
    原式=$$prod_{i=1}^nfb_i^{sumlimits_{d=1}^{n/i}mu(d)lfloorfrac{n}{id} floorlfloorfrac{m}{id} floor}$$
    (T=id)
    原式=

    [prodlimits_{T=1}^nprodlimits_{d|T}fb_d^{mu(T/d)lfloorfrac{n}{T} floorlfloorfrac{m}{T} floor} ]

    到了这里,对于外面的(prodlimits_{T=1}^n...^{lfloorfrac{n}{T} floorlfloorfrac{m}{T} floor})可以数论分块做到(sqrt{m})
    对于里面(prodlimits_{d|T}fb_d^{mu(T/d)}),暴力预处理是(nlog)的(调和级数)。
    代码如下

    /*
    * @Author: wxyww
    * @Date:   2019-02-22 15:09:25
    * @Last Modified time: 2019-02-24 08:45:10
    */
    #include<cstdio>
    #include<iostream>
    #include<cstdlib>
    #include<cstring>
    #include<algorithm>
    #include<queue>
    #include<vector>
    #include<ctime>
    using namespace std;
    typedef long long ll;
    #define int ll
    const int N = 1000000 + 10,mod = 1e9 + 7;
    ll read() {
    	ll x=0,f=1;char c=getchar();
    	while(c<'0'||c>'9') {
    		if(c=='-') f=-1;
    		c=getchar();
    	}
    	while(c>='0'&&c<='9') {
    		x=x*10+c-'0';
    		c=getchar();
    	}
    	return x*f;
    }
    int f[N],mu[N],dis[N],vis[N],g[N];
    int qm(ll x,int y) {
    	int ret = 1;
    	if(y < 0) y += mod - 1;
    	for(;y;y >>= 1,x = 1ll * x * x % mod) 
    		if(y & 1) ret = 1ll * ret * x % mod;
    	return ret;
    }
    int tmp[5];
    void pre() {
    	int END = 1000000;
    	mu[1] = 1;
    	f[1] = 1,g[0] = 1;
    	g[1] = 1;
    	int js = 0;
    	for(int i = 2;i <= END;++i) {
    		f[i] = f[i - 1] + f[i - 2];
    		g[i] = 1;
    		f[i] >= mod ? f[i] -= mod : 0;
    		if(!vis[i]) {
    			dis[++js] = i;
    			mu[i] = -1;
    		}
    		for(int j = 1;j <= js && dis[j] * i <= END;++j) {
    			vis[i * dis[j]] = 1;
    			if(i % dis[j]) mu[i * dis[j]] = -mu[i];
    			else break;
    		}
    	}
    	for(int i = 1;i <= END;++i) {
    		tmp[0] = qm(f[i],-1);
    		tmp[1] = 1;
    		tmp[2] = f[i];
    		for(int j = i;j <= END;j += i)
    			g[j] = 1ll * g[j] * tmp[mu[j / i] + 1] % mod;
    	}
    	for(int i = 1;i <= END;++i) g[i] = 1ll * g[i - 1] * g[i] % mod;
    }
    signed main() {
    	pre();
    	int T = read();
    	while(T--) {
    		int n = read(),m = read();
    		if(n > m) swap(n,m);
    		int r;
    		ll ans = 1;
    		for(int l = 1;l <= n;l = r + 1) {
    			r = min(n / (n / l),m / (m / l));
    			ans = ans * qm(g[r] * qm(g[l - 1],mod - 2) % mod , (n / l)  * (m / l) % (mod - 1)) % mod;
    		}
    		printf("%lld
    ",ans);
    	}
    	return 0;
    }
    
  • 相关阅读:
    中国软件外包IT公司最新排名
    DJ舞曲、音乐与爱好!
    linux论坛
    IBM P 系列小型机的控件面板功能~!(转用)
    JDBC Driver 驱动 For SQL 2008 Server /2005 /2000
    年报盘点:149家公司筹码大幅集中
    公式指标—精华
    观峰雨个人空间 2010 STOCK ADVICE !
    IntelliJ IDEA提示Cannot resolve symbol
    今天天变的好冷了~
  • 原文地址:https://www.cnblogs.com/wxyww/p/mofan.html
Copyright © 2011-2022 走看看