возможности переполнения в модульном возведении в квадрат при возведении в квадрат

Question

возможности переполнения в модульном возведении в квадрат при возведении в квадрат

Я ищу для реализации маленькая теорема Ферма для первичного тестирования. Вот код, который я написал:

lld expo(lld n, lld p) //2^p mod n
{
if(p==0)
return 1;
lld exp=expo(n,p/2);
if(p%2==0)
return (exp*exp)%n;
else
return (((exp*exp)%n)*2)%n;
}

bool ifPseudoPrime(lld n)
{
if(expo(n,n)==2)
return true;
else
return false;
}

НОТА: Я взял значение a(<=n-1) как 2,

Теперь число n может доходить до 10^18, Это означает, что переменная exp может достигать значений рядом 10^18, Что дополнительно подразумевает, что выражение (exp*exp) может достигать так высоко, как 10^36 следовательно вызывая переполнение. Как мне этого избежать.

Я проверил это, и он работал нормально до 10^9, я использую C ++

2

algorithm c++exponentiation primes

Решение

Другие решения

Вы можете выполнять свои умножения в несколько этапов. Например, скажем, вы хотите вычислить X * Y mod n. Возьмите X и Y и запишите их как X = 10 ^ 9 * X_1 + X_0, Y = 10 ^ 9 * Y_1 + Y_0. Затем вычислите все четыре произведения X_i * Y_j mod n и, наконец, вычислите X = 10 ^ 18 * (X_1 * Y_1 mod n) + 10 ^ 9 * (X_0 * Y_1 + X_1 * Y_0 mod n) + X_0 * Y_0. Обратите внимание, что в этом случае вы работаете с числами, равными половине максимально допустимого.

Если разделения на две части недостаточно (я подозреваю, что это так), разделите на три части, используя одну и ту же схему. Разделение на три должно работать.

Более простой подход — просто умножить школьный путь. Он соответствует предыдущему подходу, но записывает одно число на столько частей, сколько у него есть цифр.

Удачи!

0

Источник

Accepted Answer

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

Если модуль m настолько велика, что 2*(m-1) переполнение, вещи становятся действительно суетливыми, но если 2*(m-1) не переполняет, это терпимо.

Предположим, у вас есть и используется 64-битный целочисленный тип без знака.

Вы можете рассчитать модульный продукт, разделив коэффициенты на 32 и 16 битов, а затем разделить продукт на

a = a1 + (a2 << 32)    // 0 <= a1, a2 < (1 << 32)
b = b1 + (b2 << 32)    // 0 <= b1, b2 < (1 << 32)
a*b = a1*b1 + (a1*b2 << 32) + (a2*b1 << 32) + (a2*b2 << 64)

Вычислять a*b (mod m) с m <= (1 << 63), уменьшите каждый из четырех произведений по модулю m,

p1 = (a1*b1) % m;
p2 = (a1*b2) % m;
p3 = (a2*b1) % m;
p4 = (a2*b2) % m;

и самый простой способ включить сдвиги

for(i = 0; i < 32; ++i) {
p2 *= 2;
if (p2 >= m) p2 -= m;
}

то же самое для p3 и с 64 итерациями для p4, затем

s = p1+p2;
if (s >= m) s -= m;
s += p3;
if (s >= m) s -= m;
s += p4;
if (s >= m) s -= m;
return s;

Этот способ не очень быстрый, но для нескольких умножений, необходимых здесь, он может быть достаточно быстрым. Небольшое ускорение должно быть достигнуто за счет уменьшения количества смен; сначала рассчитать (p4 << 32) % m,

for(i = 0; i < 32; ++i) {
p4 *= 2;
if (p4 >= m) p4 -= m;
}

тогда все p2, p3 и текущее значение p4 нужно умножить на 2³² по модулю m,

p4 += p3;
if (p4 >= m) p4 -= m;
p4 += p2;
if (p4 >= m) p4 -= m;
for(i = 0; i < 32; ++i) {
p4 *= 2;
if (p4 >= m) p4 -= m;
}
s = p4+p1;
if (s >= m) s -= m;
return s;

8