模拟手算，计算组合数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;
    }
}