模拟手算,计算组合数C(m,n)

在笔算组合数的时候,有个很重要的步骤是将分子分母约分。而这讨论的也就是这个约分的过程。

C(m,n) = m*(m-1)*(m-2)*…*(m-n+1) / n!

下文中的分子的因子为m-n+1~m的所有整数。

首先是要储存分子的因子,在这没考虑效率的问题,直接用vector<unsigned>做容器,将分子储存在容器变量arry中,储存顺序为由大到小。

void aryMake(unsigned m,unsigned n,vector<unsigned> &arry)
//使arry成为一个以m为最大值,有n项的连续自然数数组
{
    for(size_type i=0;i!=n;++i){
        arry.push_back(m-i);
    }
}


然后进行第一轮约分操作:从2至n中,观察每一个整数是否能整除分子的某个因子。能,则约去;不能,则储存起来。
for(size_type i=n;i!=1;--i){   
    a=m/i*i;                                  //a为不大于m中i的倍数
    while(arry[m-a]%i){
        a-=i; 
        if((m-a)>(n-1)){
            arySmall.push_back(i);            //不能完全整除的,先记录下来,后面处理
            goto Break;                       //结束内层循环及结束外层循环的当前循环
        }
    
    }
    arry[m-a]/=i;
Break:
    ;
}

其中,a为不大于m中i的倍数。在这,主动寻找符合条件的因子。如在计算C(100,10)中,了解有无3的倍数时,先检验第2项,即初始值为99的项,若该值因之前的运算而不能整除3(99因能被9整除,而约成11),则检验第5(2+3)项,以此类推。若寻找不到符合条件的因子,则储存在arrySmall中。


接下来,则是用找最大公约数的方法,进行再次约分。

unsigned gcdx;
for(size_type i=0,j=arySmall.size();i!=j;++i){
    for(size_type ix=arry.size()-1;ix!=-1;--ix){
        gcdx=gcd(arry[ix],arySmall[i]);
        arry[ix]/=gcdx;
        arySmall[i]/=gcdx;
        if(arySmall[i]==1){
            break;
        }
        
    }
}

其中,gcd为寻找两形参最大公约数的函数。


最后,对arry中的元素逐个相乘。

在这,我给出一个完整代码。其中的大数乘法也仅仅是能运行而已,毫无效率可言(知识有限,现在连类都还没有学习,算法也是一窍不通。等学习了傅立叶转换,会重新编写该部分)。

#include<iostream>
#include<vector>
#include<bitset>
using namespace std;
 
typedef vector<unsigned>::size_type size_type;
 
 
void aryMake(unsigned ,unsigned ,vector<unsigned> &);
void aryReduce(unsigned ,unsigned ,vector<unsigned> &);
void aryReduce(vector<unsigned> &,vector<unsigned> &);
void aryReduce(vector<unsigned> &);
void intCalc(vector<unsigned> &);
void tempMake(unsigned ,vector<vector<bitset<8>>> &,vector<bitset<8>> &);
inline void carry(size_type ,vector<vector<bitset<8>>> & );
inline void carry(size_type i,vector<bitset<8>> &);
void bitCalc(vector<vector<bitset<8>>> & ,vector<bitset<8>> &);
inline unsigned gcd(unsigned ,unsigned );
 
int main()
{
    vector<unsigned> arry;
    unsigned m,n;
    cout<<"Calculate C(m,n)"<<endl;
    cout<<"m=";
    cin>>m;
    cin.get();                    //刷新输入缓冲流
    cout<<"n=";
    cin>>n;
 
    if(m<0||n<0||m<n){
        cout<<"Input Error!"<<endl;
        return -1;
    }
    cout<<"C("<<m<<","<<n<<") = "<<endl;
    if(n==0||m==n){
        cout<<"1"<<endl;
        return 0;
    }
 
    n=2*n<m?n:m-n;                 // C(n,m)==C(m-n,m)
    aryMake(m,n,arry);
    aryReduce(m,n,arry);
    intCalc(arry);
 
    return 0;
 
}
 
inline unsigned gcd(unsigned m,unsigned n)
//求a,b的最大公约数
{
    while(n){
        unsigned temp = m%n;
        m=n;
        n=temp;
    }
    return m;
}
 
void aryMake(unsigned m,unsigned n,vector<unsigned> &arry)
//使arry成为一个以m为最大值,有n项的连续自然数数组
{
    for(size_type i=0;i!=n;++i){
        arry.push_back(m-i);
    }
}
 
void aryReduce(unsigned m,unsigned n,vector<unsigned> &arry)
//对arry/n!进行第一次约分操作
{
    unsigned a;
    vector<unsigned> arySmall;
 
    //寻找n!中的arry的约数,并约分
    for(size_type i=n;i!=1;--i){   //n仍为arry.size
        a=m/i*i;                                  //a为不大于m中i的倍数
        while(arry[m-a]%i){
            a-=i; 
            if((m-a)>(n-1)){
                arySmall.push_back(i);            //不能完全整除的,先记录下来,后面处理
                goto Break;                       //结束内层循环及结束外层循环的当前循环
            }
        
        }
        arry[m-a]/=i;
    Break:
        ;
    }
 
    aryReduce(arry);                             //删除值为1的元素
    
 
    if(!arySmall.empty()){                                //当n!中有未被除尽的数时,进行二次约分
        aryReduce(arySmall,arry);
    }
 
}
 
