题目64：10000以下的连分数中有多少拥有奇数周期？

欧拉计划

Odd period square roots

All square roots are periodic when written as continued fractions and can be written in the form:



For example, let us consider :



If we continue we would get the following expansion:



The process can be summarised as follows:



It can be seen that the sequence is repeating. For conciseness, we use the notation , to indicate that the block (1,3,1,8) repeats indefinitely.

The first ten continued fraction representations of (irrational) square roots are:



Exactly four continued fractions, for N ≤ 13, have an odd period.

How many continued fractions for N ≤ 10000 have an odd period?

题目：

所有的平方根写作连分数的时候都是周期性的，并且可以写成如下形式：



例如 :



如果继续下去我们得到如下展开式：



这个过程可以被总结如下：



可以看出这个序列是重复的。简明起见，我们用 , 来表示 (1,3,1,8) 这个部分不断重复出现。

前十个无理平方根的连续分数表示为：

√2=[1;(2)], 周期=1
√3=[1;(1,2)], 周期=2
√5=[2;(4)], 周期=1
√6=[2;(2,4)], 周期=2
√7=[2;(1,1,1,4)], 周期=4
√8=[2;(1,4)], 周期=2
√10=[3;(6)], 周期=1
√11=[3;(3,6)], 周期=2
√12= [3;(2,6)], 周期=2
√13=[3;(1,1,1,1,6)], 周期=5

可以看出对于 N ≤ 13, 有四个连分数的周期是奇数。

对于 N ≤ 10000，有多少连分数的周期是奇数？

笑脸对世界

刚刚开始学习，看着好厉害的样子:lol:

迷雾少年

mark下

jerryxjr1220

这题的思路就是找出除整数部分后面的无理数从第几位开始循环。
a是整数部分+无理数部分，d是整数部分，c就是无理数部分。
依次存入列表，然后比较。开始循环了就记录下来。
代码如下：

1. 
   
2. import math
   
3. def gao(n):
   
4.         m = int(math.sqrt(n))
   
5.         if m * m == n:
   
6.                 return 0
   
7.         dic = {}
   
8.         d = m
   
9.         a = n
   
10.         b = -m
    
11.         c = 1
    
12.         len = 0
    
13.         while True:
    
14.                 nc =  a - b * b
    
15.                 nc /= c
    
16.                 nd = int((math.sqrt(a) - b) / nc)
    
17.                 nb = -b - nd * nc
    
18.                 t = (nb, nc, nd)
    
19.                 b = nb
    
20.                 c = nc
    
21.                 d = nd
    
22.                 if not t in dic:
    
23.                         len += 1
    
24.                         dic[t] = 1
    
25.                 else:
    
26.                         break
    
27.                        
    
28.         return len
    
29. ans = 0
    
30. for i in range(1, 10001):
    
31.         k = gao(i)
    
32.         if k % 2 != 0:
    
33.                 ans += 1
    
34. print (ans)

复制代码

jerryxjr1220

结果1322

jerryxjr1220

另外一个暴力的的方法，用decimal库，把精度提高到250位以上（200位小数，仍有3个数有误差）
同样可以得到结果......1322......7秒多，速度也还行，简单粗暴，不用考虑复杂的数学方程式

1. import decimal as dc
   
2. dc.getcontext().prec = 250
   
3. def root(n):
   
4.     blist = []
   
5.     n0 = dc.Decimal(i).sqrt()
   
6.     a0 = int(n0)
   
7.     b0 = n0 - a0
   
8.     n1 = 1/b0
   
9.     a1 = int(n1)
   
10.     b1 = n1 - a1
    
11.     while str(b1)[:11] not in blist:
    
12.         blist.append(str(b1)[:11])
    
13.         n1 = 1/b1
    
14.         a1 = int(n1)
    
15.         b1 = n1 - a1
    
16.     return len(blist)
    
17. 
    
18. count = 0
    
19. for i in range(2,10000):
    
20.     if i**0.5 != int(i**0.5):
    
