使用原始类型进行模乘的方法

Posted 2023-02-16

【中文标题】使用原始类型进行模乘的方法

【英文标题】：Ways to do modulo multiplication with primitive types

【发布时间】：2012-08-28 22:21:39

【问题描述】：

有没有办法构建例如(853467 * 21660421200929) % 100000000000007 没有 BigInteger 库（注意每个数字都适合 64 位整数但乘法结果不适合）？

这个解决方案似乎效率低：

int64_t mulmod(int64_t a, int64_t b, int64_t m) 
    if (b < a)
        std::swap(a, b);
    int64_t res = 0;
    for (int64_t i = 0; i < a; i++) 
        res += b;
        res %= m;
    
    return res;

【问题讨论】：

一方面，我建议摆脱 Microsoft 扩展并使用 int64_t。

看起来在这种情况下你可以作弊，因为你不关心大于参数__int64 m（或uint64_t对于那些赞成它的人）的任何东西，因此你只能处理64-位类型。

你读过Montgomery reduction算法吗？

@ildjarn：不，不知道，谢谢你的链接！

有趣，这在 x64 汇编中是微不足道的。

【参考方案1】：

您应该使用Russian Peasant multiplication。它使用重复加倍来计算所有值(b*2^i)%m，如果设置了a 的ith 位，则将它们相加。

uint64_t mulmod(uint64_t a, uint64_t b, uint64_t m) 
    int64_t res = 0;
    while (a != 0) 
        if (a & 1) res = (res + b) % m;
        a >>= 1;
        b = (b << 1) % m;
    
    return res;

它改进了您的算法，因为它需要 O(log(a)) 时间，而不是 O(a) 时间。

注意事项：无符号，仅当 m 为 63 位或更少时才有效。

【讨论】：

应该将res 声明为uint64_t？

【参考方案2】：

Keith Randall's answer 很好，但正如他所说，需要注意的是，它只有在 m 为 63 位或更少时才有效。

这是一个有两个优点的修改：

即使m 是 64 位，它也可以工作。 它不需要使用模运算，这在某些处理器上可能很昂贵。

