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

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

  1. diamond

    diamond New Member

    Публикаций:
    0
    Регистрация:
    21 май 2004
    Сообщения:
    507
    Адрес:
    Russia
    Код из #58 работает не во всех случаях - взял от балды dword [q] = dword [q+4] = 9ABCDEF0h, правильный остаток 0B5765398h, а код выдаёт 0B5865397h. Код из #57 не тестил.
     
  2. persicum

    persicum New Member

    Публикаций:
    0
    Регистрация:
    2 фев 2007
    Сообщения:
    947
    А может я отбираю себе код на конкурсной основе =))) Мне сначала нужно посмотреть портфолио, типа как это быстро и круто в сравнении с другими претендентами.
     
  3. Mikl___

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

    Публикаций:
    14
    Регистрация:
    25 июн 2008
    Сообщения:
    3.709
    2persicum
    Тогда облом -- я тоже отбираю работодателей на конкурсной основе =)))
    2diamond
    Не стоит. Набегает ошибка округления из-за того, что используется 10010010000000000h, а не на 1001000FF0FE0FD11h, в случае, если edx > 0EFFFFFFFh ошибка замечена и компенсирована -1, но, оказывается, я плохо продумал алго, пойду перечитаю Уоррена...
     
  4. persicum

    persicum New Member

    Публикаций:
    0
    Регистрация:
    2 фев 2007
    Сообщения:
    947
    Mikl___
    советую подойти к решению задачи более творчески. пост 27 использовал естественное разбиение числа на 64 бит как a*2^32 + b, но можно также представить его как a*2^20 + b. Тогда для нахождения искомого остатка достаточно разделить a на FFF, тут пригодятся и остаточек и частное от этого деления, если так нравится делить по Уоррену.
     
  5. Mikl___

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

    Публикаций:
    14
    Регистрация:
    25 июн 2008
    Сообщения:
    3.709
    persicum
     
  6. persicum

    persicum New Member

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

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

    Публикаций:
    14
    Регистрация:
    25 июн 2008
    Сообщения:
    3.709
    Код (Text):
    1.     mov edx,[q+4]
    2.     mov eax,[q]
    3. ;100000000h/FFF00001h=1,0002442000112558721159827002697=1001000FFh/100000000h=
    4. ;=1 + 1/1000h + 1/1000000h - 1/100000000h
    5.     mov ecx,edx
    6.     mov esi,eax
    7.     mov ebx,eax
    8.     sub ebx,edx
    9.     sbb ecx,0
    10.     shrd esi,edx,12
    11.     shr edx,12
    12.     add ebx,esi
    13.     adc ecx,edx
    14.     shrd esi,edx,12
    15.     shr edx,12
    16.     add ebx,esi
    17.     adc ecx,edx
    18. ;частное в ecx -- чтобы получить остаток умножаем ecx на 0FFF00001h и отнимаем от q
    19. ;ecx*0FFF00001=ecx*2^32 - ecx*2^20 + ecx, но так как остаток 32-разрядный
    20. ;rem = low(q) + ecx*2^20 - ecx
    21.     sub eax,ecx
    22.     shl ecx,20
    23.     add eax,ecx; остаток в eax
    diamond
    вроде бы, правильно считает остаток на всем 64-битном диаппазоне от 0 до 0FFFF.FFFF.FFFF.FFFFh, а не только, как div до 0FFF0.0000.FFFF.FFFFh, но если можно, протестируйте еще раз, я уже ни в чем не уверен...
    persicum
    Не очень сложно переделать и для C0000001h, D0000001h, F3000001h
     
  8. diamond

    diamond New Member

    Публикаций:
    0
    Регистрация:
    21 май 2004
    Сообщения:
    507
    Адрес:
    Russia
    Просто значения, взятые от балды, посчитались правильно; тогда взял edx:eax = FFF00001 * ABCDEFAB, прогнал - получил остаток FFF00001 вместо правильного 0. Добавил единицу, получил FFF00002. Добавил вместо единицы 100000h, код выдал остаток 1.
     
  9. persicum

    persicum New Member

    Публикаций:
    0
    Регистрация:
    2 фев 2007
    Сообщения:
    947
    Mikl___
    Код поста 67 просто замечательный, совсем другое дело!
    Я пока не совсем понимаю как он работает, меня от использования шифтов отворачивала необходимость вычитать, но тут идет прибавление .
    Однако, есть ряд недочетов.
    1) может возвращать недобитый остаток, например (p-1)(p-1) выдает p+1
    2) где гарантия, что невозможно переполнение в пределах допустимых значений множителей меньше p
    3) переполнение можно вызвать явно умножив -1 на -1, однако это не является недостатком процедуры, поскольку -1 лежит вне диапазона
    Это я говорю про процедуру дополренную вначале mul edx, поскольку мне нужно edx и eax вначале перемножить а уж потом брать остаток
     
  10. Mikl___

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

    Публикаций:
    14
    Регистрация:
    25 июн 2008
    Сообщения:
    3.709
    persicum
    Мне нужно сперва устранить ошибку, обнаруженную diamond, тогда и пункты 1-3, может быть, сами уйдут... А считать должно в пределах всего 64-битного диаппазона, это я с калькулятором проверял
     
  11. persicum

    persicum New Member

    Публикаций:
    0
    Регистрация:
    2 фев 2007
    Сообщения:
    947
    в подтыерждение опасения 2) нашел пример, когда допустимые входные величины вызывают потерю значимости, а именно (p-1)(p-FFFFF) должно преспокойненько давать FFFFF, а получается угадайте что? Лажа...
     
  12. Mikl___

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

    Публикаций:
    14
    Регистрация:
    25 июн 2008
    Сообщения:
    3.709
    Код (Text):
    1.     mov edx,dword ptr q+4
    2.     mov eax,dword ptr q
    3.     mov ecx,edx
    4.     mov esi,eax
    5.     mov ebx,eax
    6.     sub ebx,edx
    7.     sbb ecx,0
    8.     shrd esi,edx,12
    9.     shr edx,12
    10.     add ebx,esi
    11.     adc ecx,edx
    12.     shrd esi,edx,12
    13.     shr edx,12
    14.     add ebx,esi
    15.     adc ecx,edx
    16.     mov edi,ecx
    17. ;частное в edi
    18. ;чтобы получить остаток умножаем ecx на 0FFF00001h и отнимаем от q
    19. ;ecx*0FFF00001=ecx*2^32 - ecx*2^20 + ecx
    20.     sub eax,ecx
    21.     shl ecx,20
    22.     add eax,ecx
    23.     mov ecx,0FFF00000h
    24.     cmp ecx,eax
    25.     sbb ecx,ecx
    26.     sub edi,ecx; коррекция частного
    27.     and ecx,0FFF00001h
    28.     sub eax,ecx; остаток в eax
     
  13. diamond

    diamond New Member

    Публикаций:
    0
    Регистрация:
    21 май 2004
    Сообщения:
    507
    Адрес:
    Russia
    MS тоже любит накладывать заплатки... как показывает практика, дыр от этого меньше не становится.
    Отчётливо видна коррекция предыдущего решения. Снова взял edx:eax = FFF00001 * ABCDEFAB + 100000, снова получил остаток 1 вместо правильного 100000h.
     
  14. Mikl___

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

    Публикаций:
    14
    Регистрация:
    25 июн 2008
    Сообщения:
    3.709
    Код (Text):
    1.     mov edx,[q+4]
    2.     mov eax,[q]
    3.     test edx,edx
    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 ищем частное через умножение на
    13. ;1,0002442000112558721159827002697=1001000FFh/100000000h
    14. a1: mov ecx,edx
    15.     mov esi,eax
    16.     mov ebx,eax
    17.     sub ebx,edx
    18.     sbb ecx,0
    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.     add ebx,80000000h;<-- "округление", которого так не хватало
    28.     adc ecx,0;<------------┘
    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
    34. a2:...
    diamond
    "Чай, теперь твоя душенька довольна?" (© А.С.Пушкин)
     
  15. persicum

    persicum New Member

    Публикаций:
    0
    Регистрация:
    2 фев 2007
    Сообщения:
    947
    Mikl___
    Ну это просто звиндец! Клюв вынул - хвост увяз. Старые дыры ты залечил, не спорю, но попробуй теперь (p-1). А где гарантия, что и через десять постов оно будет абсолютно надежно? Теперь я понял, что значит деление без восстановления =))))))
     
  16. Mikl___

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

    Публикаций:
    14
    Регистрация:
    25 июн 2008
    Сообщения:
    3.709
    persicum
    Если p=0FFF00001, то
    1) (p-1)(p-1)/0FFF00001=0FFE0010000000000/0FFF00001 частное 0FFEFFFFF остаток 1
    2) (p-1)/0FFF00001=0FFF00000/0FFF00001 частное 0 остаток FFF00000
    3) (p-1)(p-FFFFF)/0FFF00001=FFD00201FFE00000/0FFF00001 частное 0FFE00001 остаток 0FFFFF
     
  17. persicum

    persicum New Member

    Публикаций:
    0
    Регистрация:
    2 фев 2007
    Сообщения:
    947
    Для примера 2) твой код выдает FFFFFFFF
     
  18. diamond

    diamond New Member

    Публикаций:
    0
    Регистрация:
    21 май 2004
    Сообщения:
    507
    Адрес:
    Russia
    Если цель всего этого - арифметика по модулю p и над числами производится много последовательных операций по модулю p, то можно воспользоваться трюком из статьи Монтгомери http://www2.itu.edu.tr/~orssi/dersler/cryptography/Montgomery.pdf : можно так кодировать вычеты по mod p, что операции сложения/вычитания устроены так же, как и при стандартном представлении, а умножение реализуется проще. При этом все константы и все входные данные нужно перекодировать (это некоторая перестановка на множестве вычетов), а все выходные данные - перекодировать обратно. Применительно к p = 0xFFF00001 три операции кодирования/умножения/декодирования выглядят так:
    Код (Text):
    1. start:
    2.     mov eax, 0x12345678 ; первый множитель в стандартном представлении
    3.     call    input2coded
    4.     push    eax
    5.     mov eax, 0x9ABCDEF0 ; второй множитель в стандартном представлении
    6.     call    input2coded
    7.     pop edx
    8.     call    coded_mul
    9.     call    coded2output
    10. ; eax = residue
    11.     ret
    12.  
    13. coded2output:
    14. ; eax -> eax
    15.     xor edx, edx
    16.     jmp coded_mul_doit
    17.  
    18. input2coded:
    19. ; eax -> eax
    20.     mov edx, 0x0FDFFF01
    21. coded_mul:
    22. ; eax*edx -> eax
    23.     mul edx
    24. coded_mul_doit:
    25.     mov esi, eax
    26.     shl eax, 20
    27.     add eax, esi
    28.     xor esi, esi
    29.     mov ebx, eax
    30.     shld    esi, eax, 20
    31.     shl ebx, 20
    32.     mov edi, eax
    33.     sub eax, ebx
    34.     sbb edi, esi
    35.     sub edx, edi
    36.     sbb eax, eax
    37.     and eax, 0xFFF00001
    38.     add eax, edx
    39.     ret
    Если производится несколько операций по mod p подряд, то между этими операциями кодировать/декодировать ничего не нужно.
     
  19. persicum

    persicum New Member

    Публикаций:
    0
    Регистрация:
    2 фев 2007
    Сообщения:
    947
    Про Монтгомери мысль была, но останавливают такие соображения
    1) лень переводить числа в другой формат
    2) грехх юзать универсальный метод с почти-Мерсеновским числом 2^32-2^20+1 !!!
    3) код #27 работает на умножениях без деления, как и Монтгомери
    4) murder , SSE2 ускоряет его в ТРИ раза !
    5) я не верю, что в Монтгомери можно складывать без компенсации. Если обычным порядком у меня произошло переполнение при сложении то я отниму p и все дела.
    А если у меня произойдет переполнение при сложении измененных данных, я разумеется не смогу взять просто младшие 32-бит, нужна такая же компинсация как и при умножении... Мне так кажется.
     
  20. diamond

    diamond New Member

    Публикаций:
    0
    Регистрация:
    21 май 2004
    Сообщения:
    507
    Адрес:
    Russia
    1) Код уже приведён.
    2) Пока что не видно работающих алгоритмов, использующих специфику числа, а то, что есть, заметно длиннее кода из #78.
    3) И использует цикл, а цикл - это ветвление, а ветвлений рекомендуется избегать.
    4) Это говорит скорее о тормознутости исходного кода :)
    5) Даже при обычном порядке если произошло переполнение при сложении, нужно это специально обрабатывать - нельзя просто взять младшие 32 бит (нужно ещё вычесть p). В методе Монтгомери всё точно так же.