扩展GCD

扩展gcd：已知ax+by=gcd(a,b)，求x,y。

①由gcd(a,b)==gcd(b,a%b)得ax+by==bx+(a%b)y，继续化下去，化成ax1+by1的形式即为递归转移式，因为递归过程中a,b不能改变，改的是x,y。

bx+(a%b)y == bx+ (a-(a/b)*b) y == bx+ay-(a/b)*b*y == ay + b(x-a/b*y)

即exgcd(ax+by)== exgcd( ay + b(x-a/b*y)  )

然后就可以递归求解了，共用相同的x,y，可以全局，下面代码是用引用实现。

②递归终止：当b==0时，式子为ax==gcd(a,b)，那么显然a=gcd(a,b),x=1,y=0;

 //扩展gcd求解x,y同时还把gcd(a,b)求出来了，exgcd的返回值就是gcd(a,b)

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

扩展gcd求出一个解后，那ax+by=gcd(a,b)的其他解为x=x0+k*b; y=y0-k*a;

那么ax+by=c就不难了，将ax+by=gcd(a,b)两边除以gcd可得，那么代码中就正常求exgcd，解x=c/gcd*x0,y=c/gcd*y0。(x0,y0为exgcd的解)

这里注意：x=c/gcd*x0 与 x=x0/gcd*c 不等价！有解情况下gcd一定是c的因子,而x0不一定是。

如8*3/2=12而8*(3/2)!=12。  当然x也可以=c*x0/gcd，不过要小心溢出。

ax+by=c无整数解的情况：c不能被gcd(a,b)整除（也即c不是gcd(a,b)的倍数）。

理由：因为ax与by均是gcd(a,b)的倍数,所以其和ax+by也是gcd(a,b)的倍数。

乘法逆元

满足a*x≡1 (mod p)的x值

通俗理解：(1/a)mod p 和 (a的逆元)mod p一样

作用：求(a/b)%p，当a过大，导致a/b无法直接求得时可以先求b对p的逆元x，再乘a去mod p，即

(a*x)%p和(a/b)%p结果是一样的

逆元求法

ll inverse(ll num,ll mod)  
{  
    ll x,y;  
    exgcd(num,mod,x,y);  
    while(x<0) x+=mod,y-=num;  
    return x;  
} 

HDU 2699 Romantic

看懂上面的就可以直接A裸题了，注意用__int64

#include <cstdio>
#include <cmath>
#include <cstring>
#include <string>
#include <set>
#include <map>
#include <stack>
#include <queue>
#include <vector>
#include <iostream>
#include <iomanip>
#include <algorithm>
using namespace std;
#define ll __int64
#define INF 0x7FFFFFFF
#define INT_MIN -(1<<31)
#define eps 10^(-6)
#define Q_CIN ios::sync_with_stdio(false);
#define REP( i , n ) for ( int i = 0 ; i < n ; ++ i )
#define FOR( i , a , b ) for ( int i = a ; i <= b ; ++ i )
#define CLR( a , x ) memset ( a , x , sizeof (a) );
#define RE freopen("1.in","r",stdin);
#define WE freopen("1.out","w",stdout);
#define MOD 10009
#define debug(x) cout<<#x<<":"<<(x)<<endl;

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

int main()
{
    //RE
    ll a,b;
    while(cin >> a>>b)
    {
        ll x,y;
        ll gcd=exgcd(a,b,x,y);
        x/=gcd;
        y/=gcd;
        if(1%gcd){cout<<"sorry"<<endl;continue;}
        while(x<0)
        {
            x+=b;y-=a;
        }
        cout<<x<<" "<<y<<endl;
    }
    return 0;
}

HDU 1576 A/B

已知n=A%9973，我一直在化n=...，但其实是要A=k*9973+n;

又设A/B=x，则A=B*x;

所以B*x=9973*k+n，即Bx-9973k=n，套公式求出x就得到A了，再取mod

#include <cstdio>
#include <cmath>
#include <cstring>
#include <string>
#include <set>
#include <map>
#include <stack>
#include <queue>
#include <vector>
#include <iostream>
#include <iomanip>
#include <algorithm>
using namespace std;
#define ll __int64
#define INF 0x7FFFFFFF
#define INT_MIN -(1<<31)
#define eps 10^(-6)
#define Q_CIN ios::sync_with_stdio(false);
#define REP( i , n ) for ( int i = 0 ; i < n ; ++ i )
#define FOR( i , a , b ) for ( int i = a ; i <= b ; ++ i )
#define CLR( a , x ) memset ( a , x , sizeof (a) );
#define RE freopen("1.in","r",stdin);
#define WE freopen("1.out","w",stdout);
#define MOD 10009
#define debug(x) cout<<#x<<":"<<(x)<<endl;

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