(注意res -= m 和temp_b -= m 行依赖64位无符号整数溢出来给出预期的结果。这应该没问题，因为无符号整数溢出在C和C++中定义良好。为此因为它是important to use unsigned integer types。）

uint64_t mulmod(uint64_t a, uint64_t b, uint64_t m) 
    uint64_t res = 0;
    uint64_t temp_b;

    /* Only needed if b may be >= m */
    if (b >= m) 
        if (m > UINT64_MAX / 2u)
            b -= m;
        else
            b %= m;
    

    while (a != 0) 
        if (a & 1) 
            /* Add b to res, modulo m, without overflow */
            if (b >= m - res) /* Equiv to if (res + b >= m), without overflow */
                res -= m;
            res += b;
        
        a >>= 1;

        /* Double b, modulo m */
        temp_b = b;
        if (b >= m - b)       /* Equiv to if (2 * b >= m), without overflow */
            temp_b -= m;
        b += temp_b;
    
    return res;

【讨论】：

我喜欢这个，因为它处理完整的 64 位值。如果你测试b < a是否在顶部，如果是，交换a和b，它可以显着加快时间，因为它更有可能是while循环可以提前退出。

如果m>UINT64_MAX / 2u，你为什么不能做b %= m？模运算会神奇地变得不稳定吗？

你绝对可以做到b %= m。但是，模运算可能很慢（取决于处理器），因此如果可能的话值得避免。因此，if (m > UINT64_MAX / 2u) b -= m; 是在m 很大的情况下避免模运算的可能优化，因此可以将模简化为简单的减法。

此评论可能有点晚了，但 res += b 和 b += temp_b 也会溢出，即使特别提到了 -= 操作，您的回答中也没有提及。不是来自 C++ 背景，所以我不确定这是否是标准行为，但也许可以在你的答案中添加？

【参考方案3】：

这两种方法都适合我。第一个与您的相同，但我将您的数字更改为明确的 ULL。第二个使用汇编符号，它应该工作得更快。 密码学中也使用了一些算法（我猜主要是基于 RSA 和 RSA 的密码学），就像已经提到的蒙哥马利减少一样，但我认为实现它们需要时间。

#include <algorithm>
#include <iostream>

__uint64_t mulmod1(__uint64_t a, __uint64_t b, __uint64_t m) 
  if (b < a)
    std::swap(a, b);
  __uint64_t res = 0;
  for (__uint64_t i = 0; i < a; i++) 
    res += b;
    res %= m;
  
  return res;


__uint64_t mulmod2(__uint64_t a, __uint64_t b, __uint64_t m) 
  __uint64_t r;
  __asm__
  ( "mulq %2\n\t"
      "divq %3"
      : "=&d" (r), "+%a" (a)
      : "rm" (b), "rm" (m)
      : "cc"
  );
  return r;


int main() 
  using namespace std;
  __uint64_t a = 853467ULL;
  __uint64_t b = 21660421200929ULL;
  __uint64_t c = 100000000000007ULL;

  cout << mulmod1(a, b, c) << endl;
  cout << mulmod2(a, b, c) << endl;
  return 0;

【讨论】：

我不知道内联汇编器，它也使用循环吗？

@ChristianAmmer 不，它不需要一个。它使用双倍宽度乘法和除法总是双倍宽度。只有在高级语言中，乘法的高部分会突然丢失。

这个例子没问题，但是如果(a * b > (2^64 - 1) * c)会失败。但我假设 OP 意味着隐含的商也是 64 位值。

关于循环：我不知道汇编程序如何计算乘法，我的意思是在幕后，但在 C++ 上不需要循环，因为我们知道，由于 %uint64，结果最多64位。 @Brett 3 个数字是 64 位的。

@Benjamin - 我只是指出汇编实现并不通用。试试：a=8534670000000000000、b=216604212009290。两者都是 64 位的，但divq 会导致异常。

【参考方案4】：

对重复加倍算法的改进是检查一次可以计算多少位而不会溢出。可以对这两个参数进行提前退出检查——加速（不太可能？）N 不是素数的事件。

例如100000000000007 == 0x00005af3107a4007，允许每次迭代计算 16（或 17）位。示例中实际迭代次数为 3。

// just a conceptual routine
int get_leading_zeroes(uint64_t n)

   int a=0;
   while ((n & 0x8000000000000000) == 0)  a++; n<<=1; 
   return a;


uint64_t mulmod(uint64_t a, uint64_t b, uint64_t n)

     uint64_t result = 0;
     int N = get_leading_zeroes(n);
     uint64_t mask = (1<<N) - 1;
     a %= n;
     b %= n;  // Make sure all values are originally in the proper range?
     // n is not necessarily a prime -- so both a & b can end up being zero
     while (a>0 && b>0)
     
         result = (result + (b & mask) * a) % n;  // no overflow
         b>>=N;
         a = (a << N) % n;
     
     return result;

【讨论】：

+1，不错的速度改进，但建议检查 n==0 以防止 get_leading_zeroes() 出现无限循环。

【参考方案5】：

您可以尝试将乘法分解为加法：

// compute (a * b) % m:

unsigned int multmod(unsigned int a, unsigned int b, unsigned int m)

    unsigned int result = 0;

    a %= m;
    b %= m;

    while (b)
    
        if (b % 2 != 0)
        
            result = (result + a) % m;
        

        a = (a * 2) % m;
        b /= 2;
    

    return result;

【讨论】：

+1 是一个可行的解决方案，我必须考虑一下才能完全理解，但它之所以有效，是因为(a * b) == (a * 2) * (b / 2)，对吧？

某些输入实际上会失败。如果m 大于1 << 63（或1 << 31，如果int 是32 位），a * 2 可能会溢出并错误地减少。

你实际上可以在每一步中减少(~0ULL/m)。例如。对于 100000000000007，您可以使用 131072 (1<<17) 而不是 2。这也解释了 harold 的评论；对于这么大的m，步长变为 1，您没有任何进展。

@harold：你是对的：第一个因素不能设置其最高位。不过，我相信这是对当前算法的函数参数值的唯一限制。

【参考方案6】：

a * b % m 等于 a * b - (a * b / m) * m

使用浮点算法逼近a * b / m。近似值留下了一个足够小的值，用于正常的 64 位整数运算，m 最多 63 位。

此方法受double 的有效位限制，通常为52 位。

uint64_t mod_mul_52(uint64_t a, uint64_t b, uint64_t m) 
    uint64_t c = (double)a * b / m - 1;
    uint64_t d = a * b - c * m;

    return d % m;

此方法受long double 的有效位限制，通常为64 位或更大。整数运算限制为 63 位。

uint64_t mod_mul_63(uint64_t a, uint64_t b, uint64_t m) 
    uint64_t c = (long double)a * b / m - 1;
    uint64_t d = a * b - c * m;

    return d % m;

这些方法要求a 和b 小于m。要处理任意的a 和b，请在计算c 之前添加这些行。

a = a % m;
b = b % m;

在这两种方法中，最终的% 操作都可以是有条件的。

return d >= m ? d % m : d;

【参考方案7】：

我可以建议对您的算法进行改进。

您实际上是通过每次添加b 来迭代计算a * b，在每次迭代后进行取模。最好每次都加上b * x，而x是确定的，这样b * x就不会溢出了。

int64_t mulmod(int64_t a, int64_t b, int64_t m)

    a %= m;
    b %= m;

    int64_t x = 1;
    int64_t bx = b;

    while (x < a)
    
        int64_t bb = bx * 2;
        if (bb <= bx)
            break; // overflow

        x *= 2;
        bx = bb;
    

    int64_t ans = 0;

    for (; x < a; a -= x)
        ans = (ans + bx) % m;

    return (ans + a*b) % m;

【讨论】：

你不能使用x=(1<<63-m)/b 吗？这是四舍五入，所以b*x <= 1<<63 - m，它不需要循环来计算。不会改变 big-O，因为 for 循环的迭代次数减少了

@MSalters：这不是引入了一个可能很昂贵的部门吗？

@CraigMcQueen：是的，但只有一个，而且循环中已经有一个模数。

这不是改进——我不知道最后一个循环的 O()，但只取 3 个“随机”数字，循环就运行了 300000 次迭代。