组合数学总结

组合数学总结

（转发请注明出处）

1.基础知识

(1)小数据范围直接预处理求组合数

例1.hdu 1799 循环多少次?

题意：

　　中文题目..不解释.

思路：

　　由于每一层从上一层+1开始,所以对于每一种瞬时的状态,对应C(n,m)中的一种情况,所以答案就是C(n,m),由于数据量比较小,直接预处理即可

代码：

View Code

 1 #include<stdio.h>
 2 int c[2010][2010];
 3 
 4 void init()
 5 {
 6     int i,j;
 7     c[0][0]=c[1][0]=c[1][1]=1;
 8     for(i=2;i<=2000;i++)
 9     {
10         for(j=0;j<=i;j++)
11         {
12             if(j==0)
13                 c[i][j]=1;
14             else
15                 c[i][j]=(c[i-1][j]+c[i-1][j-1])%1007;
16         }
17     }
18 }
19 int main()
20 {
21     int t,m,n;
22     scanf("%d",&t);
23     init();
24     while(t--)
25     {
26         scanf("%d%d",&m,&n);
27         printf("%d\n",c[n][m]);
28     }
29     return 0;
30 }

(2)大组合数取模(C(n,m)%p）

m,n巨大,p比较小的时候（p<=10^5)，且p为素数.

例1.hdu 3037 Saving Beans

题意：

　　给n棵树分配不大于m个果子,每棵树可以没有果子，问有多少种分配方法.

思路：

　　容易想到恰好分配m个果子的情况，即多重集合的组合C(n+m-1,m),但是本题是求小于等于m的所有情况,即求C(n+0-1,0)+C(n+1-1,1)+C(n+2-1,2)+C(n+3-1,3)...+C(n+m-1,m)的和,把第一项变形得ans=C(n,0)+C(n+,1)+C(n+1,2)+C(n+2,3)...+C(n+m-1,m),有由公式C(n,m) = C(n-1m)+C(n-1,m-1)得ans=C(n+m,m).

　　知道了公式以后,本题的最大难点在于对大组合数求解,需要用到Lucas定理："C(n,m)%p的值等于将n,m化为p进制数后，各对应数对的组合数的连乘积%p.",知道这个定理后就可以每次对阶乘预处理，快速求出每一位的组合数大小,继而求得相乘后的答案.

代码：

扩展欧几里得求模逆元写法：

View Code

 1 #include<stdio.h>//From  www.cnblogs.com/SolarWings/
 2 
 3 long long factor[100010];
 4 
 5 void init(long long p)//预处理阶乘,化为p进制后,C(a,b)中a,b<p,所以处理到p就足够了
 6 {
 7     long long i;
 8     factor[0]=1;
 9     for(i=1;i<=p;i++)
10     {
11         factor[i]=(factor[i-1]*i)%p;
12     }
13 }
14 
15 long long Extended_Euclid(long long a,long long b,long long &x,long long &y)//扩展欧几里得
16 {
17     long long d;
18     if(b==0)
19     {
20         x=1;y=0;return a;
21     }
22     d=Extended_Euclid(b,a%b,y,x);
23     y-=a/b*x;
24     return d;
25 }
26 long long Inv(long long a,long long n)//求模逆元
27 {
28     long long d,x,y;
29     d=Extended_Euclid(a,n,x,y);//需要用到扩展欧几里得,或者也可以用快速幂直接得到模拟元,以后的题目中再做说明
30     if(d==1) return (x%n+n)%n;
31     else return -1;
32 }
33 
34 
35 long long get(long long a,long long b,long long p) //c(n,m)=n!/(m!*(n-m)!)=n!*Inv(m!)*Inv((m-n)!),Inv代表模逆元
36 {
37     long long ans;
38     if(b>a)
39         return 0;
40     ans=factor[a];
41     ans=(ans*Inv(factor[b],p))%p;
42     ans=(ans*Inv(factor[a-b],p))%p;
43     return ans;
44 }
45 
46 long long c(long long n,long long m,long long p) //这里用到了Lucas定理
47 {
48     long long ans,a,b; ans=1;
49     while(m) //化为p进制分别求组合数
50     {
51         a=n%p;
52         b=m%p;
53         ans=(ans*get(a,b,p))%p;
54         m/=p;
55         n/=p;
56     }
57     return ans;
58 }
59 
60 
61 int main()
62 {
63     long long t,n,m,p;
64     scanf("%I64d",&t);
65     while(t--)
66     {
67         scanf("%I64d %I64d %I64d",&n,&m,&p);
68         init(p);//由于每次的p不同,所以要分别预处理
69         printf("%I64d\n",c(n+m,n,p));
70     }
71     return 0;
72 }

快速幂求模逆元写法：

