Возведение числа в степень... (возможные варианты и их улучшение)

Не ноль, а по модулю 2^16, т.к. pmullw перемножает ворды и сохраняет только младший ворд результата (lw = low word). Для двордов юзай pmuludq (SSE2). Хотя выигрыш от использования ММХ\SSE можно получить только на младших моделях P4 (Northwood и ниже), у которых нет целочисленного умножителя. На Prescott'ах от этого ничего не выиграешь, а на P6 и атлонах с их быстрыми умножителями только проиграешь. Так что - никакого простора для оптимизации

 

leo c SSE вообще ничего не получается:

Код (Text):

1. function sse_power(x,y:longint):longint;
2. begin
3. y:=y-1;
4.  asm
5.     cvtsi2ss  xmm0,x
6.     cvtsi2ss  xmm1,x
7.     mov ecx, y
8.   @rep:
9.    pmuludq xmm0, xmm1
10.    dec ecx
11.   jnz @rep
12.    movd @result, xmm0
13.    emms
14.  end;
15. end;

ноль и все... я явно что то не правильно делаю, но справочника по SSE у меня нету...

 

Ес-но, т.к. pmuludq умножает целые беззнаковые дворды, а ты неизвестно зачем вместо movd используешь cvtsi2ss, преобразуя целые в вещественные )

Ну и как ты без него собираешься (шкодить, наугад cvtXXX тыкать, где попало ? )
Скачай IA-32 и будет тебе счастие Плюс в самой мнемонике команд заложена подсказка, например, (почти) все MMX\SSE операции с целыми числами начинаются с p, а заканчиваюся указанием типа операндов. Например, pmuludq означает целочисленное умножение unsigned dwords -> qword

 

leo да ты прав выигрыша на Прескотте-2м нет. Спс за подсказки...

 

Хорошо иметь под рукой таблички латентностей от А.Фога и IA-32 Optimization, плюс до кучи AMD Optimization'ы по Athlon XP и 64. Тогда без проб и ошибок будет видно, где можно получить выигрыш, а где игра не стоит свеч. В частности mul\imul в P6 (включая Core 2) и AMD выполняются всего за 3-4 такта и никакими ухищрениями выполнить умножение быстрее не получится, кроме тривиальных случаев shl\lea при умножении на ограниченный набор констант или при векторной обработке нескольких чисел на MMX\XMM. А затачивать код под тормозные P4 смысла не имеет, т.к. жить им все равно осталось не долго

 

Никто не занимался оптимизацией умножения _больших_ целых цисел (порядка сотен и более двоичных разрядов)? и как посоветуете, в BCD или прямо так?

 

я занимался.
Прямо так и умножал - арифметически, столбиком перемножая DWORD-ы чисел.
Причем получалось что карацуба быстрее только при числах с битностью больше чем примерно 32768 бит:

Писано на сях, посему есть еще куда оптимизировать.

Кста до сих пор не понимаю - какой смысл в BCD и в строковом представлении больших чисел при их обработке? Число в памяти можно представить как записанное по модулю 2^32 а не по модулю 10 и работать с ним в несколько раз быстрее... Сколько уже библиотек видел которые сперва переводят число в десятичную строку а потом с черепашьей скоростью начинают ее мучать.

PageFault,
Ффух, пасиба - успокоил.
А то я тут груууустно так фтыкал на скорость генерации safe primes. Ну, в таком случае мой генератор простых не уступает в скорости CryptoAPI-шному Есть чему порадоваться...

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

 

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

 

CreatorCray

Если брать 32-битную арифметику, то оптимальным будет представление по базису 2^15. Обычно берут ближайшую степень десятки, чтобы не мучаться при переводе в экранное представление, то есть 10000. Такой базис например стоит(стоял? в новых версиях не знаю..) по умолчанию в Maple.

 

Vov4ick
BCD - жуткий тормоз сам по себе (фича интела), да и красивого алгоритма на них по определению не построишь
в представлении 2^32 у тебя дворд - это цифра, комманда mul - таблица умножения после которой "eax пишем, edx "в уме"", а дальше обычное умножение в столбик, только вместо десятичных цифр дворды

 

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

 

Stiver
>> оптимальным будет представление по базису 2^15

Практическим способом проверено - быстрее работает 2^32. Знак числа не нужно хранить в нем же. Я храню отдельным битовым флагом - все равно бОльшая часть практического применения этой математики - криптография. А там знак нужен в достаточно скромном кол-ве операций. Вот и реализовано у меня все через Unsigned. А все Signed операции юзают вызов Unsigned, предварительно рассчитав знак у результата.

>> Обычно берут ближайшую степень десятки, чтобы не мучаться при переводе в экранное представление

"Обычно" для каждого свое. Для меня "обычно" большие числа вообще выводить не надо. А если надо то HEX представления хватает за глаза. И беру я "обычно" степень двойки, т.к. мне результаты вычислений надо в обозримом будущем. Как думаешь, что быстрее - перевести в десятичную форму 1 результат кучи вычислений или мееедленно вычислять всю эту кучу формул но зато не потребуется перевод результата? По мне так лучше быстро вычислить...

