Умножение целого числа на рациональное без промежуточного переполнения

У меня есть структура, представляющая неотрицательное рациональное число p/q:

struct rational {
    uint64_t p;
    uint64_t q; // invariant: always > 0
};

Я хотел бы умножить свое рациональное число на uint64 n и получить целочисленный результат, округленный в меньшую сторону. То есть я хотел бы рассчитать:

uint64_t m = (n * r.p)/r.q;

избегая при этом промежуточного переполнения в n * r.p. (Конечно, окончательный результат может выйти за пределы, что допустимо.)

Как я могу это сделать? Есть ли способ сделать это без высокого умножения?

(Я посмотрел на boost::rational, но, похоже, он не предоставляет эту функцию.)


c c++ integer-arithmetic rational-numbers
person ridiculous_fish    schedule 12.08.2016    source источник
comment
разве это не сработает с uint64_t m = (n / r.q) * r.p ?   -  person dangom    schedule 12.08.2016
comment
Вычислите рациональное число n / r.q и приведите его к наименьшей форме, а затем умножьте на r.p.   -  person Barmar    schedule 12.08.2016
comment
@DanielG, Бармар: Ни один из них не помогает, если p == n и p < q, но n * p > q.   -  person lorro    schedule 12.08.2016
comment
Эта проблема меня раздражает, это так просто на ассемблере x64, но нет разумного способа записать это на C.   -  person harold    schedule 12.08.2016
comment
Связано: codeproject.com/Articles/9670/Scaling-bit-integers   -  person Ben Voigt    schedule 12.08.2016
comment
@BenVoigt только что понял это. Спасибо за ссылку   -  person dangom    schedule 12.08.2016
comment
@harold Большинство компиляторов теперь имеют встроенную функцию с высоким множителем. Но не для разделения двойной ширины.   -  person Mysticial    schedule 12.08.2016

Ответы (2)


arrow_upward
1
arrow_downward

Вы можете использовать крестьянское умножение:

// requires 0 < c < 2^63
// returns (a * b) / c.
uint64_t multDiv(uint64_t a, uint64_t b, uint64_t c) {
  uint64_t rem = 0;
  uint64_t res = (a / c) * b;
  a = a % c;
  // invariant: a_orig * b_orig = (res * c + rem) + a * b
  // a < c, rem < c.
  while (b != 0) {
    if (b & 1) {
      rem += a;
      if (rem >= c) {
        rem -= c;
        res++;
      }
    }
    b /= 2;
    a *= 2;
    if (a >= c) {
      a -= c;
      res += b;
    }
  }
  return res;
}
person Ha.    schedule 12.08.2016

arrow_upward
0
arrow_downward

Либо 128-бит, и вы можете использовать здесь умножение Карацубы; или вы можете использовать китайскую теорему об остатках для представления (n * r.p) по модулю p1, а также по модулю p2. Второй может быть медленнее.

person lorro    schedule 12.08.2016