int main()
{
    int t;
    ll a,b,n;
    cin>>t;
    while(t--)
    {
        cin>>n>>b;
        ll x,y;
        ll gcd=exgcd(b,-9973,x,y);
        x=x/gcd*n;
        y=y/gcd*n;
        while(x<0)  x+=9973;
        cout<<x%9973<<endl;
    }
    return 0;
}

POJ 1061 青蛙的约会

#include<cstdio>
#include<cmath>
#include<cstring>
#include<set>
#include<stack>
#include<queue>
#include<vector>
#include<iostream>
#include<algorithm>
using namespace std;
#define ll long long
#define INF 0x7FFFFFFF
#define INT_MIN -(1<<31)
#define eps 10^(-6)
#define Q_CIN ios::sync_with_stdio(false)
#define REP( i , n ) for ( int i = 0 ; i < n ; ++ i )
#define REV( i , n ) for ( int i = n - 1 ; i >= 0 ; -- i )
#define FOR( i , a , b ) for ( int i = a ; i <= b ; ++ i )
#define FOV( i , a , b ) for ( int i = a ; i >= b ; -- i )
#define CLR( a , x ) memset ( a , x , sizeof (a) )
#define RE freopen("1.in","r",stdin);
#define WE freopen("1.out","w",stdout);
#define NMAX 10002
#define debug(x) cout<<#x<<":"<<(x)<<endl;
ll exgcd(ll a, ll b, ll &x, ll &y)
{
    ll t, d;
    if( b == 0)
    {
        x = 1;
        y = 0;
        return a;
    }
    d = exgcd(b, a % b, x, y);
    t = x;
    x = y;
    y = t - (a / b) * y;
    return d;
}

int main()
{
    ll x, y, t, m, n, l, a, b, d,  c;
    ll x0, y0;
    //RE;
    while(cin >> x >> y >> m >> n >> l)
    {
        a = n - m;
        b = l;
        c = x - y;
        ll gcd = exgcd(a, b, x0, y0);  //exgcd()返回gcd(a,b)
        if(c % gcd)
        {
            cout << "Impossible" << endl;
            continue;
        }

        x0=x0*(c/gcd);          //将ax+by=gcd的解转为ax+by=c
        y0=y0*(c/gcd);
//  通解为x=x0+b/gcd*t(t为常数)
//  求x最小，当t与x0正负相反时x最小，假设x0为正，则t为负数，则转为x=x0-b/gcd*t_new(t_new为正数)
//  令x=0得t_new=x0*gcd/b，最小的x=x0-b/gcd*t_new，【通项x=x0+b/d*t】防止负数就再加次b/gcd
        t = x0 * gcd / b;       
        x0 = x0 - t*b /gcd ;

        if(x0 < 0)
            x0 += b / gcd;
        cout << x0 << endl;
    }
    return 0;
}

POJ 2115 C Looooops

本题没几个结论还真不会做...

假设循环次数为x，则有方程x= [ (b-a+2^k)mod(2^k) ] /c，加2^k次方是解决b-a<0的情况，跟求循环队列的长度一样道理

化简得cx=(b-a)mod(2^k)

好，下面就是搬算法导论的结论了

    ax = b(mod n)  等价于ax-ny=b 是模线性方程

定理一：设d = gcd(a, n), 假定对整数x', y', 有d = ax' + ny'， 如果d | b, 则方程ax = b(mod n)有一个解的值为x0, 满足：

　　　　　　x0 = x'(b / d)(mod n）

定理二：假设方程ax = b(mod n)有解， x0是方程的任意一个解， 则方程对模n恰有d个不同的解， 分别为： xi = x0 + i * (n / d), 其中 i = 1,2,3......d - 1

变量代入可知

我们求一次cx+(2^k)*y=gcd(c,2^k)的解就得出定理一中的x'，进而求出x0（一个解），然后求最小的非负数解x就是答案了。如果b-a不是gcd(c,2^k)的倍数就是无解（定理一那里d|b:b%d==0)

这里看别人博客都说有最小正数解的方程，表示我是没管那么多，有个特解了再求最小非负解还会难？直接while上

#include <cstdio>
#include <cmath>
#include <cstring>
#include <string>
#include <set>
#include <map>
#include <stack>
#include <queue>
#include <vector>
#include <iostream>
#include <iomanip>
#include <algorithm>
using namespace std;
#define ll __int64
#define INF 0x7FFFFFFF
#define INT_MIN -(1<<31)
#define eps 10^(-6)
#define Q_CIN ios::sync_with_stdio(false);
#define REP( i , n ) for ( int i = 0 ; i < n ; ++ i )
#define FOR( i , a , b ) for ( int i = a ; i <= b ; ++ i )
#define CLR( a , x ) memset ( a , x , sizeof (a) );
#define RE freopen("1.in","r",stdin);
#define WE freopen("1.out","w",stdout);
#define MOD 10009
#define debug(x) cout<<#x<<":"<<(x)<<endl;

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

