组合数学总结

组合数学总结

（转发请注明出处）

1.基础知识

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

例1.hdu 1799 循环多少次?

题意：

　　中文题目..不解释.

思路：

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

代码：

#include<stdio.h>
int c[2010][2010];

void init()
{
    int i,j;
    c[0][0]=c[1][0]=c[1][1]=1;
    for(i=2;i<=2000;i++)
    {
        for(j=0;j<=i;j++)
        {
            if(j==0)
                c[i][j]=1;
            else
                c[i][j]=(c[i-1][j]+c[i-1][j-1])%1007;
        }
    }
}
int main()
{
    int t,m,n;
    scanf("%d",&t);
    init();
    while(t--)
    {
        scanf("%d%d",&m,&n);
        printf("%d\n",c[n][m]);
    }
    return 0;
}

  

(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.",知道这个定理后就可以每次对阶乘预处理，快速求出每一位的组合数大小,继而求得相乘后的答案.

代码：

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

#include<stdio.h>//From  www.cnblogs.com/SolarWings/

long long factor[100010];

void init(long long p)//预处理阶乘,化为p进制后,C(a,b)中a,b<p,所以处理到p就足够了
{
    long long i;
    factor[0]=1;
    for(i=1;i<=p;i++)
    {
        factor[i]=(factor[i-1]*i)%p;
    }
}

long long Extended_Euclid(long long a,long long b,long long &x,long long &y)//扩展欧几里得
{
    long long d;
    if(b==0)
    {
        x=1;y=0;return a;
    }
    d=Extended_Euclid(b,a%b,y,x);
    y-=a/b*x;
    return d;
}
long long Inv(long long a,long long n)//求模逆元
{
    long long d,x,y;
    d=Extended_Euclid(a,n,x,y);//需要用到扩展欧几里得,或者也可以用快速幂直接得到模拟元,以后的题目中再做说明
    if(d==1) return (x%n+n)%n;
    else return -1;
}


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代表模逆元
{
    long long ans;
    if(b>a)
        return 0;
    ans=factor[a];
    ans=(ans*Inv(factor[b],p))%p;
    ans=(ans*Inv(factor[a-b],p))%p;
    return ans;
}

long long c(long long n,long long m,long long p) //这里用到了Lucas定理
{
    long long ans,a,b; ans=1;
    while(m) //化为p进制分别求组合数
    {
        a=n%p;
        b=m%p;
        ans=(ans*get(a,b,p))%p;
        m/=p;
        n/=p;
    }
    return ans;
}


int main()
{
    long long t,n,m,p;
    scanf("%I64d",&t);
    while(t--)
    {
        scanf("%I64d %I64d %I64d",&n,&m,&p);
        init(p);//由于每次的p不同,所以要分别预处理
        printf("%I64d\n",c(n+m,n,p));
    }
    return 0;
}

快速幂求模逆元写法：

#include<stdio.h>
long long fac[100010];
void init(long long p)
{
    fac[0]=1;
    for(long long i=1;i<=p;i++)
    {
        fac[i]=((long long)fac[i-1]*i)%p;
    }
    return;
}

long long quickpow(long long m,long long n,long long p)
{
    long long ans=1;
    while(n)
    {
        if(n&1)
            ans=(ans*m)%p;
        n=n>>1;
        m=(m*m)%p;
    }
    return ans;
}

long long Inv(long long x,long long p)
{
    return quickpow(x,p-2,p);
}

long long C(long long n,long long m,long long p)
{
    if(m>n)
        return 0;
    long long ans=1;
    ans=ans*fac[n];
    ans=(ans*Inv(fac[n-m],p))%p;
    ans=(ans*Inv(fac[m],p))%p;
    return ans;
}

long long getans(long long n,long long m,long long p)
{
    long long ans=1;
    init(p);
    while(n&&m)
    {
        long long a=n%p;
        long long b=m%p;
        ans=(ans*(long long)C(a,b,p))%p;
        n/=p;
        m/=p;
    }
    return ans;
}

int main()
{
    long long t,n,m,p;
    scanf("%I64d",&t);
    while(t--)
    {
        scanf("%I64d%I64d%I64d",&n,&m,&p);
        printf("%I64d\n",getans(n+m,m,p));
    }
    return 0;
}

例2.hdu 3944 DP?

题意：

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

思路：

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

代码：

#include<stdio.h>
bool mark[10010];
int num[10010],map[10010];
int prim[2000][10010];//最大的素数要9000+
void init()
{
    int i,k,j;
    k=0;
    for(i=2;i<=10000;i++)//素数筛
    {
        if(!mark[i])
        {
            map[i]=k;
            num[k++]=i;
            for(j=i*i;j<=10000;j+=i)
            {
                mark[j]=1;
            }
        }
    }
    for(i=0;i<k;i++)
    {
        int p=num[i];
        prim[i][0]=1;
        for(j=1;j<=p;j++)
        {
            prim[i][j]=(prim[i][j-1]*j)%p;//预处理出所有素数的阶乘
        }
    }
    return;
}

int Extended_Euclid(int a,int b,int &x,int &y)
{
    int d;
    if(b==0)
    {
        x=1;y=0;return a;
    }
    d=Extended_Euclid(b,a%b,y,x);
    y-=a/b*x;
    return d;
}

int Inv(int a,int n)
{
    int d,x,y;
    d=Extended_Euclid(a,n,x,y);
    if(d==1) return (x%n+n)%n;
    else return -1;
}

int c(int a,int b,int p)
{
    if(b>a)
        return 0;
    int ans=prim[map[p]][a];
    ans=(ans*(Inv(prim[map[p]][b],p)))%p;
    ans=(ans*(Inv(prim[map[p]][a-b],p)))%p;
    return ans;
}

int get(int n,int m,int p)
{
    int a,b,ans;
    ans=1;
    while(m&&n)
    {
        a=n%p;
        b=m%p;
        ans=(ans*c(a,b,p))%p;
        m/=p;
        n/=p;
    }
    return ans;
}



int main()
{
    int n,k,p,cases;
    cases=0;
    init();
    while(scanf("%d%d%d",&n,&k,&p)!=EOF)
    {
        cases++;
        if(k>n/2)//这步转化是把右边的也转化成左边的
            k=n-k;
        printf("Case #%d: %d\n",cases,(get(n+1,k,p)+n-k)%p);
    }
    return 0;
}

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)

