Euler 66

gistfile1.txt

# http://projecteuler.net/problem=66
# this is kind of a bear- it combines all the continued fraction stuff
# but it runs in 147ms

import math

squares=[]
for i in range (1,101):
	squares.append(i*i)

def CF_of_sqrt(n):
# modified this from http://eli.thegreenplace.net/2009/06/19/project-euler-problem-66-and-continued-fractions/
# who got it from http://www.maths.surrey.ac.uk/hosted-sites/R.Knott/Fibonacci/cfINTRO.html
    if n in squares: # perfect squares don't have CFs
        return [int(math.sqrt(n))]
    ans = []
    step1_num = 0
    step1_denom = 1
    while True:
        nextn = int((math.floor(math.sqrt(n)) + step1_num) / step1_denom)
        ans.append(int(nextn))
        step2_num = step1_denom
        step2_denom = step1_num - step1_denom * nextn
        step3_denom = (n - step2_denom ** 2) / step2_num
        step3_num = -step2_denom
        if step3_denom == 1:
            ans.append(ans[0] * 2)
            break
        step1_num, step1_denom = step3_num, step3_denom
    return ans

# these two functions compute the numerator and denominator of the pth convergent
def convNum(CF,p):
	if p == 1:
		return CF[0]
	if p == 2:
		return CF[0]*CF[1]+1
	nlist = CF
	while len(nlist) < p:
		nlist = nlist+nlist[1::] 	
	start = [nlist[0],nlist[0]*nlist[1]+1]
	for i in range (2,p):
		newterm = start[i-2] + start[i-1]*nlist[i]
		start.append(newterm)
	return start

def convDen(CF,p):
	if p == 1:
		return 1
	if p == 2:
		return CF[1]
	nlist = CF
	while len(nlist) < p:
		nlist = nlist+nlist[1::] 	
	start = [1,nlist[1]]
	for i in range (2,p):
		newterm = start[i-2] + start[i-1]*nlist[i]
		start.append(newterm)
	return start

def diophantine(a,b,D):
	if a*a - D*b*b == 1:
		return True
	else:
		return False

#p is our search limit- check first p convergents- by trial and error, 77 is min value
p=78
max_d = 0
max_num=0
for d in range (1,1001):
	if d not in squares:
		num = convNum(CF_of_sqrt(d),p)
		den = convDen(CF_of_sqrt(d),p)
		for j in range (0,p):
			if diophantine(num[j],den[j],d):
#				print d, j,num[j],den[j]
				if num[j]>max_num:
					max_num=num[j]
					max_d=d
				break
			elif j==p-1:
#				print d, CF_of_sqrt(d),j,"no solution"
				break		

print max_num,max_d