zoukankan      html  css  js  c++  java
  • 组合数取模&&Lucas定理题集

    题集链接:

    https://cn.vjudge.net/contest/231988

    解题之前请先了解组合数取模和Lucas定理

    A : FZU-2020 

    输出组合数C(n, m) mod p

    (1 <= m <= n <= 10^9, m <= 10^4, m < p < 10^9, p是素数)

    由于p较大,不可以打表,直接Lucas求解

     1 #include<iostream>
     2 using namespace std;
     3 typedef long long ll;
     4 const int maxn = 1e5 + 10;
     5 ll fac[maxn];//阶乘打表
     6 void init(ll p)//此处的p应该小于1e5,这样Lucas定理才适用
     7 {
     8     fac[0] = 1;
     9     for(int i = 1; i <= p; i++)
    10         fac[i] = fac[i - 1] * i % p;
    11 }
    12 ll pow(ll a, ll b, ll m)
    13 {
    14     ll ans = 1;
    15     a %= m;
    16     while(b)
    17     {
    18         if(b & 1)ans = (ans % m) * (a % m) % m;
    19         b /= 2;
    20         a = (a % m) * (a % m) % m;
    21     }
    22     ans %= m;
    23     return ans;
    24 }
    25 ll inv(ll x, ll p)//x关于p的逆元,p为素数
    26 {
    27     return pow(x, p - 2, p);
    28 }
    29 ll C(ll n, ll m, ll p)//组合数C(n, m) % p
    30 {
    31     if(m > n)return 0;
    32     ll up = 1, down = 1;//分子分母;
    33     for(int i = n - m + 1; i <= n; i++)up = up * i % p;
    34     for(int i = 1; i <= m; i++)down = down * i % p;
    35     return up * inv(down, p) % p;
    36 }
    37 ll Lucas(ll n, ll m, ll p)
    38 {
    39     if(m == 0)return 1;
    40     return C(n % p, m % p, p) * Lucas(n / p, m / p, p) % p;
    41 }
    42 int main()
    43 {
    44     ll n, m, p;
    45     int T;
    46     cin >> T;
    47     while(T--)
    48     {
    49         cin >> n >> m >> p;
    50         //init(p);
    51         cout<<Lucas(n, m, p)<<endl;
    52     }
    53     return 0;
    54 }

    B:HDU-3944

    求从顶点到达第n行,第m列的元素的最短路径值(最开始是第0行,第0列)

    只能向下或者斜向右。

    如图,如果在右侧,可以走红色路线,比如第4行第3列,与此对称的是第4行第1列(从第0行,第0列数起)

    这样,我们就可以把左侧的转化成右侧来统一计算。中间的也可以这样计算。

    左侧转化成右侧:对于n行m列而言,如果m * 2 < n 那么m = n - m

    对于右侧的如何计算呢?

    用C(x, y)表示x行y列的值,也就是组合数的值

    那么ans = C(n , m) + C(n - 1, m) + C(n - 2, m) + ... + C(m + 1, m) + C(m, m) + C(m - 1, m - 1) + C(m - 2, m - 2) + C(0, 0)

    红色部分为m

    蓝色部分可以通过组合数公式C(n, m) + C(n, m + 1) = C(n +1, m + 1)来计算

    把蓝色部分的C(m, m)写成C(m + 1, m + 1)

    那么最后两项就可以通过上面那个公式合并成C(m + 2, m + 1)

    不断合并,最后得到C(n, m) + C(n, m + 1) = C(n +1, m + 1)

    所以答案就是C(n + 1, m + 1) + m % p;

    由于计算量大,p的范围是10000内的素数,而需要计算10^6个组合数,所以需要先打表,把10000内的素数筛选出来,打表每个阶乘模上素数的值,还有逆元,这样计算组合数才不会超时,而且数组不能开long long,会超时

     1 #include<bits/stdc++.h>
     2 using namespace std;
     3 typedef long long ll;
     4 const int maxn = 1e4 + 10;
     5 int fac[maxn][maxn], inv[maxn][maxn];//阶乘打表 逆元打表
     6 int prime[maxn], pth[maxn];
     7 int pow(int a, int b, int m)
     8 {
     9     int ans = 1;
    10     a %= m;
    11     while(b)
    12     {
    13         if(b & 1)ans = (ans % m) * (a % m) % m;
    14         b /= 2;
    15         a = (a % m) * (a % m) % m;
    16     }
    17     ans %= m;
    18     return ans;
    19 }
    20 bool is_prime[maxn];
    21 void init()
    22 {
    23     int n = 10000, cnt = 0;
    24     for(int i = 2; i <= n; i++)is_prime[i] = 1;
    25     for(int i = 2; i <= n; i++)
    26     {
    27         if(is_prime[i])
    28         {
    29             prime[++cnt] = i;
    30             pth[i] = cnt;
    31             for(int j = i * i; j <= n; j += i)is_prime[j] = 0;
    32         }
    33     }
    34     //cout<<cnt<<endl;
    35     for(int i = 1; i <= cnt; i++)
    36     {
    37         fac[i][0] = inv[i][0] = 1;
    38         for(int j = 1; j < prime[i]; j++)//大于i的阶乘模上i均为0
    39         {
    40             fac[i][j] = fac[i][j - 1] * j % prime[i];
    41             inv[i][j] = pow(fac[i][j], prime[i] - 2, prime[i]);
    42         }
    43     }
    44 }
    45 
    46 int C(ll n, int m, int p)//组合数C(n, m) % p
    47 {
    48     if(m > n)return 0;
    49     if(m == n)return 1;
    50     int t = pth[p];
    51     return fac[t][n] * (inv[t][m] * inv[t][n - m] % p) % p;
    52 }
    53 
    54 int Lucas(int n, int m, int p)
    55 {
    56     //cout<<n<<" "<<m<<" "<<p<<endl;
    57     if(m == 0)return 1;
    58     return C(n % p, m % p, p) * Lucas(n / p, m / p, p) % p;
    59 }
    60 
    61 int main()
    62 {
    63     init();
    64     int n, m, p;
    65     int cases = 0;
    66     while(scanf("%d%d%d", &n, &m, &p) != EOF)
    67     {
    68         if(m <= n / 2)
    69         {
    70             m = n - m;
    71         }
    72         ll ans = (Lucas(n + 1, m + 1, p) + m) % p;
    73         printf("Case #%d: %d
    ", ++cases, ans);
    74     }
    75     return 0;
    76 }

    C:ZOJ-3557

    在n个元素的集合中取出lm和不相邻元素的种数

    插空法,由于不相邻,那就把这m个先取出来,放在剩下的n-m个数字的两侧的空位上,那么这m个一定不相邻。

    等价于n-m+1中放m个元素

    C(n - m + 1, m) % p

     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 inv(ll x, ll p)//x关于p的逆元,p为素数
    20 {
    21     return pow(x, p - 2, p);
    22 }
    23 ll C(ll n, ll m, ll p)//组合数C(n, m) % p
    24 {
    25     if(m > n)return 0;
    26     ll up = 1, down = 1;//分子分母;
    27     for(int i = n - m + 1; i <= n; i++)up = up * i % p;
    28     for(int i = 1; i <= m; i++)down = down * i % p;
    29     return up * inv(down, 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 }
    36 int main()
    37 {
    38     ll n, m, p;
    39     while(cin >> n >> m >> p)
    40     {
    41         if(n - m + 1 < m)
    42             cout<<"0"<<endl;
    43         else cout<<Lucas(n - m + 1, m, p)<<endl;
    44     }
    45     return 0;
    46 }

    D在F介绍完之后再介绍

    E:hdu-4349

    求C(n,0),C (n,1),C (n,2)...C (n,n)中奇数的数目

    首先:组合数奇偶性判断公式:n  & m == m

    当然本题要判断的组合数很多,所以不能用上述结论,只能另辟蹊径。由于是判断奇偶性,那么就是判断  是否为1,利用Lucas定理,先把化为二进制,这样它们都是01序列了。我们又知道    。这样中为0的地方对应的中的位置只有一种可能,那就是0。

    这样我们可以不用管中为0的地方,只考虑中为1的位置,可以看出,中为1的位置对应的中为0或1,其结果都是1,这样答案就是:1<<(二进制表示中1的个数)

     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 int lowbit(int x)
     7 {
     8     return x & (-x);
     9 }
    10 int main()
    11 {
    12     ll n, m, p;
    13     while(cin >> n)
    14     {
    15         int tot = 0;
    16         while(n)
    17             n -= lowbit(n), tot++;
    18         cout<<(1<<tot)<<endl;
    19     }
    20     return 0;
    21 }

    F:Gym - 100633J 

    扩展Lucas定理模板

    这里的p不是素数,用中国剩余定理合并同余方程组。

    详细介绍请看全文最开始的传送门

     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 }

    D:HDU-4373

    个for循环嵌套,有两种形式,第一类从1开始到,第二类从上一层循环当前数开始到,第一层一定是第一种类型,求总的循环的次数对364875103取余的结果。

    首先可以看出,每一个第一类循环都是一个新的开始,与前面的状态无关,所以可以把个嵌套分为几个不

         同的部分,每一个部分由第一类循环开始,最终结果相乘即可。剩下的就是第二类循环的问题,假设一个

         层循环,最大到,分析一下得到如下结果

        

         (1)只有一层,则循环次数为

     

         (2)只有前两层,则循环次数为

     

             

     

         (3)只有前三层,则循环次数为

     

             

     

          由此得到结论:第的循环次数为

    由于364875103=97 * 3761599

    用中国剩余定理合并,但是不能直接套上述模板,会超时,上述模板是对于不同的合数p分解成质数的次方,算起来很麻烦,此处364875103就等于两个质数的乘积

    这里a1 = 97, a2 = 3761599所以M1 = 3761599, M2 = 97, 求出M1关于m1的逆元,M2关于m2的逆元

    然后之需要求出组合数mod97的值为a1,组合数mod3761599的值为a2

    ans按照上述公式就可以直接算出来了。

     1 #include<bits/stdc++.h>
     2 using namespace std;
     3 typedef long long ll;
     4 const int maxn = 1e6 + 10;
     5 ll pow(ll a, ll b, ll m)
     6 {
     7     ll ans = 1;
     8     a %= m;
     9     while(b)
    10     {
    11         if(b & 1)ans = (ans % m) * (a % m) % m;
    12         b /= 2;
    13         a = (a % m) * (a % m) % m;
    14     }
    15     ans %= m;
    16     return ans;
    17 }
    18 const ll mod1 = 97;
    19 const ll mod2 = 3761599;
    20 const ll mod = mod1 * mod2;
    21 ll fac1[mod1 + 10], fac2[mod2 + 10];
    22 ll inv1, inv2;
    23 void init()
    24 {
    25     fac1[0] = fac2[0] = 1;
    26     for(int i = 1; i < mod1; i++)fac1[i] = fac1[i - 1] * i % mod1;
    27     for(int i = 1; i < mod2; i++)fac2[i] = fac2[i - 1] * i % mod2;
    28     inv1 = pow(mod2, mod1 - 2, mod1);//mod1 模上mod2的逆元
    29     inv2 = pow(mod1, mod2 - 2, mod2);//mod2 模上mod1的逆元 求出在之后的中国剩余定理合并的时候不用计算
    30 }
    31 
    32 ll C(ll n, ll m, ll p, ll fac[])
    33 {
    34     if(m > n)return 0;
    35     return fac[n] * pow(fac[m] * fac[n - m], p - 2, p) % p;
    36 }
    37 
    38 ll Lucas(ll n, ll m, ll p, ll fac[])
    39 {
    40     if(m == 0)return 1;
    41     return C(n % p, m % p, p, fac) * Lucas(n / p, m / p, p, fac) % p;
    42 }
    43 int main()
    44 {
    45     ll n, m, k, cases = 0;
    46     ll a[20];
    47     int T;
    48     init();
    49     cin >> T;
    50     while(T--)
    51     {
    52         cin >> n >> m;
    53         cin >> k;
    54         for(int i = 0; i < k; i++)cin >> a[i];
    55         a[k] = m;
    56         ll ans = 1;
    57         for(int i = 1; i <= k; i++)
    58         {
    59             ll t = a[i] - a[i - 1];
    60             ll a1 = Lucas(n + t - 1, t, mod1, fac1);
    61             ll a2 = Lucas(n + t - 1, t, mod2, fac2);
    62             ll a3 = (a1 * mod2 % mod * inv1 + a2 * mod1 % mod * inv2) % mod;//中国剩余定理合并同余方程组
    63             ans = ans * a3 % mod;
    64         }
    65         cout<<"Case #"<<++cases<<": "<<ans<<endl;
    66     }
    67     return 0;
    68 }
  • 相关阅读:
    Python常用函数
    Mock测试&Postman mockserver详细教程
    openpyxl模块
    adb 'grep' 不是内部或外部命令,也不是可运行的程序或批处理文件
    Appium-Python-Windows环境搭建笔记
    调用类方法时报错:missing 1 required positional argument: 'self'
    RE正则表达式-元字符
    微分方程
    操作系统学习记录
    Mybatis基础配置
  • 原文地址:https://www.cnblogs.com/fzl194/p/9105734.html
Copyright © 2011-2022 走看看