박치욱 님이 트위터에 올린 문제에 대한 솔루션 정리

0. 발아점 (재활용)

이전 포스트 임의의 숫자가 제곱수인지 빠르게 판별하는 법에 이어지는 포스트임.

1. 문제에 대한 나의 접근(실패)

\(a^2+b^2+4a^2b^2 = c^2\)

을 아래와 같이 변형한 뒤

\(a^2+b^2+2ab+4a^2b^2-2ab = c^2\)

아래와 같이 정리했다.

\((a+b)^2+2ab(2ab-1) = c^2\)

여기서 모든 제곱수(\(c^2\))에 대해 이러한 관계를 만족시키는 a, b를 찾는 거다.

하지만, a, b는 Brute-Force하게 루프를 돌려야 되는데, 효율성이 낮다.

(c가 100일 때 돌렸던 루프를 c가 1000일 때도 또 돌려야 함)

따라서 실패.

2. 치욱님 솔루션

https://twitter.com/chiw00k/status/217625138395480064

a, b의 최대값은 \(\sqrt {20} \cdot 10^9\)이다.

여기서, a를 1~최대값, b를 1~i^[각주:1] 까지 이중 루프를 돌린다.

이 때 \(a^2+b^2+4a^2 b^2\)을 계산해서 그 결과가 제곱수인가 판별하면 된다.

3. knauer0x 님의 홀수 배제 아이디어

https://twitter.com/knauer0x/status/217643729593442304

이전 포스트에서 언급했던 내용인데, 제곱수를 4로 나눈 나머지는 반드시 0 또는 1이다.

그런데, a, b가 모두 홀수라면 좌변의 결과를 4로 나눈 나머지는 반드시 2가 된다.

\(a^2+b^2+4a^2b^2\)

따라서, 루프를 돌릴 때 a, b가 모두 홀수인 경우는 배제하고 돌려도 된다.

4. Zaeku 님의 최적화

대규모의 Brute-Force 연산에서 연산시간을 줄이려면 물론 루프 자체의 숫자를 줄여야 된다.

그런데, 그것 만큼이나 중요한 것이 반복되는 연산의 연산횟수를 줄이는 것이다.

(이전 포스트 역시 그런 의미에서 쓴 것이다)

b의 루프에서 계속 반복적으로 사용되지만, a 루프에서 미리 계산 가능한 값이 둘 있다.

\(a^2\)과 \(4a^2\)이다.

이 두 값을 미리 계산해두면 b 루프에선 곱하기를 단 두 번만 할 수 있다.

그리고, 이전 포스트에서 언급한 제곱수 판별을 사용하면 연산 자체의 수를 최소화 할 수 있다.

5. 기타 최적화

a, b 루프에서 계속 \(a^2\), \(b^2\)을 계산하는데, 일정하게 증가하는 값이므로 더하기로 대체할 수 있다.

즉, \((a+1)^2-a^2 = 2a+1\), \((a+2)^2-a^2 = 4a+4\)를 사용하는 것이다.

더하기와 곱하기의 클럭 차이가 있기 때문에 조금이나마 시간 절감 효과를 볼 수 있다.

6. 결과

아래는 1~\(10^6\) 범위의 c에 대해 해를 찾으면서 위의 최적화를 모두 반영한 코드이다.

#include "stdafx.h"
#include <math.h>
#include <Windows.h>

#define BASENUM 64
bool bBase[BASENUM] = {false,};

void prepareToCheckSquare() {
    for (int i=0; i<BASENUM; i++) {
        bBase[i*i % BASENUM] = true;
    }
}

inline bool IsSquare(unsigned long long int &num, unsigned long long int &ret) {
    if (bBase[(int)(num & (BASENUM-1))]) {
        ret = (unsigned int)(sqrtl((long double)num)+0.5);
        return (ret*ret == num);
    } else return false;
}

int _tmain(int argc, _TCHAR* argv[])
{
    prepareToCheckSquare();

    DWORD dw = GetTickCount();

    unsigned long long int limit = 1000000LL;
    unsigned long long int limit1 = (unsigned long long int)(sqrt(20.0)*limit/10);
    unsigned long long int limit2 = limit*limit;
    unsigned long long int i, j, area;
    unsigned long long int isq, isq4, jsq, areasq;
    int count=0;

    for (i=1, isq=1, isq4=4; i<limit1; isq+=((i<<1)+1), isq4 = isq<<2, i++) {
        if (i & 1) {
            for (j=2, jsq=4; j<i; jsq+=((j<<2)+4), j+=2) {
                areasq = isq + jsq + isq4*jsq;
                if (areasq > limit2) break;

                if (IsSquare(areasq, area)) {
                    printf("#%d. a=%llu, b=%llu, c=%llu\n", ++count, i, j, area);
                }
            }
        } else {
            for (j=1, jsq=1; j<i; jsq+=((j<<1)+1), j++) {
                areasq = isq + jsq + isq4*jsq;
                if (areasq > limit2) break;

                if (IsSquare(areasq, area)) {
                    printf("#%d. a=%llu, b=%llu, c=%llu\n", ++count, i, j, area);
                }
            }
        }
    }

    dw = GetTickCount()-dw;
    printf("%u ms (%u sec)\n", dw, dw/1000);

    return 0;
}

위에 언급한 최적화를 모두 적용했으며, 수 초 이내에 답을 출력해준다.

7. 추가 고민거리

위의 코드는 unsigned long long int와 long double을 기반으로 작성되었다.
그런데 주어진 조건인 \(10^{10}\)까지 모두 계산하려면 64비트의 unsigned long long int로는 부족하다.
물론, 별도의 클래스를 하나 설계하는 방법이 있지만, 실행 효율에 악영향을 미칠 것이다.

솔루션을 더 찾아봐야겠다.

1. b의 범위를 i로 제한하는 것은 중복 솔루션의 비효율성 제거하기 위한 것임 [본문으로]