오일러 프로젝트의 열두번째 문제. 이번에는 약수가 500개 이상인 삼각수를 찾는 문제이다. (http://euler.synap.co.kr/prob_detail.php?id=12)
1부터 n까지의 자연수를 차례로 더하여 구해진 값을 삼각수라고 합니다. 예를 들어 7번째 삼각수는 1 + 2 + 3 + 4 + 5 + 6 + 7 = 28이 됩니다. 이런 식으로 삼각수를 구해 나가면 다음과 같습니다.
1, 3, 6, 10, 15, 21, 28, 36, 45, 55, ...이 삼각수들의 약수를 구해봅시다.
1: 1 3: 1, 3 6: 1, 2, 3, 6 10: 1, 2, 5, 10 15: 1, 3, 5, 15 21: 1, 3, 7, 21 28: 1, 2, 4, 7, 14, 28위에서 보듯이, 5개 이상의 약수를 갖는 첫번째 삼각수는 28입니다. 그러면 500개 이상의 약수를 갖는 가장 작은 삼각수는 얼마입니까
삼각수
삼각수는 파스칼의 삼각형에 나오는 바로 그 수이고, “1부터 연속한 자연수들의 합”을 나타내는 수이다. 사실 1~N까지 연속한 자연수의 합은 간단한 공식으로 구할 수 있는데, 이는 곧 삼각수의 일반항을 구하는 식이 된다.
약수의 개수 세기
어떤 자연수 N이 p로 나누어 떨어진다 할 때, p는 N의 약수가 된다. 이 때, 
- i = 1부터 N의 제곱근까지 1씩 증가한다.
- N이 i로 나눠지면, 약수의 개수를 +2 증가시킨다.
- N이 완전 제곱수인 경우라면 약수 1개(제곱근)를 중복으로 센 것이기 때문에 1을 빼준다.
이를 바탕으로 약수의 개수를 세는 함수는 다음과 같이 작성될 수 있다.
def number_of_divisors(n):
c = 2 # 1, n
l = n ** 0.5
k = 2
while k <= l:
if n % k == 0:
c += 2
k += 1
if int(l) ** 2 == n:
c -= 1
return c
이제 답을 구해보자. 삼각수의 일반항을 계속 구해서, 약수의 개수를 센 다음 500개를 넘어서는 시점에 루프를 중단하고 결과를 출력한다. (참고로 N의 제곱근까지에 대해서 나눠지는 수마다 2개씩 약수를 챙길 수 있으므로, 500개 이상의 약수를 가지려면 기본적으로 250**2 보다는 커야 할 것이다. )
def main():
i = 1
while True:
if number_of_divisors(n) > 500:
print(n)
break
i += 1
%time main()
# 76576500
# Wall time: 12.6 s결과값 자체가 엄청 크기 때문에 결과를 찾는데까지 12초 이상 걸렸다. 아주 큰 수의 약수의 개수를 찾는 보다 빠른 방법은 없을까?
개선 – 약수의 개수를 계산하기
사실 약수의 개수를 세는 방법은 위의 방법을 쓰는 것이 일반적으로 너무 크지 않은 숫자 범위에서는 가장 빠르게 동작한다. 하지만 이 문제에서는 답이 7천6백만을 넘는 큰 수이고, 이렇게 큰 수들에 대해서는 제곱근 범위만큼 루프를 돌더라도 몇 천번의 루프를 돌아야 하기 때문에 속도가 나지 않는다. 따라서 여기서 성능을 높이고 싶다면, 약수의 개수를 계산하는 부분을 손봐서 개선해야 한다.
어떤 자연수의 약수의 개수를 구하는 다른 방법으로는 해당 자연수를 소인수 분해하는 방법이 있다. 자연수 N을 소인수 분해한 결과가 
(p + 1)(q + 1)(r + 1)로 계산된다. 그렇다면 소인수 분해를 빠르게 할 수 있다면 그만큼 약수의 개수를 빠르게 계산해 낼 수 있을 것이다.
소인수분해하기
임의의 자연수 N을 소인수 분해하는 방법에 대해 고민해보자. 소인수 분해는 120을 
120 ==> [(2, 3), (3, 1), (5, 1)])
- N을 2로 나누어본다.
- N이 2로 나눠진다면 나눠질 때까지 계속 나눠본다. 더이상 2로 나눠지지 않을 때까지 나눈 후 몫을 다시 N으로 두고, 2로 몇 번 나누었는지를 기록한다.
- 이번에는 2가 아닌 3으로 1~2를 시행한다.
- 이렇게 다음 소인수를 계속 찾아가면서 N을 나누고, 소인수와 그 기록을 저장한다.
- N이 1이 되면 소인수 분해가 종료된다.
따라서 다음과 같이 소인수 분해 함수를 작성할 수 있다. N을 나눌 수 있는 가장 작은 소인수를 찾는 헬퍼함수를 동원하고 있다.
def factorize(n):
def divisor_finder(x):
"""x를 나눌 수 있는 가장 작은 제수를 찾는다."""
for k in (2, 3):
if n % k == 0:
return k
k, l = 5, n ** 0.5
while k <= l:
if n % k == 0:
return k
if n % (k+2) == 0:
return k + 2
k += 6
return n
result = []
while n > 1:
p = divisor_finder(n)
e = 0
while n % p is 0:
n //= p
e += 1
result.append((p, e))
return result
# 테스트
factorize(120)
# [(2, 3), (3, 1), (5, 1)]소인수분해를 할 수 있다면 지수에 1씩 더하여 그 곱을 취한 것이 약수의 개수가 되므로 다음과 같이 약수의 개수를 구하는 함수를 재작성할 수 있다.
from functools import reduce
def product(iterable):
return reduce(lambda x,y: x*y, iterable, 1)
def number_of_divisors(n):
g = [x[1] + 1 for x in factorize(n)]
return product(g)이제 새로 작성된 함수를 이용해서 다시 답을 구해보자. 그리고 시간이 얼마나 단축되는지도 확인해보자.
def main():
n = 1
while True:
x = n + ( n + 1) // 2
if number_of_divisors(n) >= 500:
print(x)
return
n += 1
# 76576500
# Wall time: 379 ms보너스 – Swift 풀이
다음은 Swift 버전의 풀이이다. (실행: http://swift.sandbox.bluemix.net/#/repl/599f88430485e2473783df9b)
import Foundation
func timeit(_ f:() ->() ) {
let s = Date()
f()
print("time: \(s.timeIntervalSinceNow * -1000)ms")
}
func factorize(_ x: Int) -> [(Int, Int)] {
var temp: [(Int, Int)] = []
var a = 1
var n = x
func divisor(over a: Int) -> Int {
switch a {
case (1...2):
if n % (a + 1) == 0 { return a + 1 }
return divisor(over: a + 1)
case 3:
if n % 5 == 0 { return 5 }
return divisor(over: 5)
default:
var k = a + 2
let ds = (k % 3 == 1) ? [4, 2] : [2, 4]
var i = 0
let l = Int(sqrt(Double(n) + 0.5))
while k <= l { if n % k == 0 { return k } k += ds[i] i = (i + 1) % 2 } return n } } while n > 1 {
a = divisor(over: a)
var e = 0
while n % a == 0 {
e += 1
n /= a
}
if e > 0 { temp.append((a, e))}
}
return temp
}
func numberOfDivisor(_ n:Int) -> Int {
let fs = factorize(n)
return fs.map{ $0.1 + 1 }.reduce(1, *)
}
timeit{
let tn: (Int) -> Int = { n in n * (n + 1) / 2 }
var x = 1
while true {
let t = tn(x)
if numberOfDivisor(t) >= 500 {
print(t)
break
}
x += 1
}
}