Этюды для программистов [неполностью, главы 1–24]

6. (Как умножить два 32-разрядных числа?) На шаге 9 требуется умножить два 32-разрядных числа, получив 64-разрядное произведение, причем оба сомножителя обязательно положительны. На многих ЭВМ имеется аппаратная возможность такого умножения, но результат нельзя получить, пользуясь языками высокого уровня. Ну и, конечно, некоторые ЭВМ не имеют подобней аппаратуры. Поэтому для выполнения этого умножения нужно написать подпрограмму, причем она должна быть эффективной, поскольку время работы алгоритма определяется главным образом временем умножения 32-разрядных чисел. Вероятно, достаточно хорошим методом будет разбиение чисел на части и моделирование ручного способа умножения. Тем не менее если нужно получить произведение uv и число и записано в виде u₁·2¹⁶ + u₀, a v — в виде v₁·2¹⁶ + v₀, то произведением будет

(2³² + 2¹⁶)u₁v₁ + 2¹⁶(u1 − u0)(v0 − v1) + (2¹⁶ + 1)u₀v₀.

Вычисление по этой формуле выполняется при помощи только 16-битовых вычитаний и умножений, а также некоторых сдвигов и сложений. Обратите внимание, что одно умножение сэкономлено.

Что можно сказать относительно деления?

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

1. (Выбор размера обратного.) Найти наименьшее j, такое, что 2^j ≥ max (m, 2n). Присвоить к значение 2^j−1.

2. (Нормализация v.) Присвоить [k, v] значение 2^k−n [n, v]. На этом шаге мы сдвигаем v влево, чтобы оно заняло k битов, причем левый бит был бы равен 1. Присвоить [2, а] значение [2, 2].

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 добавлялся к сумме, то вычисления заканчиваются, когда члены обоих рядов станут нулевыми. В результате получим примерно тысячу битов числа π. Результат, конечно, надо будет перевести в десятичную систему.

33

В алгоритм, вероятно, необходимо внести следующие изменения:

a) на шаге 1 заменить max (m, 2n) на max (2m − 2n, 2n);

b) на шаге 4 заменить 2^3·2i на 2^3·2i−1.

Прежде чем приступать к программированию алгоритма Тоома—Кука или алгоритма деления, рекомендуем тщательно разобраться в них, ознакомившись с теорией, например по книге Кнута, неоднократно цитируемой здесь — Прим. перев.

34

На самом деле a_i+1 = (a_i/5²) · (2i − 1)/(2i + 1)). Чтобы не выполнять умножение, можно хранить кроме a_i еще одно число b_i, равное (2¹⁰⁰⁰ × 16)/5²ⁱ⁻¹. Тогда переход к следующему члену осуществляется по формулам: b_i+1 = b_i/5², a_i+1 = b_i+1/(2i + 1). — Прим. перев.