HDU 5446 Lucas + 中国剩余定理 + 快速乘法

HDU 5446

题目链接：

http://acm.hdu.edu.cn/showproblem.php?pid=5446

题意：

求，n，m的范围很大

思路：

Lucas + 中国剩余定理 + 快速乘法。

Lucas定理用来求当n和m很大时的C（n,m）%p（p是素数）。证明网上有，具体公式就是

Lucas(n,m,p) = C(n%p,m%p）*Lucas(n/p,m/p,p)

中国剩余定理用来求同余方程。对方程n的最小解。

设，，则

因为本题数据会爆long long，所以需要用快速乘法，具体见代码，类似快速幂。

源码：

include

include

include

include

include

include

include

include

using namespace std;

define LL long long

const int MAXN = 15;

LL prime[MAXN];

LL fac[100005], inv[100005];

LL lv[MAXN];

void exceed_gcd(LL a, LL b, LL &ans, LL &x, LL &y)

{

if(a < b){

exceed_gcd(b, a, ans, y, x);

return;

}

if(b == 0){

x = 1, y = 0;

ans = a;

return;

}

else{

exceed_gcd(b, a % b, ans, y, x);

y -= x * (a / b);

}

}

LL ppow(LL a, LL x, LL mod)

{

LL ans = 1;

while(x){

// printf(“mod = %I64d, a = %I64d, x = %I64d, ans = %I64d\n”, mod, a, x, ans);

if(x & 1)

ans = (ans * a) % mod;

a = (a * a) % mod;

x >>= 1;

}

return ans;

}

LL C(LL n, LL m, LL mod)

{

if(n < m || n < 0 || m < 0) return 0;

// printf(“n = %I64d, m = %I64d, fac
= %I64d, in[m] = %I64d, inv[n - m] = %I64d\n”, n, m, fac
, inv[m], inv[n - m]);

return fac
* inv[m] % mod * inv[n - m] % mod;

}

LL lucas(LL n, LL m, LL mod)

{

// printf(“n = %I64d, m = %I64d\n”, n, m);

if(n < m) {/printf(“n=%I64d,m=%I64d\n”,n,m);/return 0;}

if(n == 0 || m == 0) return 1;

// printf(“n = %I64d, m = %I64d, C = %I64d, lucas = %I64d\n”, n, m, C(n % mod, m % mod, mod), lucas(n / mod, m / mod, mod));

return C(n % mod, m % mod, mod) * lucas(n / mod, m / mod, mod) % mod;

}

LL mul(LL a, LL b, LL mod)

{

a = (a % mod + mod) % mod;

b = (b % mod + mod) % mod;

LL ans = 0;

while(b){

if(b & 1)

ans = (ans + a) % mod;

b >>= 1;

a <<= 1;

a = a % mod;

}

return ans;

}

int main()

{

// freopen(“1006.in”, “r”, stdin);

LL n, m, cnt;

int t;

scanf(“%d”, &t);

while(t–){

scanf(“%I64d%I64d%I64d”, &n, &m, &cnt);

LL M = 1;

LL d, y;

LL ans = 0;

for(int i = 0 ; i < cnt ; i++){

scanf(“%I64d”, &prime[i]);

M *= prime[i];

fac[0] = 1;

for(int j = 1 ; j < prime[i] ; j++)

fac[j] = fac[j - 1] * j % prime[i];

inv[prime[i] - 1] = ppow(fac[prime[i] - 1], prime[i] - 2, prime[i]);

for(int j = prime[i] - 2 ; j >= 0 ; j–)

inv[j] = (inv[j + 1] * (j + 1)) % prime[i];

// for(int j = 0 ; j < prime[i] ; j++)

// printf(“fac[%d] = %I64d, inv[%d] = %I64d\n”, j, fac[j], j, inv[j]);

lv[i] = lucas(n, m, prime[i]);

}

for(int i = 0 ; i < cnt ; i++){

LL M1 = M / prime[i];

LL M2;

exceed_gcd(prime[i], M1, d, y, M2);

// printf(“i = %d, M1 = %I64d, M2 = %I64d, M = %I64d lv[%d] = %I64d\n”, i, M1, M2, M, i, lv[i]);

// M2 = (M2 % M + M) % M;

ans = (ans + mul(mul(lv[i], M1, M), M2, M));

}

ans = (ans % M + M) % M;

printf(“%I64d\n”, ans);

}

return 0;

}