POJ 1150 中国剩余定理

题意

传送门 POJ 1150

题解

计算最末一位非零数字，设 $n!/(n-m)!=a{2}^{e1}{5}^{e2}$ ，则原问题等价于求

$n!/(n-m)!/10^{min(e1,e2)}mod(10)$

根据威尔逊定理，阶乘 $n!$ 模素数 $p$ 可以在 $O(lo{g}_{p}n)$ 内求解，则按照中国剩余定理的思路，求 $b_1,b_2$

$$\{\begin{array}{l}n!/(n-m)!/10^{min(e1,e2)}=b_{1}mod(2)\\ n!/(n-m)!/10^{min(e1,e2)}=b_{2}mod(5)\end{array}$$

最后求解线性同余方程组即可。

#include <algorithm>
#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <queue>
#include <vector>
using namespace std;
int N, M;
// 模 2,5 的 n! 和 n 逆元
int fac2[2] = {1, 1};
int inv2[2] = {0, 1};
int fac5[5] = {1, 1, 2, 1, 4};
int inv5[5] = {0, 1, 3, 2, 4};
// 快速幂模运算
int mod_pow(int x, int n, int mod)
{
    int res = 1;
    while (n > 0)
    {
        if (n & 1)
        {
            res = res * x % mod;
        }
        x = x * x % mod;
        n >>= 1;
    }
    return res;
}
// 分解 n! = ap^e, 返回 a mod p
int mod_fact(int n, int p, int &e, int *fact)
{
    e = 0;
    if (n == 0)
    {
        return 1;
    }
    int res = mod_fact(n / p, p, e, fact);
    e += n / p;
    if (n / p % 2 != 0)
    {
        return res * (p - fact[n % p]) % p;
    }
    return res * fact[n % p] % p;
}

int gcd(int a, int b)
{
    if (b == 0)
        return a;
    return gcd(b, a % b);
}
// 求解线性同余方程组
pair<int, int> linear_congruence(const vector<int> &A, const vector<int> &B, const vector<int> &M)
{
    int x = 0, m = 1;
    for (int i = 0; i < A.size(); i++)
    {
        int a = A[i] * m, b = B[i] - A[i] * x, d = gcd(a, M[i]);
        if (b % d != 0)
        {
            return make_pair(0, -1);
        }
        int t = b / d * (i == 0 ? inv2[a / d] : inv5[a / d]) % (M[i] / d);
        x = x + m * t;
        m *= M[i] / d;
    }
    return make_pair(x % m, m);
}

void solve()
{
    int e21, e51;
    int x21 = mod_fact(N, 2, e21, fac2), x51 = mod_fact(N, 5, e51, fac5);
    int e22, e52;
    int x22 = mod_fact(N - M, 2, e22, fac2), x52 = mod_fact(N - M, 5, e52, fac5);
    int e2 = e21 - e22, e5 = e51 - e52, e = min(e2, e5);
    int x2 = x21 * inv2[x22] % 2, x5 = x51 * inv5[x52] % 5;
    vector<int> A(2), B(2), M(2);
    A[0] = 1, B[0] = x2 * inv2[mod_pow(5, e, 2)] % 2 * (e2 > e ? 0 : 1), M[0] = 2;
    A[1] = 1, B[1] = x5 * inv5[mod_pow(2, e, 5)] % 5 * (e5 > e ? 0 : 1), M[1] = 5;
    printf("%d\n", linear_congruence(A, B, M).first);
}

int main()
{
    while (~scanf("%d%d", &N, &M))
    {
        solve();
    }
}