void aryReduce(vector<unsigned> &arySmall,vector<unsigned> &arry)
//对arry/n!的第一次约分结果进行二次约分
{
    unsigned gcdx;
    for(size_type i=0,j=arySmall.size();i!=j;++i){
        for(size_type ix=arry.size()-1;ix!=-1;--ix){
            gcdx=gcd(arry[ix],arySmall[i]);
            arry[ix]/=gcdx;
            arySmall[i]/=gcdx;
            if(arySmall[i]==1){
                break;
            }
            
        }
    }
    aryReduce(arry);
}
void aryReduce(vector<unsigned> &arry)
//删除arry中值为1的元素
{
    for(vector<unsigned>::iterator iter=arry.begin();iter!=arry.end();){
        if(*iter==1){
            iter=arry.erase(iter);
        }
        else{               //容易出错的地方
            ++iter;
        }
    }
}
void intCalc(vector<unsigned> &arry)
{
    vector<bitset<8>> result;            //8字节和6字节有没不同?
    vector<vector<bitset<8>>> temp;
    result.push_back(1);
    
 
    for(size_type i=0,size=arry.size();i!=size;i++){
        tempMake(arry[i],temp,result);
        bitCalc(temp,result);
        temp.clear();
    }
    
    for(int i=result.size()-1;i!=-1;i--){
        cout<<result[i].to_ulong();
    }
    cout<<endl;
    if(result.size() == 1){
        cout<<"≈"<<result[result.size()-1].to_ulong() <<"e+"<<result.size()-1<<endl;
    }
    else{
        cout<<"≈"<<(result[result.size()-1].to_ulong()*10 
            + result[result.size()-2].to_ulong() + 5) / 10 
            <<"e+"<<result.size()-1<<endl;                      //四舍五入
    }
}
void tempMake(unsigned v,vector<vector<bitset<8>>> &temp,vector<bitset<8>> &result)
//让temp存放计算的中间变量
{
    
    for(size_type i=0,j=1;v/j!=0;j*=10){
        temp.push_back(vector<bitset<8>>(0));          //向temp新建元素     
        int x=v/j%10;                                  //计算每一位数字
            
        //补0
        for(int n=0;n!=i;n++){
            temp[i].push_back(0);
        }
 
        //计算有效数字
        for(size_type k=0,size=result.size();k!=size;k++){
            temp[i].push_back ( result[k].to_ulong() * x );
        }
        temp[i].push_back(0);
        carry(i,temp);
        i++;
    }
}
inline void carry(size_type i,vector<vector<bitset<8>>> & temp)
//进位
{
    for(size_type k=0,size=temp[i].size();k!=size;k++){
        while(temp[i][k].to_ulong()>9){
            temp[i][k+1]=temp[i][k+1].to_ulong()+1;
            temp[i][k]=temp[i][k].to_ulong()-10;
        }
    }
}
 
 
 
void bitCalc(vector<vector<bitset<8>>> & temp,vector<bitset<8>> & result)
{
    vector<bitset<8>> addTemp( temp[temp.size()-1].size()+1,0);    //创建数组用于储存
    size_type sizeAdd=temp[0].size();                 //每次只加sizeAdd个,剩下的已为0
    size_type sizeTmp=temp.size();
    size_type sizeAdT=addTemp.size();
    
    //相加
    size_type i=0;
    if(sizeTmp>=sizeAdd){
        //先计算前sizeAdd列
        for(;i!=sizeAdd;i++){
            for(size_type j=0;j!=sizeAdd;j++){
                addTemp[i]=addTemp[i].to_ulong() + temp[j][i].to_ulong();
                carry(i,addTemp);
            }
        }
        //此时i=sizeAdd
        //计算除后sizeAdd列的部分
        for(;i!=sizeAdT-sizeAdd;i++){
            for(size_type j=i-sizeAdd+1,n=0;n!=sizeAdd;n++){
                addTemp[i]=addTemp[i].to_ulong() + temp[j][i].to_ulong();
                carry(i,addTemp);
                j++;
            }
        }
        //此时i=sizeAdT-sizeAdd=sizeTmp
        //计算剩下的sizeAdd列
        for(;i!=sizeAdT-1;i++){
            for(size_type j=i-sizeAdd+1;j!=sizeTmp;j++){
                    addTemp[i]=addTemp[i].to_ulong() + temp[j][i].to_ulong();
                carry(i,addTemp);
            }
        }
    }
    else{
        //先计算前sizeAdd列
        for(;i!=sizeAdd;i++){
            for(size_type j=0;j!=sizeTmp;j++){
                addTemp[i]=addTemp[i].to_ulong() + temp[j][i].to_ulong();
                carry(i,addTemp);
            }
        }
        //此时i=sizeAdd
        if(sizeTmp!=1){
            //计算剩下部分
            for(;i!=sizeAdT-1;i++){
                for(size_type j=i-sizeAdd+1;j!=sizeTmp;j++){
                    addTemp[i]=addTemp[i].to_ulong() + temp[j][i].to_ulong();
                    carry(i,addTemp);
                }
            }
        }
    }
 
    result.clear();
    //去头位0
    while(addTemp[addTemp.size()-1].none()){
        addTemp.pop_back();
    }
    result=addTemp;
}
 
 
inline void carry(size_type i,vector<bitset<8>> & addTemp)
//进位
{
    while(addTemp[i].to_ulong()>9){
        addTemp[i+1]=addTemp[i+1].to_ulong()+1;
        addTemp[i]=addTemp[i].to_ulong()-10;
    }
}
 
 
 
原文地址:https://www.cnblogs.com/h46incon/p/1948098.html