Д. Кнут - Искусство программирования том 2 (3-е издание) - 2001 (Часть 1) (1119452), страница 94
Текст из файла (страница 94)
Штрассеном (У. Бсгвввеп), который предложил способ умножения больших чисел, более быстрый, чем это было возможно с использованием любых известных к тому времени методов. Позже вместе с Шенхаге он уточнил метод и опубликовал модифицированные алгоритмы в Сошриппй 7 (1971), 281-292. Похожие идеи, но для произвольных целых чисел, были независимо от Ф. Штрассеиа предложены Дж. М, Поллардом (3. М. Ро11егд), Ма»Ь, Сотр. 25 (1971), 365-374. Чтобы понять значимость их вклада в решение проблемы, рассмотрим для начала механизм быстрого преобразования Фурье.
Для данной последовательности К = 2" комплексных чисел (ие,..., ик а) и данного комплексного числа (37) ы = ехр(2ха(К) последовательность (Оо,..., Ок 1), определенная выражением (35), может быть бы- стро вычислена по следующей схеме. (Параметры з и 1. в зтих формулах равны либо О, либо 1, так что на каждый "проход" приходится 2" злементарных вычисли- тельных операций.) Проход О. Пусть А(е) (1» ы..., ге) = на, где 1 = (8ь а ...
»е)». Проход 1. Присвоить А(0(зз ы 1»»,..., 1е) +- А(е)(0,»ь-»м",ео)+ьа» '" 'А(е)(1 »ь-»,".,»е). Проход 2. Присвоить Айй(зь ызз»,1» з>.-,»о) +- А%(зз,о,»з-», ° . »о) +ьа» 1 а а-')аА(~1(зг, ы1,»ь-»1 ° ° °,»е) Проход й, Присвоить А(»1(зь ы,зызе) а- А(з-В( О) + ь,(аааа .аа-О»А(з-11(зь ы...,з„1), По индукции очень просто доказать., что 16 1 ) ) ' »а а(аа-а.ааа-~)ана-а.:аа-а)аи (36) е<аа,,...,а, а<а где 1 = (1ь-» 1»»о)», так что А(»1(зь ы...,зызе) и Й„где з= (зез)...зь-а)» (39) (Важно отметить, что двоичные цифры числа з в конечном результате (39) обращены.
В разделе 4.6,4 приводится дальнейший анализ такого рода преобразований.) Прежде чем приступить к поиску обратного преобразования Фурье (ие,... ...гик 1) по значениям (йе,...,йк 1), заметим, что "двойное преобразование" имеет вцд 6 = ~у ы"Оа= о<а<к . е< л<к 1 % ' а(а+г)~ ъ. па~ ~ ьа ~ = аап(гг) зааек1 о<а<к е«к (40) так как длЯ», кРатных К, сУмма геометРической пРогРессии 2;о<а к ьааа Равна нулю, Позтому обратное преобразование может быть получено точно так, как и прямое, за исключением того, что конечный результат делится на К и немного сдвинут. Возвращаясь к проблеме умножения целых чисел, предположим', что нужно вычислить произведение двух и-битовых целых чисел н и ш Будем оперировать (как и в алгоритме Т) группами битов.
Положим 2п < 2~1 < 4п, К = 2~, Ь м 2 (41) и запишем и = ((/к — 1... У1Уо)с, о = ()гк 1 .. ° КРо)с, (42) полагая числа и и о К-разрядными числами по основанию Е, так что калсдый из разрядов (/1 Или $' есть 1-битовое целое число. Поскольку 2""'Ч > и, ведущие разряды в У и гу для всех у > К/2 равны нулю, Соответствующие значения для й и ( будут выбраны позже, когда в процессе решения задачи прояснится общая ситуация и можно будет при выборе и и 1 учесть всю имеющуюся информацию.
Следующим шагом в процедуре умножения является вычисление преобразований Фурье (йо,...,йк 1) н (Оо,...,йк 1) последовательностей (ио,...,ик 1) н (ио, о,к-1), где по определению ив = (/в/2"г', о, = К/2 +~, (43) Для удобства масштабирование выполнено так, что любые и~ н ос меньше 2 ", а это, в свою очередь, обеспечивает то, что абсолютные значения (й,( и (О,( оказываются меньше 1 для всех ж Здесь возникает очевидная проблема, связанная с тем, что комплексное число ы не может быть точно представлено в двоичной системе счисления. Как же вычислить достоверное преобразование Фурье? Провидению было угодно, чтобы результат получился правильным даже в случае, если в процессе вычислений предъявить весьма скромные требования к точности представления чисел, Оставим на время эту проблему и предположим, что вычисления выполняются с бесконечной точностью, Анализ необходимой точности будет приведен несколько ниже (через пару страниц). Для полученных й, и О, положим Ф, = й,й, (О < в < К) и определим обратное преобразование Фурье (юо,, вок ~).
Как было разъяснено выше, получим юг = ~~' %ов = ~~' ~'ю11/2 так что целые числа Иг„= 2то+мю„являются коэффициентами искомого произве- дения и о=И"к-о1 + "+Иг1~+Иго (44) Поскольку 0 < И'г < (г + 1).У < КХт„каждое И'„содержит не более й + 2( бит, поэтому нетрудно получить двоичное представление для известных величин И'"ю если только й не слишком велико по сравнению с 1. Например, пусть нужно умножить и = 1234 на о = 2341 при й = 3 и ( = 4. Вычисление изображения (йо,..., йг) по оригиналу и выполняется следующим образом (см.
(12)) . (г, в, В) = (О, О, О) (О, О, 1) (О, 1, О) (О, 1, 1) (1, О, О) (1, О, 1) (1, 1, О) (1, 1, 1) 2вА1о'(г,в,в) = 2 13 4 О О О О О 2"АО~(г,в,в) = 2 13 4 О 2 13 4 О 2" А91(г,в,в) = О 13 — 2 13 2+ 41 13 2 — 41 13 2~АОО(г,в,1) = 19 — 7 — 2+131 — 2 — 131 а+3 о.— 3 а — л а+3 Здесь а = 2+4т, 3 = 13ы и ы = (1+1)/Я. Это дает нам колонку й, в табл. 1. Данные в колонке о, получаются из о аналогично. После этого находим й„умножив й, на О,. Выполняем еще одно преобразование, учитывая соотношение (40), и получаем ю, и Иг,. Итак, мы получили свертку (19), но на этот раз используя комплексные Таблица 1 УМНОЖЕНИЕ ПРН ПОМОШ14 ДИСКРЕТНОГО ПРЕОЕРАЗОВАНИН ФУРЬЕ 2 зьь = И»» 1О 69 64 125 Зб о о о 2 9.» 16 5 + 9» + 24» -4+ 2з 5-9» -2Ф 1г 5+ 9з — 24» -4 — 2» 5 — 91+ гб» числа вместо того, чтобы оставаться верными методу, оперирующему только целыми числами.
Попробуем теперь оценить время, необходимое этому методу при больших числах, если использовать пз-битовую арифметику в формате с фиксированной точкой. Из упр. 10 видно, что в процессе преобразования все величины А(й будут меньше 1 из-за масппабирования (43). Следовательно, достаточно иметь дело только с дробными пз-битовыми частями (,а 1... а ьз)2 для действительной и мнимой частей любых промежуточных результатов.
Так как входные величины из и сз являются действительными числами, можно упростить вычислительную процедуру, а именно — вместо 2К действительных значений на каждом шаге преобразований использовать только К значений (см. упр. 4.6.4-14). Чтобы не усложнять выкладки, далее эти тонкости учитывать не будем. Прежде всего необходимо вычислить ьз и ее степени.
»Тля простоты сформируем таблиззу значений ь», ..., ы . Положим зс» = ехр(2хз/2") (45) таким, что ьзз — — -1, ьзз = з» ьзз = (1+ з)/тз»2, ..., ь»ь = ьз. Если ьз„= я„+ зр„то ИМЕЕМ Ь»».„1 — — Х»4З + ЗР„+1, ГДЕ 1+ я„п, х.+1 = ~(( — Р»+1 =— (46) 2 ' " 2я ~1 [См. Я. В.. Тате» 1ЕЕЕ Тгапбасбопб ЯР-43 (1995), 1709-171Ц На вычисление ь»1, а»2,..., ь»8 затрачивается по сравнению с остальными операциями относительно немного времени, так что вполне можно использовать любой стандартный способ извлечения квадратных корней.
После ь», можно легко вычислить все степени ь»з, учитывая, что з»»~ = з»»~~ ...з»»~з~ з 1 если 7 — (29» 1 ...2120)2. (47) Поскольку каждое значение ьзг получается как результат умножения ьз, не более чем й раз, ошибки вычислений при таком методе не накапливаются. Общее время вычисления всех значений ьзз равно 0(КМ), где М вЂ” время выполнения операций умножения зп-битовых комплексных чисел, поскольку на вычисление каждого нз следующих значений ь»1 по известному предыдущему требуется одна операция умножения.
На выполнение последующих операций необходимо больше циклов, чем 0(КМ), поэтому на вычисление степеней ьз требуется сравнительно мало времени. гзо 19 2+ 4»+ 13И -2+ 13» 2-4» — 1ЗФ -7 2+ 4» — 13И -2-1З 2- 4»з + 13а 2»449 304 -26+ 64»+69И - 12$Ф -18 — 56» -26 — 64» + 125И вЂ” 69Ф -84 -26 + 64з — 69м + 125И -18+ 56» -26 — 64» — 125м * 69а гзззз» 80 о о 0 288 1000 ыг 552 Каждое из трех преобразований Фурье состоит из Ь проходов, а каждый из проходов включает К операций вида а з- Ь + ызс, так что общее время вычисления преобразований Фурье равно 0(ЬКМ) = 0(Мпй/1).
В итоге для вычисления двоичных разрядов и в согласно (44) требуется О(К(Ь+1))= 0(п + пй/1) машинных циклов. Если просуммвровать все операции, получим, что общее время умножения и-битовых чисел н и е равняется 0(п) + 0(Мпй/1) циклам. Теперь, зная, какой должна быть величина М, посмотрим, какой должна быть точность вычисления промежуточных результатов пз. Для простоты вместо поиска наилучших оценок точности будем рассматривать оценки с запасом, гарантирующие получение при помощи этого алгоритма требуемых результатов. Достаточно вычислить все значения ыб в таком виде, чтобы приближения (м")' удовлетворяли неравенству )(ыз)') < 1. Это условие легко обеспечить, если значения не округлять, а усекать до нуля, так как в (46) хз+1 + уз+, = (1 + хз + р~~ + 2яг)/(2+ 2х,).
Все операции с пьбитовыми числами в комплекснозначной арифметике с фиксированной точкой выполняются посредством замены точных вычислений вида а <- Ь+ ызс приближенными: а' +- усечение(Ь + (ыз)'с'), (48) где Ь', (ыт)' и с' — приближения, вычисленные предварительно. Все зти комплексные числа и их приближении по абсолютному значению ограничены единицей. Если )Ь'-Ь| < бы )(ьзз)'-ыз( < бз и )с'-с) < бз, то нетрудно увидеть, что будет выполняться )а' — а! < б + б1 + бз + бз, где б = )2 ™+2 ™з(= 2~~ ~,, (49) поскольку имеем )(Ы)'с' — ызс( = ~((ыз)' — Ы)с'+му(с' — с)~ < бз + бз, а б превышает максимальную погрешяость усечения.
Вычисления приближений (ыз)' начинаются с приближений ~>,' к числам, определенным в уравнениях (46), и можно считать, что вычисления в (46) выполняются с допустимой погрешностью, такой, что выполняется )ы,', — ы„( < б. Тогда из уравнения (47) следует, что )(ыз)' — мз ) < (2Й вЂ” 1)б прн всех /ц поскольку ошибка обусловлена не более чем Ь приближениями и не более чем Ь вЂ” 1 усечениями. Если перед началом любого прохода быстрого преобразования Фурье ошибка не превышает е, то операции етого прохода имеют вид (48), где бз = бз = е и бз —— (2Й-1)б, Тогда ошибки после выполнения прохода не будут превышать 2е+2Ьб.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
