반응형

 

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

 

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

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

 

denomNminuend 형태의 항의 다음 항을 계산하는 식만 찾으면 다음은 간단하다.

 

denomNminuend=denom(N+minuend)Nminuend2=N+minuend(Nminuend2)/denom가 되는데, 이 식의 정수해가 다음 항이 된다.

 

즉, an=N+minuend(Nminuend2)/denom 가 되는데, 여기서 N의 정수 파트만 계산하면 나머지는 간단하다.

 

denom(N+minuend)Nminuend2an+denomN+denomminuendan(Nminuend2)Nminuend2=an+N(an(Nminuend2)/denomminuend)(Nminuend2)/denom로 정리된다.

 

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

 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
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;
}

 

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

 

1
2
3
4
5
6
7
8
9
10
11
12
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;
  }
}

 

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

 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
#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);
}

 

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

 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
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 정도 소요됐음....

 

 

반응형

공유하기

facebook twitter kakaoTalk kakaostory naver band