Problem66.py

from collections import Counter
from fractions import Fraction

def Solve():

    '''final answer. It turns out this specific equation is called a Pell equation. for all (x,y)
    solutions to a Pell equation the rational number x/y is an approximation of d**.5 that is generated by the
    expansion of the continued fraction. Not all expansion give solutions though, so each must be tested.'''

    limit = 1000
    maxX = 0

    def nextFraction(s, a = 0, b = 1):

        a, b = -a, int((s-a**2)/b)
        i = int((s ** .5 + a)/b)
        a -= i * b

        return [i, a, b]

    def findX(d):
        root = d ** .5
        i = int(root)
        next = [i, -i, 1]
        contFractExpand = [next]
        while True:
            next = nextFraction(d, *next[1:])
            contFractExpand.append(next)
            temp = Fraction(next[0])
            for f in contFractExpand[-2::-1]:
                temp = 1 / temp + f[0]
            if not (temp.numerator**2 - d*temp.denominator**2 - 1):
                # print(d, temp)
                return(temp.numerator)

    squares = set([x ** 2 for x in range(1, 1 + int(limit ** .5))])

    for d in range(2, limit + 1):

        if sum([d + x in squares for x in [-2, 0, 2]]):
            continue

        x = findX(d)
        if x > maxX:
            maxX = x
            maxD = d
            # print(d, x)

    return maxD

def Solve1():

    '''This was the best I could do without researching Diophantaine equations.'''

    limit = 60

    maxD = 0
    maxX = 0

    def fD(x, y):
        return int(x**2 - d * y**2 - 1)

    r = 0
    squares = set([x**2 for x in range(1, 1+int(limit**.5))])

    y = 2
    for d in range(2, limit+1):

        if sum([d + x in squares for x in [-2, 0, 2]]):
            continue

        x = .5
        y = 0
        while not x == int(x):
            y += 1
            x = (1 + d * y ** 2) ** .5

        # print(d, int(x))
        if x > maxX:
            maxX = int(x)
            maxD = d

    return(maxD)


def Solve2():

    '''My first attempt. Pure numeric solving. Took an absurd amount of time for limit = 100, so never ran for anything
    larger than that.'''

    limit = 60

    maxD = 0
    maxX = 0

    def fD(x, y):
        return int(x**2 - d * y**2 - 1)

    r = 0
    squares = set([x**2 for x in range(1, 1+int(limit**.5))])

    y = 2
    for d in range(2, limit+1):

        if sum([d + x in squares for x in [-2, 0, 2]]):
            continue

        x = 0
        y = 1
        f = fD(x, y)
        while f:
            if f < 0:
                x += 1
            else:
                y += 1
            f = fD(x, y)

        # print(d, int(x))
        if x > maxX:
            maxX = int(x)
            maxD = d

    return(maxD)


def Solve3():

    '''This was my second attempt. I realized that x and y needed to have different parity so I only tested the
    d values for the y values that had opposite partiy of the x values. it finds all the d values for each x value
    and continues until it hits the limit. This takes forever, because there are lots of x values that aren't
    correct except for very large d values. i.e. it finds thousands of solutions that are outside the parameters'''

    limit = 28
    leastX = Counter()
    stopLength = limit - int(limit ** .5)

    x = 1
    while len(leastX) < stopLength:
        x += 1
        lowY = int(max(((x**2 - 1)/limit)**.5, 1))
        if (x+lowY)%2 == 0:
            lowY += 1
        for y in range(lowY, x, 2):
            d = (x**2-1)/(y**2)
            if int(d) == d <= limit and not d in leastX.keys():
                d = int(d)
                for f in range(1, int(d**.5)+1):
                    if d % (f**2) == 0:
                        leastX[d//(f**2)] = x
                        # print(x, '^2 - ', d//(f**2), '*', y, '^2 = 1', sep='')
                break

    return leastX.most_common(1)[0][0]


def Solve4():

    '''for this algorithm, I realized that x**2 = 1 mod d * y**2. I wrote a generator that would provide increments
    for x so that we weren't testing every value of x. since we are setting x and then testing the resulting y value
    there is no gain to a parity test.

    I knew from my first couple algorithms that the solution would be much larger than 1000 and during this then
    I realized that x values for d that were a perfect square +/- 2 were necessarily close to d.
    x**2 - d * y **2 = 1
    x**2 - 1 = d * y**2
    (x-1) * (x+1) = d * y**2
    x +/- 1 = y**2 AND x-/+1 = d is a possible solution, which implies
    d +/- 2 = y**2

    Since all values of d are much lower than the solution I skipped those values entirely.'''

    limit = 60

    maxD = 0
    maxX = 0

    def stepGen(d):
        mods = [x**2 % d for x in range(d)][::-1]
        cycle = []
        step = 0
        while len(mods):
            if mods.pop() == 1:
                cycle.append(step)
                step = 0
            step += 1
        yield cycle[0]
        cycle[0] += step
        while True:
            cycle = cycle[1:] + [cycle[0]]
            yield cycle[0]

    squares = set([x**2 for x in range(1, 1+int(limit**.5))])

    for d in range(2, limit + 1):

        if sum([d + x in squares for x in [-2, 0, 2]]):
            continue

        x_step = stepGen(d)

        #burn the first step, so that y != 0
        x = x_step.__next__()

        while True:
            x += x_step.__next__()
            y = ((x**2-1)/d)**.5
            if y.is_integer():
                # print(d, x)
                if x > maxX:
                    maxX = x
                    maxD = d
                break

    return(maxD)


if __name__ == '__main__':
    from time import clock
    s = clock()
    print(0, Solve(), clock() - s)
    s = clock()
    # print(1, Solve1(), clock() - s)
    # s = clock()
    # print(2, Solve2(), clock() - s)
    # s = clock()
    # print(3, Solve3(), clock() - s)
    # s = clock()
    # print(4, Solve4(), clock() - s)
    # s = clock()


if __name__ == '__main__':
    print(Solve())