Lucas定理

    xiaoxiao2021-03-25  95

    求C((m, n) % p) (p为素数)

    表达式:C(n,m)%p=C(n/p,m/p)*C(n%p,m%p)%p

    当p<1e6时

    #include<algorithm> #include<cstdio> #include<cmath> using namespace std; typedef long long ll; const int N = 100000; ll f[N]; void init(int p) { f[0] = 1; for (int i=1; i<=p; ++i) f[i] = f[i-1] * i % p; } ll pow_mod(ll a, ll x, int p) { ll ret = 1; while (x) { if (x & 1) ret = ret * a % p; a = a * a % p; x >>= 1; } return ret; } ll Lucas(ll n, ll k, int p) { //C (n, k) % p ll ret = 1; while (n && k) { ll nn = n % p, kk = k % p; if (nn < kk) return 0; //inv (f[kk]) = f[kk] ^ (p - 2) % p ret = ret * f[nn] * pow_mod (f[kk] * f[nn-kk] % p, p - 2, p) % p; n /= p, k /= p; } return ret; } int main(void) { int p = 10007, n , m; while(~scanf("%d %d", &n, &m)) { init (p); printf ("%lld\n", Lucas (n, m, p)); } return 0; } 扩展:当p较大(1e12左右)或者p不是素数时

    #include<algorithm> #include<cstdio> #include<cmath> using namespace std; typedef long long ll; ll POW(ll a, ll b, ll mod) { ll ans=1; while(b) { if(b & 1) ans = ans * a % mod; a = a *a % mod; b >>= 1; } return ans; } ll POW(ll a, ll b) { ll ans = 1; while(b) { if(b & 1) ans = ans * a; a = a * a; b >>= 1; } return ans; } ll exGcd(ll a, ll b, ll &x, ll &y) { ll t, d; if(!b) { x=1; y=0; return a; } d = exGcd(b, a % b, x, y); t = x; x = y; y = t - a / b * y; return d; } bool modular(ll *a, ll *m, ll k) { ll d, t, c, x, y, i; for(i = 2; i <= k; i++) { d = exGcd(m[1], m[i], x, y); c = a[i]-a[1]; if(c % d) return false; t = m[i] / d; x = (c / d * x % t + t) % t; a [1] =m[1] * x + a[1]; m [1] =m[1] * m[i] / d; } return true; } ll reverse(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(m > n) return 0; ll ans = 1, i, a, b; for(i = 1; i <= m; i++) { a = (n + 1 - i) % mod; b = reverse(i % mod, mod); ans = ans * a % mod * b % mod; } return ans; } ll C1(ll n, ll m, ll mod) { if(m == 0) return 1; return C(n % mod, m % mod, mod) * C1(n / mod, m / mod, mod) % mod; } ll cal(ll n, ll p, ll t) { if(!n) return 1; ll x = POW(p, t), i, y = n / x, temp = 1; for(i = 1; i <= x; i++) if(i % p) temp = temp * i%x; ll ans=POW(temp,y,x); for(i=y*x+1;i<=n;i++) if(i%p) ans=ans*i%x; return ans*cal(n/p,p,t)%x; } ll C2(ll n, ll m, ll p, ll t) { ll x=POW(p, t); ll a, b, c, ap = 0, bp = 0, cp = 0, temp; for(temp = n; temp ;temp /= p) ap += temp / p; for(temp = m; temp; temp /= p) bp += temp / p; for(temp = n - m; temp; temp /= p) cp += temp / p; ap = ap - bp - cp; ll ans = POW(p, ap, x); a = cal(n, p, t); b = cal(m, p, t); c = cal(n - m, p, t); ans = ans * a % x * reverse(b, x) % x * reverse(c, x) % x; return ans; } //计算C(n,m)%mod ll Lucas(ll n,ll m,ll mod) { ll i, t, cnt = 0; ll A[205], M[205]; for(i = 2; i * i <= mod; i++) if(mod % i == 0) { t = 0; while(mod % i == 0) { t++; mod /= i; } M[++cnt] = POW(i,t); if(t==1) A[cnt] = C1(n, m, i); else A[cnt] = C2(n, m, i, t); } if(mod > 1) { M[++cnt] = mod; A[cnt] = C1(n, m, mod); } modular(A, M, cnt); return A[1]; } int main(void) { ll mod = 10000006, n , m; while(~scanf("%lld %lld", &n, &m)) { printf ("%lld\n", Lucas (n, m, mod)); } return 0; }

    转载请注明原文地址: https://ju.6miu.com/read-15188.html

    最新回复(0)