View Code

 1 #include<stdio.h>
 2 long long fac[100010];
 3 void init(long long p)
 4 {
 5     fac[0]=1;
 6     for(long long i=1;i<=p;i++)
 7     {
 8         fac[i]=((long long)fac[i-1]*i)%p;
 9     }
10     return;
11 }
12 
13 long long quickpow(long long m,long long n,long long p)
14 {
15     long long ans=1;
16     while(n)
17     {
18         if(n&1)
19             ans=(ans*m)%p;
20         n=n>>1;
21         m=(m*m)%p;
22     }
23     return ans;
24 }
25 
26 long long Inv(long long x,long long p)
27 {
28     return quickpow(x,p-2,p);
29 }
30 
31 long long C(long long n,long long m,long long p)
32 {
33     if(m>n)
34         return 0;
35     long long ans=1;
36     ans=ans*fac[n];
37     ans=(ans*Inv(fac[n-m],p))%p;
38     ans=(ans*Inv(fac[m],p))%p;
39     return ans;
40 }
41 
42 long long getans(long long n,long long m,long long p)
43 {
44     long long ans=1;
45     init(p);
46     while(n&&m)
47     {
48         long long a=n%p;
49         long long b=m%p;
50         ans=(ans*(long long)C(a,b,p))%p;
51         n/=p;
52         m/=p;
53     }
54     return ans;
55 }
56 
57 int main()
58 {
59     long long t,n,m,p;
60     scanf("%I64d",&t);
61     while(t--)
62     {
63         scanf("%I64d%I64d%I64d",&n,&m,&p);
64         printf("%I64d\n",getans(n+m,m,p));
65     }
66     return 0;
67 }

例2.hdu 3944 DP?

题意：

　　从杨辉三角的顶端往下走,可以走直线或者对角线,求到达n行k列的一个点时,途径所有数字之和的最小值

思路：

　　容易想到,最优的方法是先沿着最外面的1走,最后沿对角线到达目标点,由于杨辉三角的对称性,可以把右边的情况全部转化到左边,即沿着左边的1一直往下走,知道走对角线可以到达目标点为止,最推出公式为：ans=C(n+1,m)+n-m,最后求解的过程和例1类似，不过本题数据量比较大，有100000组数据，也就是说要预处理100000次，注意到p是小于10000的素数,所以预处理出10000以下所有素数的阶乘一定会省下不少时间.

代码：

View Code

 1 #include<stdio.h>
 2 bool mark[10010];
 3 int num[10010],map[10010];
 4 int prim[2000][10010];//最大的素数要9000+
 5 void init()
 6 {
 7     int i,k,j;
 8     k=0;
 9     for(i=2;i<=10000;i++)//素数筛
10     {
11         if(!mark[i])
12         {
13             map[i]=k;
14             num[k++]=i;
15             for(j=i*i;j<=10000;j+=i)
16             {
17                 mark[j]=1;
18             }
19         }
20     }
21     for(i=0;i<k;i++)
22     {
23         int p=num[i];
24         prim[i][0]=1;
25         for(j=1;j<=p;j++)
26         {
27             prim[i][j]=(prim[i][j-1]*j)%p;//预处理出所有素数的阶乘
28         }
29     }
30     return;
31 }
32 
33 int Extended_Euclid(int a,int b,int &x,int &y)
34 {
35     int d;
36     if(b==0)
37     {
38         x=1;y=0;return a;
39     }
40     d=Extended_Euclid(b,a%b,y,x);
41     y-=a/b*x;
42     return d;
43 }
44 
45 int Inv(int a,int n)
46 {
47     int d,x,y;
48     d=Extended_Euclid(a,n,x,y);
49     if(d==1) return (x%n+n)%n;
50     else return -1;
51 }
52 
53 int c(int a,int b,int p)
54 {
55     if(b>a)
56         return 0;
57     int ans=prim[map[p]][a];
58     ans=(ans*(Inv(prim[map[p]][b],p)))%p;
59     ans=(ans*(Inv(prim[map[p]][a-b],p)))%p;
60     return ans;
61 }
62 
63 int get(int n,int m,int p)
64 {
65     int a,b,ans;
66     ans=1;
67     while(m&&n)
68     {
69         a=n%p;
70         b=m%p;
71         ans=(ans*c(a,b,p))%p;
72         m/=p;
73         n/=p;
74     }
75     return ans;
76 }
77 
78 
79 
80 int main()
81 {
82     int n,k,p,cases;
83     cases=0;
84     init();
85     while(scanf("%d%d%d",&n,&k,&p)!=EOF)
86     {
87         cases++;
88         if(k>n/2)//这步转化是把右边的也转化成左边的
89             k=n-k;
90         printf("Case #%d: %d\n",cases,(get(n+1,k,p)+n-k)%p);
91     }
92     return 0;
93 }

m,n可以用数组存下(一般是m,n<=10^6),p不是素数的情况.

