Problem 66 - Python

So, problem 66...It involves solving Diophantine equations, which is arguably the scariest math word to appear in a Project Euler problem up to this point (yes, it just means we are looking for solutions over the integers, but Diophantine...big word). Trying to answer this problem by brute force gets one nowhere...many of the minimal solutions are absolutely enormous, and would take forever to find by brute force (especially as one can't exactly binary search for a minimal solution...). Instead, one must find out by some means (my source: Wikipedia), that the equation in question is Pell's equation, which has nice solutions in terms of the convergents of the continued fraction expansion of the square root of D. Now, Euler did a good job of preparing most users for this problem: Problem 64 was about the continued fraction expansion of square roots, and problem 65 was about finding convergents of continued fractions. So if you are most people, by the time you reached problem 66, you should have had most of the code written already. However, as I must solve all my problems in different languages, this was not the case for me. Now, I decided to actually reuse my code for problem 65, but that meant resolving problem 65, as explained in the posts below. I did however have to port my problem 64 solution in Tcl into python, which was not too difficult. All in all, the combination of code reuse and wikipedia reading made this problem not terribly hard, but it would certainly be hard in an arbitrary language. My solution takes advantage of Python's lists, dictionaries, classes, arbitrary precision integer arithmetic, and potential for functions to have multiple return values. Resolving this problem in another language would certainly not be an easy task, and is not one I plan on doing any time too soon - this problem makes for the 5th Python solution, so its probably not a bad thing, either that I will probably not recover Python for a while. Anyway, now that this problem is done, which required writing three solutions, onwards to problem 67: which is essentially identical to problem 18, explaining why 5x as many people have solved that problem as had solved 66. Below is my solution, which runs in about 1.3s on my machine.

import math

def gcd(a, b):
    return a if b == 0 else gcd(b, a%b)

class rational:
    def __add__(self, other):
        return rational(self.num*other.denom+other.num*self.denom,
                        self.denom*other.denom)

    def inverse(self):
        return rational(self.denom, self.num)

    def __init__(self, n, d):
        g = gcd(n, d)
        self.num   = n / g
        self.denom = d / g

coefDict = {}

def coeff(n,D):
    if n == 1:
        return int(math.sqrt(D))
    else:
        lst = coefDict[D]
        return lst[(n - 2)%len(lst)]

def nextDepth(n,N,D):
    if n == N:
        return rational(coeff(n,D),1)
    else:
        return rational(coeff(n,D),1) + nextDepth(n+1,N,D).inverse()

def nthConv(n,D):
    return nextDepth(1,n,D)

def mag(x):
    return 100 if x < 10 else 10 * mag(x / 10)

def nextTerm(N, m, b):
    d  = (N - b ** 2) / m
    k  = int(math.ceil(2.0 * b / d))
    bp = k*d - b
    if bp ** 2 > N:
        k -= 1
        bp = k*d - b
    return k,d,bp

def computeExpansion(D):
    b = int(math.sqrt(D))
    coefDict[D] = []
    m = 1
    computations = []
    while True:
        a, m, b = nextTerm(D, m, b)
        c = mag(b)*m + b
        if c not in computations:
            coefDict[D].append(a)
            computations.append(c)
        else:
            break

def minx(D):
    computeExpansion(D)
    n = 1
    while True:
        test = nthConv(n,D)
        if test.num ** 2 - D * (test.denom ** 2) == 1:
            return test.num
        n += 1

ans = 0
maxsol = 0
for D in range(2,1000):
    if int(math.sqrt(D)) ** 2 != D:
        x = minx(D)
        if x > maxsol:
            ans = D
            maxsol = x
print ans