代码：

#include<stdio.h>
//本题重点在于明确如何让推出的公式，计算只是简单的大数
int a[201][510],b[510];
void init()//预处理阶乘
{
    int i,j,len,temp1,temp;
    len=0;
    a[0][0]=1;
    a[1][0]=1;
    for(i=2;i<=200;i++)
    {
        for(j=0;j<=len;j++)
        {
            a[i][j]=a[i-1][j]*i;
        }
        temp1=0;
        for(j=0;j<=len;j++)
        {
            temp=a[i][j]+temp1;
            a[i][j]=temp%10;
            temp1=temp/10;
        }
        while(temp1)
        {
            a[i][++len]=temp1%10;
            temp1/=10;
        }
    }
}
int main()
{
    init();
    int n,m,cases,temp1,temp,i,j;
    cases=0;
    while(scanf("%d%d",&m,&n)!=EOF)
    {
        cases++;
        if(m==0&&n==0)
            break;
        printf("Test #%d:\n",cases);
        if(m<n)
        {
            printf("0\n");
                continue;
        }
        temp1=0;
        for(i=0;i<=500;i++)//大数乘法计算b=(M+N)!*(M-N+1)
        {
            temp=a[m+n][i]*(m+1-n)+temp1;
            b[i]=temp%10;
            temp1=temp/10;
        }
        while(!b[i])
        {
            i--;
        }
        temp1=0;
        for(j=i;j>=0;j--)//大数除法计算b=b/(M+1)
        {
            temp=temp1*10+b[j];
            b[j]=temp/(m+1);
            temp1=temp%(m+1);
        }
        while(!b[i])
            i--;
        for(j=i;j>=0;j--)
        {
            printf("%d",b[j]);
        }
        printf("\n");
    }
    return 0;
}

例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!

代码：

#include<stdio.h>//本题的主要思想在于分解质因子
#include<string.h>
int prim[150000],cnt[150000],num;
bool mark[2000010];

int init()//初始化，筛选素数
{
    long long i,j,kiss;
    kiss=0;
    for(i=2;i<=2000000;i++)
    {
        if(!mark[i])
        {
            prim[kiss++]=i;
            for(j=i*i;j<=2000000;j+=i)
            {
                mark[j]=1;
            }
        }
    }
    return kiss;
}

