题目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就是无理数部分。
依次存入列表，然后比较。开始循环了就记录下来。
代码如下：

import math
def gao(n):
        m = int(math.sqrt(n))
        if m * m == n:
                return 0
        dic = {}
        d = m
        a = n
        b = -m
        c = 1
        len = 0
        while True:
                nc =  a - b * b
                nc /= c
                nd = int((math.sqrt(a) - b) / nc)
                nb = -b - nd * nc
                t = (nb, nc, nd)
                b = nb
                c = nc
                d = nd
                if not t in dic:
                        len += 1
                        dic[t] = 1
                else:
                        break
                        
        return len
ans = 0
for i in range(1, 10001):
        k = gao(i)
        if k % 2 != 0:
                ans += 1
print (ans)

jerryxjr1220

结果1322

jerryxjr1220

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

import decimal as dc
dc.getcontext().prec = 250
def root(n):
    blist = []
    n0 = dc.Decimal(i).sqrt()
    a0 = int(n0)
    b0 = n0 - a0
    n1 = 1/b0
    a1 = int(n1)
    b1 = n1 - a1
    while str(b1)[:11] not in blist:
        blist.append(str(b1)[:11])
        n1 = 1/b1
        a1 = int(n1)
        b1 = n1 - a1
    return len(blist)

count = 0
for i in range(2,10000):
    if i**0.5 != int(i**0.5):
        if root(i) % 2 == 1:
            count += 1
print (count)

najin

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

时间已过 0.413237 秒。
>>

永恒的蓝色梦想

C++ 43ms

#include<iostream>
#include<unordered_map>
#include<cmath>



constexpr int gcd(int a, int b, int c)noexcept {
    int t = 0;

    while (b) {
        t = a % b;
        a = b;
        b = t;
    }

    while (c) {
        t = a % c;
        a = c;
        c = t;
    }

    return a;
}



struct Quad {
    //(a * sqrt(b) + c) / d
    const short a, b, c, d;
};



int main() {
    using namespace std;
    ios::sync_with_stdio(false);

    short temp, a, b, c, d, sqrt_of_b;
    unsigned int* ptr;
    unsigned int count = 0, index;
    double sqrt_;
    unordered_map<size_t, unsigned int> indecies;


    for (b = 1; b <= 10000; b++) {
        sqrt_ = sqrt(b);

        if (fmod(sqrt_, 1.0)) {
            sqrt_of_b = sqrt_;
            a = d = 1;
            c = -sqrt_of_b;

            for (index = 1;; index++) {
                ptr = &indecies[*(unsigned long long*)(&Quad{a, b, c, d})];

                if (*ptr) {
                    count += (index - *ptr) & 1;
                    indecies.clear();
                    break;
                }

                *ptr = index;

                temp = a * a * b - c * c;
                a *= d;
                c *= -d;
                d = temp;

                temp = gcd(a, c, d);
                a /= temp;
                c /= temp;
                d /= temp;
                c -= (sqrt_of_b * a + c) / d * d;
            }
        }
    }


    cout << count << endl;
    return 0;
}

debuggerzh

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

#include<iostream>
#include<cmath>
#include<algorithm>
using namespace std;

int sq[100] = {0};  //平方数表

void init(){    //打平方数表
    for (int i = 2; i <= 101; i++)  sq[i-2] = i*i;
}

bool issquare(int n){   //判断平方数
    if (binary_search(sq, sq+100, n)) return true;
    return false;
}

int gcd(int x, int y){
    if (x % y == 0) return y;
    return gcd(y, x % y);
}

int period(int n){
    int deno = 1;   //分母
    int cons = -(int)sqrt(n);   //小数部分的有理项系数
    
    int kount = 0;
    do  //模拟手算过程
    {
        int t = n - cons*cons;
        double revert = deno / (sqrt(n) + cons);
        
        deno = t / gcd(t, deno);
        cons = (int) floor( ((revert - (int)revert) * deno - sqrt(n) + 0.5));
        
        kount++;
    } while (deno != 1);
    return kount;
}

int main(int argc, char const *argv[])
{
    init();
    int kount = 0;
    for (int n = 2; n <= 10000; n++){
        if (issquare(n))    continue;   //平方数不循环
        if (period(n) % 2 == 1) kount++;
    }
    cout << kount << endl;
    getchar();
    return 0;
}

guosl

/*
答案：1322
耗时：0.0162126秒（单核）
      0.0035869秒（六核）
*/
#include <iostream>
#include <vector>
#include <algorithm>
#include <string>
#include <cmath>
#include <omp.h>
using namespace std;

int gcd(int x, int y)//计算x,y的最大公因数
{
  if (y == 0)
    return x;
  int z = x % y;
  while (z != 0)
  {
    x = y;
    y = z;
    z = x % y;
  }
  if (y < 0)
    return (-y);
  return y;
}

