Project Euler 64: Odd period square roots
所有平方根写成连分数时都是周期性的，连分数的形式如下：
\[
\displaystyle \quad \quad \sqrt{N}=a_0+\frac 1 {a_1+\frac 1 {a_2+ \frac 1 {a3+ \dots}}}
\]
例如，让我们看一下\(\sqrt{23}\)：
\[
\quad \quad \sqrt{23}=4+\sqrt{23}-4=4+\frac 1 {\frac 1 {\sqrt{23}-4}}=4+\frac 1 {1+\frac{\sqrt{23}-3}7}
\]
如果我们继续下去，可以得到以下展开式：
\[
\displaystyle \quad \quad \sqrt{23}=4+\frac 1 {1+\frac 1 {3+ \frac 1 {1+\frac 1 {8+ \dots}}}}
\]
这个展开的过程可以展示如下：
\[
\begin{aligned}\quad \quad a_0&=4, \frac 1 {\sqrt{23}-4}=\frac {\sqrt{23}+4} 7=1+\frac {\sqrt{23}-3} 7\\\quad \quad a_1&=1, \frac 7 {\sqrt{23}-3}=\frac {7(\sqrt{23}+3)} {14}=3+\frac {\sqrt{23}-3} 2\\\quad \quad a_2&=3, \frac 2 {\sqrt{23}-3}=\frac {2(\sqrt{23}+3)} {14}=1+\frac {\sqrt{23}-4} 7\\\quad \quad a_3&=1, \frac 7 {\sqrt{23}-4}=\frac {7(\sqrt{23}+4)} 7=8+\sqrt{23}-4\\\quad \quad a_4&=8, \frac 1 {\sqrt{23}-4}=\frac {\sqrt{23}+4} 7=1+\frac {\sqrt{23}-3} 7\\\quad \quad a_5&=1, \frac 7 {\sqrt{23}-3}=\frac {7 (\sqrt{23}+3)} {14}=3+\frac {\sqrt{23}-3} 2\\\quad \quad a_6&=3, \frac 2 {\sqrt{23}-3}=\frac {2(\sqrt{23}+3)} {14}=1+\frac {\sqrt{23}-4} 7\\\quad \quad a_7&=1, \frac 7 {\sqrt{23}-4}=\frac {7(\sqrt{23}+4)} {7}=8+\sqrt{23}-4\\\end{aligned}
\]
可以看出来这个序列出现了循环，具体来说，我们用记号\(\sqrt{23}=[4;(1,3,1,8)]\)来表示(1, 3, 1, 8)这一部分是无限循环的。前十个无理数平方根的连分数表示如下：
\[
\begin{aligned}\quad \quad \sqrt{2}&=[1;(2)],period=1\\\quad \quad \sqrt{3}&=[1;(1,2)],period=2\\\quad \quad \sqrt{5}&=[2;(4)],period=1\\\quad \quad \sqrt{6}&=[2;(2,4)],period=2\\\quad \quad \sqrt{7}&=[2;(1,1,1,4)],period=4\\\quad \quad \sqrt{8}&=[2;(1,4)],period=2\\\quad \quad \sqrt{10}&=[3;(6)],period=1\\\quad \quad \sqrt{11}&=[3;(3,6)],period=2\\\quad \quad \sqrt{12}&=[3;(2,6)],period=2\\\quad \quad \sqrt{13}&=[3;(1,1,1,1,6)],period=5\\\end{aligned}
\]
即对于\(N\le13\)有四个连分数表示的奇数的周期，那么对于\(N\le10000\)有多少个连分数表示有奇数的周期？
分析：这道题并不像它的题目看上去那么复杂，核心在于寻找\(a_0,a_1,a_2\)等项之间的递推关系，找到这个递推关系后我们就可以很方便的计算连分数表示的周期。经搜索相关文献，我们发现在这个维基词条已经非常详细的描述了各个\(a\)项之间的递推关系，简述如下：对于某个非完全平方数\(S\)，其平方根为无理数且有周期，设\(a_0=\lfloor{\sqrt{S}}\rfloor,m_0=0,d_0=1\)，则三项的递推关系如下：
\[
m_{n+1}=d_na_n-m_n,\ d_{n+1}=\frac{S-m_{n+1}^2}{d_n},\ a_{n+1}=\bigg\lfloor{\frac{a_0+m_{n+1}}{d_{n+1}}}\bigg\rfloor
\]
此外根据文献推论3.3，当\(a_i=2a_0\)时循环周期即结束。如对于\(\sqrt{13}\)，\(a_5=6=2a_0\)，此时循环周期终止，其循环周期长度为五。因此，当我们计算连分数表示时，当计算到某一项是首项的二倍时，就不需要再往下计算了，此时可以直接计算出该连分数表示的循环周期长度。
使用以上原理，我们计算小于等于一万的非完全平方数的平方根的连分数表示的循环周期长度，对于完全平方数我们则可以直接路过，加快筛选的速度。最后我们统计计算出来的循环周期长度中的奇数值，即为题目所求，代码如下：
# time cost = 285 ms ± 3.41 ms
def square_root_period(n):
a0 = int(n**0.5)
if a0**2 == n:
return 0
else:
a,m,d,arr = a0,0,1,[a0]
while arr[-1] != 2*a0:
m = d*a - m
d = (n - m**2)//d
a = int((a0+m)/d)
arr.append(a)
return len(arr)-1
def main(N=10000):
c = 0
for i in range(1,N+1):
if square_root_period(i) % 2 == 1:
c += 1
return c