int getfactor(int x,int flag)//对阶乘分解质因子，表示成2^，3^，5^...
{
    int i;
    for(i=0;i<num&&prim[i]<=x;i++)
    {
        int n=x;
        while(n)
        {
            cnt[i]+=flag*(n/prim[i]);
            n/=prim[i];
        }
    }
    return i;
}

long long quickpow(long long m,long long n,long long k)//快速幂求解
{
    long long b=1;
    while(n)
    {
        if(n&1)
            b=(b*m)%k;
        n=n>>1;
        m=(m*m)%k;
    }
    return b;
}

long long solve(int n,int m)
{
    long long ans=1;
    int maxlen=getfactor(n+m,1);//后面的参数为1，表示乘法，系数为曾加
    getfactor(n+1,-1);//后面参数为2，表示除法，系数为减少
    getfactor(m,-1);
    for(int i=0;i<maxlen;i++)
    {
        ans=(ans*(quickpow(prim[i],cnt[i],20100501)))%20100501;
    }
    return ans;
}




int main()
{
    int t,m,n,i;
    num=init();
    scanf("%d",&t);
    while(t--)
    {
        scanf("%d%d",&n,&m);
        memset(cnt,0,sizeof(cnt));
        int temp=n+1-m;
        for(i=0;i<num&&prim[i]<=temp;i++)//先单独对不用阶乘的(n+1-m)处理
        {
            while(temp%prim[i]==0)
            {
                cnt[i]++;
                temp/=prim[i];
            }
        }
        printf("%lld\n",solve(n,m));
    }
    return 0;
}

(3)容斥原理

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

例1.hdu 2824

题意：

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

思路：

　　预处理求解

代码：

#include<stdio.h>
long long euler[3000010];
void init()
{
    int i,j;
    euler[1]=1;
    for(i=2;i<=3000000;i++)
    {
        if(!euler[i])
        {
            for(j=i;j<=3000000;j+=i)
            {
                if(!euler[j])
                    euler[j]=j;
                euler[j]=euler[j]*(i-1)/i;
            }
        }
        euler[i]+=euler[i-1];
    }
}

int main()
{
    int a,b;
    init();
    while(scanf("%d%d",&a,&b)!=EOF)
    {
        printf("%I64d\n",euler[b]-euler[a-1]);
    }
    return 0;
}

容斥原理应用

例1.hdu 1796

题意：

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

思路：

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

代码：

#include<stdio.h>

int m,va[20],n;
__int64 sum;

__int64 gcd(__int64 x,__int64 y)//最大公约数
{
    __int64 temp;
    while(x!=0)
    {
        temp=x;
        x=y%x;
        y=temp;
    }
    return y;
}

void dfs(__int64 ans,int now,int count)
{
    int i;
    ans=ans*va[now]/gcd(va[now],ans);//ans表示当前各因子合成的因子
    if(count&1) //如果因子个数是奇数个就加
        sum+=n/ans;
    else //因子个数是偶数个就减
        sum-=n/ans;
    for(i=now+1;i<=m;i++)//增加一个因子
    {
        dfs(ans,i,count+1);
    }
}

int check(int x,int y)
{
    if(x>y&&x%y==0)
        return 1;
    else
        return 0;
}

int main()
{
    int i,j;
    while(scanf("%I64d%d",&n,&m)!=EOF)
    {
        for(i=1;i<=m;i++)
        {
            scanf("%d",&va[i]);
            if(va[i]==0)//排除集合中的0
            {
                i--;m--;
            }
            else
            {
                for(j=1;j<i;j++)//如果存在2,那么可以直接排除4
                {
                    if(check(va[i],va[j]))
                    {
                        i--;m--;break;
                    }
                }
            }
        }
        sum=0;
        n--;
        for(i=1;i<=m;i++)//容斥求解,分别求解能只能被va[i]整除的个数(即排除后面重复的部分)
        {
            dfs(va[i],i,1);
        }
        printf("%I64d\n",sum);
    }
    return 0;
}

 例2.hdu 2197 本原串

题意：

　　中文题目...不解释.

思路：

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

代码：

#include<stdio.h>
#include<string.h>
int a[8010];

int quickpow(int m,int n,int k)//快速幂取模
{
    int b=1;
    while(n)
    {
        if(n&1)
            b=(b*m)%k;
        n=n>>1;
        m=(m*m)%k;
    }
    return b;
}

