C 中 64 位整数的模乘法

Question

我想在 C 中将两个整数与另一个（质数）整数（长度约为 50 位）相乘。

uint64_t x = ...;
uint64_t y = ...;
uint64_t q = ...;
uint64_t z = (x*y) % q;

这会产生一个小于 q 的 64 位整数，但不幸的是

x%q * y%q

对于 64 位类型来说太大了。所以在我计算余数之前它就溢出了。

我想到了一些解决方法：

• 使用 gcc 的

  __uint128_t

  会违反 C 标准并且无法与其他编译器一起工作
• 将变量转换为

  double

  或

  long double

  ，但由于浮点不准确，结果不正确
• 使用像 GMP 这样的多精度库是可行的，但这会引入新的项目依赖性，并且可能会带来一些开销
• 为此目的编写自己的乘法函数，但据我所知，我需要汇编，因此可移植性消失了

还有其他选择或您会推荐什么吗？

编辑：我忘了提及，只有高效的情况下，自写的解决方案才值得去做。例如，执行数十个低效的 % 操作会很糟糕。

Answer 1

当然你不需要汇编来自己实现任意精度乘法。

有很多著名的算法，并且（根据定义）它们可以用 C 语言实现。

请参阅维基百科了解更多信息。根据我的经验，正确获取这样的代码可能有点棘手，因此也许您最好花时间添加依赖项。

Answer 2

在您的特定情况下，最好将数字划分为四个 13 位部分，然后使用长手方法将它们相乘，对中间结果计算模。

代码如下所示：

uint64_t x = ...;
uint64_t y = ...;
uint64_t q = ...; // at most 50 bits
uint64_t z = 0;
x %= q;
y %= q;
while (1)
{
  z += x*(y & 0b1111111111111);
  z %= q;
  y >>= 13;
  if (y==0) break;
  x <<= 13;
  x %= q;
}

这仅需要 4 次乘法和 9 次模运算。

Answer 3

也许这个功能很有用。 我还没有测试它是否有效，但它应该可以工作，并且它只使用有限数量的模运算。

uint64_t XY_modulo_D_L_N(uint64_t X, uint64_t Y, uint64_t D)
{
    // return (X * Y) modulo D
    //tested for D < 2^62
    uint64_t D_max= (uint64_t)1 << 31;
    uint64_t mod_D;
    X %= D;
    Y %= D;
    if (X <= D_max && Y <= D_max)
    {
        mod_D = (X * Y) % D;
    }
    else if (X <= D_max || Y <= D_max )
    {
        uint64_t N_t;
        uint64_t RD = D % D_max;
        uint64_t KD = D / D_max;
        if (Y <= D_max)
        {
            uint64_t X1 = Y;
            Y = X;
            X = X1;
        }
        mod_D = (X * (Y % D_max)) % D_max;
        N_t = X * (Y / D_max) + (X * (Y % D_max)) / D_max;
        mod_D = (mod_D + (D_max *(N_t % KD)) % D) % D;
        mod_D = (mod_D + D -((N_t / KD) * RD) % D) % D;
    }
    else
    {
        uint64_t RD = D % D_max;
        uint64_t KD = D / D_max;
        uint64_t Rx = X % D_max;
        uint64_t Ry = Y % D_max;
        uint64_t Kx = X / D_max;
        uint64_t Ky = Y / D_max;
        uint64_t mod_D_t,N_t,R_t;

        mod_D =(Rx * Ry) % D;
        mod_D_t = (D_max * D_max) % D;
        N_t = mod_D_t;
        R_t = Kx;
        if (N_t >= (D_max / R_t))
        {
            mod_D_t = (R_t *(N_t % D_max)) % D_max;
            N_t = R_t * (N_t / D_max)+(R_t * (N_t % D_max)) / D_max;
            mod_D_t = (mod_D_t + (D_max * (N_t % KD)) % D) % D;
            N_t = (mod_D_t + D - ((N_t / KD) * RD) % D) % D;
        }
        else
            N_t *= R_t;
        R_t = Ky;
        if (N_t >= (D_max / R_t))
        {
            mod_D_t = (R_t * (N_t % D_max)) % D_max;
            N_t = R_t * (N_t / D_max) + (R_t * (N_t % D_max)) / D_max;
            mod_D_t = (mod_D_t + (D_max * (N_t % KD)) % D) % D;
            mod_D = (mod_D + (mod_D_t + D - ((N_t / KD) * RD) % D) % D) % D;
        }
        else
            mod_D =(mod_D + N_t * R_t) % D;
        N_t = (Rx * D_max) % D;
        R_t = Ky;
        if (N_t >= (D_max / R_t))
        {
            mod_D_t = (R_t * (N_t % D_max)) % D_max;
            N_t = R_t * (N_t / D_max) + (R_t * (N_t % D_max)) / D_max;
            mod_D_t = (mod_D_t + (D_max * (N_t % KD)) % D) % D;
            mod_D = (mod_D + (mod_D_t + D - ((N_t / KD) * RD) % D) % D) % D;
        }
        else
            mod_D = (mod_D + N_t * R_t) % D;
        N_t = (Ry * D_max) % D;
        R_t = Kx;
        if (N_t >= (D_max / R_t))
        {
            mod_D_t = (R_t * (N_t % D_max)) % D_max;
            N_t=R_t * (N_t / D_max) + (R_t * (N_t % D_max)) / D_max;
            mod_D_t = (mod_D_t + (D_max * (N_t % KD)) % D) % D;
            mod_D = (mod_D + (mod_D_t + D - ((N_t / KD) * RD) % D) % D) % D;
        }
        else
            mod_D = (mod_D + N_t * R_t) % D;
    }
    return mod_D;
}