Алгоритм деления без восстановления

Это будет примерно так.

Код (Text):

1. align 16
2. ab     rq 2
3. ij     rq 2
4. mask   dd 0,-1,0,-1
5. prime  dq prime,prime           ;<-Сюда подставишь числа
6. prime2 dq 2^32-Prime,2^32-Prime ;<-
7. ...
8. movaps   xmm0,[ab]
9. movaps   xmm1,xmm0
10. psrldq   xmm1,32
11. pmuludq  xmm0,xmm1
12.  
13. @loop:
14. movaps   xmm1,xmm0
15. andps    xmm1,[mask]
16. psrldq   xmm0,32
17. pmuludq  xmm0,[prime2]
18. paddq    xmm0,xmm1
19. pmovmskb eax,xmm0
20. test     eax,1111000011110000b
21. jnz @loop
22.  
23. movaps   xmm1,xmm0
24. pcmpgtd  xmm1,[prime]
25. andps    xmm1,[prime]
26. psubq    xmm0,xmm1
27.  
28. movaps   [ij],xmm0

Проверить не могу т.к. дома нет SSE2.

 

ух ты, обалденная крутизна! Мне бы биться головой об стенку целую неделю над этим... А нельзя ли пояснить только самую малость, самую творческую часть, так сказать. Как происходит выход из цикла и и команде movmsk. Понятно, что редуцируя сразу два числа, можно выйти на ситуацию когда одно число уже кончило а другое еще не отстрелялось, как узнать что старшие субполовинки обоих половин регистра уже обнулились?

 

Блин - бред написал.

Щас подумаю

 

Ну вот

Код (Text):

1. movaps   xmm0,[ab]
2. movaps   xmm1,xmm0
3. andps    xmm0,[mask]   ;выделяем младшие 32 бита у каждого из чисел
4. psrldq   xmm1,32       ;выделяем старшие 32 бита у каждого из чисел
5. pmuludq  xmm0,xmm1     ;умножаем младшую часть на старшую
6.  
7. @loop:
8. movaps   xmm1,xmm0
9. andps    xmm1,[mask]   ;выделяем младшие 32 бита у каждого из чисел
10. psrldq   xmm0,32       ;выделяем старшие 32 бита у каждого из чисел
11. pmuludq  xmm0,[prime2] ;умножаем старшую часть на 2^32-prime
12. paddq    xmm0,xmm1     ;прибавляем младшую часть
13. pcmpeqd  xmm1,xmm0     ;сравниваем старшую часть полученого числа с 0
14.                        ;если да - в старших 32 битах будет маска из единиц
15.                        ;иначе там будет ноль
16. movmskd  eax,xmm1      ;копируем старшие биты в eax - получается 2-битовая маска
17. cmp      al,11         ;если старшие 32 бита обоих чисел равны 0 - завершаем цикл
18. jnz @loop
19.  
20. movaps   xmm1,xmm0
21. pcmpgtd  xmm1,[prime]  ;если число больше prime
22. andps    xmm1,[prime]  ;вычитаем изнего prime
23. psubq    xmm0,xmm1     ;иначе вычитаем ноль
24.  
25. movaps   [ij],xmm0     ;сохраняем результат

Это не имеет значения. Если число уже "отстрелялось" дальнейшие действия не изменят его.

 

Это верно, я имел в виду чтобы не выйти раньше времени из цикла не дай бог, поскоку оба числа могут отстреляться не одновременно.

сам по себе pmul не боиться мусора в старших разрядах, если только не преследуется иных целей, специально очищать не нужно.

если равно то тоже не мешало бы вычесть и получить ноль. pcmpgtd отрабатывает такую ситуацию?

 

murder
А как считаете, если вместо pmul вставить (d shl 20) + prime - d, оно быстрее станет? С SSE2 ситуация может быть другой чем с shld eax,edx,20

 

Да, похоже, что очищать старшую часть перед умножением не обязательно.

Тогда видимо так

Код (Text):

1. movaps   xmm1,[prime]
2. pcmpgtd  xmm1,xmm0
3. pandn    xmm1,[prime]
4. psubq    xmm0,xmm1

То есть так?

Код (Text):

1. @loop:
2. movaps   xmm1,xmm0
3. pslldq   xmm0,20
4. paddq    xmm0,[prime]
5. psubq    xmm0,xmm1
6. xorps    xmm1,xmm1
7. pcmpeqd  xmm1,xmm0
8. movmskd  eax,xmm1
9. cmp      al,11
10. jnz @loop

 

Да не совсем так... Нужно тщательно выделять половинки, но двигать сначала 32 вправо а потом 20 влево большого смысла нет, видимо можно двинуть xmm0 на 12 вправо и очистить 20 младших бит через and. Но есть одна бяка, когда d=0, то будет прибавляться чистый prime а это плохо.

 