int get(int x)
{
    int ans,i;
    if(x<=8000&&a[x])
        return a[x];
    ans=quickpow(2,x,2008);
    for(i=2;i<x;i++)
    {
        if(i*i>x)//a*b=x(a<=b)这里i代表的是a
            break;
        if((x%i)==0&&(i*i!=x))
            ans=ans-get(i)-get(x/i);
        if(i*i==x)//如果是平方数则两个因子相同,减去自己即可
            ans-=get(i);
    }
    return ans-2;
}

int main()
{
    int n,i;
    a[1]=2;
    for(i=1;i<=8000;i++)//预处理出长度8000以下的本原串个数
    {
        a[i]=get(i);
    }
    while(scanf("%d",&n)!=EOF)
    {
        int tt=get(n);
        while(tt<0)
            tt+=2008;
        printf("%d\n",tt%2008);
    }
}

例3.hdu 2204 Eddy's爱好

题意：

　　中文题目...不解释.

思路：

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

代码：

#include<stdio.h>
long long max=9223372036854775807.0;
int prim[20]={2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61};

long long _pow(long long n,long long k)
{
    if(n==1||k>=64)
        return 1;
    if(k==1)
        return n;
    long long l=0,r=1000000000000000000.0,m;
    while(l<=r)
    {
        m=(l+r)>>1;
        long long res=1,i;
        for(i=1;i<=k;i++)
        {
            if(m==0||max/m<res)
                break;
            res*=m;
            if(res>n)
                break;
        }
        if(i>k)
        {
            if(res==n)
                return m;
            long long s=m+1,r=1;
            for(i=1;i<=k;i++)
            {
                if(s==0||(max/s)<r)
                    break;
                r*=s;
                if(r>n)
                    break;
            }
            if(i<=k)
                return m;
        }
        if(i<=k)
            r=m;
        else l=m;
    }
}

long long dfs(int start,long long n)
{
    long long ans=0;
    for(int i=start;i<18;i++)
    {
        long long temp=_pow(n,prim[i]);
        ans+=temp-dfs(i+1,temp);
    }
    return ans;
}
int main()
{
    long long n;
    while(scanf("%I64d",&n)!=EOF)
    {
        printf("%I64d\n",dfs(0,n));
    }
    return 0;
}

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

例4.hdu 1796 How many integers can you find

题意：

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

思路：

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

代码：

#include<stdio.h>
int va[20],m;

int gcd(int x,int y)
{
    while(x)
    {
        int temp=x;
        x=y%x;
        y=temp;
    }
    return y;
}
int dfs(int start,int n,int now)
{
    int ans=0,i;
    for(i=start;i<m;i++)
    {
        int temp=now*(va[i]/gcd(now,va[i]));
        ans+=n/temp-dfs(i+1,n,temp);
    }
    return ans;
}

int check(int x,int y)
{
    if(x<y&&y%x==0)
        return 1;
    return 0;
}

int main()
{
    int n,i,j;
    while(scanf("%d%d",&n,&m)!=EOF)
    {
        for(i=0;i<m;i++)
        {
            scanf("%d",&va[i]);
            if(va[i]==0)//排除是0的元素
            {
                i--;
                m--;
            }
            else
            {
                for(j=0;j<i;j++)//排除没有必要的元素
                {
                    if(check(va[j],va[i]))
                    {
                        i--;
                        m--;
                        break;
                    }
                }
            }
        }
        printf("%d\n",dfs(0,n-1,1));

    }
    return 0;
}

例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搭配,所以不会出现重复,可以直接容斥.

代码：

#include<stdio.h>
int prim[100010][16],cnt[100010];
long long euler[100010];

void init()
{
    long long i,j;
    euler[1]=1;
    for(i=2;i<=100000;i++)
    {
        if(!euler[i])
        {
            for(j=i;j<=100000;j+=i)
            {
                if(!euler[j])
                    euler[j]=j;
                euler[j]=euler[j]*(i-1)/i;
                prim[j][cnt[j]++]=i;//1-10W的数预处理分解质因数
            }
        }
        euler[i]+=euler[i-1];//预处理出每个euler[1]+..+euler[n]的和
    }
}

long long dfs(int start,int b,int n)
{
    long long ans;
    ans=0;
    for(int i=start;i<cnt[n];i++)
    {
        ans+=b/prim[n][i]-dfs(i+1,b/prim[n][i],n);
    }
    return ans;
}



