프로젝트 오일러 #66 디오판투스 수식

이번 문제는 난이도 25% 문제이지만, 수학적인 부분을 모르면 매우 어렵습니다.

그냥 수식만 따라가서 풀려면, 풀기 힘든 문제예요.

 

https://projecteuler.net/problem=66

Problem 66 - Project Euler

 

 

이와 같이 해가 정해지지 않는 부정방정식에서 그 중 정수해를 찾는 것은 일일이 대입하는 방법 외에는 중고등 수학과정만을 배운 사람에게는 알 수 없습니다.

 

\(x^2 - Dy^2 = 1 \) 은 디오판투스 부정방정식으로 D 값에 따라서 해가 무수하게 많이 존재할 수도 있고, 아닐 수도 있습니다.   D가 제곱수라고 한다면, 디오판투스 부정방정식은 \( x^2 = 1^2 + Dy^2 \) 형태가 되는데, 1로 시작하는 피타고라스 식을 만족하는 것은 존재하지 않으므로 정수해는 x = 1, y = 0 이 유일합니다.  D가 제곱수가 아니라면 해가 무수하게 존재합니다.

 

수학적인 것과 관련해서는 아래의 링크를 참고하세요.

https://sdev.tistory.com/213

 

이 문제는 D를 찾는 문제이니, 답은 1,000보다는 작은 값이지만, 최대값을 찾는 부분은 엄청 큰 수가 필요할 수 있습니다.

 

제가 짠 소스는 다음과 같습니다.  소스는 참고용으로 봐주세요.  Big integer 라이브러리가 필요해요.

//---------------------------------------------
//  Project Euler #66 - Diophantine equation
//    - by Aubrey Choi
//    - created at 2015-04-10
//---------------------------------------------
#include <stdio.h>
#include <memory.h>
#include <math.h>
#include "NxBigInt.h"

#define  LIMIT  1000

NxBigInt qde(int d);

int main()
{
  NxBigInt max = 0;
  int saved = 0;
  for(int d = 2 ; d <= LIMIT ; d++ )
  {
    NxBigInt x = qde(d);
    if( max < x ) { max = x; saved = d; }
  }
  printf("ans = %d\n", saved);
}

NxBigInt qde(int d)
{
  double r = sqrt(d);
  int q = 1, s = r;
  if( s*s == d ) return 0;
  NxBigInt x(r), y(1), h(1), k(0);
  while( x*x - d*y*y != 1 )
  {
    q = (d-s*s)/q;
    int a = (r+s)/q;
    s = a*q-s;

    NxBigInt ht = x, kt = y;
    x = a*x+h; h = ht;
    y = a*y+k; k = kt;
  }
  return x;
}

 

Big integer 라이브러리가 필요한 경우에는 파이썬 작성을 고려해봄직 합니다.

"""
//---------------------------------------------
//  Project Euler #66 - Diophantine equation
//    - by Aubrey Choi
//    - created at 2015-04-10
//---------------------------------------------
"""
import math

def qde(d):
  r = math.sqrt(d)
  x = int(r)
  if x*x == d: return 0
  y = 1
  h = 1
  k = 0
  a = x
  q = 1
  s = x;
  while x*x - d*y*y != 1:
    q = (d-s*s)/q
    a = int((r+s)/q)
    s = a*q-s

    ht = x
    kt = y
    x = a*x+h
    h = ht
    y = a*y+k
    k = kt

  return x

max = 0
saved = 0
for d in range(1, 1001):
  x = qde(d)
  if x > max: 
    max = x
    saved = d
print 'ans = ', max, '(', saved, ')'