例1.hdu 1133 Buy the Ticket(本题不是组合数取模,只是有助于例2的理解）

题意：

　　一群人排队买票,每张票要50元,刚开始售票处没有零钱,排队的总共有m+n个人,其中m个带了50元,n个带了100元,当售票处没有零钱可找时停止售票,问有多少种排队的方法可以把票卖完(即排队到某个位置时,50元的人数大于100元)

思路：

　　关于这个题的解释比较神奇,当时看了好多人得题解才弄明白,现在贴一个别人的解释:这是一个Catalan数的非常经典的应用，买票问题，首先我们用"0"表示用50块买票的人，用“1”表示用100块买票的人，然而假设m=4,n=3,的一个序列是：0110100显然，它不合法，然后我们把他稍微变化一下：把第一个不合法的“1”后面的所有数0位为1， 1位为0；这样我们得到了另一个序列：0111011，显然他也不是合法的，但是在这里我们关注的不是他合不合法！只是说明每个不合法的都有一个这样的序列跟他一一对应！(其实也就是说,m和n组成的不合法序列可以由一个n-1,m+1组成的序列得到)

　　所以我们计算公式就是：合法的排列方式=所有排列方式-非法排列方式

　　我们这里非法排列方式的计算就是：( $C_{m+n}^{m}$ - $C_{m+n}^{m+1}$ )*M!*N!，然而在这题，因为每个人都是不同的，所以还要乘以 M！*N！

所以得出最终方程：

　　F（N）=( $C_{m+n}^{m}$ - $C_{m+n}^{m+1}$ )*M！*N！；

然后再化简一下：

　　F(N)=(M+N)! * (M-N+1)/(M+1)

代码：

View Code

 1 #include<stdio.h>
 2 //本题重点在于明确如何让推出的公式，计算只是简单的大数
 3 int a[201][510],b[510];
 4 void init()//预处理阶乘
 5 {
 6     int i,j,len,temp1,temp;
 7     len=0;
 8     a[0][0]=1;
 9     a[1][0]=1;
10     for(i=2;i<=200;i++)
11     {
12         for(j=0;j<=len;j++)
13         {
14             a[i][j]=a[i-1][j]*i;
15         }
16         temp1=0;
17         for(j=0;j<=len;j++)
18         {
19             temp=a[i][j]+temp1;
20             a[i][j]=temp%10;
21             temp1=temp/10;
22         }
23         while(temp1)
24         {
25             a[i][++len]=temp1%10;
26             temp1/=10;
27         }
28     }
29 }
30 int main()
31 {
32     init();
33     int n,m,cases,temp1,temp,i,j;
34     cases=0;
35     while(scanf("%d%d",&m,&n)!=EOF)
36     {
37         cases++;
38         if(m==0&&n==0)
39             break;
40         printf("Test #%d:\n",cases);
41         if(m<n)
42         {
43             printf("0\n");
44                 continue;
45         }
46         temp1=0;
47         for(i=0;i<=500;i++)//大数乘法计算b=(M+N)!*(M-N+1)
48         {
49             temp=a[m+n][i]*(m+1-n)+temp1;
50             b[i]=temp%10;
51             temp1=temp/10;
52         }
53         while(!b[i])
54         {
55             i--;
56         }
57         temp1=0;
58         for(j=i;j>=0;j--)//大数除法计算b=b/(M+1)
59         {
60             temp=temp1*10+b[j];
61             b[j]=temp/(m+1);
62             temp1=temp%(m+1);
63         }
64         while(!b[i])
65             i--;
66         for(j=i;j>=0;j--)
67         {
68             printf("%d",b[j]);
69         }
70         printf("\n");
71     }
72     return 0;
73 }

例2.hdu 3398 String

题意：

　　和例1基本相同,只不过例1里每个人是不同的,而本题中所有"0"相同,所有"1"相同.

思路：

　　由题意得最后的公式为ans=C(m+n,n)-C(m+n,n+1),化简得ans=(n+1-m)*（m+n)！/(n+1)!/m!

代码：

View Code

 1 #include<stdio.h>//本题的主要思想在于分解质因子
 2 #include<string.h>
 3 int prim[150000],cnt[150000],num;
 4 bool mark[2000010];
 5 
 6 int init()//初始化，筛选素数
 7 {
 8     long long i,j,kiss;
 9     kiss=0;
10     for(i=2;i<=2000000;i++)
11     {
12         if(!mark[i])
13         {
14             prim[kiss++]=i;
15             for(j=i*i;j<=2000000;j+=i)
16             {
17                 mark[j]=1;
18             }
19         }
20     }
21     return kiss;
22 }
23 
24 int getfactor(int x,int flag)//对阶乘分解质因子，表示成2^，3^，5^...
25 {
26     int i;
27     for(i=0;i<num&&prim[i]<=x;i++)
28     {
29         int n=x;
30         while(n)
31         {
32             cnt[i]+=flag*(n/prim[i]);
33             n/=prim[i];
34         }
35     }
36     return i;
37 }
38 
39 long long quickpow(long long m,long long n,long long k)//快速幂求解
40 {
41     long long b=1;
42     while(n)
43     {
44         if(n&1)
45             b=(b*m)%k;
46         n=n>>1;
47         m=(m*m)%k;
48     }
49     return b;
50 }
51 
52 long long solve(int n,int m)
53 {
54     long long ans=1;
55     int maxlen=getfactor(n+m,1);//后面的参数为1，表示乘法，系数为曾加
56     getfactor(n+1,-1);//后面参数为2，表示除法，系数为减少
57     getfactor(m,-1);
58     for(int i=0;i<maxlen;i++)
59     {
60         ans=(ans*(quickpow(prim[i],cnt[i],20100501)))%20100501;
61     }
62     return ans;
63 }
64 
65 
66 
67 
68 int main()
69 {
70     int t,m,n,i;
71     num=init();
72     scanf("%d",&t);
73     while(t--)
74     {
75         scanf("%d%d",&n,&m);
76         memset(cnt,0,sizeof(cnt));
77         int temp=n+1-m;
78         for(i=0;i<num&&prim[i]<=temp;i++)//先单独对不用阶乘的(n+1-m)处理
79         {
80             while(temp%prim[i]==0)
81             {
82                 cnt[i]++;
83                 temp/=prim[i];
84             }
85         }
86         printf("%lld\n",solve(n,m));
87     }
88     return 0;
89 }

