容斥原理在程序設計中的應用

本文轉載自查看原文 2019-07-22 19:56 440 組合數學

1.定義

在計數時，必須注意沒有重復，沒有遺漏。為了使重疊部分不被重復計算，人們研究出一種新的計數方法，這種方法的基本思想是：先不考慮重疊的情況，把包含於某內容中的所有對象的數目先計算出來（容），然后再把計數時重復計算的數目排斥出去（斥），使得計算的結果既無遺漏又無重復，這種計數的方法稱為容斥原理。

舉個栗子，如果被計數的事物有A、B、C三類，那么，A類和B類和C類元素個數總和= A類元素個數+ B類元素個數+C類元素個數—既是A類又是B類的元素個數—既是A類又是C類的元素個數—既是B類又是C類的元素個數+既是A類又是B類而且是C類的元素個數。

（A∪B∪C = A+B+C - A∩B - B∩C - C∩A + A∩B∩C）

推導公式：

其實對於(2)，是在(1)的基礎上運用了Demorgen定理。

2.嚴格證明

對於容斥原理我們可以利用 數學歸納法證明：

證明：當

時，等式成立(證明略)。

假設

時結論成立，則當

時，

所以當

時，結論仍成立。因此對任意

，均可使所證等式成立。

注：以上證明摘自曹汝成．組合數學：華南理工大學出版社，2000年

3.應用

　　容斥原理的目的：找到不重復的總區間。

　　| A_i ∪ A_j ... ∪ A_n | = ∑ | A_i| - ∑ | A_i ∩ A_j | + ∑ | A_i ∩ A_j ∩ A_k | - ... + (-1)ⁿ × ∑ | A_i ∩ ... ∩ A_n |

1.How many integers can you find HDU - 1796

題目大意：給定一個n和一個有m個元素的序列，求小於n並且能被這m個元素中整除的整數有多少個。

分析：在區間[1,n)中，找到能被每個元素整除的元素個數，用容斥原理找到不重復元素的個數。（其實就是題2省去分解質因數的簡化版）

實現：求最大公約數，最小公倍數：

ll gcd(ll a,ll b){
    if(!b)return a;
    return gcd(b,a % b);
}
ll lcm(ll a,ll b){
    ll p = a * b / gcd(a,b);
    return p;
}

　　DFS容斥原理：按照圖示層次遞歸

ll DFS(ll i,ll w,int k,ll m){
    ll p;
    for(;i<=top;i++){
        p = lcm(a[i],w);
        ans += k * (m/p);
        DFS(i+1,p,-k,m);
    }
    return ans;
}

2.Co-prime HDU - 4135

題目大意：給定一個區間[A,B],找出在這個區間內與給定的n互質的整數的個數。

分析：直接求與n互質的數不好找，我們選擇先找不與n互質的數，將n分解質因數，再通過容斥原理找到給定區間中所有質因數的倍數。

舉個栗子，n為30，那么質因數為2、3、5, 2在[1,30]中的倍數與3在[1,30]中的倍數難免有重復，這也不難解釋我們為什么會用容斥原理了。

設每個質因數在區間中的倍數分別有A₁、A₂、A₃ ... A_n個，那么∑ | A_i ∩ A_j |就是i、j的最小公倍數在區間中的倍數個數。

| S | = ∑ | A_i| - ∑ | A_i ∩ A_j | + ∑ | A_i ∩ A_j ∩ A_k | - ... + (-1)ⁿ × ∑ | A_i ∩ ... ∩ A_n |

質因數分解

void div(ll n){
    for(int i=2; i*i<=n;i++){
        if(n % i == 0){
            a[++top] = i;
            while(n % i == 0)n /= i;
        }
        if(n == 1)break;
    }    
    if(n > 1)a[++top] = n;
}

3.Card Collector HDU - 4336

題目大意：N個物品，每次得到第i個物品的概率為pi，而且有可能什么也得不到，問期望多少次能收集到全部N個物品。

分析：這里牽涉到概率論上的一些基本概念——期望。舉個栗子，兩個元素的期望值，實際就是A1∩A2的期望值，而這里每個事件顯然是互相獨立的，那么A1∩A2顯然等於A1 + A2，因此這兩個元素的期望值應該是1/(A1 + A2),

