Рассмотрим бесконечную последовательность функций $$\begin{array}{l}{\phi }_1\left(\lambda \right)=\frac{x_1}{y_1+\lambda }\text{,}\\ {\phi }_2\left(\lambda \right)=\frac{x_2}{y_2+\lambda }\text{,}\\ \dots \text{,}\\ {\phi }_s\left(\lambda \right)=\frac{x_s}{y_s+\lambda }\text{,}\\ \dots \end{array}$$ Эти функции и есть наша главная хитрость. Заметим, что $$\begin{array}{l}{\phi }_1\left(0\right)=\frac{\overline{)x_1}}{\overline{)y_1}}\text{,}\\ {\phi }_1\left({\phi }_2\left(0\right)\right)=\frac{\overline{)x_1}}{\overline{)y_1}}+\frac{\overline{)x_2}}{\overline{)y_2}}\text{,}\\ {\phi }_1\left({\phi }_2\left({\phi }_3\left(0\right)\right)\right)=\frac{\overline{)x_1}}{\overline{)y_1}}+\frac{\overline{)x_2}}{\overline{)y_2}}+\frac{\overline{)x_3}}{\overline{)y_3}}\text{,}\\ \dots \text{,}\\ {\phi }_1\left({\phi }_2\left({\phi }_3\left(\dots {\phi }_n\left(0\right)\dots \right)\right)\right)=\frac{\overline{)x_1}}{\overline{)y_1}}+\frac{\overline{)x_2}}{\overline{)y_2}}+\frac{\overline{)x_3}}{\overline{)y_3}}+\dots +\frac{\overline{)x_n}}{\overline{)y_n}}\text{.}\end{array}$$

Пусть имеются две функции $f\left(\lambda \right)$ и $g\left(\lambda \right)$. Новую функцию $f\left(g\left(\lambda \right)\right)$ можно представлять себе как результат некоторой операции над функциями $f$ и $g$ — операции композиции. Операцию принято обозначать кружком: $$f\left(g\left(\lambda \right)\right)=\left(f\circ g\right)\left(\lambda \right)\text{.}$$ Операция композиции является ассоциативной, то есть для неё выполняется сочетательный закон: равенство $$\left(f\circ g\right)\circ h=f\circ \left(g\circ h\right)$$ выполняется для любых трёх функций $f$, $g$, $h$. Этот факт вытекает прямо из определения композиции.

Так что наша задача состоит в последовательном нахождении композиций ${\phi }_1\circ {\phi }_2\circ \dots \circ {\phi }_n$. Подстановка ноля вместо $\lambda$ в полученные функции даёт подходящие дроби для заданной непрерывной дроби.

Пора обратить внимание на то, что все функции ${\phi }_s$ принадлежат замечательному семейству дробно-линейных функций, то есть является отношением двух линейных функций. Общий вид дробно-линейной функции таков: $$f\left(x\right)=\frac{ax+b}{cx+d}$$ (числа $c$ и $d$ не должны обращаться в ноль одновременно). Дробно-линейная функция полностью задаётся четырьмя числами — своими коэффициентами.

Семейство замечательно, помимо прочего, ещё и тем, что композиция дробно-линейных функций снова является дробно-линейной. Математики говорят, что множество дробно-линейных функций замкнуто относительно операции композиции. В справедливости этого утверждения легко убедиться прямым вычислением: $$\frac{a\lambda +b}{c\lambda +d}\circ \frac{\alpha \lambda +\beta }{\gamma \lambda +\delta }=\frac{a\cdot \frac{\alpha \lambda +\beta }{\gamma \lambda +\delta }+b}{c\cdot \frac{\alpha \lambda +\beta }{\gamma \lambda +\delta }+d}=\frac{a\left(\alpha \lambda +\beta \right)+b\left(\gamma \lambda +\delta \right)}{c\left(\alpha \lambda +\beta \right)+d\left(\gamma \lambda +\delta \right)}=\frac{\left(a\alpha +b\gamma \right)\lambda +\left(a\beta +b\delta \right)}{\left(c\alpha +d\gamma \right)\lambda +\left(c\beta +d\delta \right)}\text{.}$$

Это наблюдение позволяет, зная наборы коэффициентов двух функций, легко найти коэффициенты их композиции.

Удобно (а главное полезно) поставить в соответствие дробно-линейной функции квадратную матрицу с четырьмя коэффициентами: $$\frac{a\lambda +b}{c\lambda +d}\rightleftarrows \left(\begin{array}{cc}a& b\\ c& d\end{array}\right)\text{.}$$ (о матрицах мы уже кое-что писали в разделе «Матричная версия» главы 9. «Числа Фибоначчи»). Это соответствие не является взаимно-однозначным. В то время как матрица вполне задаёт дробно-линейную функцию (кроме матриц с нулевой нижней строкой), каждая такая функция может быть задана многими разными матрицами, отличающимся друг от друга умножением всех коэффициентов на одно и то же ненулевое число.

Особенно нас должно порадовать то, что композиции двух дробно-линейных функций соответствует произведение их матриц: $$\frac{a\lambda +b}{c\lambda +d}\circ \frac{\alpha \lambda +\beta }{\gamma \lambda +\delta }\rightleftarrows \left(\begin{array}{cc}a\alpha +b\gamma & a\beta +b\delta \\ c\alpha +d\gamma & c\beta +d\delta \end{array}\right)=\left(\begin{array}{cc}a& b\\ c& d\end{array}\right)\left(\begin{array}{cc}\alpha & \beta \\ \gamma & \delta \end{array}\right)\text{.}$$

