Exponencialización modular para números altos en C++

Question

Así que he estado trabajando recientemente en una implementación de la prueba de primalidad Miller-Rabin. Lo estoy limitando al alcance de todos los números de 32 bits, porque este es un proyecto de diversión justa que estoy haciendo para familiarizarme con C++, y no quiero tener que trabajar con nada de 64 bits para Un rato. Una ventaja adicional es que el algoritmo es determinista para todos los números de 32 bits, por lo que puedo aumentar significativamente la eficiencia porque sé exactamente qué testigos probar.

Por lo tanto, para números bajos, el algoritmo funciona excepcionalmente bien. Sin embargo, parte del proceso se basa en la exponenciación modular, es decir (num^pow)% mod. así, por ejemplo,

3^2 % 5 = 
9 % 5 = 
4 

Este es el código que he estado utilizando para este exponenciación modular:

unsigned mod_pow(unsigned num, unsigned pow, unsigned mod) 
{ 
    unsigned test; 
    for(test = 1; pow; pow >>= 1) 
    { 
     if (pow & 1) 
      test = (test * num) % mod; 
     num = (num * num) % mod; 
    } 

    return test; 

} 

Como es posible que ya haya adivinado, los problemas surgen cuando los argumentos son excepcionalmente grandes números. Por ejemplo, si quiero probar el número 673109 de primalidad, yo en un momento dado tiene que encontrar:

(2^168277)% 673109

ahora 2^168277 es un número excepcionalmente grande, y en alguna parte en el proceso, desborda la prueba, lo que da como resultado una evaluación incorrecta.

en el reverso, argumentos como

4000111222^3% 1608

también evalúan de forma incorrecta, por la misma razón.

¿Alguien tiene sugerencias para la exponenciación modular de una manera que pueda prevenir este desbordamiento y/o manipularlo para producir el resultado correcto? (La forma en que lo veo, desbordamiento es más que otra forma de módulo, es decir num% (UINT_MAX + 1))

Answer 1

Exponentiation by squaring todavía "funciona" para la exponenciación de módulo. Su problema no es que 2^168277 es un número excepcionalmente grande, es que uno de sus resultados intermedios es un número bastante grande (más grande que 2^32), porque 673109 es más grande que 2^16.

Así que creo que lo siguiente será suficiente. Es posible que haya omitido un detalle, pero la idea básica funciona, y así es como el código criptográfico "real" podría hacer grandes mod-exponenciaciones (aunque no con números de 32 y 64 bits, sino con bignums que nunca deben ser más grandes que 2 * registro (módulo)):

• Comience con exponenciación cuadrando, como lo ha hecho usted.
• Realiza la cuadratura real en un entero sin signo de 64 bits.
• Reduce el módulo 673109 en cada paso para volver dentro del rango de 32 bits, como lo haces.

Obviamente, es un poco incómodo si su implementación en C++ no tiene un entero de 64 bits, aunque siempre puede falsificar uno.

Hay un ejemplo en la diapositiva 22 aquí: http://www.cs.princeton.edu/courses/archive/spr05/cos126/lectures/22.pdf, aunque usa números muy pequeños (menos de 2^16), por lo que puede no ilustrar nada que usted no sepa.

Su otro ejemplo, 4000111222^3 % 1608 funcionaría en su código actual si solo reduce 4000111222 módulo 1608 antes de comenzar. 1608 es lo suficientemente pequeño como para poder multiplicar con seguridad dos números mod-1608 en un int de 32 bits.

Answer 2

dos cosas:

• ¿Está utilizando el tipo de datos adecuado? En otras palabras, ¿UINT_MAX le permite tener 673109 como argumento?

No, no, ya que en un momento dado tiene su código no funciona porque en un momento tiene num = 2^16 y la num = ... provoca desbordamiento. Use un tipo de datos más grande para mantener este valor intermedio.