int main()
{
//    RE
    int t;
    ll a,b,n,c,k;
    while(cin>>a>>b>>c>>k,a||b||c||k)
    {
        ll x,y;
        ll gcd=exgcd(c,1LL<<k,x,y);
        if((b-a)%gcd)
        {
            cout<<"FOREVER"<<endl;
            continue;
        }
        ll x1=x*((b-a)/gcd)%(1LL<<k);
        while(x1>0)       //因为上面求出的x1不一定是最下正解，所以转求得最大负解，下面再求的第一个正解就是答案
            x1=x1-(1LL<<k)/gcd;
        while(x1<0)
            x1=x1+(1LL<<k)/gcd;
        cout<<x1<<endl;
    }
    return 0;
}

顺便贴个gcd的水题

HDU 2104 hide handkerchief

大意是n个人围圈，每次跳过m-1个人，看能否找到。

这不就是哈希表的二次（平方）探测法 存数时一定能找到位置的情况么，对应题目就是m,n互质就能找到。

#include <cstdio>
#include <cmath>
#include <cstring>
#include <string>
#include <set>
#include <map>
#include <stack>
#include <queue>
#include <vector>
#include <iostream>
#include <iomanip>
#include <algorithm>
using namespace std;
#define ll __int64
#define INF 0x7FFFFFFF
#define INT_MIN -(1<<31)
#define eps 10^(-6)
#define Q_CIN ios::sync_with_stdio(false);
#define REP( i , n ) for ( int i = 0 ; i < n ; ++ i )
#define FOR( i , a , b ) for ( int i = a ; i <= b ; ++ i )
#define CLR( a , x ) memset ( a , x , sizeof (a) );
#define RE freopen("1.in","r",stdin);
#define WE freopen("1.out","w",stdout);
#define MOD 10009
#define debug(x) cout<<#x<<":"<<(x)<<endl;

int gcd(int m,int n)
{
    return n?gcd(n,m%n):m;
}

int main()
{
//    RE
    int n,m;
    while(cin>>n>>m)
    {
        if(n==-1&&m==-1)
            break;
        if(gcd(n,m)==1)
            cout<<"YES"<<endl;
        else
            cout<<"POOR Haha"<<endl;
    }
    return 0;
}

POJ 2891 Strange Way to Express Integers

给n条x mod ai= ri，求最小正解或无解。即中国剩余定理

已知

x ≡ r1 (mod a1)

x ≡ r2 (mod a2)

两条式的情况下对于方程①解出x0后，所有解为x0+k*a1，则方程②的解必须是x0+k*a1里的一个

则等价于a1*k1+r1==a2*k2+r2，即a1k1-a2k2=r2-r1，扩展gcd解出k1,那么原x=r1+k1*a1。

求解过程中不断让r1=r1+k1*a1，a1为LCM(a1,...ai-1)。

#include <cstdio>
#include <cmath>
#include <cstring>
#include <string>
#include <set>
#include <map>
#include <stack>
#include <queue>
#include <vector>
#include <iostream>
#include <iomanip>
#include <algorithm>
using namespace std;
#define ll __int64
#define INF 0x7FFFFFFF
#define INT_MIN -(1<<31)
#define eps 10^(-6)
#define Q_CIN ios::sync_with_stdio(false);
#define REP( i , n ) for ( int i = 0 ; i < n ; ++ i )
#define FOR( i , a , b ) for ( int i = a ; i <= b ; ++ i )
#define CLR( a , x ) memset ( a , x , sizeof (a) );
#define RE freopen("1.in","r",stdin);
#define WE freopen("1.out","w",stdout);
#define MOD 10009
#define debug(x) cout<<#x<<":"<<(x)<<endl;

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

int main()
{
//    RE
    ll n,c,a1,r1,a2,r2;
    while(cin>>n)
    {
        cin>>a1>>r1;
        ll Lcm=a1;
        ll x,y;
        ll no=0;
        if(n==1&&r1>a1)    r1=a1+r1;   //POJ缺少n为1的数据，如3,4与4,3，有的话估计很多人都得WA
        for(ll i=1;i<n;i++)
        {
            cin>>a2>>r2;
            if(no)  continue;
            c=r2-r1;
            ll d=exgcd(a1,a2,x,y);  //xa1-ya2=r2-r1
            if(c%d)
            {
                no=1;
                continue;
            }
            ll s=a2/d;
//
//debug(c)debug(d)debug(x)
//debug(c/d*x)
//debug(x/d*c)
//debug(c*x/d)

            x=c/d*x;               //ax+by=c通解，血的经历c/d*x与x/d*c不等价
            x=(x%s+s)%s;          //最小正解
            r1=r1+x*a1;             //迭代
            a1=a1/d*a2;             //更新为LCM
        }
        if(!no)
            cout<<r1<<endl;
        else
            cout<<"-1\n";
    }
    return 0;
}