连分数及Pell方程的解法

　　本文代码是我为了解决[url=http://projecteuler.net/]Project Euler[/url]上的问题而写的数学工具，之前的见：
　　　　[url=http://eastsun.iteye.com/blog/207614]按字典顺序生成所有的排列[/url]
　　　　[url=http://eastsun.iteye.com/blog/240339]筛法求素数[/url]


　　所谓一个实数的[color=olive]连分数表示[/color]，是指将一个实数x写成以下形式：
　　　　[img]http://www.iteye.com/upload/picture/pic/22501/26578e33-e037-3e65-b4e3-888e38b5e83a.gif?1222869669[/img]
　　其中a0,a1,...,b1,b2,..都是自然数。
　　当其中b1,b2,..都取为1时，我们称之为[color=olive]简单连分数表示[/color](Simple Continued Fraction)
　　可以证明：每一个不是完全平方数的自然数N，其平方根可以写成[color=olive]简单连分数表示[/color]，并且其中a1,a2,..呈周期出现。
　　比如23的平方根的连分数表示为：
　　　　[img]http://www.iteye.com/upload/picture/pic/22503/923d2bb2-28c4-3e8b-99e2-141a43f7a12e.jpg?1222871288[/img]
　　并且其中[color=green][1,3,1,8][/color]就是其一个周期，就是接下来的表示是由[1,3,1,8]循环出现。
　　下面是得到一个自然数平方根的简单连分数表示的Scala代码：

/**
   &#Util.scala
   utils for mathematical algorithm,include:
   # get all primes below bound in order
   # generate all permutations in lexicographical order 
   # get simple continued fraction representation of the sqrt of n
   @author Eastsun
*/
package eastsun.math

object Util {
    /**
      Get simple continued fraction representation of the sqrt of n
    */
    def continuedFractionOfSqrt(n: Int,buf: Array[Int]):Int = {
        val sq = Math.sqrt(n)
        var (p,q) = (sq,n - sq*sq)
        if(q == 0) 0
        else{
            var idx = 0
            var an  = 0
            do {
                an = (sq + p)/q
                if(buf != null) buf(idx) = an
                idx += 1
                p = an*q - p
                q = (n - p*p)/q
            }while(an != 2*sq)
            idx
        }
    }
}

　　简单解释一下函数[color=violet]def continuedFractionOfSqrt(n: Int,buf: Array[Int]):Int[/color]的功能：该函数有两个参数，其中n表示要求其平方根连分数表示的自然数n；buf用来保存其连分数表示中a1,a2,...的一个周期(注意，没有包括a0)，该参数可以为null。函数返回一个Int，表示a1,a2,..周期的大小，也就是buf中保存数据的长度。
　　下面是对使用该函数的一个演示：
[quote]scala> var buf = new Array[Int](4)
buf: Array[Int] = Array(0, 0, 0, 0)

scala> continuedFractionOfSqrt(23,buf)
res8: Int = 4

scala> buf.mkString(",")
res9: String = 1,3,1,8

scala>[/quote]

　　OK，现在我们可以来看看Project Euler上的[url=http://projecteuler.net/index.php?section=problems&id=64]第64题[/url]了：
[quote]
64 How many continued fractions for N ≤ 10000 have an odd period?[/quote]
　　题目很简单：[color=blue]求10000以内的自然数中，平方根的连分数表示的周期长度为奇数的有多少个。[/color]
　　下面是该题的Scala解法(使用了上面的函数)：

import eastsun.math.Util._
object Euler064 extends Application {

    val res = 1.to(10000).filter{ continuedFractionOfSqrt(_,null) % 2 ==1 }.length
    println(res)
}

　　在将Project Euler66题前，先介绍一个数学名词：[url=http://mathworld.wolfram.com/PellEquation.html]佩尔方程[/url]：形如 x^2 - D×y^2 = 1的不定方程称为佩尔方程。其中D为非完全平方数的自然数。并且称其所有正整数解(x,y)中使得x最小的那个解为[color=olive]最小解[/color]。
　　佩尔方程求解与平方根的连分数表示有着很大的关联，这里我就不细说了，对数学细节干兴趣的可以参考Math World上的[url=http://mathworld.wolfram.com/PellEquation.html]Pell Equation[/url]。下面我直接给出[url=http://projecteuler.net/index.php?section=problems&id=66]Project Euler66题[/url]的叙述及其Scala代码：
[quote]Find the value of D ≤1000 in minimal solutions of x for which the largest value of x is obtained.[/quote]
　　[color=blue]就是说对D≤1000，求D使得佩尔方程x^2 - D×y^2 = 1的最小解中x的值最大。[/color]
　　下面上代码：

import eastsun.math.Util._
object Euler066 extends Application {
    val buf = new Array[Int](1000)
    var (res,max,d) = (2,3:BigInt,1)
    while(d <= 1000){
        val pd = continuedFractionOfSqrt(d,buf)
        if(pd > 0){
            val sq = Math.sqrt(d)
            var (x0,y0) = (sq:BigInt,1:BigInt)
            var (x1,y1) = ((buf(0)*sq+1):BigInt,buf(0):BigInt)
            val cnt = if(pd%2 == 1) 2*pd-1 else pd-1
            var idx = 1
            while(idx < cnt){
                var t = x1
                var a = buf(idx%pd)
                x1 = x1*a + x0
                x0 = t
                t  = y1
                y1 = y1*a + y0
                y0 = t
                idx += 1
            }
            if(x1 > max){
                max = x1
                res = d
            }
        }
        d += 1
    }
    println(res)
}