快速求解组合数（筛法+快速幂）

a quick way to calculate combinatorial number

1.Why we need this method?

Anyone who know how to print helloworld is also able to write a code to calculate a combinatorial number. However, the calculation way may be different. You can calculate this number directly, just like we used to do manually. But that’s not enough. When you need to calculate a number like C(20,13). It will probably make some error though you don’t want ,e.g. integer precision overflow. So let’s look at a new way to deal with this problem.

2.quick cloth algorithm

from some math knowledge, we can know that any number can be express as k1^a1*k2^a2*...kn^an. ki(i equals to 1,2,3…) is prime.So the first thing to do is calculate the prime we need(we talked about it yesterday).like this

bool isNotPrime[N];
vector<int>isPrime;
void make_prime(){
    isNotPrime[0]=isNotPrime[1]=1;
    for(int i=2;i<N;i++){
        if(!isNotPrime[i]){
            isPrime.push_back(i);
        }
        for(int j=0;j<isPrime.size()&&isPrime[j]*i<=N-1;j++){
            isNotPrime[i*isPrime[j]]=1;
        }
    }
}

the second step is we should get n! can be expressed as what.We use a for loop to try different prime and then get x(prime^x) in n!,m!,(n-m)!.We can use steps function here:

int steps(int x,int p){
    int ans=0,rec=p;
    while(x>=rec){
        ans+=x/rec;
        rec*=p;
    }
    return ans;
}

Then, the last step we need to do is multiply these pow numbers.This is the entire algorithm:

#include <iostream>
#include <vector>
#include <cmath>
#define N 50000
using namespace std;
bool isNotPrime[N];
vector<int>isPrime;
void make_prime(){
    isNotPrime[0]=isNotPrime[1]=1;
    for(int i=2;i<N;i++){
        if(!isNotPrime[i]){
            isPrime.push_back(i);
        }
        for(int j=0;j<isPrime.size()&&isPrime[j]*i<=N-1;j++){
            isNotPrime[i*isPrime[j]]=1;
        }
    }
}
int steps(int x,int p){
    int ans=0,rec=p;
    while(x>=rec){
        ans+=x/rec;
        rec*=p;
    }
    return ans;
}
long long C(int n,int m){
    long long ans=1;
    for(int i=0;i<isPrime.size()&&isPrime[i]<=n;i++){
        int num=steps(n,isPrime[i])-steps(m,isPrime[i])-steps(n-m,isPrime[i]);
        ans*=pow((double)isPrime[i],(double)num);
    }
    return ans;
}
int main(void){
    make_prime();
    int n,m;
    cin>>n>>m;
    cout<<C(n,m)<<endl;
}