int main()
{
    int t,i,b,d,k,cases=0;
    long long ans;
    init();
    scanf("%d",&t);
    while(t--)
    {
        cases++;
        scanf("%*d%d%*d%d%d",&b,&d,&k);
        if(k==0)
        {
            printf("Case %d: 0\n",cases);
            continue;
        }
        if(b>d)//保证b是小的那一边
        {
            int temp=b;
            b=d;
            d=temp;
        }
        b/=k;
        d/=k;
        ans=euler[b];
        for(i=b+1;i<=d;i++)
        {
            ans+=b-dfs(0,b,i);//容斥
        }
        printf("Case %d: %I64d\n",cases,ans);
    }
    return 0;
}

例6.hdu 2841 Visible Trees

题意：

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

思路:

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

代码:

#include<iostream>
using namespace std;
int cnt[100010],prim[100010][16];
bool mark[100010];

void init()
{
    int i,j;
    for(i=2;i<=100000;i++)
    {
        if(!mark[i]){
            for(j=i;j<=100000;j+=i)
            {
                mark[j]=1;
                prim[j][cnt[j]++]=i;
            }
        }
    }
}

long long dfs(int start,int va,int n)
{
    long long ans=0;
    for(int i=start;i<cnt[va];i++)
    {
        int temp=n/prim[va][i];
        ans+=temp-dfs(i+1,va,temp);
    }
    return ans;
}

int main()
{
    int t,m,n,i;init();
    long long ans;
    scanf("%d",&t);
    while(t--)
    {
        scanf("%d%d",&m,&n);
        if(m>n)
            swap(m,n);
        ans=0;
        for(i=1;i<=m;i++)
        {
            ans+=n-dfs(0,i,n);
        }
        printf("%I64d\n",ans);
    }
}

例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.

代码：

//为了加快速度,先把1-3个质因数可能的组合全部打出来了
#include<stdio.h>
long long max=9223372036854775807.0;
int prim[20]={2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61};
int prim2[20]={6,10,14,15,21,22,26,33,34,35,38,39,46,51,55,57,58};
int prim3[3]={30,42};
long long _pow(long long n,long long k)//高精度开K次方
{
    if(n==1||k>=64) return 1;
    if(k==1) return n;
    long long l=0,r=1000000000000000000.0,m;
    while(l<=r){
        m=(l+r)>>1;
        long long res=1,i;
        for(i=1;i<=k;i++){
            if(m==0||max/m<res) break;
            res*=m;
            if(res>n) break;
        }
        if(i>k){
            if(res==n) return m;
            long long s=m+1,r=1;
            for(i=1;i<=k;i++){
                if(s==0||(max/s)<r) break;
                r*=s;
                if(r>n) break;
            }
            if(i<=k) return m;
        }
        if(i<=k) r=m;
        else l=m;
    }
}

long long getsum(long long n)
{
    long long ans=0,i;
    for(i=0;i<18;i++)//一个质因子次方
    {
        long long temp=_pow(n,prim[i])-1;//1不符合,所以要减掉
        if(temp==0)
            break;
        ans+=temp;
    }
    for(i=0;i<17;i++)//两个质因子次方
    {
        long long temp=_pow(n,prim2[i])-1;
        if(temp==0)
            break;
        ans-=temp;
    }
    for(i=0;i<2;i++)//三个质因子次方
    {
        long long temp=_pow(n,prim3[i])-1;
        if(temp==0)
            break;
        ans+=temp;
    }
    return ans+1;//因为1在这道题中不算数,不符合的还要加上1
}

long long getnum(long long x)
{
    long long renum,ans,now=2;
    renum=ans=x-getsum(x);
    while(renum!=0)
    {
        long long temp=_pow(x,now);
        renum=temp-getsum(temp);
        ans+=renum*now;//最优解是now的个数为renum
        now++;
    }
    return ans;
}
    
int main()
{
    long long a,b;
    while(scanf("%I64d%I64d",&a,&b)!=EOF)
    {
        if(a==0&&b==0)
            break;
        printf("%I64d\n",getnum(b)-getnum(a-1));
    }
    return 0;
}

例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.

代码：

#include<stdio.h>
bool mark[10010];
int primnum,prim[3000],factor[66];
long long LL;

int init()    //素数筛
{    
    int i,j,kiss=0;
    for(i=2;i<=10000;i++)
    {
        if(!mark[i])
        {
            prim[kiss++]=i;
            for(j=i*i;j<=10000;j+=i)
            {
                mark[j]=1;
            }
        }
    }
    return kiss;
}

long long _pow(long long m,long long n)
{
    long long ans=1;
    while(n){
        if(n&1)    ans*=m;
        n=n>>1;    m*=m;
    }
    return ans;
}

