2

我正在用 gmp 管理一些大（128~256 位）整数。如果我想将它们乘以接近 1的双精度（0.1 < 双精度 < 10），结果仍然是一个近似整数。我需要做的一个很好的操作示例如下：

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

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

4

3 回答 3

0

这就是你想要的：

// 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（最低有效字节）函数不处理符号，但它只有一个，如果和一些要添加的 xor/inc。

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

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

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

于 2013-08-04T17:50:08.433 回答
0

我对 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 */

现在下一步取决于您希望的舍入类型。如果您希望乘以dbyi得到一个向 -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：我通过随机浏览找到了您的问题，但是如果您将其标记为“浮点数”，那么比我更有能力的人可以快速回答。

于 2013-08-05T07:54:41.847 回答
0

试试这个策略：

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

  4. 将结果转换为整数

    mpf_set_z(...)
mpf_set_d(...)
mpf_mul(...)
mpz_set_f(...)
于 2018-01-25T08:34:50.510 回答