21.         if root(i) % 2 == 1:
    
22.             count += 1
    
23. print (count)

复制代码

najin

用的matlab
浮点数真是头疼，
    1322      

时间已过 0.413237 秒。
>>

永恒的蓝色梦想

C++ 43ms

1. #include<iostream>
   
2. #include<unordered_map>
   
3. #include<cmath>
   
4. 
   
5. 
   
6. 
   
7. constexpr int gcd(int a, int b, int c)noexcept {
   
8.     int t = 0;
   
9. 
   
10.     while (b) {
    
11.         t = a % b;
    
12.         a = b;
    
13.         b = t;
    
14.     }
    
15. 
    
16.     while (c) {
    
17.         t = a % c;
    
18.         a = c;
    
19.         c = t;
    
20.     }
    
21. 
    
22.     return a;
    
23. }
    
24. 
    
25. 
    
26. 
    
27. struct Quad {
    
28.     //(a * sqrt(b) + c) / d
    
29.     const short a, b, c, d;
    
30. };
    
31. 
    
32. 
    
33. 
    
34. int main() {
    
35.     using namespace std;
    
36.     ios::sync_with_stdio(false);
    
37. 
    
38.     short temp, a, b, c, d, sqrt_of_b;
    
39.     unsigned int* ptr;
    
40.     unsigned int count = 0, index;
    
41.     double sqrt_;
    
42.     unordered_map<size_t, unsigned int> indecies;
    
43. 
    
44. 
    
45.     for (b = 1; b <= 10000; b++) {
    
46.         sqrt_ = sqrt(b);
    
47. 
    
48.         if (fmod(sqrt_, 1.0)) {
    
49.             sqrt_of_b = sqrt_;
    
50.             a = d = 1;
    
51.             c = -sqrt_of_b;
    
52. 
    
53.             for (index = 1;; index++) {
    
54.                 ptr = &indecies[*(unsigned long long*)(&Quad{a, b, c, d})];
    
55. 
    
56.                 if (*ptr) {
    
57.                     count += (index - *ptr) & 1;
    
58.                     indecies.clear();
    
59.                     break;
    
60.                 }
    
61. 
    
62.                 *ptr = index;
    
63. 
    
64.                 temp = a * a * b - c * c;
    
65.                 a *= d;
    
66.                 c *= -d;
    
67.                 d = temp;
    
68. 
    
69.                 temp = gcd(a, c, d);
    
70.                 a /= temp;
    
71.                 c /= temp;
    
72.                 d /= temp;
    
73.                 c -= (sqrt_of_b * a + c) / d * d;
    
74.             }
    
75.         }
    
76.     }
    
77. 
    
78. 
    
79.     cout << count << endl;
    
80.     return 0;
    
81. }

复制代码

debuggerzh

和求无限小数的循环节问题很像
1322

1. #include<iostream>
   
2. #include<cmath>
   
3. #include<algorithm>
   
4. using namespace std;
   
5. 
   
6. int sq[100] = {0};  //平方数表
   
7. 
   
8. void init(){    //打平方数表
   
9.     for (int i = 2; i <= 101; i++)  sq[i-2] = i*i;
   
10. }
    
11. 
    
12. bool issquare(int n){   //判断平方数
    
13.     if (binary_search(sq, sq+100, n)) return true;
    
14.     return false;
    
15. }
    
16. 
    
17. int gcd(int x, int y){
    
18.     if (x % y == 0) return y;
    
19.     return gcd(y, x % y);
    
20. }
    
21. 
    
22. int period(int n){
    
23.     int deno = 1;   //分母
    
24.     int cons = -(int)sqrt(n);   //小数部分的有理项系数
    
25.    
    
26.     int kount = 0;
    
27.     do  //模拟手算过程
    
28.     {
    
29.         int t = n - cons*cons;
    
30.         double revert = deno / (sqrt(n) + cons);
    
31.         
    
32.         deno = t / gcd(t, deno);
    
33.         cons = (int) floor( ((revert - (int)revert) * deno - sqrt(n) + 0.5));
    
34.         
    
35.         kount++;
    
36.     } while (deno != 1);
    
37.     return kount;
    
38. }
    
