连分数法解佩尔方程特解

转载自https://blog.csdn.net/wh2124335/article/details/8871535?locationNum=14&fps=1

一、佩尔方程的形式：

$$x^2-Dy^2=1,  D为正整数$$

二、关于佩尔方程的特解

特解是指佩尔方程的最小整数解，容易发现当x最小的时候y也同样达到最小。

在一般情况下，佩尔方程的特解是通过暴利枚举方法得到的，本文将介绍如何用应用连分数法求特解。

本文将不涉及证明，只介绍方法。

三、连分数法

一个实数的简单连分数表示，是指将一个实数用以下方法表示：

$$x = a_0+frac{1}{a_1 + frac{1}{a_2+frac{1}{a_3+...}}}$$

可以把连分数简记为：$x = [a_0;a_1, a_2, a_3...]$.

有理数的连分数有两种表示形式：

$$frac{p}{q} = [a_0;a_1, a_2, a_3, ...,a_n,1] 或 [a_0;a_1, a_2, a_3, ..., a_n+1]$$

所有无限连分数都是无理数，而所有无理数都可以用一种精确的方式表示成无限连分数，可以用这种方法逼近，无理数的值。

三、关于一个非完全平方数的平方根的连分数表示

可以证明：一个非完全平方数的平方根是以周期性呈现的：

比如：

简写为：

$$sqrt 22 = [4;1,2,4,2,1,8]$$

在之后就会循环出现1,2,4,2,1,8

我们不妨这样记这种连分数的形式：

$$sqrt22 = [4;<1,2,4,2,1,8>]$$

显然循环节的长度是6

并且还有个重要的特点：这个循环节一定是 $a_1$ 开始，且最后一个数 $a_n$ 一定是 $a_0$ 的2倍。

五、且佩尔方程的最小特解：

我们将 $sqrt D$ 写成连分数的形式：（相当于用连分数无线逼近平方根）

 $$sqrt D =  [a_0;<a_1, a_2, ..., a_{n-1},2a_0>]$$

并且我们记：

$$frac{p}{q} = [a_0; a_1, a_2, ..., a_{n-1}]$$

（关于计算p,q：只要按照连分数的展开形式，迭代计算即可）

其中如果记循环节长度为s

那么有如下结论：

1、如果s为偶数时。最小特解为：

$$x_0 = p, y_0 = q$$

2、如果s为奇数时，最小特解为：

$$x_0 = 2p^2+1, y_0 = 2pq$$

六、计算 $sqrt D$ 的连分数

我们希望得到准确的连分数展开，那么关键在于不用浮点型计算。

接下来以 $sqrt {22}$ 为例：

 按照这种方式，我们计算出了的连分数：$sqrt {22} = [4;<1,2,4,2,1,8>]$.

然后可以计算出来：

$$frac{197}{42} = [4;1,2,4,2,1]$$

由于循环节长度6是偶数，那么佩尔方程 $x^2 -22y^2 = 1$ 的最小特解是：$x_0 = 197, y_0 = 42$.

之后我们参照上面的例子，很容易设计出计算连分数的算法。

代码

结果很大1000以内好多结果都超long long了。。。要改成大数才行。。。

#include<cstdio>
#include<cstring>
#include<cmath>
#include<cstdlib>
using namespace std;
typedef long long ll;
ll a[20000];
bool pell_minimum_solution(ll n,ll &x0,ll &y0){
    ll m=(ll)sqrt((double)n);
    double sq=sqrt(n);
    int i=0;
    if(m*m==n)return false;//当n是完全平方数则佩尔方程无解
    a[i++]=m;
    ll b=m,c=1;
    double tmp;
    do{
        c=(n-b*b)/c;
        tmp=(sq+b)/c;
        a[i++]=(ll)(floor(tmp));
        b=a[i-1]*c-b;
        //printf("%lld %lld %lld
",a[i-1],b,c);
    }while(a[i-1]!=2*a[0]);
    ll p=1,q=0;
    for(int j=i-2;j>=0;j--){
        ll t=p;
        p=q+p*a[j];
        q=t;
        //printf("a[%d]=%lld %lld %lld
",j,a[j],p,q);
    }
    if((i-1)%2==0){x0=p;y0=q;}
    else{x0=2*p*p+1;y0=2*p*q;}
    return true;
}
 
int main(){
    ll n,x,y;
    while(~scanf("%lld",&n)){
        if(pell_minimum_solution(n,x,y)){
            printf("%lld^2-%lld*%lld^2=1	",x,n,y);
            printf("%lld-%lld=1
",x*x,n*y*y);
        }
    }

一个java打表的代码：

import java.io.*;
import java.math.*;
import java.util.*;
 
public class main {
    static long [] a = new long [1000];
    static BigInteger x;
    static BigInteger y;
    public static void main(String [] args) throws FileNotFoundException{
        FileReader fin = new FileReader ("data.in");
        File fout = new File("data.out");
        PrintStream pw = new PrintStream(fout);
        Scanner cin = new Scanner(fin);
        System.setOut(pw);
        while(cin.hasNext()){
            long n=cin.nextLong();
            if(pell_solution(n)){
                
                System.out.print("""+x+"",");
            }else{
                System.out.print(""no solution",");
            }
        }
    }
    public static boolean pell_solution(long D){
        double sq=Math.sqrt((double)D);
        long m=(long) Math.floor(sq);
        int i=0;
        if(m*m==D)return false;
        a[i++]=m;
        long b=m,c=1;
        double tmp;
        do{
            c=(D-b*b)/c;
            tmp=(sq+b)/c;
            a[i++]=(long)(Math.floor(tmp));
            b=a[i-1]*c-b;
        }while(a[i-1]!=2*a[0]);
        BigInteger p=new BigInteger("1");
        BigInteger q=new BigInteger("0");
        BigInteger t;
        for(int j=i-2;j>=0;j--){
            t=p;
            p=q.add(p.multiply(BigInteger.valueOf(a[j])));
            q=t;
        }
        if((i-1)%2==0){
            x=p;y=q;
        }else{
            x=BigInteger.valueOf(2).multiply(p).multiply(p).add(BigInteger.ONE);
            y=BigInteger.valueOf(2).multiply(p).multiply(q);
        }
        return true;
    }
}