NOIP中的数学 -- 第7课组合数

信奥教练Andy

已于 2022-08-22 16:10:56 修改

阅读量491

点赞数

分类专栏： NOIP中的数学文章标签：算法线性代数

于 2021-02-27 12:28:10 首次发布

本文链接：https://blog.csdn.net/qq_32431299/article/details/114168434

版权

NOIP中的数学专栏收录该内容

13 篇文章 19 订阅

订阅专栏

一、当n,m都很小的时候可以利用杨辉三角直接求。

C(n,m)=C(n-1,m)+C(n-1,m-1)；

注释：C（ n,m) 表示从n个苹果中选择m个苹果，假设有一个苹果分开，这个苹果就有选和不选两种方案，分别那么从剩下的n-1个苹果中选择m-1和m个，即 C( n-1, m-1) 和C( n-1, m)

给定 n 组询问，每组询问给定两个整数 a，b，请你输出 C(a,b)mod(10⁹+7) 的值。

输入格式第一行包含整数 n。

接下来 n 行，每行包含一组 a 和 b。

输出格式共 n 行，每行输出一个询问的解。

数据范围 1≤n≤10000, 1≤b≤a≤2000
输入样例：
3
3 1
5 3
2 2
输出样例：
3
10
1

#include<iostream>
using namespace std;
const int N= 2010,mod=1e9+7;
int f[N][N];
int n;
void init()
{
    for(int i = 0;i<N;i++)
        for(int j=0;j<=i;j++)
        {
            if(!j) f[i][j]=1;
            else f[i][j]=(f[i-1][j]+f[i-1][j-1])%mod;
        }
}
int main()
{
    init();
    cin >> n;
    int a,b;
    while(n--)
    {
        cin >> a>> b;
        cout<<f[a][b]<<endl;
        
    }
}

二、n和m较大，但是p为素数的时候

   还是上面那个题目，只不过数据规模稍有一些变化

1<=n<=10000;
1<=a<=b<=10⁵
在这里插入图片描述
以上是本题相关的证明，下面我们一起来看看怎么做这道题

#include<iostream>
using namespace std;
const int mod=1e9+7,N=1e5+10;
typedef long long LL;
long long fac[N],infac[N];
int quick_pow(int a, int k, int p)
{
    int res = 1;
    while (k)
    {
        if (k & 1) res = (LL)res * a % p;
        a = (LL)a * a % p;
        k >>= 1;
    }
    return res;
}
int main()
{
    int n;
    fac[0]=infac[0]=1;
    for(int i=1;i<=1e5;i++)
    {
        fac[i]=fac[i-1]*i%mod;
        infac[i]=(LL)infac[i - 1] * quick_pow(i,mod-2,mod)%mod;
    }
    cin>>n;
    while(n--)
    {
        int a,b;
        cin>>a>>b;
        cout<<(LL)fac[a] * infac[b] % mod * infac[a - b] % mod<<endl;
    }
}

三、Lucas定理是用来求 c(n,m) mod p，p为素数的值。

