Project Euler #66 - 디오판틴 방정식

프로젝트 오일러 66번 - Diophantine equation

1000보다 작은 D 중에서(단, D는 제곱수가 아니다)

x^2 - D y^2 = 1 의 최소정수해 중 x를 가장 크게 만드는 것은?

모든 D에 대해서 x의 최소해를 구해 봐야 한다.

D를 고정하고 x와 y를 늘려가며 찾아보는 건 쉬운데.. 답이 안 나온다. 한참 고민하다가 구글링..

D = 61일 때의 최소해는 x=1766319049, y = 226153980 이라고 한다. brute force는 포기. 인터넷 글을 읽다 보니, 페르마도 풀려다 못 풀었다고 한다. 페르마가 못 푼 문제를 내가 풀 수 있을 것 같지는 않고.. 풀이 방법까지 찾아보게 됐다.

디오판틴 방정식 중에서 특별히 이런 모양의 방정식을 펠방정식 - Pell's equation이라고 부른다는 것도 알게 됐다.

http://en.wikipedia.org/wiki/Chakravala_method#Examples

여기에 있는 풀이 방법을 따랐다. 우변이 ±1, ±2, or ±4 일 때는 좀 더 쉬운 특별한 방법이 있다고 하지만, 그 '특별한 방법'을 쓰지 않아도 iteration을 조금 더 진행하면 답이 나온다고 위키피디아에 친절하게 나와 있다.

#http://en.wikipedia.org/wiki/Chakravala_method#Examples
from math import sqrt

def dio_init(D):
    Di = int(sqrt(D))
    tmpdic = {Di:abs(Di**2-D),Di+1:abs((Di+1)**2-D)}
    a = min(tmpdic,key=tmpdic.get)
    return a,1,a*a-D

def dio_update(a,b,k,D):
    tmpdic = {}
    for m in range(int(sqrt(D))-abs(k),int(sqrt(D))+abs(k)+1):
        if (a + b*m) % abs(k) == 0:
            tmpdic[m] = abs(m*m - D)
    m = min(tmpdic,key=tmpdic.get)
    return (a*m+D*b)/abs(k),(a+b*m)/abs(k),(m*m-D)/k

def dio(D):
    a,b,k = dio_init(D)
    while True:
        a,b,k = dio_update(a,b,k,D)
        if k == 1: return a, b

mx, argmx = 0,0
for D in [x for x in range(2,1001) if x not in [y*y for y in range(1,32)]]:
    if dio(D)[0] > mx: mx = dio(D)[0]; argmx = D

print argmx

ps. x**.5 와 sqrt(x)는 같다. 후자는 from math import sqrt가 필요해서 보통은 0.5제곱하는 쪽을 선호한다. 하지만 가끔, sqrt로 쓰는 게 더 잘 읽히는 날은 sqrt로 쓴다.

psps. 수학이나 알고리즘 공부가 아니라 python 코딩공부를 목적으로 프로젝트 오일러 문제를 풀고 있었다. 이 문제를 접하고 나서 동기부여가 안 된다. 66번에서 그만 두기는 거시기하니까 100번까지 하고 그만 두는 걸로...