此数绵绵无绝期——欧拉计划第66题

# -*- 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))

# 计算连分数数列的下一个数字
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))

(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))

(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)))