仅以此文膜拜八年前的本人

序言

欧拉打算（Project Euler）就像 LeetCode，是一个编程答题的网站。不同于 LeetCode 的是，欧拉打算只要求用户提交最终答案即可（个别是一个数字），而不须要残缺代码。因而，能够纵情地应用本人喜爱的编程语言——不少题目甚至光靠笔和纸便能解决。

欧拉打算的第 66 题十分有意思，它的题目很简略，就是要求找出在不大于 1000 的整数中，以哪一个数字为丢番图方程的系数，能够失去所有最小解中的最大值。

能够很容易地看出方程有一个直观的暴力算法：让 y 从 1 开始递增，对于每一个 y，计算公式 Dy^2+1 的值。如果该值为平方数，那么它的平方根就是最小的 x 解。再按照这个算法求解所有 D 不大于 1000 的方程，便能够求出题目标答案。很容易用 Python 写出这个算法

# -*- coding: utf8 -*-
import math

def is_square(num: int) -> bool:
    return math.isqrt(num) ** 2 == num

def find_x(D: int) -> int:
    """求出给定 D 时，满足题目所给的丢番图方程的最小的 x。"""
    assert not is_square(D)
    y = 1
    while True:
        candidate = D * y * y + 1
        if is_square(candidate):
            return math.isqrt(candidate)
        y += 1

def solve_66(limit):
    """找出不大于 limi 的 D 中，使 find_x 的返回值最大的那一个数字。"""
    max_D = None
    max_x = None
    D = 2
    while D <= limit:
        if is_square(D):
            D += 1
            continue
        x = find_x(D)
        if max_x is None or x > max_x:
            max_D = D
            max_x = x
        D += 1
    return max_D, max_x

if __name__ == '__main__':
    D, x = solve_66(7)
    print('D is {} and x is {}'.format(D, x))

但如果将下限 limit 晋升为 1000，这个算法在有生之年是算不出后果的。

要想解决这一题，须要借助数学的力量。

佩尔方程

八年前第一次做这一题的时候，通过一番搜寻，我从这篇文章中晓得了题目中的方程叫做佩尔方程。它有规范的解法，但须要用到连分数。那么什么是连分数呢？

连分数不是一种新的数系，只是小数的另一种写法。例如能够把分数 45 除以 16 写成上面的模式

就像定义递归的数据结构一样，能够给连分数一个递归的定义。连分数要么是一个整数，要么是一个整数加上另一个连分数的倒数。除了下面的模式，连分数也能够写成更节俭篇幅的样子。比方把 45 除以 16 写成[2;1,4,3]，即把本来的式子中所有的整数局部按程序写在一对方括号之间。这种记法，看起来就像是编程语言中的数组个别。

如果用数组 [2;1,4,3] 的不同前缀来结构分式，那么后果顺次为2/1、3/1、14/5。它们是这个连分数的渐进连分数，而佩尔方程的一组解，就来自于渐进连分数的分子和分母。

以系数为 7 的佩尔方程为例，先计算出根号 7 的连分数，而后顺次尝试它的渐进连分数。前三个别离为 2/1、3/1、5/2，都不是方程的解。第四个渐进连分数8/3 才是方程的解。如果持续进步连分数的精度，还会找到第二个解127/48。持续找，还有更多，而 8 则是其中最小的 x。

所以，想要疾速算出佩尔方程的解，最重要的是找到计算一个数的平方根的连分数的算法。

计算平方根的连分数的错误方法

要计算一个数字的连分数，最重要的便是要算出所有的整数局部（a0、a2、a2等）。它们都能够根据定义间接计算

推广到一半状况，如果用变量 n 存储开平方的数字，用 numbers 存储所有已知的整数，那么用 Python 能够写出上面的算法来计算出下一个整数

# 计算连分数数列的下一个数字
import math

def compute_next_integer_part(n, numbers):
    v = math.sqrt(n)
    for a in numbers:
        v = 1 / (v - a)
    return int(v)

if __name__ == '__main__':
    n = 14
    numbers = [3, 1, 2, 1]
    v = compute_next_integer_part(n, numbers)
    print('下一个数字为{}'.format(v))

