容斥原理在程序设计中的应用

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);
}