• ¿Qué tal llevar módulo en cada posible ágape desbordamiento tales como:

  test = ((test % mod) * (num % mod)) % mod;

Editar:

unsigned mod_pow(unsigned num, unsigned pow, unsigned mod) 
{ 
    unsigned long long test; 
    unsigned long long n = num; 
    for(test = 1; pow; pow >>= 1) 
    { 
     if (pow & 1) 
      test = ((test % mod) * (n % mod)) % mod; 
     n = ((n % mod) * (n % mod)) % mod; 
    } 

    return test; /* note this is potentially lossy */ 
} 

int main(int argc, char* argv[]) 
{ 

    /* (2^168277) % 673109 */ 
    printf("%u\n", mod_pow(2, 168277, 673109)); 
    return 0; 
} 

Answer 3

Usted podría utilizar siguiente identidad:

(a * b) (mod m) === (a (mod m)) * (b (mod m)) (mod m)

Trate de usarlo de manera directa y mejorar incrementalmente.

if (pow & 1) 
     test = ((test % mod) * (num % mod)) % mod; 
    num = ((num % mod) * (num % mod)) % mod; 

Answer 4

Escribí algo para esto recientemente para RSA en C++, aunque un poco desordenado.

#include "BigInteger.h" 
#include <iostream> 
#include <sstream> 
#include <stack> 

BigInteger::BigInteger() { 
    digits.push_back(0); 
    negative = false; 
} 

BigInteger::~BigInteger() { 
} 

void BigInteger::addWithoutSign(BigInteger& c, const BigInteger& a, const BigInteger& b) { 
    int sum_n_carry = 0; 
    int n = (int)a.digits.size(); 
    if (n < (int)b.digits.size()) { 
     n = b.digits.size(); 
    } 
    c.digits.resize(n); 
    for (int i = 0; i < n; ++i) { 
     unsigned short a_digit = 0; 
     unsigned short b_digit = 0; 
     if (i < (int)a.digits.size()) { 
      a_digit = a.digits[i]; 
     } 
     if (i < (int)b.digits.size()) { 
      b_digit = b.digits[i]; 
     } 
     sum_n_carry += a_digit + b_digit; 
     c.digits[i] = (sum_n_carry & 0xFFFF); 
     sum_n_carry >>= 16; 
    } 
    if (sum_n_carry != 0) { 
     putCarryInfront(c, sum_n_carry); 
    } 
    while (c.digits.size() > 1 && c.digits.back() == 0) { 
     c.digits.pop_back(); 
    } 
    //std::cout << a.toString() << " + " << b.toString() << " == " << c.toString() << std::endl; 
} 

void BigInteger::subWithoutSign(BigInteger& c, const BigInteger& a, const BigInteger& b) { 
    int sub_n_borrow = 0; 
    int n = a.digits.size(); 
    if (n < (int)b.digits.size()) 
     n = (int)b.digits.size(); 
    c.digits.resize(n); 
    for (int i = 0; i < n; ++i) { 
     unsigned short a_digit = 0; 
     unsigned short b_digit = 0; 
     if (i < (int)a.digits.size()) 
      a_digit = a.digits[i]; 
     if (i < (int)b.digits.size()) 
      b_digit = b.digits[i]; 
     sub_n_borrow += a_digit - b_digit; 
     if (sub_n_borrow >= 0) { 
      c.digits[i] = sub_n_borrow; 
      sub_n_borrow = 0; 
     } else { 
      c.digits[i] = 0x10000 + sub_n_borrow; 
      sub_n_borrow = -1; 
     } 
    } 
    while (c.digits.size() > 1 && c.digits.back() == 0) { 
     c.digits.pop_back(); 
    } 
    //std::cout << a.toString() << " - " << b.toString() << " == " << c.toString() << std::endl; 
} 

