连分数-ProjectEuler 946

wayne

今天Project Euler发布了一道题,连分数.翻译过来就是:
已知 连分数的定义为\[a_0+ \cfrac 1{a_1+ \cfrac 1{a_2+\cfrac 1{a_3+\ddots }}}= [a_0;a_1,a_2,a_3,\ldots]\]
假设$ \alpha = [2;1,1,2,1,1,1,2,1,1,1,1,1,2,1,1,1,1,1,1,1,2,1,1,1,1,1,1,1,1,1,1,1,2,\ldots]$, 其中2中间的1的个数是连续的素数.
然后我们定义一个数$\beta = \frac{2\alpha+3}{3\alpha+2}$ , 其连分数展开的前$10$项是$[0;1,5,6,16,9,1,10,16,11]$, 和为$75$, 求$\beta$的连分数展开的前$10^8$项的和

mathe

10^8应该还好，计算若干项后，留下的还是一个alpha的线性分式，所以关键在于判断它的整数部分是多少，相信大部分情况使用alpha的近似值即可判断

mathe

所以可以维护一个alpha两个相邻的近似分数表示，后面每个分式中alpha用这两个分数替换如果整数部分相同即确定其整数部分，不然，继续寻找某个更高精度表示，确定整数部分，而后续项退回用第一个继续判断

mathe

比如可以提前算出连分数前10/11项，20/21项，30/31项...等分式

mathe

发现另外一种方法可能更好，我们知道alpha=2+1/alpha1,alpha1=1+1/alpha2
等，每个表达式可以让我们对左边的字母给出了自然的范围估计，也给出了将左边字母替换为右边的方案。
比如第一个说明2<alpha<3
利用这个就可以判断出
beta=(2alpha+3)/(3alpha+2)<1。由此得出第一项为0，表达式变为(3alpha+2)/(2alpha+3),并且范围在8/7和11/9之间，所以下一项为1，余下表达式为(2alpha+3)/(alpha-1),所以范围在7和9/2之间，无法确定范围，将第一个表达式代入，得到
(7+2/alpha1)/(1+1/alpha1)=(2+7alpha1)/(1+alpha1)
根据第二式我们知道1<alpha1<2,所以上式在9/2和16/3之间，还是不能确定整数部分，继续代入第二式得到
(9+7/alpha2)/(2+1/alpha2)=(7+9alpha2)/(1+2alpha2)
然后根据1<alpha2<2判断出上面整数部分为5，输出5并且将分式改为(1+2alpha2)/(2-alpha2).
范围3到无穷，所以将alpha2=1+1/alpha3再代入，得到(3+2/alpha3)/(1-1/alpha3)=(2+3alpha3)/(-1+alpha3)
下一步2<alpha3<3,得出上式8和11/2之间，继续代入变为
(8+3/alpha4)/(1+1/alpha4)=(3+8alpha4)/(1+alpha4)
范围11/2和19/3之间还得继续替换

wayne

设$\alpha $的第n项的渐进分数是$\frac{A_n}{B_n}$,那么存在递推公式\[ (A_n,B_n) = (1,a_n). \left(
\begin{array}{cc}
A_{n-2} & B_{n-2} \\
A_{n-1} & B_{n-1} \\
\end{array}
\right) \]
逆推表达式应该是\[ (1,a_n) =  (A_n,B_n). (
\begin{array}{cc}
A_{n-2} & B_{n-2} \\
A_{n-1} & B_{n-1} \\
\end{array}
) ^{-1}\]

而$\beta$的渐进展开应该是 $ \frac{A'_n}{B'_n}= \frac{2\alpha+3}{3\alpha+2} = \frac{2A_n+3B_n}{3A_n+2B_n}$, 于是,需要根据已知的渐进展开,反推出$\beta$连分数对应项的$a'_n$

mathe

最终转化为二阶阵的计算是对的，关键是[0,1;1,1]^p [1,0;2,1]的计算，不如速度会太慢

wayne

根据wikipedia,  $A_{n-1}B_n-A_{n}B_{n-1}=(-1)^na_1a_2...a_n = \prod_{i=1}^n (-a_i)$

mathe

我们如果把无穷连分数到某一项后记成t,那么再加上前面若干项，可以写成
[1,1,1,...,1,2+t]
这个数最后一项可以写成\(\frac{1+0t}{2+1t}\), 我们可以用矩阵\(\begin{bmatrix}1&0\\2&1\end{bmatrix}\)表示。 假设某个连分数为\(\frac{a+bt}{c+dt}\)，前面再加上一个1，值就变成\(\frac1{1+\frac{a+bt}{c+dt}}=\frac{c+dt}{(a+c)+(b+d)t}\), 相等于矩阵\(\begin{bmatrix}a&b\\c&d\end{bmatrix}\)变换为
\(\begin{bmatrix}c&d\\a+c&b+d\end{bmatrix}\), 相当于左乘了矩阵\(\begin{bmatrix}0&1\\1&1\end{bmatrix}\).
所以上面[1,1,...,1,2+t]如果有p个1,对应变成了\(M_p=\begin{bmatrix}0&1\\1&1\end{bmatrix}^p \begin{bmatrix}1&0\\2&1\end{bmatrix}\)

由此，对于本题，初始\(\alpha\)值我们可以看成一个\(2+t\), 对应\(\beta=\frac{2\alpha+3}{3\alpha+2}=\frac{7+2t}{8+3t}\), 其中\(0\lt t \lt 1\)
于是可以判断处整数部分为0，输出0， 然后变换为\(\frac{8+3t}{7+2t}\),可以判断出整数部分为1，输出1，然后变化为\(\frac{7+2t}{1+t}\).
于是我们可以判断出整数部分在7和9/2之间，但是无法确定，于是需要做第一次替换，将t替换为\(M_2=\begin{bmatrix}3&1\\5&2\end{bmatrix}\)即\(\frac{3+t}{5+2t}\)
替换结果为\(\frac{41+16t}{8+3t}\) （相当于矩阵乘法\(\begin{bmatrix}2&7\\1&1\end{bmatrix}\begin{bmatrix}3&1\\5&2\end{bmatrix}\)）同样\(0\lt t\lt 1\), 所以整数部分是5，得到下一个分数为\(\frac{8+3t}{1+t}\),所以整数部分在8和11/2之间，又无法确定, 所以我们这次需要将t替换为\(M_3\)， 然后可以确定整数部分为6.

wayne

$\beta$的连分数, 前3000项, 感觉随机性很差.

1. {0,1,5,6,16,9,1,10,16,11,11,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,11,11,11,9,1,10,11,11,11,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,11,11,11,11,11,10,1,9,11,11,11,11,11,9,1,10,11,11,11,11,11,11,16,11,11,11,11,11,11,11,11,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,11,11,11,11,11,11,11,11,11,9,1,10,11,11,11,11,11,11,11,11,11,9,1,10,11,11,11,11,11,11,11,11,11,11,11,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,11,11,11,11,11,11,11,11,11,11,11,11,11,9,1,10,11,11,11,11,11,11,11,11,11,11,11,11,11,6,11,11,11,11,11,11,11,11,11,11,11,11,11,11,16,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,16,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,6,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,9,1,10,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,16,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,9,1,10,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,9,1,10,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,6,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,6,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,9,1,10,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,9,1,10,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,9,1,10,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,16,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,6,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,9,1,10,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,6,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,10,1,9,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,10,1,9,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,16,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,6,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,6,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,9,1,10,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,9,1,10,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,16,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,10,1,9,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,10,1,9,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,6,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,10,1,9,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,11,9}

复制代码