Python: ускорить pow(base,exp,mod) для фиксированного опыта и мода или с векторизацией

Хорошие новости, плохие новости. Хорошая новость заключается в том, что когда модуль m фиксирован, есть способы ускорить вычисление a*b % m. Найдите редукцию Барретта и редукцию Монтгомери. Они работают по-разному, предварительно вычисляя константы, связанные с m, так что % m можно вычислить с помощью умножения и сдвига без необходимости деления.

Плохая новость: чтобы найти остаток, оба способа требуют (в дополнение к более дешевым операциям) двух умножений. Таким образом, они не платят в целом, если только умножение не намного дешевле, чем деление.

По этой причине они, как правило, медленнее, если только модуль не является действительно большим - от 100 до 256 бит все еще мало по современным стандартам, всего в несколько раз шире, чем собственные 64-битные машинные целые числа. Такие вещи, как быстрое умножение на основе БПФ, требуют гораздо больших целых чисел, прежде чем они окупятся.

Встроенный модульный pow CPython уже использует двоичную схему, аналогичную тому, что вы закодировали в Python, но более причудливую (если показатель степени достаточно велик, встроенный pow рассматривает его как находящийся в базе 32, потребляя 5 битов показателя степени за итерацию цикла).

При быстром внедрении редукции Монтгомери в Python и замене модульных умножений в вашем коде написанием Монтгомери modular_pow() не стало быстрее, чем встроенное, до того, как модуль вырос до десятков тысяч бит. Для входных данных около 256 бит это было примерно в 3 раза медленнее.

Это смешанная ситуация: код Python не использует трюки с основанием 32, которые могут дать существенные преимущества. Но для достаточно больших входных данных CPython использует более быстрое, чем наивное, умножение Карацубы, от которого может выиграть правописание Монтгомери без деления (целое деление CPython не имеет приемов ускорения независимо от входных размеров, а встроенный модульный pow CPython всегда использует деление для найти остатки).

Итак, краткий курс: я не знаю ничего очевидного, что вы можете сделать в CPython для ускорения одного экземпляра pow(a, b, c). Возможно, в какой-то криптографической библиотеке с кодом C есть что-то подходящее, но я не знаю.

Но другая хорошая новость заключается в том, что ваша проблема смущающе параллельна. Если у вас есть N процессоров, вы можете дать каждому из них 100000000/N ваших входов, и все они могут работать на полной скорости параллельно. Это дало бы ускорение примерно в N раз.

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

СЛЕДОВАТЬ ЗА

Справочник по прикладной криптографии (HAC), глава 14, в основном разъясняет современное состояние алгоритмов гонзо-модульного возведения в степень.

Глядя на код, GMP уже реализует все свои трюки. Это включает в себя вещи, которые я упомянул (сокращение Монтгомери и использование основания степени 2 выше 2, чтобы пережевывать больше битов экспоненты за итерацию цикла). И другие, о которых я не упомянул (например, GMP имеет специальную внутреннюю процедуру для возведения в квадрат, которая экономит циклы по общему произведению, возможно, неравных целых чисел). В целом, это небольшая гора кода реализации.

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

Поэтому, если вам нужно добиться этого, использование GMP, вероятно, будет самым быстрым способом. Как уже отмечалось, многопроцессорность — это очевидный способ получить теоретическое N-кратное ускорение с N процессорами, но, как также отмечалось, вы ничего не сказали о контексте (откуда поступают эти входные данные или что вам нужно делать с выходными данными). Так что невозможно предположить, может ли это окупиться для вас. Чем больше межпроцессного взаимодействия вам нужно, тем больше это вредит потенциальному ускорению многопроцессорной обработки.

Примечание: вы делаете именно то, что делают, например, криптосистемы с открытым ключом RSA, хотя они обычно используют большие целые числа. То есть ваша база — это их сообщение, а открытый (или закрытый) ключ RSA состоит из фиксированного показателя степени и фиксированного модуля. Только база (сообщение или зашифрованные биты) различается в зависимости от экземпляра шифрования/дешифрования. Для данного ключа показатель степени и модуль всегда одинаковы.

Многие математики мирового класса изучали эту проблему, а хакеры мирового класса закодировали алгоритмы для максимальной скорости. Вот почему вы должны отказаться от надежды на то, что есть более быстрый способ, который HAC просто забыл упомянуть ;-)

Спекулятивный

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

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

Однако не так много для простого модуля. Единственный весьма потенциально ценный трюк заключается в том, что для p простых и a не кратных p

pow(a, b, p) == pow(a, b % (p-1), p)

Это может сэкономить неограниченное время, если b может быть намного больше, чем p. Это работает, потому что по малой теореме Ферма

pow(a, p-1, p) == 1

для p простых и a не кратных p. Например,

>>> p
2347
>>> assert all(p % i != 0 for i in range(2, p))  # that is, p is prime
>>> pow(345, 1000000000000000000000000000000, p)
301
>>> 1000000000000000000000000000000 % (p-1)
1198
>>> pow(345, 1198, p) # same thing, but much faster
301

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

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

Глядя на вики-страницу. Не похоже, что ваша реализация правильная. Перемещение этих двух операторов из else значительно улучшило производительность.

Это то, что у меня есть из Википедии

def modular_pow(base, exponent, modulus):
    if modulus == 1:
        return 0
    else:
        result = 1
        base = base % modulus
        while exponent > 0:
            if exponent % 2 == 1:
                result = (result * base) % modulus
            exponent = exponent >> 1
            base = (base * base) % modulus
        return result

Вывод:

print(modular_pow(4, 13, 497))

445

Вы можете использовать оконный метод NAF для предварительного вычисления a^2, a^3,...,a^(2^w-1). Теперь вместо n продуктов и возведения в квадрат у вас есть n/w раундов продуктов. Например, в 256-битном modexp при w=4 мы делаем 6 предварительных вычислений. Но вместо 256 квадратов и произведений мы имеем 256/4=64 произведения. При стоимости 14 предварительных вычислений это серьезная экономия. Теперь 4 бита — это 2^4=16 возможных значений. Но NAF представляет их в диапазоне -w+1..w-1. Обратите внимание, что для отрицательных показателей требуется модульная инверсия a^(-1). Таким образом, простое использование кодирования для масштабирования положительных значений является более оптимальным, чем дополнительное умножение или необходимость вычисления модульной инверсии. Обратите внимание, что a^0 и a^1 не требуют вычислений.

Некоторые предварительные вычисления могут быть оптимизированы на основе показателя степени в форме NAF, поскольку вы заранее знаете, какие именно значения потребуются.

Значение w должно быть скорректировано на основе возведения в степень, но размера. Оптимальное значение можно рассчитать на основе отношения n/w к 2^w-1 или определить эмпирически.

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