Идеи реализации

Рассмотрим бесконечную последовательность функций φ 1 λ = x 1 y 1 + λ , φ 2 λ = x 2 y 2 + λ , , φ s λ = x s y s + λ , Эти функции и есть наша главная хитрость. Заметим, что φ 1 0 = x 1 y 1 , φ 1 φ 2 0 = x 1 y 1 + x 2 y 2 , φ 1 φ 2 φ 3 0 = x 1 y 1 + x 2 y 2 + x 3 y 3 , , φ 1 φ 2 φ 3 φ n 0 = x 1 y 1 + x 2 y 2 + x 3 y 3 + + x n y n .

Пусть имеются две функции f λ и g λ . Новую функцию f g λ можно представлять себе как результат некоторой операции над функциями f и g — операции композиции. Операцию принято обозначать кружком: f g λ = f g λ . Операция композиции является ассоциативной, то есть для неё выполняется сочетательный закон: равенство f g h = f g h выполняется для любых трёх функций f, g, h. Этот факт вытекает прямо из определения композиции.

Так что наша задача состоит в последовательном нахождении композиций φ 1 φ 2 φ n . Подстановка ноля вместо λ в полученные функции даёт подходящие дроби для заданной непрерывной дроби.

Пора обратить внимание на то, что все функции φ s принадлежат замечательному семейству дробно-линейных функций, то есть является отношением двух линейных функций. Общий вид дробно-линейной функции таков: f x = a x + b c x + d (числа c и d не должны обращаться в ноль одновременно). Дробно-линейная функция полностью задаётся четырьмя числами — своими коэффициентами.

Семейство замечательно, помимо прочего, ещё и тем, что композиция дробно-линейных функций снова является дробно-линейной. Математики говорят, что множество дробно-линейных функций замкнуто относительно операции композиции. В справедливости этого утверждения легко убедиться прямым вычислением: a λ + b c λ + d α λ + β γ λ + δ = a α λ + β γ λ + δ + b c α λ + β γ λ + δ + d = a α λ + β + b γ λ + δ c α λ + β + d γ λ + δ = a α + b γ λ + a β + b δ c α + d γ λ + c β + d δ .

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

Удобно (а главное полезно) поставить в соответствие дробно-линейной функции квадратную матрицу с четырьмя коэффициентами: a λ + b c λ + d a b c d . (о матрицах мы уже кое-что писали в разделе «Матричная версия» главы 9. «Числа Фибоначчи»). Это соответствие не является взаимно-однозначным. В то время как матрица вполне задаёт дробно-линейную функцию (кроме матриц с нулевой нижней строкой), каждая такая функция может быть задана многими разными матрицами, отличающимся друг от друга умножением всех коэффициентов на одно и то же ненулевое число.

Особенно нас должно порадовать то, что композиции двух дробно-линейных функций соответствует произведение их матриц: a λ + b c λ + d α λ + β γ λ + δ a α + b γ a β + b δ c α + d γ c β + d δ = a b c d α β γ δ .

В нашем случае φ s 0 x s 1 y s , следовательно, φ 1 φ 2 φ n 0 x 1 1 y 1 0 x 2 1 y 2 0 x n 1 y n = a n b n c n d n . Вычислив элементы a n , b n , c n , d n , найдём n-ю подходящую дробь. Для этого в функцию a n λ + b n c n λ + d n нужно подставить λ = 0 , тогда получим число b n d n .

Таким образом, задача вычисления n-й подходящей дроби свелась к вычислению произведения матриц 0 x 1 1 y 1 0 x 2 1 y 2 0 x n 1 y n = s = 1 n 0 x s 1 y s . Это произведение можно вычислять индуктивно, точно так же, как вычислялось произведение элементов числовой последовательности. Рассмотрим функцию на пространстве последовательностей числовых пар P x 1 y 1 x n y n = s = 1 n 0 x s 1 y s (она принимает матричные значения). При удлинении последовательности функция перевычисляется по формуле P ω x y = P ω 0 x 1 y . В качестве значения функции P для пустой последовательности следует взять единичную матрицу: P = I . Для элементов матрицы P ω = a ω b ω c ω d ω получаются такие формулы перевычисления: a ω x y = b ω , b ω x y = a ω x + b ω y , c ω x y = d ω , d ω x y = c ω x + d ω y . Итак, четыре функции a, b, c, d осуществляют индуктивное расширение функции f, причём f = b d .

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

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

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

#!/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 — «не число». Почему это «не число» ещё и с минусом, мы уже не в состоянии объяснить.

Информатика-54© А. Н. Швец