В нашем случае $${\phi }_s\rightleftarrows \left(\begin{array}{cc}0& x_s\\ 1& y_s\end{array}\right)\text{,}$$ следовательно, $${\phi }_1\circ {\phi }_2\circ \dots \circ {\phi }_n\rightleftarrows \left(\begin{array}{cc}0& x_1\\ 1& y_1\end{array}\right)\left(\begin{array}{cc}0& x_2\\ 1& y_2\end{array}\right)\dots \left(\begin{array}{cc}0& x_n\\ 1& y_n\end{array}\right)=\left(\begin{array}{cc}{a}_n& b_n\\ c_n& d_n\end{array}\right)\text{.}$$ Вычислив элементы $a_n$, $b_n$, $c_n$, $d_n$, найдём $n$-ю подходящую дробь. Для этого в функцию $\frac{a_n\lambda +b_n}{c_n\lambda +d_n}$ нужно подставить $\lambda =0$, тогда получим число $\frac{b_n}{d_n}$.

Таким образом, задача вычисления $n$-й подходящей дроби свелась к вычислению произведения матриц $$\left(\begin{array}{cc}0& x_1\\ 1& y_1\end{array}\right)\left(\begin{array}{cc}0& x_2\\ 1& y_2\end{array}\right)\dots \left(\begin{array}{cc}0& x_n\\ 1& y_n\end{array}\right)=\prod _{s=1}^n\left(\begin{array}{cc}0& x_s\\ 1& y_s\end{array}\right)\text{.}$$ Это произведение можно вычислять индуктивно, точно так же, как вычислялось произведение элементов числовой последовательности. Рассмотрим функцию на пространстве последовательностей числовых пар $$P\left(〈\left(x_1,y_1\right),\dots ,\left(x_n,y_n\right)〉\right)=\prod _{s=1}^n\left(\begin{array}{cc}0& x_s\\ 1& y_s\end{array}\right)$$ (она принимает матричные значения). При удлинении последовательности функция перевычисляется по формуле $$P\left(\omega ∟\left(x,y\right)\right)=P\left(\omega \right)\cdot \left(\begin{array}{cc}0& x\\ 1& y\end{array}\right)\text{.}$$ В качестве значения функции $P$ для пустой последовательности следует взять единичную матрицу: $P\left(〈〉\right)=I$. Для элементов матрицы $$P\left(\omega \right)=\left(\begin{array}{cc}a\left(\omega \right)& b\left(\omega \right)\\ c\left(\omega \right)& d\left(\omega \right)\end{array}\right)$$ получаются такие формулы перевычисления: $$\begin{array}{l}a\left(\omega ∟\left(x,y\right)\right)=b\left(\omega \right)\text{,}\\ b\left(\omega ∟\left(x,y\right)\right)=a\left(\omega \right)\cdot x+b\left(\omega \right)\cdot y\text{,}\\ c\left(\omega ∟\left(x,y\right)\right)=d\left(\omega \right)\text{,}\\ d\left(\omega ∟\left(x,y\right)\right)=c\left(\omega \right)\cdot x+d\left(\omega \right)\cdot y\text{.}\end{array}$$ Итак, четыре функции $a$, $b$, $c$, $d$ осуществляют индуктивное расширение функции $f$, причём $f=\frac{b}{d}$.

Приведём программу, вычисляющую последовательные приближения к числу $\pi$ по формуле Броункера. Её алгоритм в точности следует общей схеме вычисления индуктивной функции, но с одним добавлением.

Дело в том, что элементы матрицы $P\left(\omega \right)$, а значит, числители и знаменатели подходящих дробей, имеют обыкновение очень быстро расти при удлинении последовательности. Сами же дроби при этом могут оставаться ограниченными (именно так и будет, если непрерывная дробь сходится). Очень скоро числители и знаменатели превысят допустимые пределы для числовых значений языка программирования, и произойдёт арифметическое переполнение. К счастью, с этим явлением можно бороться.

Вспомним, что матрицу дробно-линейной функции можно поэлементно умножить на любое ненулевое число, и получившаяся матрица будет по-прежнему изображать ту же самую функцию. Можно, вычислив $P\left(\omega \right)$, тут же разделить все её элементы, скажем, на $c\left(\omega \right)$, тем самым превратив $c$ в единицу, и заметно уменьшив все остальные элементы. Этому делению посвящена вторая строка в теле цикла:

examples/perl/pi-brouncker.plPerl
#!/usr/bin/perl

use warnings;

my ($a, $b, $c, $d)=(1, 0, 0, 1);
my $n=1;

while()
{
	($a, $b, $c, $d)=($b, $a*$n**2+$b*2, $d, $c*$n**2+$d*2);
	($a, $b, $c, $d)=map { $_/$c } ($a, $b, $c, $d);
	$n+=2;
	print 4/($b/$d+1), "\n";
}

Попробуйте убрать эту строчку, и вы увидите, что получится:

% ./pi-brouncker.pl
2.66666666666667
3.46666666666667
2.8952380952381
3.33968253968254
2.97604617604618
3.28373848373848
3.01707181707182
(пропущено чуть меньше сотни строк вывода)
3.14830398741449
3.13492606099309
-nan
-nan
-nan
-nan
-nan
-nan
…

Программа, не добравшись даже до сотого приближения, начинает выводить странные значения -nan. Аббревиатура NaN означает Not a Number — «не число». Почему это «не число» ещё и с минусом, мы уже не в состоянии объяснить.