Точный модулярной арифметики с двойной точностью


У меня есть функция ProdBigMod которая вычисляет произведение двух двойных точные цифры \$x_1 места\$, \$x_2\$ (оба из которых являются меньше, чем \$2^{53} - 1\$) и впоследствии находит остаток \$\Mod Р\$ (где$Р\$ - простое). Я не могу использовать внешние библиотеки и должны использовать тип данных Double.

Главная проблема возникает тогда, когда продукт \x_1 места $\$ и \x_2 $\$ превышает \$2^{53} - 1\$.

Как @TobySpeights указывает, 53 важно, потому что это число mantissa бит (также называется significand отсюда и название постоянной Significand53) для двойных типов данных (см. двойной точности). Мы можем исключить некоторые из этих вопросов, сначала убедившись, что обе \$х\$ие менее \$р\$ (это достигается путем немедленного применения std::fmod). В самом деле, когда p < sqrt(2^53 - 1)мы знаем, продукт prodX = x1 * x2 < 2^53 - 1. Бороться с случае, когда p > 2^53 - 1 и prodX > 2^53 - 1мы можем воспользоваться некоторые свойства модулярной арифметики.

А именно:

$$(x_1 места \cDOT на x_2) \PMOD позволяет п = x_1 места \cDOT на (x_ не{12} + x_ не{22} + ... + x_ не{Н2}) \PMOD позволяет п$$

Где \$(x_ не{12} + x_ не{22} + ... + x_ не{Н2}) = x_2\$.

Это дает:

$$большая \( x_1 места \cDOT на x_ не{12} \PMOD позволяет п \Большой) \место\пространство + \космос\космический \большой(x_1 места \cDOT на x_ не{22} \PMOD позволяет п\Большой) \место\пространство + ... + \космос\космический \большой(x_1 места \cDOT на x_ не{Н2} \PMOD позволяет п\Большой)$$

Теперь, учитывая, что каждый из \$x_ не{И2}\$ (chunkSize ниже) являются идентичными, за исключением, возможно, окончательным (т. е. \$x_ не{Н2}\$), мы можем вычислить следующий раз:

$x_1 места $\cDOT на x_ не{12} \PMOD позволяет п \место\пространство = \космос\космический x_1 места \cDOT на chunkSize \PMOD позволяет п \место\пространство = \космос\космический chunkMod$$

Чтобы найти окончательный \$x_ не{Н2}\$ мы сначала должны определить, сколько \$x_ не{И2}\$S будет идти в \x_2$\$:

numChunkMods = floor(x2 / chunkSize)

Теперь мы можем легко получить xn2:

xn2 = (x2 - chunkSize * numChunkMods)  ==>>  part2 = x1 * xn2 (mod p)

Поставив все это, Eq. (1) можно сократить до:

x1 * x2 (mod p) = (chunkMod * numChunkMods) + part2 (mod p)

Это хорошее начало, но не полностью вышли из леса, так как мы не уверены, что продукт numChunkMods * chunkMod < 2^53 - 1. Чтобы обойти это, мы продолжаем процесс над установкой x1 = chunkMod и x2 = numChunkMods, пока продукт не менее 2^53 - 1.

#include <iostream>
#include <cmath>

const double Significand53 = 9007199254740991.0;
const double SqrtSig53 = std::floor(std::sqrt(Significand53));

double PositiveMod(double x, double m) {
    if (x < 0)
        x = x + ceil(std::abs(x) / m) * m;
    else if (x > m)
        x = std::fmod(x, m);
    return x;
}

