HDU5446 lucas定理和模线性方程组 xingxing在努力

   这道题是求解C(n, m)%M, 其中M可以被分解为几个不同素数的乘积, 假设答案是X, 那么就有 x == C(n, m) mod(M) =》x == C(n, m) mod(pi)(具体证明详见《具体数学》)  这样就将一个式子化成了几个同余方程, 我们在求解出c=C(n,m)%pi, 方程进一步化简为 x%pi = c, 这样我们解这几个同余方程组即可得到答案。。代码如下:

 

 

#include <cstdio>
#include <cstring>
#include <algorithm>

using namespace std;
typedef long long LL;

LL n, m, k;
LL p[15];           //素数
LL c[15];           //卢卡斯取摸的结果

//卢卡斯模板

LL qk_multi(LL a, LL b, LL p)
{
    LL ans = 0;
    while(b>0)
    {
        if((b&1)==1) ans = (ans+a)%p;
        a = (a+a)%p;
        b >>= 1;
    }
    return ans;
}

LL qk_mod(LL a, LL b, LL p)
{
    LL res = 1;
    while(b>0)
    {
        if((b&1)==1) res = qk_multi(res, a, p);
        a = qk_multi(a, a, p);
        b >>= 1;
    }
    return res;
}

LL Comb(LL a, LL b, LL p)
{
    if(a < b) return 0;
    if(a == b) return 1;
    if(b > a-b) b = a-b;
    LL ans=1, ca=1, cb=1;
    for(LL i=0; i<b; i++)
    {
        ca = (ca*(a-i))%p;
        cb = (cb*(b-i))%p;
    }
    ans = (ca*qk_mod(cb, p-2, p))%p;
    return ans;
}

LL Lucas(LL n, LL m, LL p)
{
    LL ans = 1;
    while(n && m && ans)
    {
        ans = (ans*Comb(n%p, m%p, p))%p;
        n/=p; m/=p;
    }
    return ans;
}

//拓展欧几里得算法
void egcd(LL a, LL b, LL&d, LL&x, LL&y)
{
    if(b == 0) { d=a; x=1; y=0; }
    else { egcd(b, a%b, d, y, x); y-=x*(a/b); }
}

LL Jie(LL a, LL b, LL c)    //返回x0
{
    LL x0, y0, gab;
    egcd(a, b, gab, x0, y0);
    if(c%gab != 0) return -1;
    else
    {
        LL minx = (c/gab*x0)%(b/gab);
        if(minx < 0) minx += b/gab;
        return minx;
    }
}

int main()
{
    int T;
    //freopen("date.txt", "r", stdin);
    //freopen("out.txt", "w", stdout);
    scanf("%d", &T);
    while(T--)
    {
        scanf("%lld%lld%lld", &n, &m, &k);
        for(int i=0; i<k; i++)
        {
            scanf("%lld", &p[i]);
            c[i] = Lucas(n, m, p[i]);
        }
        if(k==1)
        {
            printf("%lld\n", c[0]);
            continue;
        }
        LL res;
        for(int i=1; i<k; i++)
        {
            LL k1 = Jie(p[0], p[i], c[i]-c[0]);
            res = c[0] + k1*p[0];
            p[0] = p[0]*p[i];
            c[0] = res;
        }
        printf("%lld\n", res);
    }
    return 0;
}

 

posted @ 2015-11-23 20:32  xing-xing  阅读(149)  评论(0)    收藏  举报