Каждое из трех преобразований Фурье состоит из fc проходов, а каждый из проходов включает К операций вида а <- 6 + ш-с, так что общее время вычисления преобразований Фурье равно

0{кКМ) = 0{Мпк/1).

В итоге для вычисления двоичных разрядов u-v согласно (44) требуется О(Ar(fc+Z))= 0{п +пк/1) машинных циклов. Если просуммировать все операции, получим, что общее время умножения п-битовых чисел ииу равняется 0{п) + 0{Мпк/1) циклам.

Теперь, зная, какой должна быть величина М, посмотрим, какой должна быть точность вычисления промежуточных результатов тп. Для простоты вместо поиска наилучших оценок точности будем рассматривать оценки с запасом, гарантирующие получение при помощи этого алгоритма требуемых результатов. Достаточно вычислить все значения ш- в таком виде, чтобы приближения (ш-) удовлетворяли неравенству Kw-)] < 1. Это условие легко обеспечить, если значения не округлять, а усекать до нуля, так как в (46) xj + у = {I + х + у + 2x)1 {2 -f- 2х). Все операции с т-битовыми числами в комплекснозначной арифметике с фиксированной точкой выполняются посредством замены точных вычислений вида а <- Ь + сос приближенными:

а <- усечение(6 + {иУс), (48)

где 6, (ш) и с - приближения, вычисленные предварительно. Все эти комп.лекс-ные числа и их приближения по абсолютному значению ограничены единицей. Если 6-6 < Si, \{u!y-u!\ < ($2 и с-с < дз, то нетрудно увидеть, что будет выполняться о -а\<6 + б1+б2 + 6з, где

6=\2-" + 2-"Ч\ = 2/-", (49)

поскольку имеем (ш-)с- ш-с = ((ш-)- а;)с+ а;-(с- с) < 62 + 6з, а. 6 превышает максимальную погрешность усечения. Вычисления приближений (w-) начинаются с приближений wJ. к числам, определенным в уравнениях (46), и можно считать, что вычисления в (46) выполняются с допустимой погрешностью, такой, что вьшолняется - < (5. Тогда из уравнения (47) следует, что (ш-) - w- < (2fc - 1)6 при всех j, поскольку ошибка обусловлена не более чем fc приближениями и не более чем fc - 1 усечениями.

Если перед началом любого прохода быстрого преобразования Фурье ошибка не превышает е, то операции этого прохода имеют вид (48), где 6i = 63 = е и 62 = {2к-1)6. Тогда ошибки после выполнения прохода не будут превышать 2е+2к6. Так как на нулевом проходе ошибок нет, по индукции можно показать, что после j-ro прохода ошибка будет ограничена (2- - 1) • 2к6 и вычисленные значения us будут удовлетворять неравенству \{us) ~ us\ < (2* - 1) • 2к6. Аналогичное неравенство можно получить и для (vs); тогда

\{w,y - Ws\ < 2(2* -1)-2к6 + 6< (4fc2= - 2к)6.

В процессе обратного преобразования накапливаются дополнительные ошибки, однако большинство из них подавляется при делении на АГ = 2*. По этой же причине приходим к выводу, что вычисленные значения величин Wr будут удовлетворять неравенству

{ur-y-wr <2*(4fc2*-2fc)(5-f (2*-l)2fc(5; ш; - < 4fc2*d. (50)

Для округления 22*+2«;{. до правильного целого числа Wr необходимо обеспечить достаточную точность вычислений, поэтому должно выполняться неравенство

22k+2l+2+lsk+k+l/2-m 1 (gj)

Т. е. m > 3fc + 2/ + Ig fc + 7/2. Для этого достаточно потребовать, чтобы выполнялось

fc > 7 и m > 4fc + 2/. (52)

Соотношения (41) и (52) могут быть использованы для определения таких значений параметров к,1ит, чтобы на операцию умножения затрачивалось 0{п)+0(Мпк/1) единиц времени, где М - время умножения т-битовых дробных частей.

Например, для компьютера MIX предположим, что перемножаются двоичные числа, каждое из которых представлено п = 2 = 8192 бит. Можно выбрать значения параметров fc = 11, / = 8, m = 60, так что требуемые операции с т-битовыми числами будут ни чем иным, как операциями умножения с числами в арифметике с удвоенной точностью. Поэтому необходимое время выполнения М операций умножения т-битовых комплексных чисел будет сравнительно малым. Если положить, к примеру, fc = / = 15, то придется иметь дело с операциями утроенной точности, что выходит за пределы возможностей памяти машины MIX. На компьютере с большими возможностями можно, положив fc = / = 27 и m = 144, перемножить два гигабитовых числа.

Дальнейший анализ выбора значений параметров fc, / и m приводит к удивительному выводу: для большинства практических случаев можно считать М постоянным, а время выполнения процедуры умножения Шёнхаге-Штрассена пропорционально п. Суть в том, что можно выбрать fc = / и m = 6fc; такой выбор fc приводит к тому, что время выполнения всегда меньше Ign, поэтому потребуется не более чем ушестеренная точность вычислений, даже если п больше размера машинного слова. (В частности, п может быть больше размера индексного регистра, так что, возможно, числа и и и не поместятся в основной памяти.)

Таким образом, практическая сторона быстрого умножения решена, за исключением возможного усовершенствования процедуры выбора констант. Действительно, алгоритм свертки для произвольных целых чисел, приведенный в упр. 4.6.4-59, является лучшим практическим решением проблемы умножения с высокой точностью. Наш интерес к проблеме умножения больших чисел отчасти теоретический, поскольку интересно исследовать предельные ограничения на вычислительную сложность алгоритма. Итак, оставим на время практическую сторону вопроса и предположим, что число п чрезвычайно велико, возможно, значительно больше числа атомов во Вселенной. Можно считать т приблизительно равным 6Ign и использовать для умножения т-битовых чисел тот же алгоритм рекурсивно. Время его выполнения будет удовлетворять выражению Г(п) = 0(nr(logn)). Следовательно,

Г(п) < С п(С Ig п) (С Ig Ig n) (С Ig Ig Ig n)..., (53)

где произведение продолжается до тех пор, пока множитель Ig... Ign < 2.

Шёнхаге и Штрассен в своей работе показали, как улучшить эту теоретическую верхнюю границу до 0{п log п log log п), используя целые числа ш, чтобы распространить быстрое преобразование Фурье на целые числа по модулю 2* + 1. Эта верхняя граница применима к машине Тьюринга, т. е. к компьютерам с ограниченной памятью и конечным числом лент произвольной длины.

Для более мощного компьютера со случайной выборкой любого числа слов ограниченного размера Шёнхаге обратил внимание на то, что верхняя граница снижается до О(пlogn). Теперь можно выбрать параметры fc = / и m = 6fc и появляется время для формирования полной таблицы всевозможных произведений ху при О < а;,2/ < 2"/21 (Количество таких произведений равно 2* или 2*+. После этого можно вычислить любой элемент таблицы за О (к) шагов, добавив к нему одного из предшественников.) Следовательно, 0(fc2*) = 0{п) шагов будет достаточно для вычисления. В этом случае М - время, необходимое для выполнения 12-разрядной" арифметической операции по основанию 2\ а отсюда следует, что М = 0{к) = O(logn), так как 1-разрядное умножение можно выполнить по таблицам. (Время выборки слова из памяти считается пропорциональным количеству битов, содержащемуся в адресе слова.)

Более того, в 1979 году Шёнхаге обнаружил, что машина с указателями может вьшолнять умножение п-битовых чисел за 0(п) шагов (см. упр. 12).

Такие устройства (которые называются также машинами с модификацией хранения и связными автоматами) при п - оо реализуют, как нам кажется, лучшие вычислительные модели, что было описано в разделе 2.6. Таким образом, можно сделать вывод о том, что умножение за 0(п) шагов возможно как на практике, так и теоретически.

В 1986 году Д. В. Чудновский, Г. В. Чудновский, М. М. Денно (М. М. Denneau) и С. Г. Юнис (S. G. Younis) сконструировали необычный компьютер общего назначения, названный Маленьким Ферма (Little Fermat), который мог быстро умножать большие числа. Компьютер оснащен аппаратными средствами быстрого выполнения арифметических операций по модулю 2 + 1 над 257-битовыми словами. Тогда свертка массивов, состоящих из 256 слов, может выполняться с помощью 256 умножений однословных чисел вместе с тремя дискретными преобразованиями, требующими только операций сложения, вычитания и сдвига. Это позволило, основываясь на конвейерном принципе организации цикла длительностью примерно 60 НС, перемножить два 10-битовых целых числа меньше чем за 0.1 с [Ргос. Third Int. Conf. on Supercomputing 2 (1988), 498-499; Contemporary Math. 143 (1993), 136].

D. Деление. Теперь при наличии эффективных программ для умножения рассмотрим обратную проблему. 0казывс1ется, что деление может быть выполнено так же быстро, как и умножение, с точностью до постоянного множителя.

Чтобы разделить п-битовое число и на п-битовое число v, можно сначала найти п-битовое приближение к числу 1/v, затем умножить его на и, что даст приближение q Ku/v, и наконец выполнить еще одно умножение для внесения небольшой коррекции в q, чтобы убедиться, что выполняется неравенство О <и - qv < v. Исходя из сказанного, достаточно иметь эффективный алгоритм, который формировал бы по заданному п-битовому числу приближенное значение числа, обратного п-битово.му числу. Это может быть реализовано следующим алгоритмом, который использует "метод Ньютона", рассмотренный в разделе 4.3.1.

Алгоритм R (Получение обратной величины с высокой точностью). Пусть число V имеет двоичное представление v = (Сихизз • • гд i = 1- Этот алгоритм вычисляет приближение z числа 1/v, такое, что

\z - l/v\ < 2-". (54)