(3)容斥原理

前置知识--求解欧拉函数

例1.hdu 2824

题意：

　　就是最裸的欧拉函数,只不过是要求a-b的欧拉函数之和,预处理出1-n的所有欧拉函数即可　　

思路：

　　预处理求解

代码：

View Code

 1 #include<stdio.h>
 2 long long euler[3000010];
 3 void init()
 4 {
 5     int i,j;
 6     euler[1]=1;
 7     for(i=2;i<=3000000;i++)
 8     {
 9         if(!euler[i])
10         {
11             for(j=i;j<=3000000;j+=i)
12             {
13                 if(!euler[j])
14                     euler[j]=j;
15                 euler[j]=euler[j]*(i-1)/i;
16             }
17         }
18         euler[i]+=euler[i-1];
19     }
20 }
21 
22 int main()
23 {
24     int a,b;
25     init();
26     while(scanf("%d%d",&a,&b)!=EOF)
27     {
28         printf("%I64d\n",euler[b]-euler[a-1]);
29     }
30     return 0;
31 }

容斥原理应用

例1.hdu 1796

题意：

　　给一个数n和一个集合m,求1-n中能被m中任意一个数整除的数的个数

思路：

　　最裸的容斥原理,深搜求解

代码：

View Code

 1 #include<stdio.h>
 2 
 3 int m,va[20],n;
 4 __int64 sum;
 5 
 6 __int64 gcd(__int64 x,__int64 y)//最大公约数
 7 {
 8     __int64 temp;
 9     while(x!=0)
10     {
11         temp=x;
12         x=y%x;
13         y=temp;
14     }
15     return y;
16 }
17 
18 void dfs(__int64 ans,int now,int count)
19 {
20     int i;
21     ans=ans*va[now]/gcd(va[now],ans);//ans表示当前各因子合成的因子
22     if(count&1) //如果因子个数是奇数个就加
23         sum+=n/ans;
24     else //因子个数是偶数个就减
25         sum-=n/ans;
26     for(i=now+1;i<=m;i++)//增加一个因子
27     {
28         dfs(ans,i,count+1);
29     }
30 }
31 
32 int check(int x,int y)
33 {
34     if(x>y&&x%y==0)
35         return 1;
36     else
37         return 0;
38 }
39 
40 int main()
41 {
42     int i,j;
43     while(scanf("%I64d%d",&n,&m)!=EOF)
44     {
45         for(i=1;i<=m;i++)
46         {
47             scanf("%d",&va[i]);
48             if(va[i]==0)//排除集合中的0
49             {
50                 i--;m--;
51             }
52             else
53             {
54                 for(j=1;j<i;j++)//如果存在2,那么可以直接排除4
55                 {
56                     if(check(va[i],va[j]))
57                     {
58                         i--;m--;break;
59                     }
60                 }
61             }
62         }
63         sum=0;
64         n--;
65         for(i=1;i<=m;i++)//容斥求解,分别求解能只能被va[i]整除的个数(即排除后面重复的部分)
66         {
67             dfs(va[i],i,1);
68         }
69         printf("%I64d\n",sum);
70     }
71     return 0;
72 }

例2.hdu 2197 本原串

题意：

　　中文题目...不解释.

思路：

　　如果不考虑合法性，则长度为n的串有2^n种,现在要除去不是本原串的部分,对于一个长度为n的串,它的非法串可以由长度为i的本原串组成,其中i能整除n,由于小数据可能被多次用到,所以要先预处理出一部分.

代码：

View Code

 1 #include<stdio.h>
 2 #include<string.h>
 3 int a[8010];
 4 
 5 int quickpow(int m,int n,int k)//快速幂取模
 6 {
 7     int b=1;
 8     while(n)
 9     {
10         if(n&1)
11             b=(b*m)%k;
12         n=n>>1;
13         m=(m*m)%k;
14     }
15     return b;
16 }
17 
18 int get(int x)
19 {
20     int ans,i;
21     if(x<=8000&&a[x])
22         return a[x];
23     ans=quickpow(2,x,2008);
24     for(i=2;i<x;i++)
25     {
26         if(i*i>x)//a*b=x(a<=b)这里i代表的是a
27             break;
28         if((x%i)==0&&(i*i!=x))
29             ans=ans-get(i)-get(x/i);
30         if(i*i==x)//如果是平方数则两个因子相同,减去自己即可
31             ans-=get(i);
32     }
33     return ans-2;
34 }
35 
36 int main()
37 {
38     int n,i;
39     a[1]=2;
40     for(i=1;i<=8000;i++)//预处理出长度8000以下的本原串个数
41     {
42         a[i]=get(i);
43     }
44     while(scanf("%d",&n)!=EOF)
45     {
46         int tt=get(n);
47         while(tt<0)
48             tt+=2008;
49         printf("%d\n",tt%2008);
50     }
51 }

例3.hdu 2204 Eddy's爱好

题意：

　　中文题目...不解释.

思路：

　　由于任何次方数都可以用一个素数来表示如,x^6=(x^3)^2即只要计算过2次的个数,就不用再去考虑6次的个数，所以我们只需要求出1-n中有多少个数可以表示成一个数的素数次但是,像X^6这类元素,在计算2次和3次的时候各算过一次,所以还要减去重复的情况,所以明显应该用容斥来做,不过本题最恶心的地方不是思路,而是精度问题,由于系统函数pow()对大数的精度很差,所以要自己用2分模拟.

