93 lines
2.8 KiB
Python
93 lines
2.8 KiB
Python
"""
|
|
Consider quadratic Diophantine equations of the form:
|
|
x^2 - D y^2 = 1
|
|
For example, when D = 13 , the minimal solution in x is $649^2 - 13 * 180^2 = 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 * 2^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 time
|
|
from functools import wraps
|
|
from math import isqrt
|
|
from typing import Any, Callable, TypeVar
|
|
|
|
F = TypeVar("F", bound=Callable[..., Any])
|
|
|
|
|
|
def benchmark(repeat: int = 1) -> Callable[[F], F]:
|
|
if repeat < 1:
|
|
raise ValueError("repeat >= 1")
|
|
|
|
def decorator(func: F) -> F:
|
|
@wraps(func)
|
|
def wrapper(*args: Any, **kwargs: Any) -> Any:
|
|
total = 0.0
|
|
result = None
|
|
|
|
for _ in range(repeat):
|
|
start = time.perf_counter()
|
|
result = func(*args, **kwargs)
|
|
end = time.perf_counter()
|
|
total += end - start
|
|
|
|
wrapper.avg_time = total / repeat # type: ignore[attr-defined]
|
|
wrapper.total_time = total # type: ignore[attr-defined]
|
|
|
|
print(
|
|
f"[Benchmark] {func.__name__} | repeated {repeat} times | "
|
|
f"average: {wrapper.avg_time:.6f}s | total: {wrapper.total_time:.6f}s" # type: ignore[attr-defined]
|
|
)
|
|
return result
|
|
|
|
wrapper.avg_time = 0.0 # type: ignore[attr-defined]
|
|
wrapper.total_time = 0.0 # type: ignore[attr-defined]
|
|
return wrapper # type: ignore[return-value]
|
|
|
|
return decorator
|
|
|
|
|
|
def diophantine_eq(D: int) -> int | None:
|
|
a = [int(isqrt(D))]
|
|
if a[0] * a[0] == D:
|
|
return None
|
|
p = [0, 1, int(isqrt(D))]
|
|
m = 0
|
|
d = 1
|
|
while True:
|
|
m = d * a[-1] - m
|
|
d = (D - m * m) // d
|
|
a.append((a[0] + m) // d)
|
|
p.append(p[-1] * a[-1] + p[-2])
|
|
if a[-1] == 2 * a[0]:
|
|
if len(a) % 2 == 1:
|
|
return p[-2]
|
|
else:
|
|
for ai in a[1:-1]:
|
|
p.append(p[-1] * ai + p[-2])
|
|
return p[-1]
|
|
|
|
|
|
@benchmark(repeat=10)
|
|
def find_minimal_solution(limit: int) -> int:
|
|
min_x: dict[int, int] = {}
|
|
for D in range(2, limit + 1):
|
|
x = diophantine_eq(D)
|
|
if x is not None:
|
|
min_x[D] = x
|
|
return max(min_x, key=min_x.__getitem__)
|
|
|
|
|
|
if __name__ == "__main__":
|
|
max_D = 1000
|
|
dd = find_minimal_solution(max_D)
|
|
print(f"D <= {max_D} |-> D of max x = {dd}")
|