将一个大整数缩放一倍

Scaling a big integer by a double

本文关键字：缩放 一倍 整数 一个      更新时间：2023-10-16

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

（M*因子）mod M

其中M是一个大整数，因子为一个二重。我不会使用任何库，除非你想把头文件中的十几行代码称为"库"；因此，大浮动数学在这里不是一个选择。

Knuth和GMP/MPIR源代码没有答案，在这里我只找到了大整数和二重之间的乘法，这并不适用，因为第二个答案太奇特，第一个答案失去了太多精度。

根据第一性原理，用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++引用中选择了函数。

澄清：从乘法开始的所有运算都将使用（偏）大整数数学；在这里，我只是使用limb_t来实现这个目的，以便获得一个可以发布并实际运行的小玩具程序。最后一个应用程序将使用GMP的mpn_mul_1（）和mpn_rshift（）的道德等价物。

浮点数只不过是三项的乘积。这三项是符号、有效位（有时称为尾数）和指数。这三项的值计算为

（-1）符号*有效位*基数指数

基数通常是2，尽管C++标准不能保证这一点。相应地，您的计算变成

（M*因子）modM

==（M*（-1）符号*有效位*基数指数）modM

==（（-1）符号（M）+符号*abs（M*有效位*基数指数）modM

计算结果的符号应该相当琐碎。计算X*base指数是相当直接的：如果base是2，则它是合适的位移位，或者是与base的幂相乘/除法（对于正幂，左移或乘法，对于负指数则右移或除法）。假设你的大整数表示已经支持模运算，唯一有趣的项是abs（M）*有效位的乘积，但这只是一个正常的整数乘法，尽管对于你的大整型表示。我没有仔细检查过，但我认为这是你联系到的第一个答案（你形容为"太异国情调"的答案）。

剩余位用于从double中提取符号、有效位和指数。符号可以通过与0.0的比较很容易地确定，有效位和指数可以使用frexp()获得（例如，请参见此答案）。有效位作为double返回，即您可能希望将其乘以2std::numeric_limits<double>::digits并适当调整指数（我已经有一段时间没有这样做了，即我不完全确定frexp()的外部契约）。

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

相关文章：