代码：

View Code

 1 #include<stdio.h>
 2 long long max=9223372036854775807.0;
 3 int prim[20]={2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61};
 4 
 5 long long _pow(long long n,long long k)
 6 {
 7     if(n==1||k>=64)
 8         return 1;
 9     if(k==1)
10         return n;
11     long long l=0,r=1000000000000000000.0,m;
12     while(l<=r)
13     {
14         m=(l+r)>>1;
15         long long res=1,i;
16         for(i=1;i<=k;i++)
17         {
18             if(m==0||max/m<res)
19                 break;
20             res*=m;
21             if(res>n)
22                 break;
23         }
24         if(i>k)
25         {
26             if(res==n)
27                 return m;
28             long long s=m+1,r=1;
29             for(i=1;i<=k;i++)
30             {
31                 if(s==0||(max/s)<r)
32                     break;
33                 r*=s;
34                 if(r>n)
35                     break;
36             }
37             if(i<=k)
38                 return m;
39         }
40         if(i<=k)
41             r=m;
42         else l=m;
43     }
44 }
45 
46 long long dfs(int start,long long n)
47 {
48     long long ans=0;
49     for(int i=start;i<18;i++)
50     {
51         long long temp=_pow(n,prim[i]);
52         ans+=temp-dfs(i+1,temp);
53     }
54     return ans;
55 }
56 int main()
57 {
58     long long n;
59     while(scanf("%I64d",&n)!=EOF)
60     {
61         printf("%I64d\n",dfs(0,n));
62     }
63     return 0;
64 }

(本题的代码可以优化到0ms,因为本身符合情况的2个素数和3个素数的组合的情况很少,可以全部列处理出来,而4个素数的积已经超过了64,不用考虑）

例4.hdu 1796 How many integers can you find

题意：

　　给定两个数n,m和一个集合,这个集合一共有m个数,求小于n的数中,有多少个能被该集合中任意一个元素整除(只要有一个元素满足就可以)

思路：

　　赤裸裸的容斥,不过有坑爹数据,那就是集合元素有可能是0,还有一点就是,排除会得出重复答案的元素可以增强效率,比如2已经存在了,4就不用加进来了.容斥的时候注意,没有保证每个数是质数,所以不是直接的乘法,而是求最小公倍数.

代码：

View Code

 1 #include<stdio.h>
 2 int va[20],m;
 3 
 4 int gcd(int x,int y)
 5 {
 6     while(x)
 7     {
 8         int temp=x;
 9         x=y%x;
10         y=temp;
11     }
12     return y;
13 }
14 int dfs(int start,int n,int now)
15 {
16     int ans=0,i;
17     for(i=start;i<m;i++)
18     {
19         int temp=now*(va[i]/gcd(now,va[i]));
20         ans+=n/temp-dfs(i+1,n,temp);
21     }
22     return ans;
23 }
24 
25 int check(int x,int y)
26 {
27     if(x<y&&y%x==0)
28         return 1;
29     return 0;
30 }
31 
32 int main()
33 {
34     int n,i,j;
35     while(scanf("%d%d",&n,&m)!=EOF)
36     {
37         for(i=0;i<m;i++)
38         {
39             scanf("%d",&va[i]);
40             if(va[i]==0)//排除是0的元素
41             {
42                 i--;
43                 m--;
44             }
45             else
46             {
47                 for(j=0;j<i;j++)//排除没有必要的元素
48                 {
49                     if(check(va[j],va[i]))
50                     {
51                         i--;
52                         m--;
53                         break;
54                     }
55                 }
56             }
57         }
58         printf("%d\n",dfs(0,n-1,1));
59 
60     }
61     return 0;
62 }

例5.hdu 1695 GCD

题意：

　　给定五个数a,b,c,d,k,其中保证a=1,c=1,求从a-b和c-d两个区间中分别取一个数x,y使得GCD（x,y)=k的x,y组合有多少种.

思路：

　　对于每组b,d,只要从1到b或者从1到d枚举每一种情况方可得出答案,比如枚举x分别为1到b时,分别求出对于每个x,1到d有多少种符合情况的个数,个数用容斥原理来求，显然要先分解出x的素因子,考虑到很多数都会重复分解,故可以在筛素数的时候把1到10W的质因子先预处理出来.但是如果直接算的话又会发现问题x=1,y=3和x=3,y=1实际上是同一种情况,所以,假设b<d(如果不符合把b,d互换),显然1到b符合情况的组合数=euler[1]+euler[2]+...+euler[b],剩下b+1到d的情况,由于和1-b搭配,所以不会出现重复,可以直接容斥.

代码：

