Этюды для программистов

3. (Вычисление последовательных приближений к 1/v.) Выполнить шаг 4 для i = 1, 2, …, j − 1.

4. (Вычисление приближения из 2ⁱ битов.) Присвоить [2ⁱ⁺¹, d] значение

2^3·2i [2ⁱ⁻¹ + 1, a] − [2ⁱ⁻¹ + 1, a]² ([k, v]/2^k−2i)

Деление в скобках (фактически сдвиг вправо) должно выполняться до умножения; идея состоит в том, чтобы ускорить умножение, отбросив лишние биты v, ненужные в данном приближения. Хотя кажется, что результат d может содержать больше 2ⁱ⁺¹ битов, этого никогда не произойдет. Затем присвоить [2ⁱ + 1, а] значение [2ⁱ⁺¹, d]/2ⁱ⁻¹.

5. (Улучшение окончательной оценки.) Присвоить [3k, d] значение

2^2k[k + 1, а] − [k + 1, а]²·[k, v].

Затем присвоить [k + 1, а] значение

([Зk, d] + 2^2k−2)/2^2k−1.

6. (Окончательное деление.) Выдать в качестве результата

([k + 1, a] · [m, u] + 2^k+n−2)/2^2+k−1[33].

Как использовать алгоритмы?

Для нахождения π надо будет провести вычисления по одной из формул, выписанных ранее в этом пункте, с использованием ряда для арктангенса. Фактически для страховки следует использовать две формулы и затем сравнить результаты бит за битом. Значением π будет общая начальная часть этих двух результатов.

Все еще остается открытым вопрос, как с помощью алгоритмов, работающих только с целыми числами, получить очевидно дробные значения членов ряда. Пусть мы хотим вычислить я, скажем, с точностью 1000 битов. Вычислим тогда 2¹⁰⁰⁰π, умножив все числители на 2¹⁰⁰⁰. Эта процедура делает также все делимые много больше делителей (как предполагалось выше) и позволяет прекратить вычисления, когда частные станут нулевыми.

Выберем теперь (не обязательно наилучший) ряд для вычислений, скажем

π = 16 arctg(1/5) − 4 arctg(1/239).

Мы фактически будем вычислять 2¹⁰⁰⁰π, поэтому хотелось бы вычислить 2¹⁰⁰⁰·16 arctg (1/5). Первым членом соответствующего ряда будет 2¹⁰⁰⁰·16/5; назовем его a₁ (отметим, что a₁ складывается с суммой). Теперь, чтобы получить следующий член a_i+1 из a_i, поделим a₁ на 5·5·(2i − 1)[34]. Если a_i добавлялся к сумме, то вычтем a_i+1 из суммы, если a_i вычитался, прибавим a_i+i. Будем поделим a₁ на 5·5·(2i − 1). Если a_i добавлялся к сумме, то вычисления заканчиваются, когда члены обоих рядов станут нулевыми. В результате получим примерно тысячу битов числа π. Результат, конечно, надо будет перевести в десятичную систему.

Тема. Составьте программы, реализующие описанные выше алгоритмы умножения и деления, и все необходимые им служебные подпрограммы. Используйте их для вычисления я с высокой точностью при помощи одного из выписанных рядов. Проследите, чтобы ваши программы не оказались слишком тесно привязаны к вычислению π; библиотека программ для вычислений с высокой точностью может быть полезна и для других задач. Должна быть предусмотрена возможность увеличения точности счета без изменения программ, а лишь путем расширения памяти для результатов. Выходные данные должны включать статистику по использованию каждой программы, по числу выполнений каждого шага двух центральных алгоритмов и по использованию памяти. Сбор такой информации обойдется очень дешево в сравнении со всей работой.

Указания исполнителю. Этот этюд длинный и трудный. Не последнюю роль здесь играет то, что два центральных алгоритма нужно в какой-то степени принимать на веру. Однако, как это часто бывает в реальных задачах, главной проблемой является не кодирование программы, а выбор структур данных. Как представлять длинные числа? Обозначение [m, u] наводит на мысль, что всякое длинное число представляется парой аргументов длина и значение. Часть длина легко реализуется, но значение имеет, очевидно, переменную длину, и его трудно будет непосредственно хранить в памяти. Поэтому мы сделаем значение указателем на очень длинный вектор битов; тогда каждая пара будет иметь фиксированный размер. Однако имеющийся в нашем распоряжении вектор не настолько длинен, чтобы мы могли позволить себе использовать каждую его часть только по одному разу. Таким образом, нужна программа для сбора ненужной памяти. Сейчас мы фактически описали традиционную схему размещения цепочек.

Итак, в конечном итоге нам нужны кроме алгоритма умножения и деления следующие вспомогательные подпрограммы:

Выделение памяти. Получая величину длина в качестве аргумента, эта подпрограмма возвращает указатель в вектор битов, который может использоваться как значение. Начиная с бита, на который указывает значение, расположено длина битов, которые не будут использоваться ни для каких других целей.

Возврат памяти. Исходными данными для этой подпрограммы служит пара длина и значение. Связанная с ними память освобождается для повторного использования. Эту подпрограмму необходимо вызывать всякий раз, когда длина какой-либо переменной меняется.