连分数及Pell方程的解法

　　本文代码是我为了解决Project Euler上的问题而写的数学工具，之前的见：
　　　　按字典顺序生成所有的排列
　　　　筛法求素数


　　所谓一个实数的连分数表示，是指将一个实数x写成以下形式：
　　　　
　　其中a0,a1,...,b1,b2,..都是自然数。
　　当其中b1,b2,..都取为1时，我们称之为简单连分数表示(Simple Continued Fraction)
　　可以证明：每一个不是完全平方数的自然数N，其平方根可以写成简单连分数表示，并且其中a1,a2,..呈周期出现。
　　比如23的平方根的连分数表示为：
　　　　
　　并且其中[1,3,1,8]就是其一个周期，就是接下来的表示是由[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
        }
    }
}

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

引用

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>

　　OK，现在我们可以来看看Project Euler上的第64题了：

引用

64 How many continued fractions for N ≤ 10000 have an odd period?

　　题目很简单：求10000以内的自然数中，平方根的连分数表示的周期长度为奇数的有多少个。
　　下面是该题的Scala解法(使用了上面的函数)：

import eastsun.math.Util._
object Euler064 extends Application {
    
    val res = 1.to(10000).filter{ continuedFractionOfSqrt(_,null) % 2 ==1 }.length
    println(res)
}

　　在将Project Euler66题前，先介绍一个数学名词：佩尔方程：形如 x^2 - D×y^2 = 1的不定方程称为佩尔方程。其中D为非完全平方数的自然数。并且称其所有正整数解(x,y)中使得x最小的那个解为最小解。
　　佩尔方程求解与平方根的连分数表示有着很大的关联，这里我就不细说了，对数学细节干兴趣的可以参考Math World上的Pell Equation。下面我直接给出Project Euler66题的叙述及其Scala代码：

引用

Find the value of D  ≤1000 in minimal solutions of x for which the largest value of x is obtained.

　　就是说对D≤1000，求D使得佩尔方程x^2 - D×y^2 = 1的最小解中x的值最大。
　　下面上代码：

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)
}