欢迎您访问程序员文章站本站旨在为大家提供分享程序员计算机编程知识!
您现在的位置是: 首页

连分数及Pell方程的解法

程序员文章站 2022-03-13 19:17:29
...
  本文代码是我为了解决Project Euler上的问题而写的数学工具,之前的见:
    按字典顺序生成所有的排列
    筛法求素数


  所谓一个实数的连分数表示,是指将一个实数x写成以下形式:
    连分数及Pell方程的解法
            
    
    博客分类: 数学&编程 ScalaPHP.netBlog 
  其中a0,a1,...,b1,b2,..都是自然数。
  当其中b1,b2,..都取为1时,我们称之为简单连分数表示(Simple Continued Fraction)
  可以证明:每一个不是完全平方数的自然数N,其平方根可以写成简单连分数表示,并且其中a1,a2,..呈周期出现。
  比如23的平方根的连分数表示为:
    连分数及Pell方程的解法
            
    
    博客分类: 数学&编程 ScalaPHP.netBlog 
  并且其中[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)
}