题目64：10000以下的连分数中有多少拥有奇数周期？
本帖最后由 永恒的蓝色梦想 于 2020-6-5 13:30 编辑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
√3=, 周期=2
√5=, 周期=1
√6=, 周期=2
√7=, 周期=4
√8=, 周期=2
√10=, 周期=1
√11=, 周期=2
√12= , 周期=2
√13=, 周期=5
可以看出对于 N ≤ 13, 有四个连分数的周期是奇数。
对于 N ≤ 10000，有多少连分数的周期是奇数？
刚刚开始学习，看着好厉害的样子:lol: mark下 本帖最后由 jerryxjr1220 于 2016-10-14 21:22 编辑
这题的思路就是找出除整数部分后面的无理数从第几位开始循环。
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 = 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 于 2016-10-14 21:22 编辑
结果1322 本帖最后由 jerryxjr1220 于 2016-10-16 22:17 编辑
另外一个暴力的的方法，用decimal库，把精度提高到250位以上（200位小数，仍有3个数有误差）
同样可以得到结果......1322......7秒多，速度也还行，简单粗暴，不用考虑复杂的数学方程式{:5_109:}
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) 用的matlab
浮点数真是头疼，
1322
时间已过 0.413237 秒。
>> 本帖最后由 永恒的蓝色梦想 于 2021-2-20 13:18 编辑
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;
} 和求无限小数的循环节问题很像
1322
#include<iostream>
#include<cmath>
#include<algorithm>
using namespace std;
int sq = {0};//平方数表
void init(){ //打平方数表
for (int i = 2; i <= 101; i++)sq = 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;
} /*
答案：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;
}
import time as t
import decimal as dc
start = t.perf_counter()
count = 0
dc.getcontext().prec = 250
num_list =
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 =
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 于 2023-2-27 20:20 编辑
不使用浮点数开方进行求解
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
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))
页:
[1]