Алгоритм Шенхаге

Тема в разделе "WASM.A&O", создана пользователем Perec, 13 дек 2007.

  1. persicum

    persicum New Member

    Публикаций:
    0
    Регистрация:
    2 фев 2007
    Сообщения:
    947
    А мы тут в соседней ветке как раз балуемся числом FFF00001. Там есть и монтгомери и более простая редукция. Генератор 1557 - проверил и исправил. Но с этим числом удобно работать на Асме, языки могут загнуться, даже SSE2 загибается и местами требует инверсии старшего бита. Если миллиона мало, есть еще число C0000001, но оно уже слишком далеко от 2^32 и не так удобно. Хотя если набивать 24 или 31 бит то разницы не будет. А если 32 тогда в пользу первого.

    А для языков и трех-байтовой начинки удобно число 10*2^24+1
     
  2. persicum

    persicum New Member

    Публикаций:
    0
    Регистрация:
    2 фев 2007
    Сообщения:
    947
    а в плавучем Штрассене я сильно разочаровался, там с ростом длины можно все меньше и меньше полезных бит класть в double, сначала 24, потом 16, потом 8, а потом единичные вообще
     
  3. persicum

    persicum New Member

    Публикаций:
    0
    Регистрация:
    2 фев 2007
    Сообщения:
    947
    блин, кажись бред написал, числа с полиномами перепутал... Для чисел еще переносы нужно учитывать... Тогда отличия плавающей версий и целочисленной не в полне ясны...
     
  4. maxdiver

    maxdiver Max

    Публикаций:
    0
    Регистрация:
    18 июл 2006
    Сообщения:
    308
    Адрес:
    Саратов
    persicum
    Да ладно, вродь всё правильно :)
    Нам же не генератор нужен, а примитивный корень n-ой степени из единицы. Т.е. число, n-я степень которого равна 1, а все меньшие степени не дают единицу.
    Код (Text):
    1.     int cur = 1;
    2.     int mod = 7340033;
    3.     for (int i=0; i<(1<<20); ++i) {
    4.         cur = (cur * 5) % mod;
    5.         if (cur == 1 && i!=(1<<20)-1)
    6.             throw;
    7.     }
    8.     cout << cur << endl;
    Выводит 1, как и должно.

    Ссылочку дай что ли )
     
  5. persicum

    persicum New Member

    Публикаций:
    0
    Регистрация:
    2 фев 2007
    Сообщения:
    947
    1557 - это 2^20-th root, если хочется...

    Расскажи лучше, в целочисленном штрассене тоже чтоли нужно по нескольку бит класть в ячейку чтобы избежать переполнения? Тогда вчем кайф перед плавучей версией? Никак нельзя скажем для очень длинных чисел класть 24bit в dword?
     
  6. persicum

    persicum New Member

    Публикаций:
    0
    Регистрация:
    2 фев 2007
    Сообщения:
    947
    Профессионалы простым перебором не работают =)))
    По малой теореме Ферма a^(p-1)=1, т.е. a^(7*2^20)=(a^7)^(2^20)=1.

    Берем 3^7=2187 и убеждаемся, что энто хороший генератор! =)))
     
  7. maxdiver

    maxdiver Max

    Публикаций:
    0
    Регистрация:
    18 июл 2006
    Сообщения:
    308
    Адрес:
    Саратов
    persicum
    Чёт я не знаю, что у тебя за арифметика, вот например в online calculator посчитал:
    А почему именно три, а не другое число? :)
    Этот способ же тоже не гарантирует, что получится примитивный корень, так что перебирать тоже придётся.

    Не понял, о чём ты :)
    О том, что я int64 во всех произведениях использую? Ну казалось бы, если не брать всякие хитрые методы типа Монтгомери, правильный результат по модулю не посчитаешь, из-за переполнений получится фиг знает что.
     
  8. persicum

    persicum New Member

    Публикаций:
    0
    Регистрация:
    2 фев 2007
    Сообщения:
    947
    Не гарантирует, согласен. Если взять a=2 то получим, что 2^7 является не 20-тым, а только 19-м корнем. Но его можно использовать для 2^19-точечного FFT. Иногда так и делают. Хотя конечно лучше иметь самый полный генератор.

    Опти, у тебя там деление на каждом шагу (очевидно % просто DIV-ует чтобы взять остаток). Откуда ты такой модуль откапал - 7*2^20+1? Попробуй сначала чтоли с Ферматным 65537. Там остаток - просто a - b

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

    Например, у меня модуль 5*2^25+1. Я могу класть в ячейку три байта запросто.
    Но если я попробую перемножить скажем FFFFFF на ЕЕЕЕЕЕ используя вектора (0, FFFFFF), (0,EEEEEE), то я получу в результате ifft(fft*fft) такой вектор - (0, 2F92C61) - только кусочек произведения =(((
     
  9. maxdiver

    maxdiver Max

    Публикаций:
    0
    Регистрация:
    18 июл 2006
    Сообщения:
    308
    Адрес:
    Саратов
    persicum
    Ну 7*2^20+1 простое, а числа поменьше вида k*2^20+1 не простые. А почему такое большое - ну например при длинах входных чисел 10^5 элементы результата умножения полиномов могут достигать (9^2)*(10^5), т.е. примерно 8*(10^6). Уже больше модуля. Поэтому даже для входных данных такого размера модуль нужен ещё больше...
    А про деление на каждом шагу я сразу и сказал - если его заменить чем-нибудь, хотя бы умножением, будет всё намного быстрей.

    Ну если модуль брать достаточно большим, то не будет опасности. Но как-то получается, что этот модуль даже для умножения десятичных чисел должен быть не меньше N*(9^2), а для полиномов с более-менее большими коэффициентами вообще беда получается... Да, по скорости не очень хорошо выходит, но зато можно быть уверенным в точности ответа.

    Ясно, модуль-то должен вмещать в себя коэффициенты произведения двух полиномов, здесь же далеко не так.
     
  10. persicum

    persicum New Member

    Публикаций:
    0
    Регистрация:
    2 фев 2007
    Сообщения:
    947
    Тот кто косит под Эйнштейна мог бы написать процедуру редукции регистровой пары mod 7*2^20+1 без деления и даже без Монтгомери