!](https://img-blog.csdnimg.cn/20210506150853593.png)
也就是Lucas(n,m)%p=Lucas(n/p,m/p)*C(n%p,m%p)%p

证明看这里

求上式的时候，Lucas递归出口为m=0时返回1

求C(n%p, m%p)%p的时候，此处写成C(n, m)%p(p是素数，n和m均小于p)

C(n, m)%p = n! / (m ! * (n - m )!) % p = n! * mod_inverse[m! * (n - m)!, p] % p

由于p是素数，有费马小定理可知，m! * (n - m)! 关于p的逆元就是m! * (n - m)!的p-2次方。

#include<iostream>
#include<algorithm>

using namespace std;

typedef long long LL;

int qmi(int a,int k,int p)
{
    int res = 1;
    while(k)
    {
        if(k&1)res = (LL)res*a%p;
        a = (LL)a*a%p;
        k>>=1;
    }
    return res;
}

int C(int a,int b,int p)//自变量类型int
{
    if(b>a)return 0;//漏了边界条件
    int res = 1;
    // a!/(b!(a-b)!) = (a-b+1)*...*a / b! 分子有b项
    for(int i=1,j=a;i<=b;i++,j--)//i<=b而不是<
    {
        res = (LL)res*j%p;
        res = (LL)res*qmi(i,p-2,p)%p;
    }
    return res;
}
//对公式敲
int lucas(LL a,LL b,int p)
{
    if(a<p && b<p)return C(a,b,p);//lucas递归终点是C_{b0}^{a0}
    return (LL)C(a%p,b%p,p)*lucas(a/p,b/p,p)%p;//a%p后肯定是<p的,所以可以用C(),但a/p后不一定<p 所以用lucas继续递归
}

int main()
{
    int n;
    cin >> n;
    while(n--)
    {
        LL a,b;
        int p;
        cin >> a >> b >> p;
        cout << lucas(a,b,p) << endl;
    }
    return 0;
}

p较小的时候预处理出1-p内所有阶乘%p的值，然后用快速幂求出逆元，就可以求出解。p较大的时候只能逐项求出分母和分子模上p的值，然后通过快速幂求逆元求解。

P较小，打表：

 1 const int maxn = 1e5 + 10;
 2 ll fac[maxn];//阶乘打表
 3 void init(ll p)//此处的p应该小于1e5，这样Lucas定理才适用
 4 {
 5     fac[0] = 1;
 6     for(int i = 1; i <= p; i++)
 7         fac[i] = fac[i - 1] * i % p;
 8 }
 9 ll pow(ll a, ll b, ll m)
10 {
11     ll ans = 1;
12     a %= m;
13     while(b)
14     {
15         if(b & 1)ans = (ans % m) * (a % m) % m;
16         b /= 2;
17         a = (a % m) * (a % m) % m;
18     }
19     ans %= m;
20     return ans;
21 }
22 ll inv(ll x, ll p)//x关于p的逆元，p为素数
23 {
24     return pow(x, p - 2, p);
25 }
26 ll C(ll n, ll m, ll p)//组合数C(n, m) % p
27 {
28     if(m > n)return 0;
29     return fac[n] * inv(fac[m] * fac[n - m], p) % p;
30 }
31 ll Lucas(ll n, ll m, ll p)
32 {
33     if(m == 0)return 1;
34     return C(n % p, m % p, p) * Lucas(n / p, m / p, p) % p;
35 }

四、n，m较大且p不为素数的时候

扩展Lucas定理：

在这里插入图片描述

1 #include<bits/stdc++.h>
 2 using namespace std;
 3 typedef long long ll;
 4 const int maxn = 1e6 + 10;
 5 const int mod = 1e9 + 7;
 6 ll pow(ll a, ll b, ll m)
 7 {
 8     ll ans = 1;
 9     a %= m;
10     while(b)
11     {
12         if(b & 1)ans = (ans % m) * (a % m) % m;
13         b /= 2;
14         a = (a % m) * (a % m) % m;
15     }
16     ans %= m;
17     return ans;
18 }
19 ll extgcd(ll a, ll b, ll& x, ll& y)
20 //求解ax+by=gcd(a, b)
21 //返回值为gcd(a, b)
22 {
23     ll d = a;
24     if(b)
25     {
26         d = extgcd(b, a % b, y, x);
27         y -= (a / b) * x;
28     }
29     else x = 1, y = 0;
30     return d;
31 }
32 ll mod_inverse(ll a, ll m)
33 //求解a关于模上m的逆元
34 //返回-1表示逆元不存在
35 {
36     ll x, y;
37     ll d = extgcd(a, m, x, y);
38     return d == 1 ? (m + x % m) % m : -1;
39 }
40 
41 ll Mul(ll n, ll pi, ll pk)//计算n! mod pk的部分值  pk为pi的ki次方
42 //算出的答案不包括pi的幂的那一部分
43 {
44     if(!n)return 1;
45     ll ans = 1;
46     if(n / pk)
47     {
48         for(ll i = 2; i <= pk; i++) //求出循环节乘积
49             if(i % pi)ans = ans * i % pk;
50         ans = pow(ans, n / pk, pk); //循环节次数为n / pk
51     }
52     for(ll i = 2; i <= n % pk; i++)
53         if(i % pi)ans = ans * i % pk;
54     return ans * Mul(n / pi, pi, pk) % pk;//递归求解
55 }
56 
57 ll C(ll n, ll m, ll p, ll pi, ll pk)//计算组合数C(n, m) mod pk的值 pk为pi的ki次方
58 {
59     if(m > n)return 0;
60     ll a = Mul(n, pi, pk), b = Mul(m, pi, pk), c = Mul(n - m, pi, pk);
61     ll k = 0, ans;//k为pi的幂值
62     for(ll i = n; i; i /= pi)k += i / pi;
63     for(ll i = m; i; i /= pi)k -= i / pi;
64     for(ll i = n - m; i; i /= pi)k -= i / pi;
65     ans = a * mod_inverse(b, pk) % pk * mod_inverse(c, pk) % pk * pow(pi, k, pk) % pk;//ans就是n! mod pk的值
66     ans = ans * (p / pk) % p * mod_inverse(p / pk, pk) % p;//此时用剩余定理合并解
67     return ans;
68 }
69 
70 ll Lucas(ll n, ll m, ll p)
71 {
72     ll x = p;
73     ll ans = 0;
74     for(ll i = 2; i <= p; i++)
75     {
76         if(x % i == 0)
77         {
78             ll pk = 1;
79             while(x % i == 0)pk *= i, x /= i;
80             ans = (ans + C(n, m, p, i, pk)) % p;
81         }
82     }
83     return ans;
84 }
85 
86 int main()
87 {
88     ll n, m, p;
89     while(cin >> n >> m >> p)
90     {
91         cout<<Lucas(n, m, p)<<endl;
92     }
93     return 0;
94 }

以上内容来自博客：https://www.cnblogs.com/fzl194/p/9095177.html
深表感谢

五、直接求C(a,b)

输入 a,b，求 C(a,b) 的值。
注意结果可能很大，需要使用高精度计算。
输入格式
共一行，包含两个整数 a 和 b。
输出格式
共一行，输出 C(a,b) 的值。
数据范围
1≤b≤a≤5000
输入样例：
5 3
输出样例：
10

#include<iostream>
#include<algorithm>
#include<vector>

using namespace std;

const int N=5010;

int primes[N],cnt;
int sum[N];
bool st[N];

void get_primes(int n)
{
    for(int i=2;i<=n;i++)
    {
        if(!st[i])primes[cnt++]=i;
        for(int j=0;primes[j]*i<=n;j++)
        {
            st[primes[j]*i]=true;
            if(i%primes[j]==0)break;//==0每次漏
        }
    }
}
// 对p的各个<=a的次数算整除下取整倍数
int get(int n,int p)
{
    int res =0;
    while(n)
    {
        res+=n/p;
        n/=p;
    }
    return res;
}
//高精度乘
vector<int> mul(vector<int> a, int b)
{
    vector<int> c;
    int t = 0;
    for (int i = 0; i < a.size(); i ++ )
    {
        t += a[i] * b;
        c.push_back(t % 10);
        t /= 10;
    }
    while (t)
    {
        c.push_back(t % 10);
        t /= 10;
    }
    // while(C.size()>1 && C.back()==0) C.pop_back();//考虑b==0时才有pop多余的0 b!=0不需要这行
    return c;
}

int main()
{
    int a,b;
    cin >> a >> b;
    get_primes(a);

    for(int i=0;i<cnt;i++)
    {
        int p = primes[i];
        sum[i] = get(a,p)-get(a-b,p)-get(b,p);//是a-b不是b-a
    }

    vector<int> res;
    res.push_back(1);

    for (int i = 0; i < cnt; i ++ )
        for (int j = 0; j < sum[i]; j ++ )//primes[i]的次数
            res = mul(res, primes[i]);

    for (int i = res.size() - 1; i >= 0; i -- ) printf("%d", res[i]);
    puts("");

    return 0;
}

练习题目：
3356 鸣人的影分身

【题目描述】
在火影忍者的世界里，令敌人捉摸不透是非常关键的。我们的主角漩涡鸣人所拥有的一个招数多重影分身之术就是一个很好的例子。影分身是由鸣人身体的查克拉能量制造的，使用的查克拉越多，制造出的影分身越强。针对不同的作战情况，鸣人可以选择制造出各种强度的影分身，有的用来佯攻，有的用来发起致命一击。
那么问题来了，假设鸣人的查克拉能量为M，他影分身的个数为N，那么制造影分身时有多少种（用K表示）不同的分配方法？另外还有以下条件：
①影分身最低为0点查克拉
②每个分身分到的查克拉都是整数
③查克拉可以不分配完 ④这里每个影分身被认为都是不同的【输入】
第一行是测试数据的数目t (0 <= t <=10) 以下每行均包含二个整数M和N，以空格分开。1<=M，N<=1000000。
【输出】
对输入的每组数据M和N，用一行输出相应的K，结果对1e9+7取模。
【样例输入】
2
3 2
4 3
【样例输出】
10
35

【解题思路】
在这里插入图片描述
【参考答案】

#include <bits/stdc++.h>
using namespace std;
#define ll long long
const ll mod = 1000000007;
ll Count(ll n, int p)
{
    ll cnt = 0;
    for (ll i = p; i <= n; i *= p)
        cnt += n / i;
    return cnt;
}
ll PowerMod(ll a, ll b, ll c)
{
    ll ans = 1;
    a = a % c;
    while (b > 0) {
        if (b % 2 == 1)    ans = (ans * a) % c;
        b = b / 2;
        a = (a * a) % c;
    }
    return ans;
}
int main()
{
    int t;
    cin >> t;
    while (t--) {
        ll ans = 1;
        int n, m;
        cin >> m >> n;
        m = m + n; // m + (n + 1) - 1
        n = n; // (n + 1) - 1
        bool primes[m + 1];
        memset(primes, 0, sizeof(primes));
        for (int i = 2; i <= m; i++) {
            if (!primes[i]) {
                ll cnt = Count(m, i) - Count(n, i) - Count(m - n, i);
                ans = (ans * PowerMod(i, cnt, mod)) % mod;
                for (ll j = (ll)i * i; j <= m; j += i)
                    primes[j] = true;
            }
        }
        cout << ans << endl;
    }
}

练习题目

青原樱

信奥教练Andy

关注

0
点赞
踩
3

收藏

觉得还不错? 一键收藏
打赏
0
评论
NOIP中的数学 -- 第7课组合数

1.当n,m都很小的时候可以利用杨辉三角直接求。C(n,m)=C(n-1,m)+C(n-1,m-1)；2、n和m较大，但是p为素数的时候Lucas定理是用来求 c(n,m) mod p，p为素数的值。C(n,m)%p=C(n/p,m/p)*C(n%p,m%p)%p也就是Lucas(n,m)%p=Lucas(n/p,m/p)*C(n%p,m%p)%p求上式的时候，Lucas递归出口为m=0时返回1求C(n%p, m%p)%p的时候，此处写成C(n, m)%p(p是素数，n和m均小于p)C(n,
复制链接

扫一扫