C++ 中大数的模块化乘法

modular multiplication of large numbers in c++

本文关键字：模块化 C++      更新时间：2023-10-16

我有三个整数A，B（小于10^12）和C（小于10^15）。我想计算 （A * B） % C。我知道那件事

(A * B) % C = ((A % C) * (B % C)) % C

但是假设如果A = B = 10^11，那么上面的表达式将导致整数溢出。对于上述情况，是否有任何简单的解决方案，或者我必须使用快速乘法算法。

如果我必须使用快速乘法算法，那么我应该使用哪种算法。

编辑：我已经尝试了上述问题C++（不会导致溢出，不知道为什么），但答案不应该是零吗？

提前谢谢。

你可以使用Schrage的方法解决这个问题。这允许您将两个有符号数相乘a，并以一定的模数m z两者，而不会生成大于该模数的中间数。

它基于模m的近似分解，

m = aq + r

即

q = [m / a]

和

r = m mod a

其中[]表示整数部分。如果 r < q 和 0 < z < m − 1 ，则 a(z mod q) 和 r[z / q] 都在 0,...,m − 1 和

az mod m = a(z mod q) − r[z / q]

如果这是负数，则添加m。

[这种技术经常用于线性同余随机数生成器]。

给定您的公式和以下变体：

(A + B) mod C = ((A mod C) + (B mod C)) mod C

您可以使用分而治之的方法开发一种既简单又快速的算法：

#include <iostream>
long bigMod(long  a, long  b, long c) {
    if (a == 0 || b == 0) {
        return 0;
    }
    if (a == 1) {
        return b;
    }
    if (b == 1) {
        return a;
    } 
    // Returns: (a * b/2) mod c
    long a2 = bigMod(a, b / 2, c);
    // Even factor
    if ((b & 1) == 0) {
        // [((a * b/2) mod c) + ((a * b/2) mod c)] mod c
        return (a2 + a2) % c;
    } else {
        // Odd exponent
        // [(a mod c) + ((a * b/2) mod c) + ((a * b/2) mod c)] mod c
        return ((a % c) + (a2 + a2)) % c;
    }
}
int main() { 
    // Use the min(a, b) as the second parameter
    // This prints: 27
    std::cout << bigMod(64545, 58971, 144) << std::endl;
    return 0;
}

这是O(log N)

已更新：修复了设置高位a % c时的错误。（帽子提示：凯文·霍普斯）

如果您正在寻找简单而不是快速，那么您可以使用以下内容：

typedef unsigned long long u64;
u64 multiplyModulo(u64 a, u64 b, u64 c)
{
    u64 result = 0;
    a %= c;
    b %= c;
    while(b) {
        if(b & 0x1) {
            result += a;
            result %= c;
        }
        b >>= 1;
        if(a < c - a) {
            a <<= 1;
        } else {
            a -= (c - a);
        }
    }
    return result;
}

抱歉，当变量"a"保留设置了高位的值时，godel9 的算法将产生不正确的结果。这是因为"a <<= 1"会丢失信息。这是一个更正的算法，适用于任何整数类型，有符号或无符号。

template <typename IntType>
IntType add(IntType a, IntType b, IntType c)
    {
    assert(c > 0 && 0 <= a && a < c && 0 <= b && b < c);
    IntType room = (c - 1) - a;
    if (b <= room)
        a += b;
    else
        a = b - room - 1;
    return a;
    }
template <typename IntType>
IntType mod(IntType a, IntType c)
    {
    assert(c > 0);
    IntType q = a / c; // q may be negative
    a -= q * c; // now -c < a && a < c
    if (a < 0)
        a += c;
    return a;
    }
template <typename IntType>
IntType multiplyModulo(IntType a, IntType b, IntType c)
    {
    IntType result = 0;
    a = mod(a, c);
    b = mod(b, c);
    if (b > a)
        std::swap(a, b);
    while (b)
        {
        if (b & 0x1)
            result = add(result, a, c);
        a = add(a, a, c);
        b >>= 1;
        }
    return result;
    }

在这种情况下，A 和 B 是 40 位数字，C 是 50 位数字，这在 64 位模式下不是问题，如果您有内部函数或可以编写汇编代码以使用 64 位乘以 64 位乘法产生 128 位结果（乘积实际上是 80 位），之后您将 128 位除以 50 位除数以产生 50 位余数（模）。

根据处理器的不同，通过乘以 81 位（或更少）常量来实现除以 50 位常量可能会更快。再次假设 64 位处理器，它将需要 4 次乘法和一些加法，然后移动 4 乘积的上位以产生商。然后使用商乘以 50 位模数并减去（从 80 位乘积）来产生 50 位余数。

相关文章：