Модульное возведение в степень для больших чисел в С++
Итак, я недавно работал над реализацией теста примитивности Миллера-Рабина. Я ограничиваю его в объеме всех 32-битных чисел, потому что это проект для забавы, который я делаю, чтобы познакомиться с С++, и я не хочу работать с чем-либо 64-битным для какое-то время. Дополнительный бонус заключается в том, что алгоритм детерминирован для всех 32-битных чисел, поэтому я могу значительно повысить эффективность, потому что я точно знаю, какие свидетели должны тестировать.
Итак, для низких чисел алгоритм работает исключительно хорошо. Однако часть процесса зависит от модульного возведения в степень, то есть (num ^ pow)% mod. так, например,
3 ^ 2 % 5 =
9 % 5 =
4
вот код, который я использовал для этого модульного возведения в степень:
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;
}
Как вы, возможно, уже догадались, проблемы возникают, когда аргументы являются исключительно большими числами. Например, если я хочу проверить число 673109 для простоты, мне нужно в какой-то момент найти:
(2 ^ 168277)% 673109
теперь 2 ^ 168277 является исключительно большим числом, и где-то в процессе он переполняет тест, что приводит к неправильной оценке.
на обратной стороне, такие аргументы, как
4000111222 ^ 3% 1608
также неверно оценивают по той же причине.
Есть ли у кого-нибудь предложения по модульной возведения в степень, которые могут предотвратить это переполнение и/или манипулировать им для получения правильного результата? (как я его вижу, переполнение - это еще одна форма модуляции, то есть num% (UINT_MAX + 1))
Ответы
Ответ 1
Экспоненциальное возведение в квадрат по-прежнему "работает" для модуляции экспоненциальности. Ваша проблема заключается не в том, что 2 ^ 168277 является исключительно большим числом, так что один из ваших промежуточных результатов - довольно большое число (больше 2 ^ 32), потому что 673109 больше, чем 2 ^ 16.
Итак, я думаю, что это будет сделано. Возможно, я пропустил деталь, но основная идея работает, и вот как "реальный" криптокод может сделать большое mod-exponentiation (хотя и не с 32 и 64-битными номерами, а с бонусами, которые никогда не должны быть больше, чем 2 * log (модуль)):
- Начните с возведения в степень по квадрату, как и у вас.
- Выполните фактическое возведение в квадрат в 64-битном беззнаковом целое.
- Уменьшите modulo 673109 на каждом шаге, чтобы вернуться в 32-битный диапазон, как и вы.
Очевидно, что немного неудобно, если ваша реализация на С++ не имеет 64-битного целого числа, хотя вы всегда можете подделать ее.
Вот пример на слайде 22 здесь: http://www.cs.princeton.edu/courses/archive/spr05/cos126/lectures/22.pdf, хотя он использует очень маленькие числа (менее 2 ^ 16), поэтому он может не иллюстрировать все, что вы еще не знаете.
Другой пример, 4000111222 ^ 3 % 1608 будет работать в вашем текущем коде, если вы просто уменьшите 4000111222 modulo 1608 до начала. 1608 достаточно мал, чтобы вы могли безопасно размножать любые два номера мод-1608 в 32-битном int.
Ответ 2
Я написал что-то для этого недавно для RSA на С++, но немного беспорядочно.
#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();
И пример использования.
BigInteger a("87682374682734687"), b("435983748957348957349857345"), c("2348927349872344")
// Will Calculate pow(87682374682734687, 435983748957348957349857345) % 2348927349872344
a.powAssignUnderMod(b, c);
Его быстрый тоже, и имеет неограниченное количество цифр.
Ответ 3
Две вещи:
- Используете ли вы соответствующий тип данных? Другими словами, разрешает ли UINT_MAX иметь 673109 в качестве аргумента?
Нет, это не так, поскольку в какой-то момент у вас есть код, который не работает, потому что в какой-то момент у вас есть num = 2^16, а num = ... вызывает переполнение. Для хранения этого промежуточного значения используйте больший тип данных.
-
Как насчет приема по модулю при каждом возможном опережении переполнения, например:
test = ((test % mod) * (num % mod)) % mod;
Edit:
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;
}
Ответ 4
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);
}
}
Ответ 5
LL для 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;
}
Используйте приведенную выше рекурсивную функцию для нахождения mod exp числа. Это не приведет к переполнению, потому что оно вычисляется снизу вверх.
Пример тестового прогона для:
a = 2 и k = 168277 показывает, что результат равен 518358, и функция работает в O(log(k)) время;
Ответ 6
Вы можете использовать следующее удостоверение:
(a * b) (mod m) === (a (mod m)) * (b (mod m)) (mod m)
Попробуйте использовать его простым способом и постепенно улучшайте.
if (pow & 1)
test = ((test % mod) * (num % mod)) % mod;
num = ((num % mod) * (num % mod)) % mod;