double ProdBigMod(double x1, double x2, double p) {
    double result = 0, prodX;

    x1 = PositiveMod(x1, p);
    x2 = PositiveMod(x2, p);
    prodX = x1 * x2;

    if (prodX < p) {
        result = prodX;
    } else if (p < SqrtSig53 || prodX < Significand53) {
        result = std::fmod(prodX, p);
    } else {
        double numChunkMods, part1 = Significand53;
        double chunkSize, chunkMod, part2;

        while (part1 >= Significand53) {
            chunkSize = std::floor(Significand53 / x1);  // Ensures chunkMod < 2^53 - 1
            chunkMod = std::fmod(x1 * chunkSize, p);
            numChunkMods = std::floor(x2 / chunkSize);
            part2 = std::fmod((x2 - chunkSize * numChunkMods) * x1, p);
            part1 = numChunkMods * chunkMod;
            x1 = chunkMod;
            x2 = numChunkMods;
            result = std::fmod(result + part2, p);
        }

        result = std::fmod(part1 + result, p);
    }

    return result;
}

Вот пример вызова функции:

int main() {
    double test = ProdBigMod(914806066069, 497967734853, 732164213243);
    std::cout  << std::fixed;
    std::cout << test << "\n";
    return 0;
}

Output: 85635829849.000000

Я уже сравнил выше на выборку случайных чисел в диапазоне 10^12 С gmp аналоговый и, кажется, чтобы дать правильные результаты. Однако, я не совсем уверен, если моя логика или реализации пуленепробиваемые. Я также интересно, если это может быть более эффективным.

Здесь расположена онлайн "большое количество калькулятор", который вы можете проверить результаты в отношении: https://defuse.ca/big-number-calculator.htm

Например, ввод (914806066069 * 497967734853) % 732164213243 в выражение поле и нажмите кнопку "Рассчитать".



Комментарии
1 ответ

Целое число-Несс

Это должно быть сделано более ясным, что ProdBigMod() чтобы работать только с целым числом значений и не полный спектр double включая значения с дробной частей, Нэн и непостижимых.

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

OFF-на-один:

"... когда p < sqrt(2^53 - 1) ..." должно быть "... когда p <= sqrt(2^53) ...".
Это позволяет немного больше p.

Напомним, что something mod p будет в большинстве, p - 1.

Это отражается на неправильный код: x > m --> x >= m

double PositiveMod(double x, double m) {
if (x < 0)
...
// else if (x > m)
else if (x >= m)
x = std::fmod(x, m);
return x;
}

и

// } else if (p < SqrtSig53 || ...
} else if (p <= SqrtSig53 || ...

Удобоносимость

Использование Significand53 и другой код полагается на double как binary64. Простой тест позволит предотвратить ряд ошибочных сборники, даже если не повысить мобильность.

#if DBL_MANT_DIG != 53
#error TBD code
#endif

Альтернативная точность теста.

Код prodX = x1 * x2; и различные тесты, когда код может проверить FE_INEXACT после умножения.

feclearexcept(FE_INEXACT);
prodX = x1 * x2
if (fetestexcept(FE_INEXACT)) Handle_inexact_product();

Сомнительный код

Я не уверен x = x + ceil(std::abs(x) / m) * m; работать, как ожидалось для всех x, mособенно когда в математике частное это только за целое число, но округление до ceil(). Четкую альтернативу:

double PositiveMod(double x, double m) {
double y = std::fmod(x, m);
if (y < 0.0) {
y = (m < 0) ? y - m : y + m;
}
return y;
}

Альтернативный код:

Использовать 64-разрядные или более целочисленных вычислениях. См mulmodmax() в модульное возведение в степень диапазон без ограничений. По меньшей мере 19 десятичной цифры числа и 12 Здесь.

библиотеки FMOD() правильности.

Примечание: спецификация std::fmod() не требует результат будет самый лучший из возможных ответов, пока с стандартом IEEE 754 присоединения, это. Хорошая библиотека будет осуществлять точный результат. Реф:
Это FMOD() точное когда Y является целым числом?
.


Заключение

В целом, используя ФП математике для решения целочисленной задачи имеет различные неожиданные углу проблемы и так трудно застраховаться код является правильным для всех x1,x2,p.

5
ответ дан 10 февраля 2018 в 08:02 Источник Поделиться