39. 
    
40. int main(int argc, char const *argv[])
    
41. {
    
42.     init();
    
43.     int kount = 0;
    
44.     for (int n = 2; n <= 10000; n++){
    
45.         if (issquare(n))    continue;   //平方数不循环
    
46.         if (period(n) % 2 == 1) kount++;
    
47.     }
    
48.     cout << kount << endl;
    
49.     getchar();
    
50.     return 0;
    
51. }

复制代码

guosl

1. /*
   
2. 答案：1322
   
3. 耗时：0.0162126秒（单核）
   
4.       0.0035869秒（六核）
   
5. */
   
6. #include <iostream>
   
7. #include <vector>
   
8. #include <algorithm>
   
9. #include <string>
   
10. #include <cmath>
    
11. #include <omp.h>
    
12. using namespace std;
    
13. 
    
14. int gcd(int x, int y)//计算x,y的最大公因数
    
15. {
    
16.   if (y == 0)
    
17.     return x;
    
18.   int z = x % y;
    
19.   while (z != 0)
    
20.   {
    
21.     x = y;
    
22.     y = z;
    
23.     z = x % y;
    
24.   }
    
25.   if (y < 0)
    
26.     return (-y);
    
27.   return y;
    
28. }
    
29. 
    
30. struct stItem //对应分式 (a*sqrt(D)+b)/c,d=sqrt(D)
    
31. {
    
32.   int a, b, c;
    
33.   stItem(void)
    
34.   {
    
35.     a = 0;
    
36.     b = 0;
    
37.     c = 1;
    
38.   }
    
39.   stItem(int x, int y, int z)
    
40.   {
    
41.     c = gcd(x, y);
    
42.     c = gcd(c, z);
    
43.     a = x / c;
    
44.     b = y / c;
    
45.     c = z / c;
    
46.   }
    
47.   stItem(const stItem& x)
    
48.   {
    
49.     memcpy_s(this, sizeof(stItem), &x, sizeof(stItem));
    
50.   }
    
51.   const stItem& operator=(const stItem& x)
    
52.   {
    
53.     memcpy_s(this, sizeof(stItem), &x, sizeof(stItem));
    
54.     return *this;
    
55.   }
    
56.   bool operator==(const stItem& x) const
    
57.   {
    
58.     return (a == x.a && b == x.b && c == x.c);
    
59.   }
    
60. };
    
61. 
    
62. /*
    
63. 返回0还没有找到周期，否则就是周期的长度
    
64. a,b,c进入时为本次的三元组，返回时是下次的三元组
    
65. vItems 记录三元组，用于查找周期
    
66. */
    
67. int f(int nD, int rt, int& a, int& b, int& c, vector<stItem>& vItems)
    
68. {
    
69.   int q = int((a * rt + b) / c);
    
70.   int a1 = a * c;
    
71.   int b1 = c * (c * q - b);
    
72.   int c1 = a * a * nD - (c * q - b) * (c * q - b);
    
73.   stItem it(a1, b1, c1);
    
74.   a = it.a;
    
75.   b = it.b;
    
76.   c = it.c;
    
77.   vector<stItem>::iterator itr = find(vItems.begin(), vItems.end(), it);
    
78.   if (itr == vItems.end())
    
79.   {
    
80.     vItems.push_back(it);
    
81.     return 0;
    
82.   }
    
83.   return int(vItems.end() - itr);
    
84. }
    
85. 
    
86. int main(void)
    
