Термист Термомеханическое упрочнение арматурного проката технология, средства, разработка

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

Арифметические вычисления с высокой точностью

Чарлз Уэзерелл

Предыдущие разделы:
Арифметические вычисления числа π с высокой точностью. Общие соображения.
Как можно быстро умножать?
Алгоритм Тоома - Кука
Алгоритм быстрого умножения Тоома - Кука
Комментарии к алгоритму Тоома - Кука
Что можно сказать относительно деления?

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

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

Выберем теперь (не обязательно наилучший) ряд для вычислений, скажем
π = 16arctg(1/5) - 4arctg(1/239).

Мы фактически будем вычислять 2¹⁰⁰⁰π, поэтому хотелось бы вычислить 2¹⁰⁰⁰·16arctg(1/5). Первым членом соответствующего ряда будет 2¹⁰⁰⁰·16/5; назовем его a₁ (отметим, что a₁ складывается с суммой). Теперь, чтобы получить следующий член a_i+1 из a_i, поделим a₁на 5·5·(2i-1).
Примечание переводчика:
На самом деле a_i+1 = (a_i/5²)·((2i-1)/(2i+1)). Чтобы не выполнять умножение, можно хранить кроме a_i еще одно число b_i, равное (2¹⁰⁰⁰·16)/5^2i-1. Тогда переход к следующему члену осуществляется по формулам: b_i+1 = b_i/5², a_i+1 = b_i+1/(2i+1).
Если a_i добавлялся к сумме, то вычтем a_i+1 из суммы, если a_i вычитался, прибавим a_i+1. Затем поделим a_i на 5·5·(2i-1) (см. Примечание выше). Если a_i добавлялся к сумме, то вычисления заканчиваются, когда члены обоих рядов станут нулевыми. В результате получим примерно тысячу битов числа π. Результат, конечно, надо будет перевести в десятичную систему.

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

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

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

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

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

Уплотнение памяти. Эта подпрограмма должна просмотреть всю используемую память и попытаться объединить неиспользуемые отрезки вектора битов в более длинные отрезки. Обычно подпрограмма уплотнения вызывается в тех случаях, когда подпрограмма выделения памяти не смогла найти достаточно длинную цепочку последовательных битов. Поскольку такая возможность может не потребоваться при решении коротких задач, эту подпрограмму можно запрограммировать позднее. Существует много способов хранения информации о неиспользуемой памяти.

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

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

Вычитание. Эта подпрограмма аналогична подпрограмме сложения и выдает разность двух длинных чисел.

Подавление нулей. Исходными данными этой подпрограммы служит длинное число, а результатом должно быть более короткое длинное число, имеющее то же значение, но без старших нулей. Если окажется, что исходное число равно нулю, то результатом должно быть [1, 0].

Короткое умножение. Исходными данными служат два длинных числа длиной точно 32 бита; результатом должно быть их 64-разрядное произведение. Эту операцию можно выполнять справа налево, как в ручном методе.

Умножение длинного на короткое. Исходными данными служат длинное число и обычное число, по величине равное 64 или меньше; результатом должно быть их произведение в виде длинного числа. Эту операцию можно выполнять справа налево, как вручную.

Деление длинного на короткое. Исходными данными служат длинное число и обычное число, не превосходящее 64, а результатом должно быть длинное частное от деления длинного числа на короткое. Эту операцию можно выполнять слева направо, как это делается вручную.

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

Инструментовка. В качестве языка реализации сразу же приходит на ум Паскаль: в этом языке хорошие структуры данных и управляющие структуры. Однако Паскаль не позволяет легко переводить внутреннее битовое представление в битовые цепочки, доступные программисту, и обратно. Языки более низкого уровня - BLISS и XPL - обеспечивают более прямой доступ к ЭВМ за счет некоторой потери выразительности и надежности. Хорошая защищенность языков высокого уровня и доступ к машинному представлению сочетаются в PL/1, но обычно за это приходится расплачиваться временем выполнения. Для данного этюда важно также время, которое вы потеряете, пытаясь постичь некоторые весьма изощренные возможности PL/1. Интересной представляется реализация на Траке, поскольку в этом случае автоматически решается задача распределения памяти для цепочек.

Длительность исполнения. Одному исполнителю на 5 недель или двум на 3 недели.

Развитие темы. Как только у нас появляется арифметика высокой точности, сразу же возникает много интересных задач. Одна из них - точное вычисление числа e. Ряд для e особенно прост:
e = sum(1/i!) для i=0:∞,
где 0! = 1. Любой студент, изучающий математический анализ, может придумать еще очень много рядов и констант.

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

Кроме того, имеет смысл хранить длинные числа не в двоичной системе счисления, а в десятичной (конечно, не по одной в элементе массива, а по столько цифр, сколько помещается в обычном целом числе). При этом отпадает необходимость в программе перевода. Теперь арифметические программы могут работать несколько медленнее (но вовсе не обязательно; далеко не все компиляторы используют команды сдвига для умножения и деления на степени двойки), тем не менее в целом можно получить выигрыш в скорости, поскольку время работы алгоритма перевода длинных чисел из двоичной системы в десятичную (описанного у Кнута) пропорционально n², т. е. того же порядка, что и время всех остальных вычислений.

С помощью лишь этих программ сложения и деления можно вычислить многие математические константы: π, e, √2, дт(2) и т. д. Реализация такого усеченного варианта потребует, вероятно, не более одной человеко-недели. Сложные динамические структуры данных уже не потребуются; у нас будет всего два-три длинных числа известного размера, для представления которых вполне подойдут массивы Фортрана.

Читать дальше:
Литература

Подготовлено по публикации: Уэзерелл Ч. Этюды для программистов: Пер. с англ. - М.: Мир, 1982. - 288 с. Стр. 125 - 143.

К началу страницы

Web-сайт “Термист” (termist.com)
Термомеханическое упрочнение арматурного проката

Отсутствие ссылки на использованный материал является нарушением заповеди "Не укради"

Редактор сайта: Гунькин И.А. (termist.com@gmail.com)