프로젝트 오일러 141번

A positive integer, \(n\), is divided by \(d\) and the quotient and remainder are \(q\) and \(r\) respectively. In addition \(d, \;q\) and \(r\) are consecutive positive integer terms in a geometric sequence, but not necessarily in that order.

For example, \(58\) divided by \(6\) has quotient \(9\) and remainder \(4\). It can also be seen that \(4, \;6,\;9\) are consecutive terms in a geometric sequence \( \left ( {\rm common \;ratio}\; \dfrac{3}{2} \right )\).

We will call such numbers, \(n\), progressive.

Some progressive number, such as \(9\) and \(10404 = 102^2\), happen to also be perfect squares.

The sum of all progressive perfect squares below one hundred thousand is \(124657\).

Find the sum of all progressive perfect squares below one trillion \(\left ( 10 ^{12} \right )\).

---

자연수 \(n\) 을 자연수 \(d\) 로 나누었을 때의 몫과 나머지를 각각 \(q, \; r\) 이라고 하자. 어떤 경우에는 \(d, \; q,\; r\) 이 등비수열을 이루는 것이 가능하다. (순서가 꼭 \(d, \;q, \;r\) 일 필요는 없다.)

예를 들어, \(58\) 을 \(6\) 으로 나누면 몫이 \(9\) 이고 나머지가 \(4\) 인데, \(4, \; 6, \; 9\) 는 공비가 \(\dfrac{3}{2}\) 인 등비수열을 이룬다. 이런 경우에 \(n\) (여기서는 \(58\)) 을 progressive 하다고 한다. 

또한 \(9\) 나 \(10404=102^2\) 과 같은 progressive들은 동시에 완전제곱수이기도 한데 이런 수들을 progressive 완전제곱수라고 한다.

\(10^5\) 미만의 수 중 progressive 완전제곱수들의 총합은 \(124657\) 이다.

그렇다면 \(10^12\) 미만의 수 중 progressive 완전제곱수들의 총합은 얼마일까? 

등비수열을 이루는 세 수를 다음과 같이 생각할 있다. \[a, \;\; ar, \;\; ar^2\] 이때, 이 세 수를 이용하여 \(n\) 을 만드는 방법은 \(ar^2 \times ar + a\) 와 \(ar^2 \times a + ar\) 이 가능한데, 이 중 \(ar^2 \times a + ar = ar(ar+1)\) 이 되므로 절대로 완전제곱수가 될 수 없다. 따라서 \(n= ar^2 \times ar + a\) 가 되어야 한다. 

또한 \(a, \; ar, \; ar^2\) 중 최소의 수를 \(a\) 라고 한다면 \(r>1\) 이 되어야 하고, 문제에서 봤듯이 \(r\) 은 정수가 아닌 유리수도 가능하므로, \(r\) 을 다음과 같이 표현할 수 있다. \[r=\dfrac{q}{p} \;\; (단, \; p, \;q\; 는\; 서로소인 \; 자연수,\;\; q>p \ge 1)\] 서로소인 조건을 두는 이유는 \(p=1\) 인 경우와 중복되게 하지 않기 위해서다.

이제 \(r= \dfrac{q}{p}\) 로 바꾸면 세 수는 \[a, \;\; a\times \dfrac{q}{p}, \;\; a \times \dfrac{q^2}{p^2}\] 이 된다. 그런데 이 세수가 모두 자연수가 되어야 하므로 \(a\) 는 \(p^2\) 의 배수가 되어야 한다. 따라서 \(a=kp^2\) (단, \(k\) 는 자연수)이라고 할 수 있다. 결국 세 수는 다음과 같은 형태를 갖게 된다. \[kp^2, \;\; kpq, \;\; kq^2\] 위 세 수로 만들 수 있는 \(n\) 은 \[n=kpq^3 + kp^2\] 이 된다. 

이제 다음과 세 개의 루프를 돌리면서 완전제곱수가 되는 경우만 더해주면 된다. 

limit = 10 ** 12
result = []
for q in range(2, 10000):
    for p in range(1, q):
        if p ** 2 * (q ** 3 + 1) > limit or gcd(q, p) > 1:
            continue
        k = 1
        n = q ** 3 * p + p ** 2
        while n < limit:
            if issquare(n):
                if n not in result:
                    result.append(n)
            k += 1
            n = q ** 3 * p * k ** 2 + k * p ** 2
print(sum(result))

\(q\) 의 범위를 \(10^4\) 으로 하는 이유는 \(p=1, \; k=1\) 인 상황에서 \(n\) 값을 기준으로 결정한 것이다. 또한 함수 gcd 는 최대공약수를 구해주는 함수, issuqare 는 완전제곱수를 판별해주는 함수이다. 

다소 시간이 오래 걸리기는 하지만 그래도 이 방법이 생각해 낸 가장 좋은 방법인듯...

혹시 더 나은 방법이 있다면 공유 !!!