欧拉计划 Project Euler64（奇周期平方根）题解

题干

奇数平方根

所有的平方根写成如下连分数表示时都是周期性重复的：
$$\sqrt{N}=a_0+\frac{1}{a_1+\frac{1}{a_2+\frac{1}{a_3+\dots }}}$$
例如，让我们考虑$\sqrt{23}$
$\sqrt{23}=4+\sqrt{23}-4=4+\frac{1}{\frac{1}{\sqrt{23}-4}}=4+\frac{1}{1+\frac{\sqrt{23}-3}{7}}$
如果我们继续展开，会得到：
$$\sqrt{23}=4+\frac{1}{1+\frac{1}{3+\frac{1}{1+\frac{1}{8+\dots }}}}$$
这个过程可以总结如下：
$a_0=4,\frac{1}{\sqrt{23}-4}=\frac{\sqrt{23}+4}{7}=1+\frac{\sqrt{23}-3}{7}$

$a_1=1,\frac{7}{\sqrt{23}-3}=\frac{7(\sqrt{23}+3)}{14}=3+\frac{\sqrt{23}-3}{2}$

$a_2=3,\frac{2}{\sqrt{23}-3}=\frac{2(\sqrt{23}+3)}{14}=1+\frac{\sqrt{23}-4}{7}$

$a_3=1,\frac{7}{\sqrt{23}-4}=\frac{7(\sqrt{23}+4)}{7}=8+\sqrt{23}-4$

$a_4=8,\frac{1}{\sqrt{23}-4}=\frac{\sqrt{23}+4}{7}=1+\frac{\sqrt{23}-3}{7}$

$a_5=1,\frac{7}{\sqrt{23}-3}=\frac{7(\sqrt{23}+3)}{14}=3+\frac{\sqrt{23}-3}{2}$

$a_6=3,\frac{2}{\sqrt{23}-3}=\frac{2(\sqrt{23}+3)}{14}=1+\frac{\sqrt{23}-4}{7}$

$a_7=1,\frac{7}{\sqrt{23}-4}=\frac{7(\sqrt{23}+4)}{7}=8+\sqrt{23}-4$

可以看出序列正在重复。我们将其简记为 $\sqrt{23}=[4;(1,3,1,8)]$，表示$(1,3,1,8)$在此之后无限循环。
前$10$个（无理数）平方根的连分数表示是：
$\sqrt{2}=[1;(2)]$,周期为$1$

$\sqrt{3}=[1;(1,2)]$, 周期为$2$

$\sqrt{5}=[2;(4)]$, 周期为$1$

$\sqrt{6}=[2;(2,4)]$, 周期为$2$

$\sqrt{7}=[2;(1,1,1,4)]$, 周期为$4$

$\sqrt{8}=[2;(1,4)]$, 周期为$2$

$\sqrt{10}=[3;(6)]$,周期为$1$

$\sqrt{11}=[3;(3,6)]$, 周期为$2$

$\sqrt{12}=[3;(2,6)]$, 周期为$2$

$\sqrt{13}=[3;(1,1,1,1,6)]$, 周期为$5$
在$N\le 13$中，恰好有$4$个连分数表示的周期是奇数。
在$N\le 10000$中，有多少个连分数表示的周期是奇数？

思路

我们需要找出$N\le 10000$的范围内，有多少个整数$N$使得其平方根$\sqrt{N}$的连分数表示的周期长度是奇数。我们只考虑$N$不是完全平方数的情况，因为完全平方数的平方根是整数（有理数），其连分数表示是有限的，没有周期。

连分数展开算法

求$\sqrt{N}$的连分数表示$[a_0:(a_1,a_2,\dots ,a_p)]$的标准算法如下:

1. 令$x_0=\sqrt{N}$
2. 计算$a_0=\lfloor x_0\rfloor =\lfloor \sqrt{N}\rfloor$
3. 我们寻找一种形式${\xi }_k=\frac{\sqrt{N}+m_k}{d_k}$，其中$m_k,d_k$是整数
4. 初始状态: $m_0=0,d_0=1$
5. 迭代计算$a_k,m_{k+1},d_{k+1}(k\ge 0)$:
   • $a_k=\lfloor \frac{\sqrt{N}+m_k}{d_k}\rfloor =\lfloor \frac{a_0+m_k}{d_k}\rfloor$
   • $m_{k+1}=d_k{a}_k-m_k$
   • $d_{k+1}=\frac{N-m_{k+1}^2}{d_k}$
6. 这个序列$a_1,a_2,\dots$是周期性的，周期从$a_1$开始
7. 当$a_k=2a_0$时,达到周期的最后一个元素$a_p$，此时的$k$就是周期长度$p$

code

// 1322
#include <bits/stdc++.h>using namespace std;using ll = long long;void solve() {int li = 10000;int ans = 0;for (int N = 2; N <= li; ++N) {int a0 = static_cast<int>(sqrt(N));if (a0 * a0 == N) continue; // 完全平方数直接跳过即可int p = 0; // 周期长度// 第一轮参数  对应 m1 d1 a1int m = a0;int d = N - m * m;int a = static_cast<int>(floor((double)(a0 + m) / d));while (a != 2 * a0) {m = d * a - m;d = (N - m * m) / d;a = static_cast<int>(floor((double)(a0 + m) / d));p++;}p++;if (p % 2) {ans++;}}cout << ans << "\n";
}int main() {ios::sync_with_stdio(false);cin.tie(nullptr);int tt = 1; // cin >> tt;while (tt--) {solve();}return 0;
}