View Code

 1 #include<stdio.h>
 2 int prim[100010][16],cnt[100010];
 3 long long euler[100010];
 4 
 5 void init()
 6 {
 7     long long i,j;
 8     euler[1]=1;
 9     for(i=2;i<=100000;i++)
10     {
11         if(!euler[i])
12         {
13             for(j=i;j<=100000;j+=i)
14             {
15                 if(!euler[j])
16                     euler[j]=j;
17                 euler[j]=euler[j]*(i-1)/i;
18                 prim[j][cnt[j]++]=i;//1-10W的数预处理分解质因数
19             }
20         }
21         euler[i]+=euler[i-1];//预处理出每个euler[1]+..+euler[n]的和
22     }
23 }
24 
25 long long dfs(int start,int b,int n)
26 {
27     long long ans;
28     ans=0;
29     for(int i=start;i<cnt[n];i++)
30     {
31         ans+=b/prim[n][i]-dfs(i+1,b/prim[n][i],n);
32     }
33     return ans;
34 }
35 
36 
37 
38 int main()
39 {
40     int t,i,b,d,k,cases=0;
41     long long ans;
42     init();
43     scanf("%d",&t);
44     while(t--)
45     {
46         cases++;
47         scanf("%*d%d%*d%d%d",&b,&d,&k);
48         if(k==0)
49         {
50             printf("Case %d: 0\n",cases);
51             continue;
52         }
53         if(b>d)//保证b是小的那一边
54         {
55             int temp=b;
56             b=d;
57             d=temp;
58         }
59         b/=k;
60         d/=k;
61         ans=euler[b];
62         for(i=b+1;i<=d;i++)
63         {
64             ans+=b-dfs(0,b,i);//容斥
65         }
66         printf("Case %d: %I64d\n",cases,ans);
67     }
68     return 0;
69 }

例6.hdu 2841 Visible Trees

题意：

　　一个人站在(0,0)点观看一个m*n的矩阵上的点,如果同意直线上已经有一个点,则该点会被前面的点挡住,简化一下,就是求满足(x,y)中GCD(x,y)==1的点的个数,和例5相比要简单很多.

思路:

　　基本和例5一样,只要把X轴或Y轴枚举一遍就可以了,为了提高效率,枚举元素较少的一遍,然后对枚举的每一个元素用容斥求个数.

代码:

View Code

 1 #include<iostream>
 2 using namespace std;
 3 int cnt[100010],prim[100010][16];
 4 bool mark[100010];
 5 
 6 void init()
 7 {
 8     int i,j;
 9     for(i=2;i<=100000;i++)
10     {
11         if(!mark[i]){
12             for(j=i;j<=100000;j+=i)
13             {
14                 mark[j]=1;
15                 prim[j][cnt[j]++]=i;
16             }
17         }
18     }
19 }
20 
21 long long dfs(int start,int va,int n)
22 {
23     long long ans=0;
24     for(int i=start;i<cnt[va];i++)
25     {
26         int temp=n/prim[va][i];
27         ans+=temp-dfs(i+1,va,temp);
28     }
29     return ans;
30 }
31 
32 int main()
33 {
34     int t,m,n,i;init();
35     long long ans;
36     scanf("%d",&t);
37     while(t--)
38     {
39         scanf("%d%d",&m,&n);
40         if(m>n)
41             swap(m,n);
42         ans=0;
43         for(i=1;i<=m;i++)
44         {
45             ans+=n-dfs(0,i,n);
46         }
47         printf("%I64d\n",ans);
48     }
49 }

例7.hdu 3208 Integer’s Power

题意：

　　给一个范围a-b,对于这个范围中的每一个元素I,如果它能表示成I=K^n,那么这个数的值为n(当有多个n符合时取最大的),然后求a-b所有元素值的和.

思路:

　　刚看到这个题先感觉和例3有点像,但是这道题要求的比那道题要详细一点,WA了一个下午以后才恍然大悟,对于一个区间1-n能表示成I=K^n的元素的个数用上一题的解法就可以求出,那么假设共有m个数可以表示成n(n!=1)次方的形式,他们的底数K一定是1-m这几个数,对于区间1-m,如果一个元素不能表示成I=K^n的形式,那么这个元素当前就是最优解,如果还可以表示成I=K^n的形式,那么他一定有值更大的n.

代码：

View Code

 1 //为了加快速度,先把1-3个质因数可能的组合全部打出来了
 2 #include<stdio.h>
 3 long long max=9223372036854775807.0;
 4 int prim[20]={2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61};
 5 int prim2[20]={6,10,14,15,21,22,26,33,34,35,38,39,46,51,55,57,58};
 6 int prim3[3]={30,42};
 7 long long _pow(long long n,long long k)//高精度开K次方
 8 {
 9     if(n==1||k>=64) return 1;
10     if(k==1) return n;
11     long long l=0,r=1000000000000000000.0,m;
12     while(l<=r){
13         m=(l+r)>>1;
14         long long res=1,i;
15         for(i=1;i<=k;i++){
16             if(m==0||max/m<res) break;
17             res*=m;
18             if(res>n) break;
19         }
20         if(i>k){
21             if(res==n) return m;
22             long long s=m+1,r=1;
23             for(i=1;i<=k;i++){
24                 if(s==0||(max/s)<r) break;
25                 r*=s;
26                 if(r>n) break;
27             }
28             if(i<=k) return m;
29         }
30         if(i<=k) r=m;
31         else l=m;
32     }
33 }
34 
35 long long getsum(long long n)
36 {
37     long long ans=0,i;
38     for(i=0;i<18;i++)//一个质因子次方
39     {
40         long long temp=_pow(n,prim[i])-1;//1不符合,所以要减掉
41         if(temp==0)
42             break;
43         ans+=temp;
44     }
45     for(i=0;i<17;i++)//两个质因子次方
46     {
47         long long temp=_pow(n,prim2[i])-1;
48         if(temp==0)
49             break;
50         ans-=temp;
51     }
52     for(i=0;i<2;i++)//三个质因子次方
53     {
54         long long temp=_pow(n,prim3[i])-1;
55         if(temp==0)
56             break;
57         ans+=temp;
58     }
59     return ans+1;//因为1在这道题中不算数,不符合的还要加上1
60 }
61 
62 long long getnum(long long x)
63 {
64     long long renum,ans,now=2;
65     renum=ans=x-getsum(x);
66     while(renum!=0)
67     {
68         long long temp=_pow(x,now);
69         renum=temp-getsum(temp);
70         ans+=renum*now;//最优解是now的个数为renum
71         now++;
72     }
73     return ans;
74 }
75     
76 int main()
77 {
78     long long a,b;
79     while(scanf("%I64d%I64d",&a,&b)!=EOF)
80     {
81         if(a==0&&b==0)
82             break;
83         printf("%I64d\n",getnum(b)-getnum(a-1));
84     }
85     return 0;
86 }

