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

Итак, я недавно работал над реализацией теста простоты Миллера-Рабина. Я ограничиваю его набором всех 32-битных чисел, потому что это проект для развлечения, который я делаю, чтобы познакомиться с C ++, и я не хочу работать с чем-то 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 - это исключительно большое число, и где-то в процессе оно переполняет test, что приводит к неправильной оценке.

с обратной стороны, такие аргументы, как

4000111222 ^ 3 % 1608

также оценивают неправильно по той же причине.

Есть ли у кого-нибудь предложения по модульному возведению в степень, которые могут предотвратить это переполнение и / или манипулировать им для получения правильного результата? (как я вижу, переполнение - это просто еще одна форма по модулю, то есть num% (UINT_MAX + 1))


person Axel Magnuson    schedule 05.02.2010    source источник


Ответы (6)


Возведение в квадрат путем возведения в квадрат все еще "работает" для возведения в степень по модулю. Ваша проблема не в том, что 2 ^ 168277 - исключительно большое число, а в том, что одним из ваших промежуточных результатов является довольно большое число (больше 2 ^ 32), потому что 673109 больше, чем 2 ^ 16.

Так что я думаю, что подойдет следующее. Возможно, я упустил какую-то деталь, но основная идея работает, и именно так «настоящий» криптографический код может выполнять большое мод-возведение в степень (хотя и не с 32- и 64-битными числами, а с большими числами, которые никогда не должны становиться больше, чем 2 * лог (модуль)):

  • Начните с возведения в степень возведения в квадрат, как вы это сделали.
  • Выполните фактическое возведение в квадрат 64-битного целого числа без знака.
  • Уменьшайте по модулю 673109 на каждом шаге, чтобы вернуться в 32-битный диапазон, как вы это делаете.

Очевидно, что это немного неудобно, если в вашей реализации на C ++ нет 64-битного целого числа, хотя вы всегда можете подделать его.

На слайде 22 есть пример: http://www.cs.princeton.edu/courses/archive/spr05/cos126/lectures/22.pdf, хотя в нем используются очень маленькие числа (менее 2 ^ 16), поэтому он может не проиллюстрировать то, чего вы еще не сделали. знать.

Ваш другой пример, 4000111222 ^ 3 % 1608 будет работать в вашем текущем коде, если вы просто уменьшите 4000111222 по модулю 1608 перед тем, как начать. 1608 достаточно мал, чтобы вы могли безопасно умножить любые два числа mod-1608 на 32-битное int.

person Steve Jessop    schedule 05.02.2010
comment
спасибо, чувак, это помогло. Просто из любопытства, знаете ли вы какие-либо методы, которые не потребовали бы использования большего объема памяти? Я уверен, что они пригодятся. - person Axel Magnuson; 05.02.2010
comment
Не то, что я знаю из. Вам нужно умножить два числа до 673108, mod 673109. Очевидно, что вы можете разбить вещи и выполнить длинное умножение с меньшими цифрами, скажем, 2 ^ 10. Но как только вы реализуете умножение и деление в программном обеспечении, вы можете просто реализовать его для особого случая умножения двух 32-битных значений, чтобы получить 64-битный результат, а затем деления для извлечения 32-битного остатка. Могут быть некоторые жесткие оптимизации, которые делают необходимый минимум, но я их не знаю, и подделать 64-битное int в C ++ не так сложно. - person Steve Jessop; 05.02.2010

Я недавно написал кое-что для RSA на C ++, хотя и немного запутанно.

#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);

Это тоже быстро и имеет неограниченное количество цифр.

person clinux    schedule 15.07.2010
comment
Спасибо, что поделился! Вопрос, это цифра std :: vector ‹unsigned short›? - person darius; 23.06.2016
comment
Да, но под капотом работает база 65536, а не база 10. - person clinux; 30.09.2017

Две вещи:

  • Вы используете соответствующий тип данных? Другими словами, позволяет ли UINT_MAX использовать в качестве аргумента 673109?

Нет, это не так, поскольку в какой-то момент у вас есть Ваш код, который не работает, потому что в какой-то момент у вас есть num = 2^16, а num = ... вызывает переполнение. Используйте больший тип данных для хранения этого промежуточного значения.

  • Как насчет того, чтобы взять по модулю все возможные возможности переполнения, такие как:

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

Редактировать:

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;
}
person dirkgently    schedule 05.02.2010

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)) время;

person abkds    schedule 09.09.2014

Вы можете использовать следующую личность:

(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;
person Alexander Poluektov    schedule 05.02.2010
comment
спасибо вам обоим за ваши предложения, но по характеру алгоритма и test, и num всегда будут меньше, чем mod, поэтому: {(test% mod) = test} и {(num% mod) = test}, поэтому идентичность может это мне не поможет, потому что функция не работает, даже если num и test меньше mod. Кроме того, беззнаковые целые числа позволяют мне использовать 673109 в качестве аргумента. UINT_MAX = 4 294 967 295 для моего компьютера. - person Axel Magnuson; 05.02.2010

person    schedule
comment
это был код, который я написал на java очень быстро. Я использовал пример 11 ^ 644mod 645. = 1. мы знаем, что двоичный код 645 равен 1010000100. Я как бы обманул и жестко закодировал переменные, но он работает нормально. - person ShowLove; 04.10.2013
comment
вывод был Array (0) x (121) finalVal (1) Array (0) x (451) finalVal (1) Array (1) x (226) finalVal (451) Array (0) x (121) finalVal (451) Массив (0) x (451) finalVal (451) Массив (0) x (226) finalVal (451) Массив (0) x (121) finalVal (451) Массив (1) x (451) finalVal (391) Массив ( 0) x (226) finalVal (391) Массив (1) x (121) finalVal (1) конечное значение: 1 - person ShowLove; 04.10.2013