详解几个方程的解法

一:求解不定方程 ax+by=c

  先给出两个介绍扩展欧几里德很详细的网站:

  http://www.cnblogs.com/frog112111/archive/2012/08/19/2646012.html

  http://www.cnblogs.com/ka200812/archive/2011/09/02/2164404.html

扩展欧几里德的递归代码:

//扩展欧几里德的递归代码:
int exgcd(int a,int b,int &x,int &y)
{
    if(b==0)
    {
    //b==0,也就是说gcd(a,0)==a。
    //原式变为ax+by==a --> x==1,y==0。
        x=1;
        y=0;
        return a;
    }
    int d=exgcd(b,a%b,x,y);
    int t=x;
    x=y;   //更新当前的x
    y=t-a/b*y;   //更新当前的y
    return d; //返回最大公约数
}
View Code

  扩展欧几里德的应用:

  (一):求解不定方程:

  利用扩展欧几里德,可以求出

  ax+by=gcd(a,b)的一组解,x1、y1

  则通解为:

  x = x1+b/gcd*i

  y = y1-a/gcd*i

  那么对于一般的不定式ax+by=c,

  若gcd(a,b)整除c,则方程有解

  一组解为:x0=x1*(c/gcd(a,b)),y0=y1*(c/gcd(a,b))。

  通解为:x=x0+b/gcd*i,y=y0-a/gcd*i。

  否则无解。

  若ax+by=c有解,且gcd(a,b)=d,x的特解为x*(c/d)。

  令ans=x*(c/d),s=b/d;

  则方程的x最小整数解为:(ans%s+s)%s;

  (二):求解同余方程

  ax+by=c,实际上它等价于同余方程,ax =c mod b,这样同余方程就可以用扩展欧几里德算法求解。

  (三):求解模的逆元

  另外如果a b互质那么ax+by=gcd(a,b) = 1,转换为同余方程ax = 1 mod b这样x就变成了a的关于mod b的逆元。

  求解扩展欧几里德过后,只需x=(x%b+b)%b,求出最小正整数解x,即为a的逆元。

  也就是说求逆元也不过是扩展欧几里德算法的一个特殊情况。

  两个简单的求解不定方程的例子:

  POJ 2115  C Looooops

#include <iostream>
#include <stdio.h>
using namespace std;
long long A,B,C;
int k;
long long exgcd(long long a,long long b,long long &x,long long &y) {
    if(b==0) {
        x=1;
        y=0;
        return a;
    }
    long long d=exgcd(b,a%b,x,y);
    long long tmp=x;
    x=y;
    y=tmp-(a/b)*y;
    return d;
}
int main() {
    while(scanf("%I64d%I64d%I64d%d",&A,&B,&C,&k)!=EOF) {
        if(A==0 && B==0 && C==0 && k==0)
            break;
        long long a,b,c,x,y,d,s;
        a=C;
        b=(long long)1<<k;
        c=B-A;
        d=exgcd(a,b,x,y);
        if(c%d==0) {
            x=x*c/d;
            s=b/d;
            x=(x%s+s)%s;
            printf("%I64d
",x);
        } else
            printf("FOREVER
");

    }
    return 0;
}
View Code

  POJ 2142  The Balance

http://www.cnblogs.com/chenxiwenruo/p/3557315.html

二:解方程x^2+y^2=z^2

  设不定方程:x^2+y^2=z^2
  若正整数三元组(x,y,z)满足上述方程,则称为毕达哥拉斯三元组。
  若gcd(x,y,z)=1,则称为本原的毕达哥拉斯三元组。

  定理:
  正整数x,y,z构成一个本原的毕达哥拉斯三元组且y为偶数,当且仅当存在互素的正整数m,n(m>n),其中m,n的奇偶性不同,
  并且满足
    x=m^2-n^2,y=2*m*n, z=m^2+n^2

  例题:

  POJ 1305 Fermat vs. Pythagoras

  http://www.cnblogs.com/chenxiwenruo/p/3558297.html

三:佩尔方程x^2-dy^2=1

  形如x^2-dy^2=1 (d>1且d不为完全平方数)的不定方程称为佩尔方程。

   若d是完全平方数,那么x^2-dy^2=1是无解的。

  

 

  求解佩尔方程,可以用暴力求解法。

  还有一种更高效的解法,连分数法,有兴趣的可以自己看看。

  暴力解法代码:

void solve(int d){
    y=1;
    while(1){
        x=(long long)sqrt(double(d*y*y+1));
        if(x*x-d*y*y==1){
            break;
        }
        y++;
    }

}
View Code

  给出两个简单例题:

  POJ    1320  Street Numbers

  求解两个不相等的正整数n和m (n<M),使得1+2+3+…+n=n+(n+1)+…+m。

  思路:实际上求解(2m+1)^2-8n^2=1。令x=2m+1,y=n。可以利用迭代公式求解出前十组的(n,m)。

 

  HDU 3292  No more tricks, Mr Nanguo

  求解X^2-NY^2=1的第K大的满足式子的X,如果无解则输出相应的句子。利用矩阵迭代+快速幂即可求出。

#include <iostream>
#include <cstdio>
#include <string.h>
#include <algorithm>
#include <math.h>
using namespace std;
const int mod=8191;
int n,k;
int vis[30];
long long x,y;
void solve(int d){
    y=1;
    while(1){
        x=(long long)sqrt(double(d*y*y+1));
        if(x*x-d*y*y==1){
            break;
        }
        y++;
    }

}
struct Matrix{
    long long m[2][2];
    void init(long long x,long long y){
        m[0][0]=x;m[0][1]=n*y%mod;
        m[1][0]=y;m[1][1]=x;
    }
    void initE(){
        memset(m,0,sizeof(m));
        m[0][0]=m[1][1]=1;
    }
}A;
Matrix operator*(Matrix a,Matrix b){
    Matrix c;
    for(int i=0;i<2;i++){
        for(int j=0;j<2;j++){
            c.m[i][j]=0;
            for(int k=0;k<2;k++){
                c.m[i][j]=(c.m[i][j]+a.m[i][k]*b.m[k][j]%mod)%mod;
            }
        }
    }
    return c;
}
Matrix quickPow(Matrix a,int b){
    Matrix ret;
    ret.initE();
    while(b){
        if(b&1)
            ret=ret*a;
        a=a*a;
        b=b>>1;
    }
    return ret;
}

int main()
{
    memset(vis,0,sizeof(vis));
    for(int i=1;i*i<30;i++){
        vis[i*i]=1;  //若是完全平方数,则无解,标记为1
    }
    while(scanf("%d%d",&n,&k)!=EOF){
        if(!vis[n]){
            solve(n);
            x=x%mod;y=y%mod;
            A.init(x,y);
            Matrix B;
            B=quickPow(A,k-1);
            long long xk,yk;
            xk=(B.m[0][0]*x%mod+B.m[0][1]*y%mod)%mod;
            printf("%I64d
",xk);
        }
        else{
            printf("No answers can meet such conditions
");
        }

    }
    return 0;
}
View Code

四:解高次同余方程A^x=B(mod C)

  Baby Step Giant Step算法:复杂度O( sqrt(C) )

  请戳:

  http://www.cnblogs.com/chenxiwenruo/p/3554885.html

原文地址:https://www.cnblogs.com/chenxiwenruo/p/3639340.html