Vov4ick
>> перевод может занять больше времени, чем собственно счёт

Ну, если это сложение двух чисел - то да... может...
Ну а если это RSA 4096 то...

 

CreatorCray

Хм, было бы действительно интересно посмотреть код сравнения. Ведь результат простых операций(сложение, умножение и т.д.) над 32-битными числами не влезает снова в 32 бита, нужно обрабатывать переполнение. При базисе 2^15 переполнения не происходит.

Не совсем понял, какой знак? Все "цифры" в представлении числа всегда неотрицательны.

Конечно, для каждого отдельного случая оптимальное представление может отличаться.

 

Stiver

При базисе 2^32 тоже не происходит - младшая часть в eax, старшая в edx - это гораздо удобнее для "умножения в столбик", чем 2^15

приспособить к цепочечному умножению больших чисел знаковую команду imul гораздо сложнее, чем беззнаковую mul, потому лучше множить модули, а затем учитывать знак

А где, если не секрет, кроме криптографии (где на экран их выводить не нужно) реально нужны большие числа?
Для научных и финансовых рассчётов 80 бит девать не куда, а с ними и FPU прекрасно справляется

 

Stiver
>> результат простых операций(сложение, умножение и т.д.) над 32-битными числами не влезает снова в 32 бита

Сложение 2-х 32-битных чисел => 32 битное число + бит переноса в CF, который будет прибавлен к следующей паре чисел путем adc

Умножение 2-х 32-битных чисел => 64 битное число в edx:eax

Как видишь все в проце есть чтоб с максимальным удобством юзать именно 32 бита

 

Угу, именно поэтому я и оговорил выше: Если брать 32-битную арифметику. Если брать 64 бита, да еще и конкретную реализацию современных x86, то и базис можно взять другой. Хотя Maple тоже сильно не дураки писали..

Y_Mur

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

 

так нулевая степень и так за циклом: цикл завершиться так и не начавшись и функа выдаст 1
циклов по идее столько сколько бит в числе, так что о каком лишнем цикле речь??
хотя да при степени ноль будет лишнее умножение, но это легко исправить)

 

UbIvItS
Ну, начнем с того, что и у меня и у тебя в коде ошибка )

мой код:

Код (Text):

1. int pow (int v, int e)
2. {
3.     int res = v;
4.  
5.     int loop = _bit_scan_reverse (e); // тута я ошибся в первом посте. "-1" не надо
6.     while (loop--)
7.     {
8.         res *= res;
9.         if (e & (1<<loop))
10.             res *= v;
11.     }
12.     return res;
13. }

твой код:

Код (Text):

1. int pow (int v, int e)
2. {
3.     int res = 1;
4.  
5.     int loop = _bit_scan_reverse (e);
6.     for(int i=0; i<=loop; i++) // твоя ошибка была тут. Надо "<=" а не "<"
7.     {
8.         if (e & (1<<i))  
9.             res *= v;
10.         v *= v;
11.     }
12.     return res;
13. }

Итак, теперь они 100% рабочие - я проверял

>> так нулевая степень и так за циклом
Ну я имел в виду if (!e) return 1;
Для int это несущественно. Но для 1024 битных чисел уже весомо.
К тому-ж bsr при нулевом аргументе - это UB )) Без if ну никак не обойтись.
На всяк случай напомню - bsr(1) = 0, т.к. bsr возвращает ПОЗИЦИЮ старшего значащего бита.

>> циклов по идее столько сколько бит в числе, так что о каком лишнем цикле речь
Неа ) Циклов можно смело делать на 1 меньше, т.к. старший значащий всегда == 1

Далее касательно циклов и лишних умножений:
давай развернем твой и мой циклы при возведении в скажем 5-ю степень

вычисляем X^5
v = X
e = 5

мой код:

Код (Text):

1. res = v
2. loop = _bit_scan_reverse (e) -> e == 2
3.  
4. // 1
5. res *= res;
6. if (e & (1<<1)) res *= v;
7. // 2
8. res *= res;
9. if (e & (1<<0)) res *= v;
10.  
11. return res;

твой:

Код (Text):

1. res = 1
2. loop = _bit_scan_reverse (e) -> e == 2
3.  
4. // 1
5. if (e & (1<<0)) res *= v;
6. v *= v;
7. // 2
8. if (e & (1<<1)) res *= v;
9. v *= v;
10. // 3
11. if (e & (1<<2)) res *= v;
12. v *= v; // вот это умножение тут вообще не выполняет полезной работы
13.  
14. return res;

ну и на цикл больше получается.
У меня за счет того, что 1*Х уже вычислено (и равно Х ), и СТАРШИЙ ЗНАЧАЩИЙ БИТ числа, позицию которого возвращает bsr, ВСЕГДА установлен - то я могу сэкономить один цикл + выброшу никому не нужное умножение на самом последнем цикле, когда его результат уже не будет использован.

