Project Euler #66: Diophantine equation

  • + 0 comments

    I solved this problem using Pell's Equation fundermental solution.

    import math
    
    def is_perfect_square(n):
        """Check if a number is a perfect square."""
        root = int(math.sqrt(n))
        return root * root == n
    
    def continued_fraction_sqrt(n):
        """Compute the continued fraction representation of sqrt(n)."""
        m, d, a = 0, 1, int(math.sqrt(n))
        initial_a = a
        period = []
        
        if initial_a * initial_a == n:  # Perfect square check
            return None, []
        
        while True:
            m = d * a - m
            d = (n - m * m) // d
            a = (initial_a + m) // d
            period.append(a)
            if a == 2 * initial_a:  # Period repeats here
                break
        
        return initial_a, period
    
    def solve_pell_equation(n):
        """Solve Pell's equation x^2 - n * y^2 = 1 for the minimal solution."""
        if is_perfect_square(n):
            raise ValueError(f"The given number {n} is a perfect square, no solution exists.")
        
        a0, period = continued_fraction_sqrt(n)
        if not period:
            return None
        
        # Compute convergents using the periodic continued fraction
        p_prev, p_curr = 1, a0
        q_prev, q_curr = 0, 1
        
        for k in range(1, len(period) * 2 + 1):  # Repeat until period length is hit
            a_k = period[(k - 1) % len(period)]
            p_next = a_k * p_curr + p_prev
            q_next = a_k * q_curr + q_prev
            
            # Check if the solution satisfies Pell's equation
            if p_next * p_next - n * q_next * q_next == 1:
                return p_next, q_next
            
            # Update for the next iteration
            p_prev, p_curr = p_curr, p_next
            q_prev, q_curr = q_curr, q_next
    
    if __name__ == '__main__':
        N = int(input())
        max_x = 0
        index = 2
        for i in range(2, N + 1):
            try:
                x, y = solve_pell_equation(i)
                if x > max_x:
                    max_x = x
                    index = i
            except ValueError as e:
                pass
        print(index)