轉載自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\n",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\n",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\t",x,n,y); printf("%lld-%lld=1\n",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; } }