例8.poj 1091 跳蚤

题意：

　　中文题目不解释.

思路：

　　刚看到题目并没有什么明确的思路.......后来忍不住就看了一下Discuss( 唉~~顿时觉得自己真是弱爆了,数学知识给谁吃了- - ),对于每张卡片,假设上面的每个数字分别是x1,x2,x3....xn,那么这个跳蚤走的路线可以表示成a1*x1+a2*x2+...+an*xn,(a1,a2..an)为任意整数(可以使负的表示反方向),只要这个值可能等于1,当前的卡片就是合法的,假设gcd(x1,x2...xn)=L,那么刚才路线的公式可以变式为a1*(x1/L)*L+a2*(x2/L)*L+a3*(x3/L)*L+....+an*(xn/L)*L=1,由于a1可以表示任意整数,所以可以直接写成a1*L+a2*L+a3*L+....+an*L=1,--->可以看出L必须是0,才有可能等于1.这样本题就转化成了,求n个小于等于m的数,使得gcd(x1,x2...xn)=1.

代码：

View Code

 1 #include<stdio.h>
 2 bool mark[10010];
 3 int primnum,prim[3000],factor[66];
 4 long long LL;
 5 
 6 int init()    //素数筛
 7 {    
 8     int i,j,kiss=0;
 9     for(i=2;i<=10000;i++)
10     {
11         if(!mark[i])
12         {
13             prim[kiss++]=i;
14             for(j=i*i;j<=10000;j+=i)
15             {
16                 mark[j]=1;
17             }
18         }
19     }
20     return kiss;
21 }
22 
23 long long _pow(long long m,long long n)
24 {
25     long long ans=1;
26     while(n){
27         if(n&1)    ans*=m;
28         n=n>>1;    m*=m;
29     }
30     return ans;
31 }
32 
33 int getfactor(int x)    //分解质因子
34 {
35     int kiss=0,i;
36     for(i=0;i<primnum;i++)
37     {
38         if(x%prim[i]==0)
39         {
40             factor[kiss++]=prim[i];
41             while(x%prim[i]==0) x/=prim[i];
42         }
43     }
44     if(x!=1) factor[kiss++]=x;
45     return kiss;
46 }
47 
48 void dfs(int start,int end,int maxnum,int cnt,int now,int x,int n)
49 {    //开始元素,结束元素,目标质因子个数,当前质因子个数,m的值,n的值
50     if(cnt==maxnum)
51         LL+=_pow(x/now,n);
52     for(int i=start;i<end;i++)
53         dfs(i+1,end,maxnum,cnt+1,now*factor[i],x,n);
54 }
55 
56 long long getans(int n,int m)
57 {
58     int factornum=getfactor(m);//对m分解质因子
59     long long ans=_pow(m,n);//所有可能的情况,包括非法的情况
60     for(int i=1;i<=factornum;i++)
61     {
62         LL=0; dfs(0,factornum,i,0,1,m,n);//分别求出由i个质因子组合作为公因数时可能的情况数,再用容斥原理求出总共的不合法状态
63         if(i&1) ans-=LL;
64         else ans+=LL;
65     }
66     return ans;
67 }
68 
69 
70 int main()
71 {
72     int n,m;
73     primnum=init();
74     while(scanf("%d%d",&n,&m)!=EOF)
75         printf("%lld\n",getans(n,m));
76     return 0;
77 }

PS：这道题本来用long long是要暴掉的,但是看到说数据比较水就懒得写大数了,练练思想就好(∩_∩)~

(4)斯特灵数

第一类斯特灵数

第一类斯特灵数的定义：　　

　　第一类Stirling数是有正负的，其绝对值是n个元素的项目分作k个环排列的方法数目

递推公式：

　　stl[i][j]=stl[i-1][j-1]+(i-1)*stl[i-1][j].注意和第二类斯特灵数区分,前半部分表示i-1个元素中形成j-1个环,剩下的一个元素自成一环,后半部分表示i-1个元素已经组成了j个环,由于排列的位置不同,对应的环也不同(而对于第二类斯特灵数是相同的)所以,在每一个环里有k种方法,其中K=元素数目,因此总共有i-1种放法.

例1.hdu 3625 Examining the Rooms

题意：

