仅以此文膜拜八年前的本人
序言
欧拉打算（Project Euler）就像LeetCode，是一个编程答题的网站。不同于LeetCode的是，欧拉打算只要求用户提交最终答案即可（个别是一个数字），而不须要残缺代码。因而，能够纵情地应用本人喜爱的编程语言——不少题目甚至光靠笔和纸便能解决。
欧拉打算的第66题十分有意思，它的题目很简略，就是要求找出在不大于1000的整数中，以哪一个数字为丢番图方程的系数，能够失去所有最小解中的最大值。
能够很容易地看出方程有一个直观的暴力算法：让y从1开始递增，对于每一个y，计算公式Dy^2+1的值。如果该值为平方数，那么它的平方根就是最小的x解。再按照这个算法求解所有D不大于1000的方程，便能够求出题目标答案。很容易用Python写出这个算法
# -*- coding: utf8 -*-import mathdef is_square(num: int) -> bool: return math.isqrt(num) ** 2 == numdef 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 += 1def 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_xif __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 mathdef 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。有趣味的读者能够试一下暴力解法要花多久能力算到这个数字。
全文完。
浏览原文