int BigInteger::cmpWithoutSign(const BigInteger& a, const BigInteger& b) { 
    int n = (int)a.digits.size(); 
    if (n < (int)b.digits.size()) 
     n = (int)b.digits.size(); 
    //std::cout << "cmp(" << a.toString() << ", " << b.toString() << ") == "; 
    for (int i = n-1; i >= 0; --i) { 
     unsigned short a_digit = 0; 
     unsigned short b_digit = 0; 
     if (i < (int)a.digits.size()) 
      a_digit = a.digits[i]; 
     if (i < (int)b.digits.size()) 
      b_digit = b.digits[i]; 
     if (a_digit < b_digit) { 
      //std::cout << "-1" << std::endl; 
      return -1; 
     } else if (a_digit > b_digit) { 
      //std::cout << "+1" << std::endl; 
      return +1; 
     } 
    } 
    //std::cout << "0" << std::endl; 
    return 0; 
} 

void BigInteger::multByDigitWithoutSign(BigInteger& c, const BigInteger& a, unsigned short b) { 
    unsigned int mult_n_carry = 0; 
    c.digits.clear(); 
    c.digits.resize(a.digits.size()); 
    for (int i = 0; i < (int)a.digits.size(); ++i) { 
     unsigned short a_digit = 0; 
     unsigned short b_digit = b; 
     if (i < (int)a.digits.size()) 
      a_digit = a.digits[i]; 
     mult_n_carry += a_digit * b_digit; 
     c.digits[i] = (mult_n_carry & 0xFFFF); 
     mult_n_carry >>= 16; 
    } 
    if (mult_n_carry != 0) { 
     putCarryInfront(c, mult_n_carry); 
    } 
    //std::cout << a.toString() << " x " << b << " == " << c.toString() << std::endl; 
} 

void BigInteger::shiftLeftByBase(BigInteger& b, const BigInteger& a, int times) { 
    b.digits.resize(a.digits.size() + times); 
    for (int i = 0; i < times; ++i) { 
     b.digits[i] = 0; 
    } 
    for (int i = 0; i < (int)a.digits.size(); ++i) { 
     b.digits[i + times] = a.digits[i]; 
    } 
} 

void BigInteger::shiftRight(BigInteger& a) { 
    //std::cout << "shr " << a.toString() << " == "; 
    for (int i = 0; i < (int)a.digits.size(); ++i) { 
     a.digits[i] >>= 1; 
     if (i+1 < (int)a.digits.size()) { 
      if ((a.digits[i+1] & 0x1) != 0) { 
       a.digits[i] |= 0x8000; 
      } 
     } 
    } 
    //std::cout << a.toString() << std::endl; 
} 

void BigInteger::shiftLeft(BigInteger& a) { 
    bool lastBit = false; 
    for (int i = 0; i < (int)a.digits.size(); ++i) { 
     bool bit = (a.digits[i] & 0x8000) != 0; 
     a.digits[i] <<= 1; 
     if (lastBit) 
      a.digits[i] |= 1; 
     lastBit = bit; 
    } 
    if (lastBit) { 
     a.digits.push_back(1); 
    } 
} 

void BigInteger::putCarryInfront(BigInteger& a, unsigned short carry) { 
    BigInteger b; 
    b.negative = a.negative; 
    b.digits.resize(a.digits.size() + 1); 
    b.digits[a.digits.size()] = carry; 
    for (int i = 0; i < (int)a.digits.size(); ++i) { 
     b.digits[i] = a.digits[i]; 
    } 
    a.digits.swap(b.digits); 
} 

void BigInteger::divideWithoutSign(BigInteger& c, BigInteger& d, const BigInteger& a, const BigInteger& b) { 
    c.digits.clear(); 
    c.digits.push_back(0); 
    BigInteger two("2"); 
    BigInteger e = b; 
    BigInteger f("1"); 
    BigInteger g = a; 
    BigInteger one("1"); 
    while (cmpWithoutSign(g, e) >= 0) { 
     shiftLeft(e); 
     shiftLeft(f); 
    } 
    shiftRight(e); 
    shiftRight(f); 
    while (cmpWithoutSign(g, b) >= 0) { 
     g -= e; 
     c += f; 
     while (cmpWithoutSign(g, e) < 0) { 
      shiftRight(e); 
      shiftRight(f); 
     } 
    } 
    e = c; 
    e *= b; 
    f = a; 
    f -= e; 
    d = f; 
} 