　　一个侦探需要打开n间房子,现在有n把钥匙分别放在n个房子里(也就是说一共有n!种放法),当手里没有钥匙的时候需要撞开一扇门,从这个门里拿出钥匙再来开别的门,其中1号门中有重要人物,不能撞开,现在给一个最多撞门的次数,问有多大的概率可以打开所有的门.

分析:

　　从撞开第一个门开始,就相当于进入了一个环,所以可以形成环的数目就相当于需要撞门的次数,但是由于特殊的要求,1不能自成一环,所以每次要减掉1自成一环的情况.

代码：

View Code

 1 #include<stdio.h>
 2 long long stl[30][30],fac[30];
 3 void init()
 4 {
 5     int i,j;
 6     stl[0][0]=0;stl[1][0]=0;stl[1][1]=1;
 7     fac[1]=1;
 8     fac[0]=1;
 9     for(i=2;i<=20;i++)
10     {
11         for(j=1;j<=i;j++)
12         {
13             stl[i][j]=stl[i-1][j-1]+(i-1)*stl[i-1][j];
14         }
15         fac[i]=fac[i-1]*i;
16     }
17 }
18 int main()
19 {
20     int t,n,k,i;
21     long long ans;
22     init();
23     scanf("%d",&t);
24     while(t--)
25     {
26         ans=0;
27         scanf("%d%d",&n,&k);
28         for(i=1;i<=k;i++)
29         {
30             ans+=stl[n][i]-stl[n-1][i-1];
31         }
32         printf("%.4lf\n",ans*1.0/fac[n]);
33     }
34     return 0;
35 }

第二类斯特灵数

第二类斯特灵数的定义：

　　第二类Stirling数是n个元素的集定义k个等价类的方法数目

递推公式：

　　stl[i][j]=stl[i-1][j]*j+stl[i-1][j-1],这个很好想,i-1个元素找出j个等价类,剩下一个元素可能在这j个类中的一个,所以情况数为stl[i-1][j]*j,而从i-1个元素中找出j-1个等价类,剩下的一个元素要自成1类,所以情况数为stl[i-1][j-1].实现的时候可以利用数组的初始值是0不对每个stl[i][0]赋值.

例1.hdu 2643

题意：

　　给你2种符号"<","=",和n个元素,问他们之间可能有几种组合.

思路：

　　把相等的元素看做一组,然后每一组从小到大排列假设有x组,则共有X!种排列方法.

代码：

View Code

 1 #include<stdio.h>
 2 long long stl[102][102],fac[102];
 3 
 4 void init()//预处理斯特灵数和阶乘
 5 {
 6     int i,j;
 7     stl[1][1]=1;stl[1][0]=0;
 8     for(i=2;i<=100;i++)
 9         for(j=1;j<=i;j++)
10             stl[i][j]=(stl[i-1][j]*j+stl[i-1][j-1])%20090126;
11     fac[0]=1;
12     for(i=1;i<=100;i++)
13         fac[i]=(fac[i-1]*i)%20090126;
14 }
15 
16 int main()
17 {
18     long long ans;
19     int t,n,i;
20     init();
21     scanf("%d",&t);
22     while(t--)
23     {
24         ans=0;
25         scanf("%d",&n);
26         for(i=1;i<=n;i++)
27             ans=(ans+stl[n][i]*fac[i])%20090126;
28         printf("%I64d\n",ans);
29     }
30     return 0;
31 }

(5)DP！？

　　有一类题目,一般是让求满足条件的状态数,这时候一般都是从前面满足条件的小状态转移到大一点的状态,本来和组合数学关系并不大,更倾向于DP的范畴,但是鉴于这些题目是做组合数学专题的时候遇到的,又比较有思辨性,就记录下来吧.

例1.hdu 2431 Counting Problem

题意：

　　给一个n*n的棋盘,往上面放2*n个棋子,使得每行每列都有且只有2个棋子,问有多少种放法.

思路：

　　每一个n*n的棋盘都可以分成a*a和b*b的棋盘,其中(a+b)=n且a>1,b>1,所以可以直接从前面的合法情况推后面的.

代码：

View Code

 1 #include<stdio.h>
 2 int ans[510];
 3 void init()
 4 {
 5     int i,j;
 6     ans[0]=1;
 7     ans[1]=0;
 8     for(i=2;i<=500;i++)
 9     {
10         for(j=i;j<=500;j++)
11             ans[j]=(ans[j]+ans[j-i])%1000007;
12     }
13 }
14 int main()
15 {
16     int t,n;
17     init();
18     scanf("%d",&t);
19     while(t--)
20     {
21         scanf("%d",&n);
22         printf("%d\n",ans[n]);
23     }
24     return 0;
25 }

例2.hdu 1246 自共轭Ferrers图
题意:

　　中文题目,不解释..

思路：

　　对于一个大于3的奇数K,可以套在最外层小于K个的合法状态上形成一个新的状态

代码:

View Code

 1 #include<stdio.h>
 2 int ans[310];
 3 void init()
 4 {
 5     int i,j;
 6     ans[0]=1;ans[1]=1;ans[2]=0;
 7     for(i=3;i<=300;i+=2)
 8     {
 9         for(j=300;j>=i;j--)
10         {
11             ans[j]+=ans[j-i];
12         }
13     }
14 }
15 
16 
17 int main()
18 {
19     int n;
20     init();
21     while(scanf("%d",&n)!=EOF)
22     {
23         printf("%d\n",ans[n]);
24     }
25     return 0;
26 }