제곱근의 연분수 전개시 반복되는 숫자가 홀수 개인 경우는?

 

오랜만에 풀어보는 오일러 프로젝트 하나.

 

제곱근을 연분수로 전개하면 나오는 값은 순환하게 되는데, 이 값이 홀수 개인 경우를 찾는 문제다.

이 문제는 프로그래밍 난이도 보다는 수학적인 의미를 이해하는 쪽이 더 어려운 문제.

 

\( \frac { denom } { \sqrt {N} - minuend } \) 형태의 항의 다음 항을 계산하는 식만 찾으면 다음은 간단하다.

 

\( \frac {denom} { \sqrt {N} - minuend } = \frac {denom ( \sqrt{N} + minuend ) } { N - minuend ^ 2 } = \frac { \sqrt {N} + minuend } {(N - minuend ^ 2) / denom} \)가 되는데, 이 식의 정수해가 다음 항이 된다.

 

즉, \( a_{n} = \lfloor \frac { \sqrt {N} + minuend } {(N - minuend ^ 2) / denom} \rfloor \) 가 되는데, 여기서 \( \sqrt {N} \)의 정수 파트만 계산하면 나머지는 간단하다.

 

\( \frac {denom ( \sqrt{N} + minuend ) } { N - minuend ^ 2 } \)는 \( a_{n} + \frac {denom \sqrt{N}+denom \cdot minuend - a_{n}(N - minuend ^ 2 )}{ N - minuend ^ 2 } = a_{n} + \frac {\sqrt{N} - ( a_{n} (N-minuend ^ 2 )/denom - minuend)}{( N - minuend ^ 2 ) / denom } \)로 정리된다.

 

해를 찾는 과정에서 \( \sqrt{N} \)는 정수 파트만 필요하므로 아래 함수를 사용하면 더 효율적이다.

 

unsigned isqrt(unsigned num) noexcept {
  unsigned res = 0;
  unsigned bit = 1 << (sizeof(unsigned) * 8 - 2); // The second-to-top bit is set: 1 << 30 for 32 bits

  // "bit" starts at the highest power of four <= the argument.
  while (bit > num) { bit >>= 2; }

  while (bit != 0) {
    if (num >= res + bit) {
      num -= res + bit;
      res = (res >> 1) + bit;
    }
    else {
      res >>= 1;
    }
    bit >>= 2;
  }
  return res;
}

 

각 항의 다음 값을 계산하는 함수는 다음과 같이 간략하게 쓸 수 있다.

 

void GetSquareRootVals(const unsigned N, const unsigned isqrtN, const unsigned denom, const unsigned minuend, unsigned& aNew, unsigned& denomNew, unsigned& minuendNew)
{
  if (isqrtN * isqrtN == N) {
    denomNew = 1;
    aNew = minuendNew = 0;
  }
  else {
    denomNew = (N - minuend * minuend) / denom;
    aNew = (isqrtN + minuend) / denomNew;
    minuendNew = aNew * denomNew - minuend;
  }
}

 

모든 내용을 다 갈아넣은(?) 결과는 아래와 같다.

 

#include <stdio.h>
#include <vector>

using namespace std;

unsigned isqrt(unsigned num) noexcept {
  unsigned res = 0;
  unsigned bit = 1 << (sizeof(unsigned) * 8 - 2); // The second-to-top bit is set: 1 << 30 for 32 bits

  // "bit" starts at the highest power of four <= the argument.
  while (bit > num) { bit >>= 2; }

  while (bit != 0) {
    if (num >= res + bit) {
      num -= res + bit;
      res = (res >> 1) + bit;
    }
    else {
      res >>= 1;
    }
    bit >>= 2;
  }
  return res;
}

void GetSquareRootVals(const unsigned N, const unsigned isqrtN, const unsigned denom, const unsigned minuend, unsigned& aNew, unsigned& denomNew, unsigned& minuendNew)
{
  if (isqrtN * isqrtN == N) {
    denomNew = 1;
    aNew = minuendNew = 0;
  }
  else {
    denomNew = (N - minuend * minuend) / denom;
    aNew = (isqrtN + minuend) / denomNew;
    minuendNew = aNew * denomNew - minuend;
  }
}

struct result {
  unsigned a, d, m;
};

bool IsPeriodOdd(const unsigned N, const bool printOddOnly) noexcept {
  unsigned sqrN = isqrt(N);

  if (sqrN * sqrN == N) {
    if (!printOddOnly) {
      printf("sqrt(%u) = [%u], period = 0\n", N, sqrN);
    }
    return false;
  }

  unsigned a, d, m;
  vector<struct result> vResult;
  struct result resultN = { sqrN, 1, sqrN };

  vResult.push_back(resultN);
  for (;;) {
    GetSquareRootVals(N, sqrN, resultN.d, resultN.m, a, d, m);
    struct result resultT = { a, d, m };
    if (vResult.size() > 1 && !memcmp(&vResult[1], &resultT, sizeof(resultT))) {
      break;
    }

    vResult.push_back(resultT);
    resultN = resultT;
  }

  int sz = (int)vResult.size();

  if (!printOddOnly || ((sz - 1) & 1)) {
    printf("sqrt(%u) = [%u;(", N, vResult[0].a);

    for (int i = 1; i < sz; ++i) {
      printf("%u", vResult[i].a);
      printf((i == (sz - 1)) ? ")" : ",");
    }
    printf("], period = %d\n", sz - 1);
  }

  return ((sz - 1) & 1);
}

int main()
{
  int count = 0;
  for (unsigned i = 1; i <= 10000; ++i) {
    if (IsPeriodOdd(i, true)) {
      ++count;
    }
  }

  printf("odd count: %d\n", count);
}

 

결과의 일부는 아래와 같다.

 

sqrt(2) = [1;(2)], period = 1
sqrt(5) = [2;(4)], period = 1
sqrt(10) = [3;(6)], period = 1
sqrt(13) = [3;(1,1,1,1,6)], period = 5
sqrt(17) = [4;(8)], period = 1
sqrt(26) = [5;(10)], period = 1
sqrt(29) = [5;(2,1,1,2,10)], period = 5
sqrt(37) = [6;(12)], period = 1
sqrt(41) = [6;(2,2,12)], period = 3
sqrt(50) = [7;(14)], period = 1
sqrt(53) = [7;(3,1,1,3,14)], period = 5
sqrt(58) = [7;(1,1,1,1,1,1,14)], period = 7
sqrt(61) = [7;(1,4,3,1,2,2,1,3,4,1,14)], period = 11
sqrt(65) = [8;(16)], period = 1
sqrt(73) = [8;(1,1,5,5,1,1,16)], period = 7
sqrt(74) = [8;(1,1,1,1,16)], period = 5
sqrt(82) = [9;(18)], period = 1
sqrt(85) = [9;(4,1,1,4,18)], period = 5
sqrt(89) = [9;(2,3,3,2,18)], period = 5
sqrt(97) = [9;(1,5,1,1,1,1,1,1,5,1,18)], period = 11
...하략...

 

코드 자체는 단순하기 때문에 굳이 더 최적화할 필요는 없어보인다.

내 PC에서 화면 출력 없이 10,000까지 계산하는데 15ms 정도 소요됐음....