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)