I admire myself eight years ago with this article
The preface
Project Euler, like LeetCode, is a programming quiz site. Unlike LeetCode, Euler’s plan simply requires the user to submit a final answer (usually a number), not the full code. So, feel free to use your favorite programming language — many problems can be solved with pen and paper alone.
Question 66 of euler’s plan is very interesting, its title is very simple, is required to find in the integer not greater than 1000, which number is the coefficient of diophantine equation, can get the maximum of all the smallest solutions.
It is easy to see that the equation has an intuitive brute force algorithm: let y increase from 1, and for each y, evaluate the formula Dy^2+1. If the value is square, then its square root is the smallest x solution. And then you solve all the equations with D not greater than 1,000, and you get the answer. This algorithm is easy to write in Python
# -*- coding: utf8 -*-
import math
def is_square(num: int) - >bool:
return math.isqrt(num) ** 2 == num
def find_x(D: int) - >int:
"" Find the smallest x that satisfies the Diophantine equation given D. "" "
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) :
Find the number in D not greater than limi that maximizes the return value of 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))
Copy the codeBut if you raise the limit to 1000, the algorithm won’t work in a lifetime.
To solve this problem, we need the power of mathematics.
Pell equation
When I first did this problem eight years ago, after some searching, I learned from this article that the equation in question was called the Pell equation. It has a standard solution, but it uses continued fractions. So what is a continued fraction?
Continued fraction is not a new number system, just another way of writing decimals. For example, the fraction 45 divided by 16 can be written as follows
Just as we define recursive data structures, we can give continued fractions a recursive definition. A continued fraction is either an integer or the reciprocal of an integer plus another continued fraction. In addition to the above form, even fractions can be written in a more efficient form. For example, to write 45 divided by 16 as [2;1,4,3], all the integer parts of the original formula are written in order between square brackets. This notation looks like an array in a programming language.
If the fractions are constructed with different prefixes of the array [2;1,4,3], the results are 2/1, 3/1, and 14/5, respectively. They are the asymptotic continued fraction of this continued fraction, and a set of solutions to the Pell equation comes from the numerator and denominator of the asymptotic continued fraction.
Take the Pell equation with a coefficient of 7 as an example, first calculate the continued fraction of the square root of 7, and then try its asymptotic continued fraction in turn. The first three, 2/1, 3/1 and 5/2, are not solutions to the equation. The fourth progressive continued fraction, 8/3, is the solution. If we continue to improve the accuracy of continued fractions, a second solution, 127/48, will be found. Keep looking. There’s more, and 8 is the smallest x.
So, to quickly solve the Pell equation, the most important thing is to find an algorithm for calculating the continued fraction of the square root of a number.
The wrong way to calculate the continued fraction of square root
To compute the continued fraction of a number, the most important thing is to compute all the whole number parts (A0, A2, A2, etc.). They can all be computed directly by definition
Generalizing to half the case, if you use the variable n to store the square root of the number, and use numbers to store all known integers, then in Python you can write the following algorithm to compute the next integer
# Calculate the next number in the continuous fraction sequence
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('The next number is {}'.format(v))
Copy the codeUnfortunately, the results of this algorithm can be very inaccurate due to the accuracy of the calculation.
The correct way to calculate the continued fraction of square roots
To get the correct result, you need to eliminate as much as possible of the error introduced in the calculation of 1 / (v-a), so you have to remove the floating point number from the denominator.
In this site, the author presents a table for calculating the continued fraction of the square root of 14 as an example
You can see that x1, x2, and x3 are of the form SQRT (n)+a /b, which is much better for controlling errors. So does every x that we compute fit this format? The answer is yes, and it can be proved by mathematical induction (written in LaTeX with a screenshot to make writing the formula easier).
In the process of proving this, we also have the recursive formula for a in the numerator and b in the denominator, and now we can write the correct code for calculating the integer part of the continued fraction.
Common Lisp is used to implement the above algorithm
To implement this algorithm while writing elegant code, I’ll use the object-oriented features of Common Lisp. The first is to define a class to represent a continued fraction with increasing precision
(defpackage #:com.liutos.cf
(:use #:cl(a))in-package #:com.liutos.cf)
(defclass <cf> ()
((a
:documentation "The number added to the square root of a numerator or numerator in mathematical induction."
:initform 0)
(b
:documentation "Denominator in mathematical induction."
:initform 1)
(numbers
:documentation "An array of the integral parts of a continued fraction."
:initform nil)
(origin
:documentation "The number being squared."
:initarg :origin(a)):documentation "The continued fraction representing the square root of the integer ORIGIN."))
Copy the codeThen define the “interface” that the class needs to implement.
(defgeneric advance (cf)
(:documentation "Let the continued fraction CF be raised to the next accuracy."(a))defgeneric into-rational (cf)
(:documentation "Converts continued fraction CF to a value of rational type."))
Copy the codeFinally, the two interfaces are implemented
(defmethod advance ((cf <cf>))
"Calculate the next batch of a's and B's, and the whole part of the continued fraction, according to the recursive formula."
(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))
Copy the codeCommon Lisp’s rational numeric types have been a huge benefit to me in implementing the Into-Rational approach. It allows me to write code that is straightforward without worrying about introducing errors into the calculation of (/ 1 v).
The problem solving
Follow up on your success and use Common Lisp to solve 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(a))let ((target-d)
(max-x 0(a))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)))
Copy the codeAnswer D what is not, but as the answer is 16421658242965910275055840472270471049 x. Interested readers can try the brute force solution to figure out how long it takes to get to that number.
The full text.
Read the original