将大整数按双倍缩放

Question

我需要将大整数（几百位）放大一倍。特别是，我需要计算

(M * 因子) mod M

其中 M 是一个大整数，factor 是一个双精度数。除非您想将头文件中的十几行代码称为“库”，否则我不会使用任何库；因此大浮点数学在这里不是一个选项。

Knuth 和 GMP/MPIR 源代码没有答案，在这里我发现只有 Multiplication between big integers and doubles 并不真正适用，因为第二个答案太奇特而且第一个失去了太多精度。

从第一原理开始并使用 uint64_t 模拟大整数我想出了这个（可使用 64 位 VC++ 或 gcc/MinGW64 运行）：

#include <cassert>
#include <cfloat>
#include <climits>
#include <cmath>
#include <cstdint>
#include <cstdio>

#include <intrin.h>   // VC++, MinGW

#define PX(format,expression)  std::printf("\n%35s  ==  " format, #expression, expression);

typedef std::uint64_t limb_t;

// precision will be the lower of LIMB_BITS and DBL_MANT_DIG
enum {  LIMB_BITS = sizeof(limb_t) * CHAR_BIT  };

// simulate (M * factor) mod M with a 'big integer' M consisting of a single limb
void test_mod_mul (limb_t modulus, double factor)
{
   assert( factor >= 0 );

   // extract the fractional part of the factor and discard the integer portion

   double ignored_integer_part;
   double fraction = std::modf(factor, &ignored_integer_part);

   // extract the significand (aligned at the upper end of the limb) and the exponent

   int exponent;
   limb_t significand = limb_t(std::ldexp(std::frexp(fraction, &exponent), LIMB_BITS));

   // multiply modulus and single-limb significand; the product will have (n + 1) limbs

   limb_t hi;
/* limb_t lo = */_umul128(modulus, significand, &hi);

   // The result comprises at most n upper limbs of the product; the lowest limb will be 
   // discarded in any case, and potentially more. Factors >= 1 could be handled as well,
   // by dropping the modf() and handling exponents > 0 via left shift.

   limb_t result = hi;

   if (exponent)
   {
      assert( exponent < 0 );

      result >>= -exponent;
   }

   PX("%014llX", result);
   PX("%014llX", limb_t(double(modulus) * fraction));
}

int main ()
{
   limb_t const M = 0x123456789ABCDEull;  // <= 53 bits (for checking with doubles)

   test_mod_mul(M, 1 - DBL_EPSILON);
   test_mod_mul(M, 0.005615234375);
   test_mod_mul(M, 9.005615234375);
   test_mod_mul(M, std::ldexp(1, -16));
   test_mod_mul(M, std::ldexp(1, -32));
   test_mod_mul(M, std::ldexp(1, -52));
}

乘法和移位将在我的应用程序中使用大整数数学完成，但原理应该相同。

基本方法是正确的还是仅因为我在这里使用玩具整数进行测试才有效？我对浮点数学一无所知，我从C++ reference 中挑选了函数。

澄清：从乘法开始的所有事情都将使用（部分）大整数数学来完成；在这里，我只使用limb_t 来获得一个可以发布并实际运行的小型玩具程序。最终的应用程序将使用 GMP 的 mpn_mul_1() 和 mpn_rshift() 的道德等效项。

Answer 1

浮点数不过是三项的乘积。 sign、significand（有时称为尾数）和指数这三个项。这三个项的值计算为

(-1)^符号 * significand * base^指数

基数通常为 2，尽管 C++ 标准不保证这一点。相应地，你的计算就变成了

(M * 因子) mod M

== (M * (-1)^sign * significand * base^指数) mod M

== ((-1)^{sign(M) + sign} * abs(M) * significand * base^exponent) mod M

计算结果的符号应该相当简单。计算 X * base^exponent 相当简单：如果基数为 2，它要么是适当的位移位，要么是与 / 的乘法除以基数（正指数左移或乘法，负指数右移或除法）。假设你的大整数表示已经支持模运算，唯一有趣的术语是 abs(M) * significand 的乘法，但这只是一个普通的整数乘法，尽管对于你的大整数表示。我没有仔细检查，但我认为这是你链接到的第一个答案（你描述为“太异国情调”的答案）。

剩下的部分是从double 中提取 sign、significand 和 exponent。通过与0.0 的比较可以轻松确定符号，并且可以使用frexp() 获得significand 和exponent（例如参见this answer）。 significand 以 double 的形式返回，也就是说，您可能希望将其乘以 2^{std::numeric_limits<double>::digits} 并适当调整指数（我在一段时间，即我不完全确定frexp()的确切合同。

回答您的问题：我不熟悉您使用的 GMP 操作，但我认为您所做的操作确实执行了上述计算。