BigInteger::BigInteger(const BigInteger& other) { 
    digits = other.digits; 
    negative = other.negative; 
} 

BigInteger::BigInteger(const char* other) { 
    digits.push_back(0); 
    negative = false; 
    BigInteger ten; 
    ten.digits[0] = 10; 
    const char* c = other; 
    bool make_negative = false; 
    if (*c == '-') { 
     make_negative = true; 
     ++c; 
    } 
    while (*c != 0) { 
     BigInteger digit; 
     digit.digits[0] = *c - '0'; 
     *this *= ten; 
     *this += digit; 
     ++c; 
    } 
    negative = make_negative; 
} 

bool BigInteger::isOdd() const { 
    return (digits[0] & 0x1) != 0; 
} 

BigInteger& BigInteger::operator=(const BigInteger& other) { 
    if (this == &other) // handle self assignment 
     return *this; 
    digits = other.digits; 
    negative = other.negative; 
    return *this; 
} 

BigInteger& BigInteger::operator+=(const BigInteger& other) { 
    BigInteger result; 
    if (negative) { 
     if (other.negative) { 
      result.negative = true; 
      addWithoutSign(result, *this, other); 
     } else { 
      int a = cmpWithoutSign(*this, other); 
      if (a < 0) { 
       result.negative = false; 
       subWithoutSign(result, other, *this); 
      } else if (a > 0) { 
       result.negative = true; 
       subWithoutSign(result, *this, other); 
      } else { 
       result.negative = false; 
       result.digits.clear(); 
       result.digits.push_back(0); 
      } 
     } 
    } else { 
     if (other.negative) { 
      int a = cmpWithoutSign(*this, other); 
      if (a < 0) { 
       result.negative = true; 
       subWithoutSign(result, other, *this); 
      } else if (a > 0) { 
       result.negative = false; 
       subWithoutSign(result, *this, other); 
      } else { 
       result.negative = false; 
       result.digits.clear(); 
       result.digits.push_back(0); 
      } 
     } else { 
      result.negative = false; 
      addWithoutSign(result, *this, other); 
     } 
    } 
    negative = result.negative; 
    digits.swap(result.digits); 
    return *this; 
} 

BigInteger& BigInteger::operator-=(const BigInteger& other) { 
    BigInteger neg_other = other; 
    neg_other.negative = !neg_other.negative; 
    return *this += neg_other; 
} 

BigInteger& BigInteger::operator*=(const BigInteger& other) { 
    BigInteger result; 
    for (int i = 0; i < (int)digits.size(); ++i) { 
     BigInteger mult; 
     multByDigitWithoutSign(mult, other, digits[i]); 
     BigInteger shift; 
     shiftLeftByBase(shift, mult, i); 
     BigInteger add; 
     addWithoutSign(add, result, shift); 
     result = add; 
    } 
    if (negative != other.negative) { 
     result.negative = true; 
    } else { 
     result.negative = false; 
    } 
    //std::cout << toString() << " x " << other.toString() << " == " << result.toString() << std::endl; 
    negative = result.negative; 
    digits.swap(result.digits); 
    return *this; 
} 

BigInteger& BigInteger::operator/=(const BigInteger& other) { 
    BigInteger result, tmp; 
    divideWithoutSign(result, tmp, *this, other); 
    result.negative = (negative != other.negative); 
    negative = result.negative; 
    digits.swap(result.digits); 
    return *this; 
} 

BigInteger& BigInteger::operator%=(const BigInteger& other) { 
    BigInteger c, d; 
    divideWithoutSign(c, d, *this, other); 
    *this = d; 
    return *this; 
} 

