1. Если вы только начинаете программировать на ассемблере и не знаете с чего начать, тогда попробуйте среду разработки ASM Visual IDE
    (c) на правах рекламы
    Скрыть объявление

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

Тема в разделе "WASM.A&O", создана пользователем Mikl___, 16 янв 2009.

  1. murder

    murder Member

    Публикаций:
    0
    Регистрация:
    3 июн 2007
    Сообщения:
    630
    Это будет примерно так.
    Код (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.
     
  2. persicum

    persicum New Member

    Публикаций:
    0
    Регистрация:
    2 фев 2007
    Сообщения:
    947
    ух ты, обалденная крутизна! Мне бы биться головой об стенку целую неделю над этим... А нельзя ли пояснить только самую малость, самую творческую часть, так сказать. Как происходит выход из цикла и и команде movmsk. Понятно, что редуцируя сразу два числа, можно выйти на ситуацию когда одно число уже кончило а другое еще не отстрелялось, как узнать что старшие субполовинки обоих половин регистра уже обнулились?
     
  3. murder

    murder Member

    Публикаций:
    0
    Регистрация:
    3 июн 2007
    Сообщения:
    630
    Блин - бред написал.

    Щас подумаю
     
  4. murder

    murder Member

    Публикаций:
    0
    Регистрация:
    3 июн 2007
    Сообщения:
    630
    Ну вот
    Код (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     ;сохраняем результат
    Это не имеет значения. Если число уже "отстрелялось" дальнейшие действия не изменят его.
     
  5. persicum

    persicum New Member

    Публикаций:
    0
    Регистрация:
    2 фев 2007
    Сообщения:
    947
    Это верно, я имел в виду чтобы не выйти раньше времени из цикла не дай бог, поскоку оба числа могут отстреляться не одновременно.

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

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

    persicum New Member

    Публикаций:
    0
    Регистрация:
    2 фев 2007
    Сообщения:
    947
    murder
    А как считаете, если вместо pmul вставить (d shl 20) + prime - d, оно быстрее станет? С SSE2 ситуация может быть другой чем с shld eax,edx,20
     
  7. murder

    murder Member

    Публикаций:
    0
    Регистрация:
    3 июн 2007
    Сообщения:
    630
    Да, похоже, что очищать старшую часть перед умножением не обязательно.

    Тогда видимо так
    Код (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
     
  8. persicum

    persicum New Member

    Публикаций:
    0
    Регистрация:
    2 фев 2007
    Сообщения:
    947
    Да не совсем так... Нужно тщательно выделять половинки, но двигать сначала 32 вправо а потом 20 влево большого смысла нет, видимо можно двинуть xmm0 на 12 вправо и очистить 20 младших бит через and. Но есть одна бяка, когда d=0, то будет прибавляться чистый prime а это плохо.
     
  9. persicum

    persicum New Member

    Публикаций:
    0
    Регистрация:
    2 фев 2007
    Сообщения:
    947
    А что там, на SSE2 нету great or equal = not less?
     
  10. murder

    murder Member

    Публикаций:
    0
    Регистрация:
    3 июн 2007
    Сообщения:
    630
    Для сравнения больше или равно применяется pcmpgtd+pandn. В SSE присутствует инструкция cmpnltps, но она работает с вещественными числами. SSE2 - это просто расширенный MMX.

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

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

    persicum New Member

    Публикаций:
    0
    Регистрация:
    2 фев 2007
    Сообщения:
    947
    Угу, так и сделаю... Уж больно с шифтами много гимора. Надеюсь PMUL тормозить не будет плюс досрочный выход если много ведущих нулей надеюсь там реализован в железе.

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

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

    persicum New Member

    Публикаций:
    0
    Регистрация:
    2 фев 2007
    Сообщения:
    947
    murder, похоже наш код работать не будет, так как pcmpgtd производит знаковое сравнение. У нее есть вариант для квадров или для даблов без знака? Второе было бы еще лучше...

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

    Mikl___ Супермодератор Команда форума

    Публикаций:
    14
    Регистрация:
    25 июн 2008
    Сообщения:
    2.993
    Извините за молчание, у меня вот такой вариант
    Код (Text):
    1. .686P
    2. .model flat
    3. include windows.inc
    4. include masm32.inc
    5. includelib user32.lib
    6. includelib masm32.lib
    7. extern _imp__wsprintfA:dword
    8. extern _imp__MessageBoxA@16:dword
    9. .code
    10. start:  mov edx,dword ptr q+4
    11.     mov eax,dword ptr q
    12.     test edx,edx
    13.     jnz a1
    14. ;если edx = 0 и eax < 0FFF00001  тогда остаток = eax
    15. ;если edx = 0 и eax >= 0FFF00001 тогда остаток = eax - 0FFF00001
    16.     sub eax,0FFF00001h
    17.     setnc dl;<-- это только для сравнения с div вообще-то не нужно
    18.     push edx;<-- это только для сравнения с div вообще-то не нужно
    19.     sbb esi,esi
    20.     and esi,0FFF00001h
    21.     add eax,esi
    22.     jmp a2
    23. ;если edx > 0 ищем частное через умножение на 1001h и сдвиг вправо на 44 бита
    24. ;100000000h/FFF00001h=1,0002442000112558721159827002697=1001h/1000h
    25. a1: mov ecx,edx
    26.     shld ecx,eax,12
    27.     add ecx,edx
    28.     shr ecx,12
    29.     push ecx;<-- это только для сравнения с div вообще-то не нужно
    30. ;частное в ecx
    31. ;чтобы получить остаток умножаем ecx на 0FFF00001h и отнимаем от q
    32. ;ecx*0FFF00001=ecx*2^32 - ecx*2^20 + ecx
    33.     sub eax,ecx
    34.     shl ecx,20
    35.     add eax,ecx
    36. a2: push eax
    37.     mov edx,dword ptr q+4
    38.     mov eax,dword ptr q
    39.         mov ecx,0FFF00001h
    40.     div ecx
    41.     push eax
    42.     push edx
    43.         push offset format
    44.         push offset buffer 
    45.     call _imp__wsprintfA
    46.     add esp,4*6;корректируем стек
    47.     push MB_OK
    48.     push offset MsgCaption
    49.     push offset buffer 
    50.     push 0
    51.     call _imp__MessageBoxA@16
    52.     ret
    53. MsgCaption      db "persicum",0
    54. buffer db 20 dup (0)
    55. q dq ?;<--подставляем любое число из диаппазона от 0 до 0F.FFFF.FFFF.FFFFh для проверки
    56. format db 'остаток через div  %08Xh',0Ah,'частное через div %08Xh',0Ah
    57. db '------------------------',0Ah
    58. db 'остаток через сдвиги %08Xh',0Ah,'частное через сдвиги %08Xh',0
    59. end start
    К сожалению, без ошибок работает, если q в диапазоне от 0 до 4503599627370495. Для распространения на весь 64-битный диапазон, может быть, murder, leo или diamond помогут...
     
  14. persicum

    persicum New Member

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

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

    murder Member

    Публикаций:
    0
    Регистрация:
    3 июн 2007
    Сообщения:
    630
    Можно ещё так
    Код (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]
     
  16. diamond

    diamond New Member

    Публикаций:
    0
    Регистрация:
    21 май 2004
    Сообщения:
    507
    Адрес:
    Russia
    Замена деления 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
     
  17. Mikl___

    Mikl___ Супермодератор Команда форума

    Публикаций:
    14
    Регистрация:
    25 июн 2008
    Сообщения:
    2.993
    persicum
    Еще один вариант -- считает правильно в диаппазоне от 0 до 0EFFF.FFFF.FFFF.FFFFh, обходимся без умножения одними сложениями и вычитаниями (на этот вариант потенции должно хватить :) )
    Код (Text):
    1. start:  mov edx,[q+4]
    2.     mov eax,[q]
    3.     test edx,edx
    4.     jnz a1
    5. ;если edx = 0 и eax < 0xFFF00001,  тогда остаток = eax
    6. ;если edx = 0 и eax >= 0xFFF00001, тогда остаток = eax - 0xFFF00001
    7.     sub eax,0FFF00001h
    8.     setnc dl;<-- это только для сравнения с div, вообще-то не нужно
    9.     push edx;<-- это только для сравнения с div, вообще-то не нужно
    10.     sbb esi,esi
    11.     and esi,0FFF00001h
    12.     add eax,esi
    13.     jmp a2
    14. ;если edx > 0 ищем частное через умножение на 1001001h
    15. ;100000000h/FFF00001h=1,0002442000112558721159827002697=1001001h/1000000h
    16. a1: mov ecx,edx
    17.     mov esi,eax
    18.     mov ebx,eax
    19.     shrd esi,edx,12
    20.     shr edx,12
    21.     add ebx,esi
    22.     adc ecx,edx
    23.     shrd esi,edx,12
    24.     shr edx,12
    25.     add ebx,esi
    26.     adc ecx,edx
    27.     push ecx
    28. ;частное в ecx
    29. ;чтобы получить остаток умножаем ecx на 0FFF00001h и отнимаем от q
    30. ;ecx*0FFF00001=ecx*2^32 - ecx*2^20 + ecx
    31.     sub eax,ecx
    32.     shl ecx,20
    33.     add eax,ecx
    34. a2: push eax;<-- это только для сравнения с div, вообще-то не нужно
    35.     mov edx,[q+4]
    36.     mov eax,[q]
    37.         mov ecx,0FFF00001h
    38.     div ecx
    39.     push eax
    40.     push edx
    41.         push offset format
    42.         push offset buffer 
    43.     call _imp__wsprintfA
    44.     add esp,4*6;корректируем стек
    45.     push MB_OK
    46.     push offset MsgCaption
    47.     push offset buffer 
    48.     push 0
    49.     call _imp__MessageBoxA@16
    50. . . .
    51. MsgCaption      db "persicum",0
    52. buffer db 20 dup (0)
    53. q dq ?;<--подставляем любое число из диаппазона от 0 до 0EFFFFFFFFFFFFFFFh для проверки
    54. format db 'остаток через div  %08Xh',0Ah,'частное через div %08Xh',0Ah
    55. db '------------------------',0Ah
    56. db 'остаток через сдвиги %08Xh',0Ah,'частное через сдвиги %08Xh',0
     
  18. Mikl___

    Mikl___ Супермодератор Команда форума

    Публикаций:
    14
    Регистрация:
    25 июн 2008
    Сообщения:
    2.993
    persicum
    Третий и надеюсь окончательный вариант -- считаем в диаппазоне от 0 до 0FFF00000FFFFFFFF,
    если содержимое edx больше или равно 0FFF00001h вычисление остатка через деление на FFF00001h должно прерывать программу через переполнение
    Код (Text):
    1.     mov edx,[q+4]
    2.     mov eax,[q]
    3.     test edx,edx;<-- этот кусок до метки a1 можно убрать, просто вычислений меньше
    4.     jnz a1; если edx равен нулю
    5. ;если edx = 0 и eax < 0FFF00001  тогда остаток = eax
    6. ;если edx = 0 и eax >= 0FFF00001 тогда остаток = eax - 0FFF00001
    7.     sub eax,0FFF00001h 
    8.     sbb esi,esi
    9.     and esi,0FFF00001h
    10.     add eax,esi
    11.     jmp a2; в eax остаток
    12. ;если edx > 0 ищем частное через умножение на 1001h и сдвиг вправо на 44 бита
    13. ;100000000h/FFF00001h=1,0002442000112558721159827002697=1001001h/1000000h
    14. a1: mov ecx,edx
    15.     cmp edx,0F0000000h
    16.     sbb edi,edi; придется компенсировать -1
    17.     not edi; так как умножаем на 1,001001h, а не на 1,001000FF0FE0FD11h
    18.     mov esi,eax
    19.     mov ebx,eax
    20.     shrd esi,edx,12
    21.     shr edx,12
    22.     add ebx,esi
    23.     adc ecx,edx
    24.     shrd esi,edx,12
    25.     shr edx,12
    26.     add ebx,esi
    27.     adc ecx,edx
    28.     add ecx,edi; компенсация равна 0 или -1, если edx больше или равно 0F0000000h
    29. ;частное в ecx, чтобы получить остаток умножаем ecx на 0FFF00001h и отнимаем от q
    30. ;ecx*0FFF00001=ecx*2^32 - ecx*2^20 + ecx
    31.     sub eax,ecx
    32.     shl ecx,20
    33.     add eax,ecx; в eax остаток
    34. a2:
     
  19. persicum

    persicum New Member

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

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

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

    Mikl___ Супермодератор Команда форума

    Публикаций:
    14
    Регистрация:
    25 июн 2008
    Сообщения:
    2.993
    сейчас диапазон от 0 до 0FFF00000FFFFFFFFh, так же как и через div, иначе div даст выход по int0, а замер скорости, все-таки корректнее делать вам самому, я же не знаю конечной цели и оборудование, на котром будет работать находится на вашей стороне