我們可以看到，A1的期望顯然和A1∩A2的期望是有重疊部分，那么這里就解釋了為什么會用到容斥原理了。

void DFS(int i,double w,int k){
    double p;
    for(;i<=n;i++){
        p = a[i] + w;
        ans += k * (1 / p);
        DFS(i+1,p,-k);
    }
}

4.Visible Trees HDU - 2841 、儀仗隊 Luogu P2158 [SDOI2008]

HDU——2841 一個n*m的方格紙，格點上有樹，求從(0,0)處看，能看到幾棵樹。

P2158 [SDOI2008] 作為體育委員，C君負責這次運動會儀仗隊的訓練。儀仗隊是由學生組成的N * N的方陣，為了保證隊伍在行進中整齊划一，C君會跟在儀仗隊的左后方，根據其視線所及的學生人數來判斷隊伍是否整齊(如下圖) 在，C君希望你告訴他隊伍整齊時能看到的學生人數。

分析：如何判斷一個人是否被遮擋？它與觀測點連線的斜率，和離觀測點更近的一個人與觀測點連線的斜率相等。

舉個栗子，觀察上圖，后面的人顯然就被前面的人遮擋了。若觀測點為（0,0），前一個數為（2,1），后一個數為（4,2）。那么后一個人就被前一個人遮擋了。

如果一個一個點求斜率再比較，數據范圍 m,n<=100,000，O(mn)算法是超時算法。

抽象一點說，若一個人的坐標有不為1的公因數，即坐標為（x*q,y*q），q ∈ N^* ,那么勢必會被坐標為（x,y）的人遮擋。

所以就是求解x∈[1,m],y∈[1,n]中互質的數的個數。以列來枚舉，列數為i,i∈[1,m],則就是求[1,n]中與i互質的數。

如何求該區間內與i互質的數呢？還是找到i的所有質因數，利用容斥定理求解。

其實第二題還有更巧妙的解法。在找到i的質因數后，利用歐拉函數求解。

由於第二題圖是n×n的正方形，所以將圖分為對稱的兩半，則是求[1,i]中與i互質的數的個數。我們知道歐拉函數的定義φ(x)為在小於等於x的數中與x互質的數的個數，那么通過歐拉函數計算公式，就可以更快，更直接地求解了。

附蒟蒻代碼一篇：

// P2158
#include<cstdio>
#include<cstring>
#include<iostream>
#include<cmath>
using namespace std;
// 歐拉函數模板 
int main(){
    int k,n,t = 0;
    long long ans = 0;
    cin>>k;// k*k 的正方形矩陣 
    for(int j = 2; j < k; j++)// 高從2到k-1,為什么是k-1不是k?
    //因為最后一列是k個節點，但我們判斷是否互質是根據實際長度，即節點之間的距離判斷的 
    {
        n = j;//一個小細節，由於找質因數時n被除值會改變，所以不直接用 j,以免影響外循環 
        long long int anss = 1;
        for(int i=2; i*i<=j; i++){
            if(n % i == 0){
                t = 0;
                while(n % i == 0)n /= i,t++;
                anss *= (long long int)(i-1) * (int)pow(i,t-1);
            }
            if(n <= 1)break;
        }
        if(n > 1)anss *= (n-1) ;
        ans += anss;
    }
    if(k == 1)printf("0");// 1*1的矩陣要特判，此時只能看到0個人 
    else
        printf("%lld",ans*2 + 3);// +3是考慮到直線x = 0,y = 0,y = x都只能看到一個點 
    
    return 0;
}

5.被破壞的玉米地

輸入一個整數n（0≤n≤200）,表示圓圈的個數。以下n行每行有兩個整數x和y，由空格分開，代表圓心坐標。

當兩個圓的圓心坐標為（0，0）和（1，0）時圖3.1兩個圓所覆蓋的面積為5.0548。編寫程序需要輸出統計的總面積，四舍五入到小數點后四位。

分析：

題意要求求出覆蓋的總面積，考慮到圓之間可能有相交的情況，我們不能僅僅求所有圓的面積之和。圓之間可能有2個圓相交的部分、3個圓相交的部分、4個圓相交的部分......

