大整数和双精度数之间的乘法答案

【问题标题】：Multiplication between big integers and doubles大整数和双精度数之间的乘法
【发布时间】：2013-05-27 12:39:13
【问题描述】：

我正在使用 gmp 管理一些大（128~256 位）整数。如果我想将它们相乘成一个 接近 1 (0.1

int i = 1000000000000000000 * 1.23456789

我在 gmp 文档中进行了搜索，但没有找到用于此功能的函数，因此我最终编写了这段似乎运行良好的代码：

mpz_mult_d(mpz_class & r, const mpz_class & i, double d, int prec=10) {
  if (prec > 15) prec=15; //avoids overflows
  uint_fast64_t m = (uint_fast64_t) floor(d);
  r = i * m;
  uint_fast64_t pos=1;
  for (uint_fast8_t j=0; j<prec; j++) {
    const double posd = (double) pos;
    m = ((uint_fast64_t) floor(d * posd * 10.)) -
        ((uint_fast64_t) floor(d * posd)) * 10;
    pos*=10;
    r += (i * m) /pos;
  }
}

你能告诉我你的想法吗？你有什么建议让它更健壮或更快？

【问题讨论】：

这是Code Review 的问题，不是 StackOverflow 的问题 :)
抱歉，我不知道那个分支。然而，这只是我对一个非常精确的问题的解决方案。请考虑回答一般问题，并最终仅对代码发表评论。谢谢！
如果你需要一个近似值，你为什么不把大整数转换成双精度数呢？
GMP 支持任意精度的有理数和浮点数。将您的两个值转换为其中一种格式，然后让 GMP 进行乘法运算。
您应该查看 MPFR (mpfr.org)，它相当于浮点数的 GMP。它与 GMP 一样容易/简单。

标签： c++ c++11 integer double gmp

【解决方案1】：

这就是你想要的：

// BYTE lint[_N]   ... lint[0]=MSB, lint[_N-1]=LSB
void mul(BYTE *c,BYTE *a,double b)  // c[_N]=a[_N]*b
    {
    int i; DWORD cc;
    double q[_N+1],aa,bb;
    for (q[0]=0.0,i=0;i<_N;)        // mul,carry down
        {
        bb=double(a[i])*b; aa=floor(bb); bb-=aa;
        q[i]+=aa; i++;
        q[i]=bb*256.0;
        }
    cc=0; if (q[_N]>127.0) cc=1.0;  // round
    for (i=_N-1;i>=0;i--)           // carry up
        {
        double aa,bb;
        cc+=q[i];
        c[i]=cc&255;
        cc>>=8;
        }
    }

_N 是每个大整数的位数/8，大整数是 _N 个字节的数组，其中第一个字节是 MSB（最高有效字节），最后一个字节是 LSB（最低有效字节）函数不处理符号，但它只是一个 if 和一些要添加的 xor/inc。

问题是即使对于您的号码 1.23456789，double 的精度也很低！！！由于精度损失，结果并不准确（1234387129122386944 而不是 1234567890000000000）我认为我的代码比你的更快，甚至更精确，因为我不需要将 mul/mod/div 数字乘以 10，而是我使用在可能的地方进行位移，而不是位移 10 位，而是位移 256 位（8 位）。如果您需要比使用长算术更高的精度。您可以通过使用更大的数字（16,32，...位）来加速此代码

我用于精确 astro 计算的长算法通常是定点 256.256 位数字，由 2*8 DWORDs + 符号组成，但当然要慢得多，而且一些测角函数很难实现，但如果你只想要基本函数而不是编写自己的 lon 算法并不难。

如果您希望数字经常以可读的形式出现，最好在速度/大小之间进行折衷，并考虑不使用二进制编码的数字，而是使用 BCD 编码的数字

【讨论】：

【解决方案2】：

我对 C++ 或 GMP 都不太熟悉，我可以在没有语法错误的情况下建议源代码，但你所做的比它应该做的更复杂，并且可能引入不必要的近似。

相反，我建议你像这样编写函数mpz_mult_d()：

mpz_mult_d(mpz_class & r, const mpz_class & i, double d) {
  d = ldexp(d, 52); /* exact, no overflow because 1 <= d <= 10 */
  unsigned long long l = d; /* exact because d is an integer */
  p = l * i; /* exact, in GMP */
  (quotient, remainder) = p / 2^52; /* in GMP */

现在下一步取决于您希望的舍入类型。如果您希望将d 与i 相乘以得到向-inf 舍入的结果，只需返回quotient 作为函数的结果。如果您希望将结果四舍五入到最接近的整数，您必须查看remainder：

 assert(0 <= remainder);  /* proper Euclidean division */
 assert(remainder < 2^52);
 if (remainder < 2^51) return quotient;
 if (remainder > 2^51) return quotient + 1; /* in GMP */
 if (remainder == 2^51) return quotient + (quotient & 1); /* in GMP, round to “even” */

PS：我通过随机浏览找到了您的问题，但如果您将其标记为“浮点”，那么比我更有能力的人可以快速回答。

【讨论】：

【解决方案3】：

试试这个策略：

将整数值转换为大浮点数
将双精度值转换为大浮点数
做产品

将结果转换为整数

mpf_set_z(...)
mpf_set_d(...)
mpf_mul(...)
mpz_set_f(...)

【讨论】：