POJ 2427 Smith's Problem Pell方程

题目链接 :  http://poj.org/problem?id=2427

PELL方程几个学习的网址: 

http://mathworld.wolfram.com/PellEquation.html     wolfram的讲解

http://hi.baidu.com/aekdycoin/item/a45f7c37850e5b9db80c03d1     AC神的博客 

http://blog.csdn.net/acdreamers/article/details/8529686    acdreamer的博客  (从这里知道的思路...

Pell方程 :  形如 X² - D*Y² = 1 的式子我们称作Pell方程 (D为正整数)

Pell方程的推广形式 :  形如A*X² - B*Y² = C 的式子我们称作Pell方程的推广 (其中 A,B,C均为正整数)   

本题是Pell方程的最小根

按照Pell方程连分数解法的定义 , 只需要求出sqrt(N)的连分数即可 

于是我苦翻了一天数论书看懂了连分数的性质...公式在初等数论及其应用 P370

所以我们只需要一直求连分数sqrt(N)的 收敛子p/q p,q就是最后我们要求的答案

但是这题不需要化成连分数的向量形式即 [a1;a2,a3...] , 因为double精度误差很大,BigDecimal很不方便 而且 找循环节很复杂

这题只需要用下面的定理即可 ...本题中 P0 = 0,Q0 = 1  ,这里因为求出Pk,Qk会非常大 , 所以用Java的BigInteger实现更方便

import java.util.*;
import java.io.BufferedInputStream;
import java.math.*;
public class Main {
    public static void main(String[] args) {
        Scanner cin = new Scanner(new BufferedInputStream(System.in)); 
        while(cin.hasNext()){
            int n = cin.nextInt();
            double p = Math.sqrt(n);
            int k = (int)p;
            if(k == p) {
                System.out.println("No solution!");
                continue;
            }else {
                /*
                 *  逐项求 sqrt(D) 的连分数 用p/q表示 
                 *  公式在初等数论及其应用P370
                 */
                BigInteger x = BigInteger.ONE; //p
                BigInteger y = BigInteger.ONE; //q 
                BigInteger a,N,P1,Q1,P2,Q2,ak,p1,q1,p2,q2;
                q1 = p2 = P1 = BigInteger.ZERO;
                p1 = q2 = Q1 = BigInteger.ONE;
                N = BigInteger.valueOf(n); // N = n;
                a = BigInteger.valueOf(k); // a = [sqrt(n)]
                ak = a; //ak
                while(!x.multiply(x).subtract(N.multiply(y).multiply(y)).equals(BigInteger.ONE)){
                    x = ak.multiply(p1).add(p2); // p[k] = ak * p[k-1] + p[k-2]
                    y = ak.multiply(q1).add(q2); // q[k] = ak * q[k-1] + q[k-2]
                    P2 = ak.multiply(Q1).subtract(P1);  //P2=P[k+1];
                    Q2 = N.subtract(P2.multiply(P2)).divide(Q1); //Q2=Q[k+1];
                    ak = P2.add(a).divide(Q2); //ak 
                    
                    P1 = P2;  
                    Q1 = Q2;  
                    
                    p2 = p1; 
                    p1 = x; 
                    q2 = q1;
                    q1 = y;
                }
                System.out.println(x+" "+y);
            }
        }
    }
}