87. {
    
88.   double t = omp_get_wtime();
    
89.   int nCount = 0;
    
90. #pragma omp parallel for reduction(+:nCount)
    
91.   for (int i = 2; i <= 10000; ++i)
    
92.   {
    
93.     int j = int(sqrt(double(i)));
    
94.     if (j * j == i)
    
95.       continue;
    
96.     vector<stItem> vItems;
    
97.     int a = 1, b = 0, c = 1;//刚开始只有sqrt(i)
    
98.     stItem it(a, b, c);
    
99.     vItems.push_back(it);
    
100.     while (true)//展开sqrt(i)
     
101.     {
     
102.       int k = f(i, j, a, b, c, vItems);
     
103.       if (k > 0)//直到出现周期
     
104.       {
     
105.         nCount += (k & 1);//计数奇数周期
     
106.         break;
     
107.       }
     
108.     }
     
109.   }
     
110.   t = omp_get_wtime() - t;
     
111.   cout << nCount << endl << t << endl;
     
112.   return 0;
     
113. }
     

复制代码

B1tetheDust

1. import time as t
   
2. import decimal as dc
   
3. 
   
4. start = t.perf_counter()
   
5. count = 0
   
6. dc.getcontext().prec = 250
   
7. num_list = [n for n in range(10001)]
   
8. nums = []
   
9. for n in range(101):
   
10.     num_list.remove(n ** 2)
    
11. for n in num_list:
    
12.     radicand = dc.Decimal(n).sqrt()
    
13.     first_num = int(radicand)
    
14.     decimal_list = [first_num]
    
15.     while True:
    
16.         radicand = dc.Decimal(1) / dc.Decimal(radicand - int(radicand))
    
17.         if str(radicand)[:11] not in decimal_list:
    
18.             decimal_list.append(str(radicand)[:11])
    
19.         else:
    
20.             len_decimal = len(decimal_list)
    
21.             if (len_decimal - decimal_list.index(str(radicand)[:11])) % 2:
    
22.                 count += 1
    
23.             break
    
24. 
    
25. print(count)
    
26. print("It costs %f s" % (t.perf_counter() - start))
    

复制代码

1322
It costs 2.698963 s

TLM

不使用浮点数开方进行求解
python
1322
Time:1.935354232788086s
但问题的关键是为什么一次取整之后，就可以循环了

1. import time
   
2. st=time.time()
   
3. def getGYS(u,v):
   
4.     # 获取公约数
   
5.     if u<v:
   
6.         u,v=v,u
   
7.     while v>0:
   
8.         u,v=v,u%v
   
9.     return u
   
10. def get(N,a0=1):
    
11.     ans=[]
    
12.     # 第一次取整
    
13.     while (a0+1)**2<N:
    
14.         a0=a0+1
    
15.     if (a0+1)**2==N:
    
16.         return [a0+1],a0
    
17.     ans.append(a0)
    
18.     a1=-a0
    
19.     a2=N-a0**2
    
20.     a0=1
    
21.     chushi=(a0,a1,a2)
    
22.     while True:
    
23.         # 循环取整
    
24.         '''
    
25.         a3<(a0√N-a1)/a2<a3+1
    
26.         转化为(a3*a2+a1)**2<a0**2*N<((a3+1)*a2+a1)**2
    
27.         获取a3
    
28.         '''
    
29.         a3=0
    
30.         while ((a3+1)*a2+a1)**2<a0**2*N:
    
31.             a3=a3+1
    
32.         ans.append(a3)
    
33.         # 重新转化为(a0√N-a1)/a2的形式
    
34.         a0,a1,a2=a2*a0,-a2*a0*(a1+a2*a3),a0**2*N-(a1+a2*a3)**2
    
35.         gys=getGYS(a2,a0)
    
36.         a0,a1,a2=a0//gys,a1//gys,a2//gys
    
37.         if (a0,a1,a2)==chushi:
    
38.             return ans,a0
    
39. # N=13
    
40. # print(get(N))
    
41. n=0
    
42. N=10000
    
43. a0=1
    
44. for i in range(2,N+1):
    
45.     ans,a0=get(i,a0)
    
46.     if len(ans)%2==0:
    
47.         n=n+1
    
48. print(n)
    
49. print('Time:{}s'.format(time.time()-st))

复制代码