А что там, на SSE2 нету great or equal = not less?

 

Для сравнения больше или равно применяется pcmpgtd+pandn. В SSE присутствует инструкция cmpnltps, но она работает с вещественными числами. SSE2 - это просто расширенный MMX.

Я тоже об этом подумал. Тогда лучше оставить умножение.

Что-то Mikl___ сюда не заглядывает. Отправь ему сообщение в личку.

 

Угу, так и сделаю... Уж больно с шифтами много гимора. Надеюсь PMUL тормозить не будет плюс досрочный выход если много ведущих нулей надеюсь там реализован в железе.

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

еще часто используют остатки при вычислении хешей=ключей в базах данных.
Читаем: Number of tables per database 2^32 - 2^20 - 1 = 4293918719
Похожее число, но там прощще - d*2^20 + d, проблемы с нулем нет.

 

murder, похоже наш код работать не будет, так как pcmpgtd производит знаковое сравнение. У нее есть вариант для квадров или для даблов без знака? Второе было бы еще лучше...

Придумал! Инвертировать старшие биты перед сравнением... Елыпалы, опять гимморр...

 

murder, спасиьо за поддержку, числодробилку ака имитацию сдвигового регистра отладил, все работает хорошо, напрягает только необходимость дополнительных шагов по инверсии старших битов для команды pcmdgtd, поскоку она знаковая. Не ожидал такой глупости. В общем , прогоняю мельницу по очереди два раза, потом объединяю все в один xmm и добиваю остаток.

На заценку варианта Майкла___ уже потенции не осталось, еще нужно много работы чтобы интегрировать наш SSE2 код ко мне в прогу

 

Можно ещё так

Код (Text):

1. value    rq 1
2. cw       rw 1
3. rdivider dq 1/4293918721.0
4. divider  dq 4293918721.0 ;0FFF00001h
5. ...
6. fstcw   [cw]
7. or      [cw],111100000000b ;точность 64 бита, округление с отбрасыванием дробной части
8. fldcw   [cw]
9. ...
10. fild    [value]
11. fmul    [rdivider]
12. mov     eax,dword[value]
13. frndint
14. fmul    [divider]
15. fistp   [value]
16. sub     eax,dword[value]

 

Замена деления 64-битного числа на 0xFFF00001 умножениями для полного диапазона выглядит так: умножение на 0x1001000FF0FE0FD11 и сдвиг вправо на 96. Если считать, что частное влезает в 32 бита (а если делится число, меньшее 0xFFF00001^2, то частное заведомо влезает), и расписать по методу Карацубы, то один div можно заменить на следующий эквивалентный код:

Код (Text):

1.     mov edx, [q+4]
2.     mov eax, [q]
3. ; edx:eax mod 0xFFF00001 -> edi
4.     mov ebp, edx
5.     mov edi, eax
6.     mul [const1]
7.     mov ecx, eax
8.     mov esi, edx
9.     mov eax, edi
10.     add eax, ebp
11.     sbb ebx, ebx
12.     mul [const2]
13.     and ebx, [const2]
14.     add edx, ebx
15.     sub eax, ecx
16.     sbb edx, esi
17.     add eax, esi
18.     adc edx, 0
19.     mov ecx, eax
20.     mov esi, edx
21.     mov eax, ebp
22.     mul [const3]
23.     sub ecx, eax
24.     sbb esi, edx
25.     add eax, esi
26.     adc edx, 0
27.     add eax, edi
28.     adc edx, ebp
29. ; edx = частное
30.     sub edi, edx
31.     shl edx, 20
32.     add edi, edx
33. ; edi = остаток
34. ...
35. ; используемые константы
36. align 4
37. const1  dd  0FE0FD11h
38. const2  dd  1000FFh + 0FE0FD11h
39. const3  dd  1000FFh

 

Ну этого явно маловато будет =))) Процедура должна уметь редуцировать цисла до FFE0 0100 0000 0000 хотя бы, т.е. (p-1)^2. Иначе она практически бесполезна

А зачем козе баян? Процедура поста №27 работает как в АЛУшном варианте, так и практически лоб-в-лоб в SSE2 варианте. Требует только три укороченных умножения, в редчайших случаях четыре. И ващще, деление - это только один из методов получения остатка, самый наивный. Но есть и другие - например сумма цифр числа в десятичной системе есть остаток от деления на 9. Мне то частное совсем не нужно, а тока остаток и нужен. Кто работал с RSA тот знает что это такое, хотя 32-бит модуль будет коротковат для RSA, но хорош для других целей =)))

Mikl___
А нельзя ли провести замер скорости кода Вашего кода и из поста №27?