т.е. мой код до свертки был примерно такой:

Код (Text):

1. res = 1
2. loop = _bit_scan_reverse (e) -> e == 2
3.  
4. if (e & (1<<2)) res *= v; // я свернул тут, т.к. этот бит всегда == 1
5. // 1
6. res *= res;
7. if (e & (1<<1)) res *= v;
8. // 2
9. res *= res;
10. if (e & (1<<0)) res *= v;
11.  
12. return res;

+ при обратном просмотре бит экспоненты "лишнее" возведение в квадрат (res *= res не нужно по алгоритму.

В принципе после доработки напильником и твою функцию можно ускорить:

Код (Text):

1. int pow (int v, int e)
2. {
3.     int res = 1;
4.  
5.     int loop = _bit_scan_reverse (e);
6.     for(int i=0; i<loop; i++) // твоя ошибка была тут. Надо "<=" а не "<"
7.     {
8.         if (e & (1<<i))  
9.             res *= v;
10.         v *= v;
11.     }
12.     return res * v;
13. }

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

PS. Sorry за столь длинное и наверное путанное объяснение.

 

Один сек. - не понятно как может быть на цикл больше, если на конце всегда 1.
в твоем примере: X^5 по правилу разложения числа на степень 2 -ки должно быть 3 итерации: X, X^2, X^4 и мой код это делает, а ты радуешься, что у тебя 2-е итерации. Теперь давай посмотрим на строчку в твоем коде "res*=v" - из -за нее ты сможешь получить , по идее, только степени вида 2^i & 2^i+1.
и еще почему ты e просматривашь со старших бит к младшим - это вообще будет работать некорректно, а то, что с 5 усе вышло так ето из-за симметрии binary(5)=101

 

UbIvItS
>> с 5 усе вышло так ето из-за симметрии
От жеж неверующий. )) число 5 было взято с потолка! Вместо него можно взять вообще любое другое число. 13 например.

>> это вообще будет работать некорректно
Хе хе... Вот уже больше чем полгода как оно работает корректно для чисел начиная от 4096 бит. И почему то ни разу не отработало некорректно. Что я делаю не так?

ВНИМАТЕЛЬНО прочитай алгоритм по шагам. Ну или возьми да прогони тестом оба алгоритма с рандомными значениями.

Или тебе надо доказательства что мой алгоритм работает? )))
Вот тебе расклад для возведения в 13-ю степень (она хоть не симметрична?)

Код (Text):

1. loop = 3;// 1101
2. res = v;
3. // в этой точке res = v = v^1
4.  
5. // 1
6. res *= res;
7. if (e & (1<<2)) res *= v;
8. // в этой точке res = (v ^ 2) * v = v^3
9.  
10. // 2
11. res *= res;
12. if (e & (1<<1)) res *= v;
13. // в этой точке res = (v ^ 3) ^ 2 = v^6
14.  
15. // 3
16. res *= res;
17. if (e & (1<<0)) res *= v;
18. // в этой точке res = ((v ^ 6) ^ 2) * v = v^13

давай теперь посмотрим как это же возведение осуществляет твой алгоритм:

Код (Text):

1. оригинальное значение v буду обозначать v' чтоб не путаццо
2. loop = 3
3. res = 1
4. // в этой точке res = 1 = v' ^ 0
5.  
6. // 1
7. if (e & (1<<0)) res *= v;
8. v *= v;
9. // в этой точке res = 1 * v' = v' ^ 1
10. // v = v' ^ 2
11.  
12. // 2
13. if (e & (1<<1)) res *= v;
14. v *= v;
15. // в этой точке res не изменился
16. // v = (v' ^2 ) ^ 2 = v' ^ 4
17.  
18. // 3
19. if (e & (1<<2)) res *= v;
20. v *= v;
21. // в этой точке res = v' * (v' ^ 4) = v' ^ 5
22. // v = (v' ^ 4) ^ 2 = v' ^ 8
23.  
24. // 4
25. if (e & (1<<3)) res *= v;
26. v *= v;
27. // в этой точке res = (v' ^ 5) * (v' ^ 8) = v' ^ 13
28. // v = (v' ^ 8) ^ 2 = v' ^ 16

Теперь понятнее стало?

>> Один сек. - не понятно как может быть на цикл больше, если на конце всегда 1
я цикл для бита, который "на конце всегда 1" просто не делаю. У меня циклов на 1 меньше чем кол-во бит в числе. Т.е. е=255 - 8 бит, а у меня 7 циклов. е = 1 - 1 бит, у меня 0 циклов.
И достижимо это именно из-за обратного порядка просмотра т.к. умножение для этого бита, что всегда == 1, можно заменить присвоением.

Ну елы палы, пройдись по алго пошагово дебаггером - сам поймешь всю его математику...