lucas+CRT

lucas+CRT

#10229. 「一本通 6.6 例 4」古代猪文

题目是要求(g^{sum{C_{n}^{n/i}}}) % 999911659;(i能整除n)

模数是质数，phi[999911659] = 999911658, 根据欧拉降幂，g的幂数可以%phi[999911659],由于此时的模数不是一个质数，然后就想着用exlucas来求，但是亲测这样会TLE，所以就只能把999911658拆成4个质数，用lucas结合CRT来求。

假设n是(C_{n}^{m}) ,a[i]是999911658的素数因子，n % a[i] = b[j],根据方程组用CRT就可以求出(C_{n}^{m}) 了

#include <cstdio>
#define LL long long
using namespace std;
void exgcd(LL a, LL b, LL &x, LL &y) {
    if (!b) {
        x = 1, y = 0;
        return;
    }
    exgcd(b, a % b, x, y);
    LL t = x;
    x = y;
    y = t - (a / b) * y;
}
LL inv(LL a, LL b) {
    LL x, y;
    exgcd(a, b, x, y);
    return (x % b + b) % b;
}
LL C(LL n, LL m, LL mod) {
    if(n == m) return 1;
   	LL inv[36000];
   	inv[1] = 1;
   	inv[0] = 1;
    LL f = 1;
    for(int i = n - m + 1; i <= n; i++) f = f * i % mod;//(n-m+1)*(n-m+2)*...*n
    for(int i = 2; i <= m; i++) inv[i] = (mod - mod / i) * 1LL * inv[mod % i] % mod;//打表1~m关于mod的逆元
    for(int i = 2; i <= m; i++) inv[i] = inv[i] * inv[i-1] % mod;//m!关于mod的逆元
    return (f * inv[m]) % mod;
}
LL lucas(LL n, LL m, LL mod) {
    if(m == 0) return 1LL;
    return (C(n % mod, m % mod, mod) * lucas(n / mod, m / mod, mod)) % mod;
}
LL power_mod(LL a, LL b, LL p) {
    LL ret = 1;
    while(b) {
        if(b & 1) ret = ret * a % p;
        a = a * a % p;
        b >>= 1;
    }
    return ret;
}
int main() {
    LL mmod = 999911659;
    LL nmod[5] = {2, 3, 4679, 35617};
    LL n, g, ans = 1, sum = 999911658;
    scanf("%lld %lld", &n, &g);
    for (LL i = 1; i * i <= n; i++) {
        if (n % i == 0) {
            LL a[5], b[5];
            for(int j = 0; j < 4; j++) {
                b[j] = lucas(n, i, nmod[j]);
                a[j] = nmod[j];
            }
            LL num = 0;
            for(int j = 0; j < 4; j++) {
                num = (num + (sum / a[j]) * b[j] * inv(sum / a[j], a[j])) % sum;
            }
            num = num % sum;
            ans = ans * power_mod(g, num, mmod) % mmod;
            if (i * i != n) {
                for(int j = 0; j < 4; j++) {
                    b[j] = lucas(n, n/i, nmod[j]);
                    a[j] = nmod[j];
                    //printf("b %lld
", b[j]);
                }
                //printf("i sum %lld %lld
", i, num);
                num = 0;
                for(int j = 0; j < 4; j++) {
                    num = (num + (sum / a[j]) * b[j] * inv(sum / a[j], a[j])) % sum;
                    //printf("%d %lld
",j, num);
                }
                num = num % sum;
                ans = ans * power_mod(g, num, mmod) % mmod;
            }
        }
    }
    printf("%lld
", ans % mmod);
    return 0;
}