但如何求圓之間相交的面積呢？我們分類討論來看一看。

一、2個圓相交

　　我們稱這種為α型交

　　　　稱這種為β型交

二、3個圓相交

三、4個圓相交

四、幸運的是，由於四個圓大小相同，沒有五個圓互不重疊並且有共同相交面積的情況。

接下來就是判斷圓之間的相交類型，並通過容斥原理求覆蓋總面積。

由於代碼太過繁雜，蒟蒻摘抄了大佬代碼：

#include <stdlib.h>
#include <stdio.h>
#include <iostream>
#include <math.h>
using namespace std; 
#define MAX 80
const double pi=3.14159265358979324;
/**原點坐標**/
struct coordinate
{
    int x;
    int y;
};
/**計算任意兩個坐標之間的距離”,這里的距離定義為橫坐標之差與縱坐標之差的和**/
int f1(coordinate k1,coordinate k2)
{
    int l;
    if(abs(k1.x-k2.x)==0&&abs(k1.y-k2.y)==2){
        return 0;
    }
    if(abs(k1.x-k2.x)==2&&abs(k1.y-k2.y)==0){
        return 0;
    }
    l=abs(k1.x-k2.x)+abs(k1.y-k2.y);
    if (l>2)    return 0;
    else return l;
}
/**快速排序算法的比較函數**/
int compare( const void *a, const void *b )
{
    coordinate *arg1=(coordinate*)a;
    coordinate *arg2=(coordinate*)b;
    if (arg1->x<arg2->x)
        return -1;
    else if(arg1->x>arg2->x)
        return 1;
    else if (arg1->y<arg2->y)
        return -1;
    else if(arg1->y>arg2->y)
        return 1;
    else return 0;        
};

int main()
{
    int i,n;
    double area;
    coordinate *a;
 
    double duration;
    
    //輸入圓心的個數
    cin>>n;
    
    a=(coordinate*)malloc(sizeof(coordinate)*n);
    qsort(a,n,sizeof(coordinate),compare);

    int *List;
    List=(int*)malloc(sizeof(int)*n*n);
    
    int j,k,l,count1=0,count2=0,count3=0,count4=0;
    area=n*pi;
    
    //四種情況的相交面積 
    double s1=(2.0/3.0)*pi-sqrt(3.0)/2.0;
    double s2=pi/2.0-1.0;
    double s3=5.0/12.0*pi-sqrt(3.0)/2.0;
    double s4=pi/3.0-1+0.25/(sin(75*pi/180)*sin(75*pi/180));
    
    //分析兩個圓相交的情況 
    for (i=0;i<n-1;i++)
        for (j=i+1;j<n;j++)
        {
            *(List+n*i+j)=f1(a[i],a[j]);
            *(List+n*j+i)=*(List+n*i+j);
            if (*(List+n*i+j)==1)
                count1++;
            else if ((*(List+n*i+j)==2)&&(abs(a[i].x-a[j].x)==1))
                count2++;
        }
        
    //分析三個圓或者四個圓相交的情況 
    for (i=0;i<n-2;i++)
        for (j=i+1;j<n-1;j++)
            for (k=j+1;k<n;k++)
            {
                bool check=true;
                int ans=*(List+n*i+j)+*(List+n*j+k)+*(List+n*i+k);
                if(*(List+n*i+j)==0||*(List+n*j+k)==0||*(List+n*i+k)!=0)
                if(ans==4&&check)
                {
                    count3++;
                }
            }
 
    for (i=0;i<n;i++)
        for (j=0;j<n;j++)
            for (k=0;k<n;k++)
                for(l=0;l<n;l++){
                        if ((*(List+n*i+j)==1)&&(*(List+n*j+k)==1)&&(*(List+n*k+l)==1)&&(*(List+n*i+l)==1)&&(i!=j)&&(i!=k)&&(i!=l)&&(j!=k)&&(j!=l)&&(k!=l))
                            count4++;
                }        
    
    area=area-s1*count1-s2*count2+s3*count3-s4*count4/8;
 
    free(a);
    free(List);
}