【拓展Lucas】模板

\(C_n^m \mod p\),写得太丑了qwq。
第一次写拓展Lucas竟然是在胡策的时候qwq写了两个半小时啊_(:з」∠)还写挂了一个地方qwq
当然今天胡策我也是第一次写中国剩余定理(ˇˍˇ)
↑平时懒得动手的后果→

#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;

int ipow2(int a, int b) {
    int ret = 1, w = a;
    while (b) {
        if (b & 1) ret = ret * w;
        w = w * w;
        b >>= 1;
    }
    return ret;
}

int ipow(int a, int b, int c) {
    int ret = 1, w = a;
    while (b) {
        if (b & 1) ret = 1ll * ret * w % c;
        w = 1ll * w * w % c;
        b >>= 1;
    }
    return ret;
}

int pic;

const int N = 100003;

struct data {int num, zhi;};

data work(int n, int p, int c) {
    if (n == 1 || n == 0) return (data) {1, 0};
    int ret1, las = n % pic, ret2 = 1;
    
    for (int i = 1; i <= las; ++i) {
        if (i % p == 0) continue;
        ret2 = 1ll * ret2 * i % pic;
    }
    
    if (n >= pic) {
        ret1 = ret2;
        for (int i = las + 1; i < pic; ++i) {
            if (i % p == 0) continue;
            ret1 = 1ll * ret1 * i % pic;
        }
        ret1 = ipow(ret1, n / pic, pic);
    } else
        ret1 = 1;
    
    
    int tim = n / p;
    data s = work(n / p, p, c);
    return (data) {1ll * ret1 * ret2 % pic * s.num % pic, tim + s.zhi};
}

void exgcd(int a, int b, ll &x, ll &y) {
    if (b == 0) {x = 1; y = 0;}
    else {
        exgcd(b, a % b, x, y);
        ll t = x;
        x = y;
        y = t - a / b * y;
    }
}

int inv(int a, int b) {
    ll x, y;
    exgcd(a, b, x, y);
    return (x % b + b) % b;
}

int Lucas(int n, int m, int p, int c) {
    data s1, s2, s3;
    pic = ipow2(p, c);
    s1 = work(n, p, c);
    s2 = work(m, p, c);
    s3 = work(n - m, p, c);
    
    int tim = s1.zhi - s2.zhi - s3.zhi;
    int num1 = s1.num, num2 = s2.num, num3 = s3.num;
    num2 = inv(num2, pic); num3 = inv(num3, pic);
    num1 = 1ll * num1 * num2 % pic * num3 % pic;
    for (int i = 1; i <= tim; ++i)
        num1 = 1ll * num1 * p % pic;
    return num1;
}

int prime[N], cnt = 0, cc[N], pc[N];

void Div(int p) {
    int sq = ceil(sqrt(p));
    for (int i = 2; i <= sq; ++i) {
        if (p % i == 0) {
            prime[++cnt] = i;
            pc[cnt] = 1;
            cc[cnt] = 0;
            while (p % i == 0) {++cc[cnt]; p /= i; pc[cnt] *= i;}
        }
    }
    if (p > 1) prime[++cnt] = p, cc[cnt] = 1, pc[cnt] = p;
}

int n, m, k, p, x;
int num[N];

int main() {
    int T; scanf("%d", &T);
    while (T--) {
        scanf("%d%d%d%d", &n, &m, &k, &p);
        
        cnt = 0;
        Div(p);
        
        for (int i = 1; i <= cnt; ++i)
            num[i] = Lucas(n, m, prime[i], cc[i]);
        
        x = 0;
        for (int i = 1; i <= cnt; ++i)
            (x += (1ll * (p / pc[i]) * inv((p / pc[i]) % pc[i], pc[i]) % p * num[i] % p)) %= p;
        printf("%d\n", x);
    }
    
    return 0;
}

转载于:https://www.cnblogs.com/abclzr/p/6413218.html

  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值