bool BigInteger::operator>(const BigInteger& other) const { 
    if (negative) { 
     if (other.negative) { 
      return cmpWithoutSign(*this, other) < 0; 
     } else { 
      return false; 
     } 
    } else { 
     if (other.negative) { 
      return true; 
     } else { 
      return cmpWithoutSign(*this, other) > 0; 
     } 
    } 
} 

BigInteger& BigInteger::powAssignUnderMod(const BigInteger& exponent, const BigInteger& modulus) { 
    BigInteger zero("0"); 
    BigInteger one("1"); 
    BigInteger e = exponent; 
    BigInteger base = *this; 
    *this = one; 
    while (cmpWithoutSign(e, zero) != 0) { 
     //std::cout << e.toString() << " : " << toString() << " : " << base.toString() << std::endl; 
     if (e.isOdd()) { 
      *this *= base; 
      *this %= modulus; 
     } 
     shiftRight(e); 
     base *= BigInteger(base); 
     base %= modulus; 
    } 
    return *this; 
} 

std::string BigInteger::toString() const { 
    std::ostringstream os; 
    if (negative) 
     os << "-"; 
    BigInteger tmp = *this; 
    BigInteger zero("0"); 
    BigInteger ten("10"); 
    tmp.negative = false; 
    std::stack<char> s; 
    while (cmpWithoutSign(tmp, zero) != 0) { 
     BigInteger tmp2, tmp3; 
     divideWithoutSign(tmp2, tmp3, tmp, ten); 
     s.push((char)(tmp3.digits[0] + '0')); 
     tmp = tmp2; 
    } 
    while (!s.empty()) { 
     os << s.top(); 
     s.pop(); 
    } 
    /* 
    for (int i = digits.size()-1; i >= 0; --i) { 
     os << digits[i]; 
     if (i != 0) { 
      os << ","; 
     } 
    } 
    */ 
    return os.str(); 

Y un ejemplo de uso.

BigInteger a("87682374682734687"), b("435983748957348957349857345"), c("2348927349872344") 

// Will Calculate pow(87682374682734687, 435983748957348957349857345) % 2348927349872344 
a.powAssignUnderMod(b, c); 

También es rápido y tiene un número ilimitado de dígitos.

Answer 5

package playTime; 

    public class play { 

     public static long count = 0; 
     public static long binSlots = 10; 
     public static long y = 645; 
     public static long finalValue = 1; 
     public static long x = 11; 

     public static void main(String[] args){ 

      int[] binArray = new int[]{0,0,1,0,0,0,0,1,0,1}; 

      x = BME(x, count, binArray); 

      System.out.print("\nfinal value:"+finalValue); 

     } 

     public static long BME(long x, long count, int[] binArray){ 

      if(count == binSlots){ 
       return finalValue; 
      } 

      if(binArray[(int) count] == 1){ 
       finalValue = finalValue*x%y; 
      } 

      x = (x*x)%y; 
      System.out.print("Array("+binArray[(int) count]+") " 
          +"x("+x+")" +" finalVal("+    finalValue + ")\n"); 

      count++; 


      return BME(x, count,binArray); 
     } 

    } 

Answer 6

LL es para long long int

LL power_mod(LL a, LL k) { 
    if (k == 0) 
     return 1; 
    LL temp = power(a, k/2); 
    LL res; 

    res = ((temp % P) * (temp % P)) % P; 
    if (k % 2 == 1) 
     res = ((a % P) * (res % P)) % P; 
    return res; 
} 

utilizar la función recursiva de arriba para encontrar la exp mod del número. Esto no dará lugar a un desbordamiento porque se calcula de abajo hacia arriba.

Muestra de prueba de funcionamiento para: a = 2 y k = 168277 muestra la salida para ser 518.358 que es correcta y la función se ejecuta en O(log(k)) tiempo;