struct stItem //对应分式 (a*sqrt(D)+b)/c,d=sqrt(D)
{
  int a, b, c;
  stItem(void)
  {
    a = 0;
    b = 0;
    c = 1;
  }
  stItem(int x, int y, int z)
  {
    c = gcd(x, y);
    c = gcd(c, z);
    a = x / c;
    b = y / c;
    c = z / c;
  }
  stItem(const stItem& x)
  {
    memcpy_s(this, sizeof(stItem), &x, sizeof(stItem));
  }
  const stItem& operator=(const stItem& x)
  {
    memcpy_s(this, sizeof(stItem), &x, sizeof(stItem));
    return *this;
  }
  bool operator==(const stItem& x) const
  {
    return (a == x.a && b == x.b && c == x.c);
  }
};

/*
返回0还没有找到周期，否则就是周期的长度
a,b,c进入时为本次的三元组，返回时是下次的三元组
vItems 记录三元组，用于查找周期
*/
int f(int nD, int rt, int& a, int& b, int& c, vector<stItem>& vItems)
{
  int q = int((a * rt + b) / c);
  int a1 = a * c;
  int b1 = c * (c * q - b);
  int c1 = a * a * nD - (c * q - b) * (c * q - b);
  stItem it(a1, b1, c1);
  a = it.a;
  b = it.b;
  c = it.c;
  vector<stItem>::iterator itr = find(vItems.begin(), vItems.end(), it);
  if (itr == vItems.end())
  {
    vItems.push_back(it);
    return 0;
  }
  return int(vItems.end() - itr);
}

int main(void)
{
  double t = omp_get_wtime();
  int nCount = 0;
#pragma omp parallel for reduction(+:nCount)
  for (int i = 2; i <= 10000; ++i)
  {
    int j = int(sqrt(double(i)));
    if (j * j == i)
      continue;
    vector<stItem> vItems;
    int a = 1, b = 0, c = 1;//刚开始只有sqrt(i)
    stItem it(a, b, c);
    vItems.push_back(it);
    while (true)//展开sqrt(i)
    {
      int k = f(i, j, a, b, c, vItems);
      if (k > 0)//直到出现周期
      {
        nCount += (k & 1);//计数奇数周期
        break;
      }
    }
  }
  t = omp_get_wtime() - t;
  cout << nCount << endl << t << endl;
  return 0;
}

B1tetheDust

import time as t
import decimal as dc

start = t.perf_counter()
count = 0
dc.getcontext().prec = 250
num_list = [n for n in range(10001)]
nums = []
for n in range(101):
    num_list.remove(n ** 2)
for n in num_list:
    radicand = dc.Decimal(n).sqrt()
    first_num = int(radicand)
    decimal_list = [first_num]
    while True:
        radicand = dc.Decimal(1) / dc.Decimal(radicand - int(radicand))
        if str(radicand)[:11] not in decimal_list:
            decimal_list.append(str(radicand)[:11])
        else:
            len_decimal = len(decimal_list)
            if (len_decimal - decimal_list.index(str(radicand)[:11])) % 2:
                count += 1
            break

print(count)
print("It costs %f s" % (t.perf_counter() - start))

1322
It costs 2.698963 s

TLM

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

import time
st=time.time()
def getGYS(u,v):
    # 获取公约数
    if u<v:
        u,v=v,u
    while v>0:
        u,v=v,u%v
    return u
def get(N,a0=1):
    ans=[]
    # 第一次取整
    while (a0+1)**2<N:
        a0=a0+1
    if (a0+1)**2==N:
        return [a0+1],a0
    ans.append(a0)
    a1=-a0
    a2=N-a0**2
    a0=1
    chushi=(a0,a1,a2)
    while True:
        # 循环取整
        '''
        a3<(a0√N-a1)/a2<a3+1
        转化为(a3*a2+a1)**2<a0**2*N<((a3+1)*a2+a1)**2
        获取a3
        '''
        a3=0
        while ((a3+1)*a2+a1)**2<a0**2*N:
            a3=a3+1
        ans.append(a3)
        # 重新转化为(a0√N-a1)/a2的形式
        a0,a1,a2=a2*a0,-a2*a0*(a1+a2*a3),a0**2*N-(a1+a2*a3)**2
        gys=getGYS(a2,a0)
        a0,a1,a2=a0//gys,a1//gys,a2//gys
        if (a0,a1,a2)==chushi:
            return ans,a0
# N=13
# print(get(N))
n=0
N=10000
a0=1
for i in range(2,N+1):
    ans,a0=get(i,a0)
    if len(ans)%2==0:
        n=n+1
print(n)
print('Time:{}s'.format(time.time()-st))