本文深入讲解组合数学的基础概念与核心技巧，包括小数据范围的组合数预处理、大组合数取模、容斥原理的应用及斯特灵数的计算。通过多个实例解析，如求解组合数、计算本原字符串数量、解决复杂计数问题等，帮助读者掌握组合数学在算法竞赛中的实际应用。
组合数学总结
（转发请注明出处）
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组成的序列得到)
所以我们计算公式就是：合法的排列方式=所有排列方式-非法排列方式
我们这里非法排列方式的计算 就是：(
- 
)*M!*N!，然而在这题，因为每个人都是不同的，所以还要乘以 M！*N！
所以得出最终方程：
F（N）=(-)*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 }
扫一扫
416
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
