大整数和双精度之间的乘法 [英] Multiplication between big integers and doubles

问题描述

我使用gmp管理一些大（128〜256位）整数。有一点是我想将它们乘以两个接近1 （0.1

I am managing some big (128~256bits) integers with gmp. It has come a point were I would like to multiply them for a double close to 1 (0.1 < double < 10), the result being still an approximated integer. A good example of the operation I need to do is the following:

int i = 1000000000000000000 * 1.23456789

我在gmp文档中搜索，但我没有找到一个函数，所以我最后写这个代码，似乎工作良好：

I searched in the gmp documentation but I didn't find a function for this, so I ended up writing this code which seems to work well:

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;
  }
}

你能告诉我你觉得怎么样？你有什么建议让它更强壮或更快？

Can you please tell me what do you think? Do you have any suggestion to make it more robust or faster?

推荐答案

我不是那么熟悉C ++或GMP什么可能会建议源代码没有语法错误，但你所做的是比它应该更复杂，可以引入不必要的近似。

I am not so familiar with either C++ or GMP what I could suggest source code without syntax errors, but what you are doing is more complicated than it should and can introduce unnecessary approximation.

相反，我建议你写函数 mpz_mult_d（）如下：

Instead, I suggest you write function mpz_mult_d() like this:

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舍入的结果，只需返回商作为函数的结果。如果您希望将结果四舍五入为最接近的整数，则必须查看 remaining ：

And now the next step depends on the kind of rounding you wish. If you wish the multiplication of d by i to give a result rounded toward -inf, just return quotient as result of the function. If you wish a result rounded to the nearest integer, you must look at 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：我通过随机浏览找到你的问题，但如果你标记为浮点，比我更能胜任的人可以很快地回答。

---

PS: I found your question by random browsing but if you had tagged it "floating-point", people more competent than me could have answered it quickly.

这篇关于大整数和双精度之间的乘法的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！