Consider quadratic Diophantine equations of the form:
+$x^2 – Dy^2 = 1$
+For example, when D=13, the minimal solution in x is 6492 – 13×1802 = 1.
+It can be assumed that there are no solutions in positive integers when D is square.
+By finding minimal solutions in x for D = {2, 3, 5, 6, 7}, we obtain the following:
+$3^2 – 2×2^2 = 1$
+$2^2 – 3×1^2 = 1$
+$9^2 – 5×4^2 = 1$
+$5^2 – 6×2^2 = 1$
+$8^2 – 7×3^2 = 1$
+Hence, by considering minimal solutions in x for D ≤ 7, the largest x is obtained when D=5.
+Find the value of D ≤ 1000 in minimal solutions of x for which the largest value of x is obtained.
+import math
+
+def get_floor_sqrt(n):
+ return math.floor(math.sqrt(n))
+
+assert(get_floor_sqrt(5) == 2)
+assert(get_floor_sqrt(23) == 4)
+
+def next_expansion_1(current_a, current_nominator, current_denominator, original_number):
+ cn = current_nominator
+ cd = current_denominator
+ cn, cd = cd, (original_number - cd * cd) // cn
+ next_a = math.floor((math.sqrt(original_number) + cn) // cd)
+ cn = cn - next_a * cd
+ cn *= -1
+
+ return next_a, cn, cd
+
+def get_continued_fraction_sequence(n):
+
+ # If number is a square number we return it.
+ floor_sqrt = get_floor_sqrt(n)
+ if n == floor_sqrt * floor_sqrt:
+ return ((floor_sqrt), [])
+
+ # Otherwise, we calculate the next expansion till we
+ # encounter a step a second time.
+ a = floor_sqrt
+ cn = a
+ cd = 1
+ sequence = []
+ previous_steps = []
+
+ while not (a, cn, cd) in previous_steps :
+ #print("a: {} cn: {} cd: {}".format(a, cn, cd))
+ previous_steps.append((a, cn, cd))
+ a, cn, cd = next_expansion_1(a, cd, cn, n)
+ sequence.append(a)
+ sequence.pop()
+ return ((floor_sqrt), sequence)
+
+def get_continued_fraction_sequence(n):
+
+ # If number is a square number we return it.
+ floor_sqrt = get_floor_sqrt(n)
+ if n == floor_sqrt * floor_sqrt:
+ return ((floor_sqrt), [])
+
+ # Otherwise, we calculate the next expansion till we
+ # encounter a step a second time.
+ a = floor_sqrt
+ cn = a
+ cd = 1
+ sequence = []
+ previous_steps = []
+
+ while not (a, cn, cd) in previous_steps :
+ #print("a: {} cn: {} cd: {}".format(a, cn, cd))
+ previous_steps.append((a, cn, cd))
+ a, cn, cd = next_expansion_1(a, cd, cn, n)
+ sequence.append(a)
+ sequence.pop()
+ return ((floor_sqrt), sequence)
+
+assert(get_continued_fraction_sequence(1) == ((1), []))
+assert(get_continued_fraction_sequence(4) == ((2), []))
+assert(get_continued_fraction_sequence(25) == ((5), []))
+assert(get_continued_fraction_sequence(2) == ((1), [2]))
+assert(get_continued_fraction_sequence(3) == ((1), [1,2]))
+assert(get_continued_fraction_sequence(5) == ((2), [4]))
+assert(get_continued_fraction_sequence(13) == ((3), [1,1,1,1,6]))
+def gcd(a, b):
+ if b > a:
+ a, b = b, a
+ while a % b != 0:
+ a, b = b, a % b
+ return b
+
+def add_fractions(n1, d1, n2, d2):
+ d = d1 * d2
+ n1 = n1 * (d // d1)
+ n2 = n2 * (d // d2)
+ n = n1 + n2
+ p = gcd(n, d)
+ return (n // p, d // p)
+
+def next_expansion_2(previous_numerator, previous_denumerator, value):
+ if previous_numerator == 0:
+ return (value, 1)
+ return add_fractions(previous_denumerator, previous_numerator, value, 1)
+
+def get_fractions(n, x):
+ #print("get_fractions(n={}, x={})".format(n, x))
+ # Get sequence of a_x
+ first_value, sequence = get_continued_fraction_sequence(n)
+ sequence = [first_value] + math.ceil(x / len(sequence)) * sequence
+ sequence = sequence[:x + 1]
+ sequence = sequence[::-1]
+
+ n, d = 0, 1
+ for s in sequence:
+ n, d = next_expansion_2(n, d, s)
+ return (n, d)
+
+assert(get_fractions(7, 0) == (2, 1))
+assert(get_fractions(7, 1) == (3, 1))
+assert(get_fractions(7, 2) == (5, 2))
+assert(get_fractions(7, 3) == (8, 3))
+def get_minimal_solution(d):
+ for i in range(0, 100):
+ x, y = get_fractions(d, i)
+ if x * x - d * y * y == 1:
+ return((x, y))
+
+assert(get_minimal_solution(2) == (3, 2))
+assert(get_minimal_solution(3) == (2, 1))
+assert(get_minimal_solution(85) == (285769, 30996))
+assert(get_minimal_solution(109) == (158070671986249, 15140424455100))
+from math import sqrt
+
+def is_square(n):
+ return sqrt(n).is_integer()
+
+print(is_square(4))
+x_max = 0
+d_max = 0
+
+for d in range(2, 1001):
+ if is_square(d):
+ continue
+
+ x, y = get_minimal_solution(d)
+ if x > x_max:
+ x_max = x
+ d_max = d
+
+print("d: {} x: {}".format(d_max, x_max))
+s = d_max
+assert(s == 661)
+
+