Умножение между большими целыми и двойными

Question

Я управляю некоторыми большими (128~256бит) целыми числами с помощью gmp. Настал момент, когда я хотел бы умножить их на двойное число, близкое к 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;
  }
}

Подскажите, пожалуйста, что вы думаете? У вас есть предложение сделать его более надежным или быстрым?

Answer 1

Я не очень знаком с 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: я нашел ваш вопрос по случайному просмотру, но если бы вы пометили его "с плавающей точкой", люди, более компетентные, чем я, могли бы ответить на него быстро.

Answer 2

Это то, что вы хотели:

// 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 на большое целое, большое int - это массив байтов _N, где первый байт - это MSB (старший байт), а последний - байт LSB (младший байт), функция не обрабатывает сигнатуру, но она только одна, если и немного xor/inc, чтобы добавить.

Беда в том, что double имеет низкую точность даже для вашего числа 1.23456789!!! из-за потери точности результат не является точным, как это должно быть (1234387129122386944 вместо 1234567890000000000) Я думаю, что мой код намного быстрее и даже точнее, чем ваш, потому что мне не нужно умножать числа mul / mod / div на 10, вместо этого я использую сдвиг бит там, где это возможно, и не на 10 цифр, а на 256 цифр (8 бит). если вам нужно больше точности, чем использовать длинную арифметику. Вы можете ускорить этот код, используя большие цифры (16,32, ... бит)

Моя длинная арифметика для точных астро-вычислений обычно состоит из чисел с фиксированной точкой 256,256 бит, состоящих из 2*8 DWORD + знак, но, конечно, намного медленнее, и некоторые гониометрические функции очень сложно реализовать, но если вам нужны только базовые функции, а не свой собственный код Арифметика не так сложна.

Кроме того, если вы хотите, чтобы числа часто читались в удобной для чтения форме, было бы неплохо пойти на компромисс между скоростью и размером и подумать не о том, чтобы использовать двоичные числа, а числа в формате BCD.

Answer 3

Попробуйте эту стратегию:

1. Преобразовать целочисленное значение в большое число
2. Конвертировать двойное значение в большое число
3. Сделать продукт

4. Преобразовать результат в целое число

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