int getfactor(int x)    //分解质因子
{
    int kiss=0,i;
    for(i=0;i<primnum;i++)
    {
        if(x%prim[i]==0)
        {
            factor[kiss++]=prim[i];
            while(x%prim[i]==0) x/=prim[i];
        }
    }
    if(x!=1) factor[kiss++]=x;
    return kiss;
}

void dfs(int start,int end,int maxnum,int cnt,int now,int x,int n)
{    //开始元素,结束元素,目标质因子个数,当前质因子个数,m的值,n的值
    if(cnt==maxnum)
        LL+=_pow(x/now,n);
    for(int i=start;i<end;i++)
        dfs(i+1,end,maxnum,cnt+1,now*factor[i],x,n);
}

long long getans(int n,int m)
{
    int factornum=getfactor(m);//对m分解质因子
    long long ans=_pow(m,n);//所有可能的情况,包括非法的情况
    for(int i=1;i<=factornum;i++)
    {
        LL=0; dfs(0,factornum,i,0,1,m,n);//分别求出由i个质因子组合作为公因数时可能的情况数,再用容斥原理求出总共的不合法状态
        if(i&1) ans-=LL;
        else ans+=LL;
    }
    return ans;
}


int main()
{
    int n,m;
    primnum=init();
    while(scanf("%d%d",&n,&m)!=EOF)
        printf("%lld\n",getans(n,m));
    return 0;
}

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自成一环的情况.

代码：

#include<stdio.h>
long long stl[30][30],fac[30];
void init()
{
    int i,j;
    stl[0][0]=0;stl[1][0]=0;stl[1][1]=1;
    fac[1]=1;
    fac[0]=1;
    for(i=2;i<=20;i++)
    {
        for(j=1;j<=i;j++)
        {
            stl[i][j]=stl[i-1][j-1]+(i-1)*stl[i-1][j];
        }
        fac[i]=fac[i-1]*i;
    }
}
int main()
{
    int t,n,k,i;
    long long ans;
    init();
    scanf("%d",&t);
    while(t--)
    {
        ans=0;
        scanf("%d%d",&n,&k);
        for(i=1;i<=k;i++)
        {
            ans+=stl[n][i]-stl[n-1][i-1];
        }
        printf("%.4lf\n",ans*1.0/fac[n]);
    }
    return 0;
}

第二类斯特灵数

第二类斯特灵数的定义：

　　第二类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!种排列方法.

代码：

#include<stdio.h>
long long stl[102][102],fac[102];

void init()//预处理斯特灵数和阶乘
{
    int i,j;
    stl[1][1]=1;stl[1][0]=0;
    for(i=2;i<=100;i++)
        for(j=1;j<=i;j++)
            stl[i][j]=(stl[i-1][j]*j+stl[i-1][j-1])%20090126;
    fac[0]=1;
    for(i=1;i<=100;i++)
        fac[i]=(fac[i-1]*i)%20090126;
}

int main()
{
    long long ans;
    int t,n,i;
    init();
    scanf("%d",&t);
    while(t--)
    {
        ans=0;
        scanf("%d",&n);
        for(i=1;i<=n;i++)
            ans=(ans+stl[n][i]*fac[i])%20090126;
        printf("%I64d\n",ans);
    }
    return 0;
}

(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,所以可以直接从前面的合法情况推后面的.

代码：

#include<stdio.h>
int ans[510];
void init()
{
    int i,j;
    ans[0]=1;
    ans[1]=0;
    for(i=2;i<=500;i++)
    {
        for(j=i;j<=500;j++)
            ans[j]=(ans[j]+ans[j-i])%1000007;
    }
}
int main()
{
    int t,n;
    init();
    scanf("%d",&t);
    while(t--)
    {
        scanf("%d",&n);
        printf("%d\n",ans[n]);
    }
    return 0;
}

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

　　中文题目,不解释..

思路：

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

代码:

#include<stdio.h>
int ans[310];
void init()
{
    int i,j;
    ans[0]=1;ans[1]=1;ans[2]=0;
    for(i=3;i<=300;i+=2)
    {
        for(j=300;j>=i;j--)
        {
            ans[j]+=ans[j-i];
        }
    }
}


int main()
{
    int n;
    init();
    while(scanf("%d",&n)!=EOF)
    {
        printf("%d\n",ans[n]);
    }
    return 0;
}