遗憾的是，这个算法算进去的数字会因为计算上的精度误差而导致失之毫厘谬以千里。

计算平方根的连分数的正确办法

要想计算出正确的后果，就须要尽可能地打消在计算 1 / (v - a) 的时候引入的误差，因而必须把浮点数从分母中除去。

在这个网站中，作者以计算根号 14 的连分数为例，列出了一个表格

能够看到 x1、x2，以及x3 都是形如 (sqrt(n)+a)/b 这样的格局，这样的式子更利于管制误差。那么是否每一个待计算的 x 都合乎这种格局呢？答案是必定的，能够用数学归纳法予以证实（为了不便写公式，用 LaTeX 写好后截了图）

在这个证实过程中，还失去了分子中的 a 以及分母中的 b 的递推公式，当初能够写出正确的计算连分数整数局部的代码了。

用 Common Lisp 实现上述算法

为了在实现这个算法的同时还要写出优雅的代码，我会用上 Common Lisp 的面向对象个性。首先是定义一个类来示意一个能够一直进步精度的连分数

(defpackage #:com.liutos.cf
  (:use #:cl))

(in-package #:com.liutos.cf)

(defclass <cf> ()
  ((a
    :documentation "数学归纳法中、分子中与平方根相加的数"
    :initform 0)
   (b
    :documentation "数学归纳法中的分母"
    :initform 1)
   (numbers
    :documentation "连分数中的整数局部顺次组成的数组。"
    :initform nil)
   (origin
    :documentation "被开平方的数字"
    :initarg :origin))
  (:documentation "示意整数 ORIGIN 的平方根的连分数。"))

接着再定义这个类须要实现的“接口”

(defgeneric advance (cf)
  (:documentation "让连分数 CF 进步到下一个精度。"))

(defgeneric into-rational (cf)
  (:documentation "将连分数 CF 转换为有理数类型的值。"))

最初来实现上述两个接口

(defmethod advance ((cf <cf>))
  "依据递推公式计算出下一批 a、b，以及连分数的整数局部。"
  (let* ((a (slot-value cf 'a))
         (b (slot-value cf 'b))
         (n (slot-value cf 'origin))
         (m (truncate (+ (sqrt n) a) b)))
    (let ((a (- (* b m) a))
          (b (/ (- n (expt (- a (* b m)) 2)) b)))
      (setf (slot-value cf 'a) a
            (slot-value cf 'b) b
            (slot-value cf 'numbers) (append (slot-value cf'numbers) (list m))))
    (values)))

(defmethod into-rational ((cf <cf>))
  (let* ((numbers (reverse (slot-value cf 'numbers)))
         (v (first numbers)))
    (dolist (n (rest numbers))
      (setf v
            (+ n (/ 1 v))))
    v))

在实现 into-rational 办法上，Common Lisp 的有理数数值类型给我带来了极大的便当，它使我不用放心计算 (/ 1 v) 的时候会引入误差，代码写起来简略直白。

解题

乘胜追击，用 Common Lisp 解答第 66 题

(defun find-min-x (D)
  (let ((cf (make-instance '<cf> :origin D)))
    (loop
       (advance cf)
       (let* ((ratio (into-rational cf))
              (x (numerator ratio))
              (y (denominator ratio)))
         (when (= (- (* x x) (* D y y)) 1)
           (return-from find-min-x x))))))

(defun square-p (n)
  (let ((rt (sqrt n)))
    (= rt (truncate rt))))

(defun pro66 (&optional (bnd 1000))
  (let ((target-d)
    (max-x 0))
    (loop :for i :from 2 :to bnd
       :when (not (square-p i))
       :do (let ((x (find-min-x i)))
         (if (> x max-x)
         (setf target-d i
               max-x x))))
    (values target-d max-x)))

答案的 D 是多少就不说了，不过作为答案的 x 是 16421658242965910275055840472270471049。有趣味的读者能够试一下暴力解法要花多久能力算到这个数字。

全文完。

浏览原文