【总结】扩展欧几里得算法

来源:互联网 发布:java json传值到前台 编辑:程序博客网 时间:2024/06/07 07:32

前言


    最近的一连串考试中不幸考到了自己学习得最不扎实的部分——数学,还是想一点点把那些无论基础还是高深的知识点捡起来 ……

概念

    扩展欧几里德算法,在已知a, b求解一组x,y,使它们满足贝祖等式: ax+by = gcd(a, b) =d,常用于求解模线性方程及方程组。

                                                ——百度百科
    贴上觉得简洁易懂的讲解链接http://blog.csdn.net/acmore_xiong/article/details/47694909
    代码

void Exgcd(int a,int b,int &x,int &y,int &g) {    if(b == 0) {        g = a; x = 1, y = 0; return;    }    Exgcd(b,a%b,y,x,g);    y -= (a/b)*x;}

    然后来看求解不定整数方程的应用吧
    对于 ax+by = c (a,b不全为0)
    g = gcd(a,b)
    若方程无解 c % g != 0
    有解时: c’ = c/g
    然后我们可以用上面的代码解 ax + by = g
    得到一组特殊解 x0,y0, 于是 x1=x0c/g, y1=y0c/g就是原方程的一组特殊解啦
    通解可以表示为 xi=x1+tb/g,yi=y1ta/g (t∈Z)

例题


    青蛙的约会这种经典题不多说
    上昨天考试题Math

题意

    T 组数据,每组问二项式(x^(a/b)+x^(c/d))^n 展开后有多少有理项,又有多少方案使得有理项不相邻,方案数%p.

数据范围

    |a|,|b|,|c|,|d|,T<=105,n<=109,0=<p<=107

思路:

    60 分算法:
    暴力计算有理项数量 x,则 ans = x!*(N+1-x)!*C(N+2-x,x),费马小定理求逆元即可再搞组合数.注意快速幂要先将底数%Mod 否则易炸.
    对C(N+2-x,x)的理解:画图,有N+1-x个无理项,为满足要求,x个有理项只能放在两端或无理项之间的共N+2-x个位置上.
    满分算法:
    将(x^(a/b)+x^(c/d))^n 展开共 n+1 项,每一项可以表示为 x^((a/b)k+(c/d)(n-k))(k∈[0,n])= x^(((ad-bc)*k+bcn)/bd).
    于是问题按照如下转化
    (adbc)k+bcn0(mod(bd))>
    (bcad)k+bdybcn(mod(bd))>
    (bcad)k+bdybcn
    可以使用扩展欧几里得求解 k,计算 k∈[0,n]的可行解个数.
    记可行解个数为 A,则 Ans = A! * (n+1-A)! * C(n+2-A,A). 组合数部分费马小定理+逆元计算,但是要注意 C 应单独写函数判断两种情况,对于 C(m,n): 1)m < n 2)n==0 ,否则在遇到模数P==n+2-A 而 A==0 的情况答案错误.

代码

60分暴力

#include<cstdio>using namespace std;typedef long long LL;const int sm = 8e3 + 50;int T,Mod;LL Ans;LL Fac[sm];LL QMod(LL a,int b) {    LL ret = 1, x = a % Mod;    while(b) {        if(b&1) ret = 1ll*ret*x%Mod;        x = 1ll*x*x%Mod;        b >>= 1;    }    return ret;}int main() {    freopen("math.in","r",stdin);    freopen("math.out","w",stdout);    LL x,y;     LL N,A,B,C,D,Up,Down,Tmp;    scanf("%d%d",&T,&Mod);    Fac[0] = 1;    for(int i = 1; i <= 8010; ++ i)        Fac[i] = 1ll* Fac[i-1] * i % Mod;    for(int i = 1; i <= T; ++ i) {        scanf("%lld%lld%lld%lld%lld",&N,&A,&B,&C,&D);        x = 0;        for(int j = 0; j <= N; ++ j) {            Up = 1ll*A*D*(N-j)+1ll*C*B*j;            Down = 1ll*B*D;            if(Up % Down == 0)++x;        }        y = N + 1 - x;        if(N-2*x+2 <= 0) Ans = 0;        else {            Tmp = 1ll*Fac[N+2-x]*QMod(1ll*Fac[x]*Fac[N+2-2*x],Mod-2)%Mod;            Ans = 1ll*Fac[x]*Fac[y]%Mod;            Ans = Ans*Tmp%Mod;        }        printf("%lld %lld\n",x,Ans);    }    return 0;}

AC代码

#include<cstdio>#include<cmath>#include<algorithm>using namespace std;typedef long long LL;const int sm = 1e7 + 10;const int B = 15, sn = sm/B + 10;int T,P;int A;LL n, Ans;LL gcd(LL u,LL v) { return (v==0)?u:gcd(v,u%v); }LL ret;LL Qmod(LL a,LL b) {    ret = 1; a %= P;    while(b) {        if(b&1) ret = ret*a%P;        a = a*a%P;        b>>=1;    }    return ret;}LL x,y;void Exgcd(LL a,LL b,LL &x,LL &y,LL &g) {    if(b == 0) {        g = a; x = 1, y = 0; return;    }    Exgcd(b,a%b,y,x,g);    y -= (a/b)*x;}LL Abs,g;LL a,b,c,d;bool flag;int Calc(LL aa,LL bb) {     Exgcd(aa,bb,x,y,g);    if(b*c*n%g) return ret = 0;     x%=(b*d);    x = x*b % abs(b*d) *c % abs(b*d) *n % abs(b*d)/g;    x = (x % abs(b*d/g) + abs(b*d/g) ) % abs(b*d/g);    ret = ((x<=n+1)?1:0) + (n-x)/(abs(b*d/g));    return ret;}LL fac[sn];void Init() {    LL tmp = 1, Ct = 0, cnt = 0;    fac[0] = 1;    for(int i = 1; i <= 1e7; ++i) {        if(cnt%B == 0) fac[++Ct] = fac[Ct-1] * i % P, cnt = 1;        else fac[Ct] = fac[Ct] * i % P, ++cnt;    } }LL t,q,l;LL Fac(LL u) {    t = u/B; q = 1;    if(u%B!=0) {        l = u%B;        for(int i = 1; i <= l; ++ i)                q = q * (t*B+i) % P;     }    return fac[t] * q % P; }LL C(LL m,LL n) {    if(n==0)  return 1;    if(n> m)  return 0;    return Fac(m)*Qmod(Fac(n)*Fac(m-n)%P,P-2)%P;}int main() {    freopen("math.in","r",stdin);    freopen("math.out","w",stdout);    scanf("%d%d",&T,&P);    Init();    for(int i = 1; i <= T; ++i){        scanf("%lld%lld%lld%lld%lld",&n,&a,&b,&c,&d);        A = Calc(b*c-a*d,b*d);        if(A >= P || n + 1 - A >= P || A > n + 2 - A) Ans = 0;        else Ans = Fac(A) * Fac(n+1-A) % P * C(n+2-A,A) % P;        printf("%d %lld\n",A,Ans);    }    return 0;}
原创粉丝点击