diff --git a/project_euler/problem_066/__init__.py b/project_euler/problem_066/__init__.py new file mode 100644 index 000000000000..e69de29bb2d1 diff --git a/project_euler/problem_066/sol1.py b/project_euler/problem_066/sol1.py new file mode 100644 index 000000000000..7e619d10ae60 --- /dev/null +++ b/project_euler/problem_066/sol1.py @@ -0,0 +1,191 @@ +""" +Project Euler Problem 61: https://projecteuler.net/problem=66 + +Problem statement- +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 * 1^2 = 1 +9^2 - 5 * 4^2 = 1 (highest value of x) +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 is obtained. + + +Solution explanation - +The above given equation is known as Pell's equation. The solution uses the standard +method for solving pell's equation. Which is available on the wikipedia page. +The solution uses the continued fraction method for finding close approximations +of sqrt(D) which are of the form h/k. Because for these solutions the h and k are +solutions for x and y respectively. But we dont want an exact approximation because +there is a 1 on the RHS. We consider upto 2 periodic cycles, which guarantee a +solution. We check in the continued fraction considering newer fractions and checking +if numerator(x) and denominator(y) satisfy the diophantine equation. + + +References- +Wikipedia - https://en.wikipedia.org/wiki/Pell%27s_equation + + +""" + +import math + + +def is_perfect_square(num: int) -> bool: + """ + Checks if the number inputted is a perfect square or not. + Returns True if perfect square else False. + + >>> is_perfect_square(25) + True + >>> is_perfect_square(3) + False + + Time complexity - O(log2(num)) + + """ + + if num < 0: + return False + root = math.isqrt(num) + return root * root == num + + +def continued_fraction(num: int) -> list[int]: + """ + Computes and returns the continued fraction of the square root of the given integer. + Note for irrational numbers, the continued fraction stops when there is repetition. + The returned list is of the form -> [a0;a1,a2,a3...an] + where a1,a2,a3...an represents the minimumperiodic part of the continued fraction. + + Time complexity - O(sqrt(num)) + + >>> continued_fraction(2) + [1, 2] + >>> continued_fraction(13) + [3, 1, 1, 1, 1, 6] + """ + + ans = [] + + # Initialize variables for continued fraction + m = 0 + d = 1 + a0 = math.isqrt(num) + + # this is the end point from where the sequence starts to repeat + terminate_point = 2 * a0 + + ans.append(a0) + an = a0 + + # loop till termination point is reached, storing each an value + while an != terminate_point: + m = d * an - m + d = (num - m**2) / d + an = math.floor((a0 + m) / d) + ans.append(an) + + return ans + + +def diophantine_equation(var_x: int, var_y: int, var_d: int) -> int: + """ + Returns the LHS of the diophantine equation using x, y and D. + + >>> diophantine_equation(3, 2, 2) + 1 + >>> diophantine_equation(8, 3, 7) + 1 + >>> diophantine_equation(3, 4, 5) + -71 + """ + + return var_x**2 - var_d * var_y**2 + + +def solution(max_d: int = 1000) -> int: + """ + This function provides a solution to the problem 66 from project Euler. + The original problem is for D<= 1000. But there is a optional parameter + max_d, by changing it you can find solutions for any D theoretically. + + Time complexity - O(max_d * (log2(max_d) + sqrt(max_d))) + But since sqrt(max_d) dominates log2(max_d) + Final time complexity is - O(max_d^3/2) + + >>> solution() + 661 + >>> solution(10000) + 9949 + """ + + # Initialise answer which stored max value of D. + answer = -1 + # Initialise max_x which stores max value of x. + max_x = -1 + + # loop over all values of D + for d in range(2, max_d + 1): + # skip if d is prefect square + if is_perfect_square(d): + continue + + # Get continued fraction and then + # Add the periodic part of the fraction to itself because we are + # guaranteed to find the solution in 2 cycles of the periodic part. + # example - + # fraction = [3,1,1,2] + # First element is the integer part, rest is periodic fractional part. + # Add the fractional part on itself. + # fraction becomes [3,1,1,2,1,1,2]. + # This contains our solution for sure. + fraction = continued_fraction(d) + fraction += fraction[1:] + + num_0 = 0 + num_1 = 1 + den_0 = 1 + den_1 = 0 + + # This is a recursive relation to find the next numerator and denominator + num_2 = fraction[0] * num_1 + num_0 + den_2 = fraction[0] * den_1 + den_0 + + n = 1 + # Use the recursive relation to get new, more accurate numerator and + # denominator, then check if they satisfy the given diophantine equation. + while diophantine_equation(num_2, den_2, d) != 1: + num_0 = num_1 + num_1 = num_2 + num_2 = fraction[n] * num_1 + num_0 + + den_0 = den_1 + den_1 = den_2 + den_2 = fraction[n] * den_1 + den_0 + + n += 1 + + # If we see a new larger value of x, we store the corresponding D + # value as our answer. + if num_2 > max_x: + answer = d + max_x = num_2 + + return answer + + +if __name__ == "__main__": + print(f"{solution(10000) = }")