Модульное возведение в степень для большого мод в C++

Question

Это код, который я использую для расчета (n^p)%mod. К сожалению, он не работает при больших значениях mod (в моем случае mod = 10000000000ULL), когда я вызываю его из метода main(). Есть идеи; Зачем?

ull powMod(ull n, ull p, ull mod) { 
    ull ans = 1; 
    n = n%mod; 
    while(p) { 
     if(p%2 == 1) { 
      ans = (ans*n)%mod; 
     } 
     n = (n*n)%mod; 
     p /= 2; 
    } 
    return ans; 
} 

Здесь ull является ЬурейиМ для unsigned long long.

Answer 1

Да, вы можете сделать это на C++. Как указывали другие, вы не можете это сделать напрямую. Используя небольшое падение теории чисел, можно разложить проблему на две управляемые подпрограммы.

Сначала подумайте, что 10^10 = 2^10 * 5^10. Оба фактора являются взаимными, поэтому вы можете использовать Chinese remainder theorem, чтобы найти мощность по модулю 10^10, используя полномочия modulo 2^10 и modulo 5^10.

Обратите внимание, что в следующем коде были найдены магия значения u2 и u5 с помощью Extended Euclidean Algorithm. Вам не нужно программировать этот алгоритм самостоятельно, потому что эти значения являются константами. Я использую maxima и его функцию gcdex, чтобы вычислить их.

Вот модифицированная версия:

typedef unsigned long long ull; 

ull const M = 10000000000ull; 

ull pow_mod10_10(ull n, ull p) { 
    ull const m2 = 1024; // 2^10 
    ull const m5 = 9765625; // 5^10 
    ull const M2 = 9765625; // 5^10 = M/m2 
    ull const M5 = 1024; // 2^10 = M/m5 
    ull const u2 = 841;  // u2*M2 = 1 mod m2 
    ull const u5 = 1745224; // u5*M5 = 1 mod m5 

    ull b2 = 1; 
    ull b5 = 1; 
    ull n2 = n % m2; 
    ull n5 = n % m5; 

    while(p) { 
    if(p%2 == 1) { 
     b2 = (b2*n2)%m2; 
     b5 = (b5*n5)%m5; 
    } 
    n2 = (n2*n2)%m2; 
    n5 = (n5*n5)%m5; 
    p /= 2; 
    } 

    ull np = (((b2*u2)%M)*M2)%M; 
    np += (((b5*u5)%M)*M5)%M; 
    np %= M; 
    return np; 
} 

Answer 2

Возможно, одна из возможных проблем заключается в том, что когда вы делаете (a*b)%c, сама часть a*b может переполняться, что приводит к неправильному ответу. Один из способов обойти это использовать тождество,

(a*b)%c 

эквивалентно

(a%c * b%c)%c 

Это позволит предотвратить переполнение в промежуточных умножений, а также.

Answer 3

Кажется, что вы не можете избежать этого.

Если mod является 10000000000ULL, в (a*b)%c в вашей программе, как a и b меньше, чем мода, поэтому мы относимся к ним, как 9999999999ULL, a*b будет 99999999980000000001, но unsigned long long может выразить только 2^64-1=18446744073709551615 < 99999999980000000001 так что ваш метод будет переполнение.

Answer 4

Код линии

n = (n*n)%mod; 

выполняется повторно. Пока n меньше, чем мода, это может привести к оценке (mod-1) * (mod-1) в определенный момент времени.

На входе n может быть не так много, но указанная строка кода увеличивает n в цикле.