只有lisp能做——欧拉计划第66题

仅以此文膜拜八年前的自己

序言

欧拉计划（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。有兴趣的读者可以试一下暴力解法要花多久才能算到这个数字。

全文完。

阅读原文