В этом упражнении точно вычисляются матрицы, обратные к матрицам Гильберта, и обсуждается, как их использовать для проверки точности алгоритма обращения.
22.
π эр Квадрат,
или Арифметические вычисления с высокой точностью
Математику часто считают сухой наукой, однако и математику творили люди. Одной из самых печальных была судьба Уильямса Шенкса, жившего в девятнадцатом веке и посвятившего себя вычислению числа π с высокой точностью. Закончив многолетний труд, Шенкс в 1837 г. опубликовал значение π до 707-го десятичного знака, впоследствии исправив некоторые знаки. Может быть, надо счесть за благо, что Шенкс умер в 1882 г., поскольку в 1946 г. было показано, что его вычисления ошибочны начиная с 528-го десятичного знака. Фактически Шенкс не продвинулся дальше своих предшественников.
Полученное Шенксом значение было проверено, вероятно, с помощью механических устройств, а компьютер был впервые использован для вычисления π только в 1949 г.; это была машина ENIAC. Даже тогда проект был монументальным. Джордж У. Рейтуиснер писал: «Поскольку получить машину в рабочее время было практически невозможно, мы воспользовались разрешением выполнить эту работу за 4 выходных дня в период летних отпусков, когда ENIAC стоял без дела». Собственно вычисления (не программирование!) заняли 70 часов: было получено несколько больше 2 000 цифр. Все это время приходилось постоянно обслуживать компьютер, поскольку из-за ограниченности его возможностей требовались постоянная перфорация и ввод промежуточных результатов. Те первые программисты так же далеки от нынешних, как Шенкс далек от них.
Как бы мы стали вычислять π? Во-первых, необходимо выражение, которое можно вычислять. Ряд
π/4 = 1 − 1/3 + 1/5 − 1/7 + 1/9 − …
довольно прост для понимания, но он ужасно медленно сходится. Гораздо лучше ряд для арктангенса
arctg x = x − х3/3 + х5/5 − х7/7 + …, |x| ≤ 1
Объединим его с формулой сложения для тангенса
tg (a + b) = (tg a + tg b)/(1 − tg a · tg b)
и выберем а и b так, чтобы tg (a + b) = 1 = tg π/4. (Учитывая, что tg (arctg x) = х для −π/2 < х < π/2, можно взять, например, а = arctg (1/2), b = arctg (1/3).)
Тогда
arctg (tg (a + b)) = a + b = arctg 1 = π/4,
и теперь можно использовать приведенный выше ряд для нахождения a и b. На практике чаще всего используются следующие суммы:
π/4 = 4 arctg (1/5) − arctg (1/239),
π/4 = 8 arctg (1/10) − 4 arctg (1/515) − arctg (1/239),
π/4 = 3 arctg (1/4) + arctg (1/20) + arctg (1/1985).
Теперь мы собираемся просуммировать эти ряды на ЭВМ. Как известно, все, что нужно для суммирования, — это простой итерационный цикл, но тут возникает одна проблема. Точность вычислений на ЭВМ ограничена, а весь смысл этого упражнения в том, чтобы найти много-много цифр числа π, значительно превзойдя обычную точность. Первое, что приходит в голову, — промоделировать ручные методы выполнения арифметических действий. Будем представлять числа очень большими целочисленными массивами (по одной десятичной цифре в каждом элементе), тогда ясно, как составить программы сложения, вычитания и умножения. Запрограммировать ручной метод деления несколько сложнее, но все же возможно. Неприемлемым, однако, оказывается время выполнения алгоритмов. Хотя на это редко обращают внимание, но при ручных методах для умножения или деления n-значных чисел требуется время, пропорциональное n². Если речь идет об операциях над числами из тысяч цифр, то такие расходы будут нам не по карману. К счастью, имеются лучшие алгоритмы.
Алгоритм быстрого умножения Тоома—Кука, описываемый Кнутом, зиждется на четырех основных идеях[31]. Вот первая из них. Пусть нам известен способ выполнения некоторой операции над исходными данными размера n за время T(n). Если эту операцию удастся разбить на r частей, выполнение каждой из которых займет менее чем T(n)/r шагов, то такое разбиение позволит улучшить общее время, если, конечно, считать, что вспомогательные организационные расходы не сведут экономию на нет. Пусть, далее, каждая из r частей есть применение того же алгоритма к исходным данным длины n/r и каждая часть может быть разбита аналогичным образом. Тогда можно продолжать это разбиение, пока мы не получим столь короткие исходные данные, что вычисления для них станут тривиальными и займут лишь небольшой фиксированный отрезок времени. Этот принцип разделяй и властвуй обычно дает выигрыш во времени работы алгоритма по крайней мере в log n раз; так, классический метод умножения требует времени n², и его можно свести к 
, что существенно лучше при больших n (не забывайте, что у обеих функций стоимости имеются постоянные множители).
31
Как здесь не отметить, что самую плодотворную идею по части быстрого умножения вы можете позаимствовать у кроликов. Их многочисленное потомство — тому порука.