Project Euler #100: Arranged Probability.py

import time


start_time = time.time()
# b(b-1)/(b+r)(b+r-1)... b= (1+sqrt(8r^2+1))//2 + r
# k^2-8r^2=1 (3,1)
# (xx'+yy'd)+(xy'+x'y)sqrt(d) => 3*x+8*y, x+3*y


def solution(limit):
    x, y = 3, 1

    while 2*y+(x+1)//2 < limit:
        x, y = 3*x+8*y, x+3*y

        yield y+(x+1)//2


print([x for x in solution(10**12)][-1])


print("--- %s seconds ---" % (time.time() - start_time))