2

gmp でいくつかの大きな (128 ~ 256 ビット) 整数を管理しています。1 に近いdouble (0.1 < double < 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

私は構文エラーのないソース コードを提案できる 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 */

次のステップは、希望する丸めの種類によって異なります。dbyの乗算で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: たまたまブラウジングしてあなたの質問を見つけましたが、「floating-point」というタグを付けていれば、私より有能な人がすぐに回答できたはずです。

于 2013-08-05T07:54:41.847 に答える