오일러 프로젝트 46

오일러 프로젝트 46 번은 은근히 시간이 많이 걸리는 문제이다. 모든 소수가 아닌 홀수가 소수와 제곱수의 두 배의 합으로 나타낼 수 있다는 추측의 반례를 찾는 것이다. 

크리스티안 골드바흐는 모든 홀수인 합성수를 (소수 + 2×제곱수)로 나타낼 수 있다고 주장했습니다.

9 = 7 + 2×12

15 = 7 + 2×22

21 = 3 + 2×32

25 = 7 + 2×32

27 = 19 + 2×22

33 = 31 + 2×12

이 추측은 잘못되었음이 밝혀졌습니다. 위와 같은 방법으로 나타낼 수 없는 가장 작은 홀수 합성수는 얼마입니까?(http://euler.synap.co.kr/prob_detail.php?id=46)

접근

임의의 자연수 N에 대해서 위 추측의 반례인지를 확인하는 코드를 작성하는 것은 간단하다.

  1. N은 홀수여야 함
  2. N은 소수가 아니어야 함
  3. N의 절반보다 작은 제곱수의 목록을 만들고
  4. N보다 작은 소수체를 만든다.
  5. N보다 작은 제곱수  S에 대해서 N – S * 2 의 값을 계산하고, 이 값을 소수체에서 검사해서 있으면 추측을 만족, 없으면 반례가 된다.

그런데 매 N을 변경할 때마다 소수체를 만드는 것이 좀 비효율적으로 보인다. 따라서 다음과 같이 개선한다.

  1. 소수 판별은 소수체를 사용하지 않고 효율적인 소수 검사 함수를 사용한다.[^1]
  2. 제곱수는 2배하기로 하였으므로 제곱근의 절반까지의 범위에서만 구한다.

이에 대한 Swift 풀이는 아래와 같다.

import Foundation

func isPrime(_ n: Int) -> Bool {
  if n < 2 { return false }
  if (2...3) ~= n { return true }
  if n % 2 == 0 || n % 3 == 0 { return false }
  let l = Int(sqrt(Double(n)) + 0.5)
  var k = 5
  while k <= l {
    if n % k == 0 || n % (k + 2) == 0 { return false }
    k += 6
  }
  return true
}

func squared(upto n: Int) -> [Int] {
  let l = Int(sqrt(Double(n)) / 2 + 0.5)
  return (1...l).map{ $0 * $0 }
}

func test(_ n: Int) -> Bool {
  if isPrime(n) { return false }
  let sq = squared(upto: n)
  for s in sq where s * 2 < n - 1 {
    let p = n - s * 2
    if isPrime(p) { return false }
  }
  return false
}

main: do {
  var a = 9
  while !test(a) { a += 2 }
  print(a)
}

파이썬 풀이

귀찮아서 빨리 타이핑하고, 가독성은 그냥 무시했다. 이 포스트의 수정 전 방식은 소수체를 매번 만드는 거 였는데 그 때는 3초가량, 지금은 0.1초 미만이 걸린다.

def memoize(f):
  memo = dict()
  def wrapped(n):
    if n in memo: returm memo[n]
    r = f(n)
    memo[n] = r
    return r
  return wrapped

@memoize
def is_prime(n):
  if n < 2: return False
  if n in (2, 3): return True
  if n % 2 is 0 or n % 3 is 0 : return False
  if n < 10 : return True
  k, l = 5, n **0.5
  while k <= l:
    if n % k is 0 or n % (k+2) is 0: return False
    a += 6
  return True

def squared(n):
  l = int((n/2)**0.5) + 1
  return [x*x for x in range(1, l)]

def test(n):
  if is_prime(n): return False
  for s in squared(n):
    p = n - s * 2
    if is_prime(p): return False
  return True

def solve():
  a = 9
  while not test(a): a += 2
  print(a)

%time solve()