Д. Кнут - Искусство программирования том 2 (3-е издание) - 2001 (Часть 1) (1119452), страница 106
Текст из файла (страница 106)
Поэтому вычисления с однократной точностью можно выполнять на первых пяти шагах, а операции с многократной точностью— только при вычислении значений -4ио + 11ео и 7ио — 19во. В этом случае получаем и = 1268728, в = 279726; дальнейшие вычисления можно продолжить подобным же образом с числами и' = 1268, е' = 280, и" = 1269, в" = 279 и т, д. При наличии большего накопителя можно было бы увеличить количество шагов, на которых вычисления выполняются с однократной точностью. Из примера видно, что в один сэожный шаг объединяются только пять циклов алгоритма Евклида, а если бы, скажем, размер слова равнялся десяти разря~шм, можно было бы объединить в один шаг до двенадцати циклов. Из результатов, доказанных в разделе 4.5.3, следует, что на каждой ктерации число циклов с многократной точностью, которые можно заменить циклами с однократной точностью, пропорционально размеру слова, используемому в вычислениях с однократной точностью. Метод Лемера можно сформулировать следующим образом.
Алгоритм 1 (Ало»репьи Евклида длл баоьщит чисел). Пусть и и е — представляемые с однократной точностью неотрицательные целые числа, такие, что и > е. Этот алгоритм вычисляет наибольший общий делитель чисел и и е, используя вспомогательные р.разрядиые переменные й, О, А, В, С, Р, Т, о, которые представлены с однократной точностью, и вспомогательные переменные г и в, которые представлены с многократной точностью. Ь1.
[Начальная установка.) Если число е достаточно мало, чтобы быть представленным в формате с однократной точностью, то йсй(и, е) вычисляется по алгоритму А, и на этом вычисления заканчиваются. В противном случае обозначим р ведущих разрядов числа и через й, а соответствующие разряды числа е— через е. Другими словами, если используется представление чисел в системе счисления по основанию Ь, то й +- [и/Ь"] и й +- [е/ЬЬ], где Ь вЂ” наименьшее возможное число, удовлетворяющее условию й с Ь". Присвоить А +- 1„В +- О„С +- О, Р о- 1.
(Эти переменные представляют коэффициенты в (20), где и=Аио+Вео и е=Сио+Рео в равносильных операциях алгоритма А над числамн с многократной точностью. Кроме того, и' = й+ В, е' ж й+ Р, ии = й+ А, е" = й+ С (30) в обозначениях рассмотренного выше примера.) Ь2. [Проверить частное.) Присвоить о+- [(й+А)/(в+С)]. Если й~ [(й+В)/(е+Р)], то перейти к шагу Ь4.
(На этом шаге проверяется выполнение условия д' р о" в обозначениях того же примера. В процессе вычислений на этом шаге при некоторых обстоятельствах, когда й = Ь» — 1 и А = 1 илн когда 6 = Ь» — 1 и Р = 1, может возникнуть переполнение при выполнении операций в формате с однократной точностью. В силу равенств (30) всегда будут выполняться условия О<0+А<Ь, 0<й+ССЬ, 0 < й + В < Ь», О < е + Р < Ь».
(21) Может оказаться, что выполнится одно нз равенств й + С = 0 и й + Р = О, но не оба одновременно. Поэтому попытка деления на куль на этом шаге означает "Перейти непосредственно к шагу 1А".) ЬЗ. [Имитация алгоритма Евклида.] Присвоить Т +- А — оС, А +- С, С +- Т, Т <-  — йР, В +- Р, Р +- Т, Т +- й — ве, й +- е, е +- Т и возвратиться к шагу Ь2, (Эти вычисления с однократной точностью равносильны операциям с многократной точностью в процедуре (28) с учетом (29).) Ь4. [Ш)а, иа котором выполняются вычисления с многократной точностью.] Если В = О, то, используя деление с многократной точностью, присвоить г +- и шоп е, и +- е, е +- Г. (Это может случиться только тогда, когда с помощью операции с однократной точностью нельзя моделировать операцию с многократной точностью.
Отсюда следует, что алгоритму Евклида требуется очень большое частное, что может произойти крайне редко.) В противном случае присвоить г +- Аи, Ф +- Ф+ Ве, щ <- Си, щ +- щ + Ре, и +- 1, е +- ю (выполняя непосредственно операции с многократной точностью). Возвратиться к шагу 1,1. ! С учетом неравенств (31) в процессе вычислений значения величин А, В, С, О представляются с однократной точностью. Для реалнзацнн алгоритма 1 требуется несколько более сложная программа, чем для алгоритма В, но прн опернрованнн болыпнмн числами этот алгоритм на многих компьютерах выполняется быстрее. Подобным образом можно ускорить выполнение бинарного алгоритма В в завершающей стадии (см. упр. 38).
Пренмущество алгоритма 1 заключается в следующем: он определяет последовательность частных, получаемых прн выполнении алгорнтма Евклида, что используется в многочнсленных приложениях (см., например, упр. 43, 47, 49, упр. 51 в разделе 4.5.3„ а также упр. 4.5,3-46). вАпалмэ бинарного алгоритма. В заключенне этого раздела для обоснованна установленных ранее формул проанализируем время выполнения алгоритма В.
Выясняется, что точно описать поведение алгоритма В крайне затруднительно, но можно начать аналнз этого алгорнтма с его приближенной модели. Предположнм, что числа и н о нечетны, и > о н (1би) ж т, (18о) = и. (32) (Таням образом, и является (т+1)-битовым числом, а о — (и+ 1)-бнтовым числом.) Рассмотрнм выполнение в алгоритме В цикла "вычнтанне н сдвиг", т.
е. операцню, которая начинается на шаге Вб, а прекращается после окончания выполнения шага В5. Каждый цикл "вычитание н сдвиг" прн и > о вычисляет разность и — о н сдвигает эту величину вправо до тех пор, пока не будет получено нечетное число и', которое замещает число и. Если входные числа случайны, можно ожидать, что прнмерно в половине случаев и' = (и — о)/2, примерно в четверти случаев и' = (и — о)/4, прнмерно в одной восьмой случаев и' = (и — о)/8 н т.
д. Получаем (ЗЗ) (18и') = т — й — г, где й — число разрядов., на которые было сдвинуто вправо число и — о, а г есть (18 и) — (18(и — о)~ — количество битов, потерянных слева во время вычктання числа о нз числа и. Заметим, что г < 1 прн т > и+ 2, а г > 1 прн т = и. Взаимосвязь между й н г довольно беспорядочная (см. упр.
20), однако Ричард Брент (Испагб Вгепс) нашел изящный способ анаанза поведения приближенной модели алгоритма, положнв и н о достаточно большими, такими, чтобы отношенне о/и имело непрерывное распределение прн дискретном изменении к. [См. А)йог11)ппв апс( Сошр!ех11у, ед1ьед Ьу Л, Г. Тгац(з (Хе»" Уог1с Асабеш(с Ргевэ, 1976), 321-355.) Предположим, что и н о — большие целые числа, возможно, случайные, но обязательно нечетные н нх отношение подчинено определенному закону распределения. Тогда на шаге Вб младшие значашне биты велнчнны 1 = и — о могут быть случайными, но велнчнна 1 будет четной.
Следовательно, 1 будет нечетно кратным 2ь с вероятностью 2 "; это приближенная вероятность того, что в цикле "вычитание н сдвиг" потребуется выполнить Й сдвигов вправо. Другими словамн, получено подходящее приближение, опнсывающее поведение алгоритма В прн сделанных предположениях о том, что переход от шага В4 к шагу ВЗ всегда будет происходить с вероятностью 1/2. Пусть 6'в(х) — вероятность того, что ппп(и, о) /гпах(и, о) будет > х после выполнення с учетом этих предположений и циклов "вычнтанне н сдвиг".
Если и > о н если выполнено точно й сдвигов вправо, то отношение Х = и/и изменится на Х' = пйп(2«и/(и — е), (ц — е)/2«в) = п1!п(2«Х/(1- Х), (1 — Х)/2«Х). Таким образом, неравенство Х' > х будет справедливо тогда н только тогда, когда 2«Х/(1 — Х) > х н (1 — Х)/2«Х > х, а это то же самое, что н 1 1 — < Х < —. 1+ 2"/х 1+ 2«х (34) Поэтому С„(х) удовлетворяет интересному рекуррентному соотношению С +1(х) = ~2 «(«"'„( — ) -6' ( — )), (33) 1«(х) =~~',2 (б( —,) -С( — „)) прнО<*< 1; (36) «>1 С(О) ю 1; 0(1) = О.
(37) Пусть 5( ) 1 Ц ( 1 ) 1 ~ ( 1 ) 1 ( 1 ) 2 «с ( „ ); (38) тогда 0(х) = б(1/х) — йх). (39) Естественно определить 6(1/х) = -~(х), (40) так что уравнение (39) справедливо для всех х > О, Поскольку х изменяется от 0 до оз, Я(х) увеличивается от 0 до 1. Следовательно, 0(х) уменьшается от +1 до -1. Конечно, прн т > 1 функция 1«(х) перестает быть вероятностью, тем не менее она имеет смысл (см, упр. 23). допустим, что имеются степенные ряды а(х), )1(х), ъ (х), б (х) Л(х) и(х) и (х), г (х) и р(х), такие, что С'(х) = а(х) 18х+ Р(х) + ~ (Уп(х)сое2нп«1$х+ Ач(х) сбп2тп«)бх) Я(х) = Л(х))бх+д(х)+ ~~ (о (х)ссм2ягп18х+т (х)в(п2ягп)бх), ~аж1 р(х) = «'(1+х) = р1х+ртх +рзх + р«х + р«х + рех + (41) (42) (43) где Се(х) = 1- х прн 0 < х < 1, Проведенные вычислительные эксперименты показали, что П„(х) быстро стремятся к предельному распределению С (х) = С(х) несмотря на то, что формальное доказательство скоднмости представляется неочевидным.
Будем полагать, что существует функция распределения ««(х). Тогда она удовлетворяет уравнению 4,0 т з.о Рис, 10. Предельное распределение отношений в бинарном алгоритме вычисления на4~- болынего общего делителя. поскольку можно показать, что этим свойством при и > 1 обладает решение ст„(х) уравнения (35) (см., например, упр. 30). Эти степенные ряды сходятся при 14 с 1. Какие выводы можно сделать относительно а(х), ..., р(х) в уравнениях (36)-(43)7 Прежде всего, из уравнений (38), (40) н (43) имеем Соответственно равенство (42) выполняется тогда и только тогда, когда (45) (46) (47) 2Л(х) = Л(2х); 2д(х) =д( )+Л(2х) -42 ')1 2а„,(х) = ам(2х), 2тщ(х) = т„,(2х) прн гп > 1. Из (45) следует, что Л(х) есть просто константа, кратная х; так как она отрицательна, можно записать Л(х) = — Лх.
(48) (соответствую1ций коэФФициент приобретает вид Л = 0,39792 26811 88316 64407 67071 61142 65498 23098+, (49) но способ, которым его можно вычислить, неизвестен.) Из уравнения (46) следует, что р1 = -Л и 2да = 24144 — 2арь прн Й > 1; другими словами, Из уравнения (47) следует также, что оба степенных ряда (51) т (х) = т х иш(Х) = 4гшХ„ являются просто линейными функциями, (Но это не распространяется на функции Ъ (х) и А (х).) Если в уравнении (44) заменить х на 1/2х, то получим (52) 28(1/2х) = 8(1/х) + 6(х/(1 + х)), 1.0 О.З ои б/т/ О.О -ой -ол -од -о.з -1.0 О.о 28(х) = С(1/(1+ 2х)) + Я(2х) = 8(2х) — р(2х). (44) дь = ра/(1 — 2 ) при й > 2.