Алгоритм деления 2W-разрядных чисел с использованием операций с числами разрядностью W

На примере 32-битных целых чисел рассматривается масштабируемый алгоритм деления, использующий числа с двукратно меньшей (16 бит) разрядностью. Для иллюстрации работоспособности алгоритма приведен код тестового приложения на языке С++.

Описание алгоритма

Представим некоторое -битное число в виде:

$X = (A, B) = A \cdot {2}^{W} + B,$

где , — старшая и младшая части -битных компонент числа, соответственно.

Тогда деление двух чисел, и , можно записать в виде дроби:

$Z = {{X} \over {Y}} = {{A \cdot {2}^{W} + B} \over {C \cdot {2}^{W} + D}}.$

Заметим, что если -битное число. Ограничимся этим случаем. Для в примере ниже на С++ реализован итеративный алгоритм деления «широкого» числа на «узкое», см. Приложение А.

Считаем далее, что $A \geq C$ , иначе — результат деления равен нулю. Представим число в виде:

$A = \left \lfloor {{A} \over {C}} \right \rfloor \cdot C + A \mod C = q \cdot C + r.$

Здесь $q = \left \lfloor {{A} \over {C}} \right \rfloor$ — целая часть от деления, а $r = A \mod C$ — остаток от деления на .

Перепишем дробь:

$Z = {{q \cdot C \cdot {2}^{W} + B + r \cdot {2}^{W} } \over {C \cdot {2}^{W} + D}}.$

Пренебрежем слагаемым , чтобы упростить формулу:

$Z_1 = {(q \cdot C + r)}{{{2}^{W}} \over {C \cdot {2}^{W} + D}} \leq Z.$

Выделим слагаемое в качестве главной компоненты:

$Z_1 = q + {{r \cdot {2}^{W} - D \cdot q} \over {C \cdot {2}^{W} + D}}.$

Сделаем замену переменных $\Delta=2^W-D.$ Дело в том, что «тяжелые» случаи соответствуют параметру достаточно большому, поэтому введенный параметр «дельта» при этом будет мал и формула быстрее приведет к результату:

$Z_1 = q + {{(r-q) \cdot {2}^{W} + \Delta \cdot q} \over {(C+1) \cdot {2}^{W} - \Delta}}.$

В числителе и знаменателе дроби стоят «широкие» числа (разрядностью ). Так как допускается использовать алгоритм деления широкого числа на узкое, то добьемся «узости» числа в знаменателе, вынеся множитель за скобки и выполняя деление последовательно:

$Z_1 = q + { { {(r-q) \cdot {2}^{W} + \Delta \cdot q} \over {{2}^{W} - {{\Delta} \over {C+1}}} } \over {C + 1} }.$

Таким образом, «широкий» числитель последовательно делим на «узкие» знаменатели. Первый знаменатель иногда может равняться множителю , что необходимо отслеживать в алгоритме: в регистрах ЦПУ число при этом обнулится и возникнет исключение. Альтернативно, можно заранее вычесть единицу и не заботиться о граничных условиях, так как для данного алгоритма $\Delta$ , что в итоге даст окончательный вариант:

$\Delta = 2^W - 1 - D,$ $Z_1 = q + {{ {(r-q) \cdot {2}^{W} + \Delta \cdot q} \over {{2}^{W} - 1 - {{\Delta} \over {C+1}}} } \over {C + 1} },$ $q = \left \lfloor {{A} \over {C}} \right \rfloor,$ $r = A \mod C.$

Численный эксперимент показал, что достаточно одной вышеописанной итерации. Физически это объясняется тем, что алгоритм разработан специально для

Hidden text

#include <limits.h> #include <algorithm> #include <cassert> #include <cmath> #include <cstdint> #include <iostream> #include <list> #include <random> #include <string>  using u64 = uint64_t; using i64 = int64_t; using u32 = uint32_t; using u16 = uint16_t;  using ULOW = u16; // Тип половинок: старшей и младшей частей составного числа.  static_assert(CHAR_BIT == 8);  struct Quadrupole { // Структура для задания дроби (A*M + B) / (C*M + D).                     // M - множитель системы счисления, 2^W, W - битовая ширина половинок.     ULOW A;     ULOW B;     ULOW C;     ULOW D; };  struct Signess { // Структура для задания знаков двух чисел.     int s1;     int s2; };  static auto const seed = std::random_device{}();  /***  * Генератор случайных чисел.  */ auto roll_ulow = [urbg = std::mt19937{seed},                 distr = std::uniform_int_distribution<ULOW>{}]() mutable -> ULOW {     return distr(urbg); };  static constexpr char DIGITS[10]{'0', '1', '2', '3', '4',                                  '5', '6', '7', '8', '9'};  // High/Low структура 32-битного числа со знаком и флагом переполнения. // Для иллюстрации алгоритма деления двух U32 чисел реализованы основные // арифметические операторы, кроме умножения двух U32 чисел. // Если дополнительно реализовать полноценное умножение и оператор побитового сдвига,  // то данная структура позволяет себя масштабировать: реализовывать 128-битные числа,  // 256-битные и т.д. struct U32 {     // Битовая полуширина половинок.     static constexpr int mHalfWidth = (sizeof(ULOW) * CHAR_BIT) / 2;     // Наибольшее значение половинок, M-1.     static constexpr ULOW mMaxULOW = ULOW(-1);     ULOW mHigh = 0;     ULOW mLow = 0;     int mSign = 0;     int mOverflow = 0;  bool is_negative() const {     return mSign != 0; }  bool is_overflow() const {     return mOverflow != 0; }  bool is_zero() const {     return (mLow | mHigh | mOverflow) == 0; }  constexpr static U32 get_zero() {     return U32 {.mHigh = 0, .mLow = 0, .mSign = 0, .mOverflow = 0}; }  constexpr static U32 get_unit() {     return U32 {.mHigh = 0, .mLow = 1, .mSign = 0, .mOverflow = 0}; }  constexpr static U32 get_unit_neg() {     return U32 {.mHigh = 0, .mLow = 1, .mSign = 1, .mOverflow = 0}; }  U32 operator+(U32 rhs) const {     U32 result{};     U32 X = *this;     if (X.is_negative() && !rhs.is_negative()) {         X.mSign = 0;         result = rhs - X;         return result;     }     if (!X.is_negative() && rhs.is_negative()) {         rhs.mSign = 0;         result = X - rhs;         return result;     }     result.mLow = X.mLow + rhs.mLow;     const ULOW c1 = result.mLow < std::min(X.mLow, rhs.mLow);     result.mHigh = X.mHigh + rhs.mHigh;     const int c2 = result.mHigh < std::min(X.mHigh, rhs.mHigh);     ULOW tmp = result.mHigh;     result.mHigh = tmp + c1;     const int c3 = result.mHigh < std::min(tmp, c1);     result.mOverflow = c2 || c3;     if (X.mSign && rhs.mSign) {         result.mSign = 1;     }     return result; }  U32& operator+=(U32 other) {     *this = *this + other;     return *this; }  U32 operator-(U32 rhs) const {     U32 result{};     U32 X = *this;     if (X.is_negative() && !rhs.is_negative()) {         rhs.mSign = 1;         result = rhs + X;         return result;     }     if (!X.is_negative() && rhs.is_negative()) {         rhs.mSign = 0;         result = X + rhs;         return result;     }     if (X.is_negative() && rhs.is_negative()) {         rhs.mSign = 0;         X.mSign = 0;         result = rhs - X;         return result;     }     if (X.is_zero()) {         result = rhs;         result.mSign = rhs.mSign ^ 1;         return result;     }     result.mLow = X.mLow - rhs.mLow;     result.mHigh = X.mHigh - rhs.mHigh;     const bool borrow = X.mLow < rhs.mLow;     const bool hasUnit = X.mHigh > rhs.mHigh;     if (borrow && hasUnit) {         result.mHigh -= ULOW(1);     }     if (borrow && !hasUnit) {         result = rhs - X;         result.mSign ^= 1;         return result;     }     if (!borrow && X.mHigh < rhs.mHigh) {         result.mHigh = -result.mHigh - ULOW(result.mLow != 0);         result.mLow = -result.mLow;         result.mSign = 1;     }     return result; }  U32& operator-=(U32 other) {     *this = *this - other;     return *this; }  U32 mult16(ULOW x, ULOW y) const {     constexpr ULOW MASK = (ULOW(1) << mHalfWidth) - 1;     const ULOW x_low = x & MASK;     const ULOW y_low = y & MASK;     const ULOW x_high = x >> mHalfWidth;     const ULOW y_high = y >> mHalfWidth;     const ULOW t1 = x_low * y_low;     const ULOW t = t1 >> mHalfWidth;     const ULOW t21 = x_low * y_high;     const ULOW q = t21 >> mHalfWidth;     const ULOW p = t21 & MASK;     const ULOW t22 = x_high * y_low;     const ULOW s = t22 >> mHalfWidth;     const ULOW r = t22 & MASK;     const ULOW t3 = x_high * y_high;     U32 result{};     result.mLow = t1;     const ULOW div = (q + s) + ((p + r + t) >> mHalfWidth);     const ULOW mod = (t21 << mHalfWidth) + (t22 << mHalfWidth);     result.mLow += mod;     result.mHigh += div;     result.mHigh += t3;     result.mOverflow = result.mHigh < t3 ? 1 : 0;     return result; }  U32 operator*(ULOW rhs) const {     U32 result = mult16(mLow, rhs);     U32 tmp = mult16(mHigh, rhs);     tmp.mHigh = tmp.mLow;     tmp.mLow = 0;     result += tmp;     result.mSign = this->mSign;     return result; }  U32 operator/(ULOW y) const {     U32 X = *this;     const int sign = X.mSign;     X.mSign = 0;     ULOW Q = X.mHigh / y;     ULOW R = X.mHigh % y;     ULOW N = R * (mMaxULOW / y) + (X.mLow / y);     U32 result{.mHigh = Q, .mLow = N};     U32 E = X - result * y;     while (E.mHigh != 0 || E.mLow >= y) {         Q = E.mHigh / y;         R = E.mHigh % y;         N = R * (mMaxULOW / y) + (E.mLow / y);         const U32 tmp {.mHigh = Q, .mLow = N};         result += tmp;         E -= tmp * y;     }     result.mSign = sign;     return result; }  U32& operator/=(ULOW y) {     *this = *this / y;     return *this; }  U32 operator/(const U32 other) const {     U32 X = *this;     U32 Y = other;     constexpr U32 ZERO = get_zero();     constexpr U32 UNIT = get_unit();     constexpr U32 UNIT_NEG = get_unit_neg();     if (X.mHigh < Y.mHigh) {         return ZERO;     }     if (Y.mHigh == 0) {         U32 result = X / Y.mLow;         result.mSign ^= (Y.is_negative());         return result;     }     const int make_sign_inverse = (X.mSign != Y.mSign);     X.mSign = 0;     Y.mSign = 0;     const ULOW Q = X.mHigh / Y.mHigh;     const ULOW R = X.mHigh % Y.mHigh;     const ULOW Delta = mMaxULOW - Y.mLow;     const U32 DeltaQ = mult16(Delta, Q);     U32 W1 = U32{.mHigh = R, .mLow = 0} - U32{.mHigh = Q, .mLow = 0};     W1 = W1 + DeltaQ;     ULOW C1 = (Y.mHigh < mMaxULOW) ? Y.mHigh + 1u : mMaxULOW;     ULOW W2 = mMaxULOW - Delta / C1;     U32 Quotient = W1 / W2;     Quotient = Quotient / C1;     U32 result = U32{.mHigh = 0, .mLow = Q} + Quotient;     U32 N = Y * result.mLow;     if (N.mOverflow != 0) {         result.mLow -= 1;         N = Y * result.mLow;     }     U32 Error = (result.is_negative()) ? X + N : X - N;     U32 More = Error - Y;     bool do_inc = !More.is_negative();     bool do_dec = Error.is_negative();     while (do_dec || do_inc) {         result += (do_inc ? UNIT : (do_dec ? UNIT_NEG : ZERO));         if (do_dec) {             Error += Y;         }         if (do_inc) {             Error -= Y;         }         More = Error - Y;         do_inc = !More.is_negative();         do_dec = Error.is_negative();     }     result.mSign ^= make_sign_inverse;     return result; }  /**  * Возвращает строковое представление числа.  */ std::string value() const {     std::string result{};     if (this->is_overflow()) {         result = "Overflow";         return result;     }     U32 X = *this;     constexpr int multiplier_mod10 = mMaxULOW % 10 + 1;     while (!X.is_zero()) {         const int d =             ((X.mLow % 10) + multiplier_mod10 * (X.mHigh % 10)) % 10;         result.push_back(DIGITS[d]);         X /= 10;     }     if (this->is_negative() && !this->is_zero()) {         result.push_back('-');     }     std::reverse(result.begin(), result.end());     return result.length() != 0 ? result : "0"; } };  bool test_div(U32 z1, U32 z2) {     const U32 z3 = z1 / z2;     i64 x = ((u64)z1.mHigh << 16) + (u64)z1.mLow;     i64 y = ((u64)z2.mHigh << 16) + (u64)z2.mLow;     x = z1.is_negative() ? -x : x;     y = z2.is_negative() ? -y : y;     const i64 z_built_in = x / y;     i64 z = ((u64)z3.mHigh << 16) + (u64)z3.mLow;     z = z3.is_negative() ? -z : z;     bool is_ok = z_built_in == z;     if (!is_ok) {         std::cout << "Assertion:\n";         std::cout << "Built-in: " << x << " / " << y << " = " << z_built_in                   << '\n';         std::cout << "This algorithm: " << x << " / " << y << " = "                   << z3.value() << '\n';         std::cout << "Parameters: ABCD, s1, s2: " << z1.mHigh << ", " << z1.mLow                   << ", " << z2.mHigh << ", " << z2.mLow << ", " << z1.mSign                   << ", " << z2.mSign << '\n';     }     return is_ok; }  void test_division_quick() {     const std::list<Quadrupole> my_cases = {         {51774, 28457, 50018, 10280}, {28792, 5507, 37, 64804},         {65258, 18362, 87, 35198},    {65526, 63280, 198, 52129},         {56139, 10364, 39, 36881},    {65498, 60804, 204, 20825},         {58092, 52199, 1, 57003},     {65498, 60804, 204, 20825},         {64666, 34598, 1, 60805},     {30903, 7652, 143, 48035},         {30161, 40182, 3351, 26310},  {40824, 35384, 13, 49151},         {60215, 18033, 165, 58003},   {42499, 42189, 4, 58879},         {16384, 16384, 0, 1},         {16384, 16384, 1, 0},         {16384, 16384, 1, 65535},     {16384, 16384, 1, 1},         {16384, 16384, 1, 65535},     {16384, 16384, 65535, 1},         {16384, 16384, 1, 65535},     {16384, 16384, 65535, 65535},         {16384, 16384, 65535, 65535}, {16384, 16384, 65535, 65535}         };     auto make_test = [](const Quadrupole q) -> bool {         return test_div({.mHigh = q.A, .mLow = q.B},                         {.mHigh = q.C, .mLow = q.D});     };     bool is_ok = true;     for (const auto& element : my_cases) {         is_ok &= make_test(element);     }     assert(is_ok); }  void test_division_randomly(long long N) {     if (N < 1) {         std::cout << "Skipped!\n";         return;     }     auto make_test = [](const Quadrupole q, const Signess s) -> bool {         return test_div(U32{.mHigh = q.A, .mLow = q.B, .mSign = s.s1},                         U32{.mHigh = q.C, .mLow = q.D, .mSign = s.s2});     };     long long counter = 0;     long long ext = 0;     bool is_ok = true;     while (ext < N) {         ++counter;         const Quadrupole q{roll_ulow(), roll_ulow(), roll_ulow(), roll_ulow()};         const Signess s{roll_ulow() % 2, roll_ulow() % 2};         if (q.C == 0 && q.D == 0) {             continue;         }         is_ok &= make_test(q, s);         assert(is_ok);         if (counter % (1ll << 27) == 0) {             ext++;             std::cout << "... iterations: " << counter << ". External: " <<                  ext << " from " << N << '\n';         }     } }  int main(int argc, char* argv[]) {     long long N = 10;     if (argc > 1) {         N = std::atoi(argv[1]);         std::cout << "You set " << N << " external iterations\n";     }     std::cout << "Test division quick\n";     test_division_quick();     std::cout << "Ok\n";     std::cout << "Test division randomly...\n";     test_division_randomly(N);     std::cout << "Ok\n";     return 0; }

Приложение А Итеративный алгоритм деления «широкого» числа на «узкое»

Пусть имеем дробь

в виде частного и остатка от деления на

$A = Q \cdot C + R,$

где — частное от деления на , а — остаток от такого деления. Тогда искомую дробь можно переписать в виде:

$Z = Q \cdot 2^w + {{R \cdot 2^w + B} \over {C}}.$

Заменим на в дроби и пренебрежем слагаемым , тогда первое приближение к результату можно записать в виде:

$Z = Q \cdot 2^w + R \cdot { \left \lfloor {{2^w-1} \over {C}} \right \rfloor} + \left \lfloor {{B} \over {C}} \right \rfloor.$

Далее вычисляем ошибку приближения:

$E = X - Z \cdot C,$

которую подаем на следующий шаг итерации если она больше или равна знаменателю .

Выводы

Предложен и протестирован алгоритм деления чисел, состоящих из старшей и младшей половинок, масштабирующийся на произвольную разрядность кратно ${2}^{N}$ . Данный алгоритм в некотором смысле является вариантом «умного» деления в столбик: сначала вычисляется первое приближение, равное делению старших половинок числа, а затем — второе, равное скорректированному первому. Корректор равен последовательному делению некоторого широкого числа на два узких.

Предложенный алгоритм легко масштабируется на 128-битный вариант с использованием встроенной 64-битной арифметики. Однако, вариант с масштабированием, например, на 256-бит, требует реализации в структуре U128 полноценного умножения, что можно сделать, масштабируя реализованный оператор «половинчатого» умножения: U32 на u16. Также потребуется реализация оператора побитового сдвига. В конечном итоге, при реализации всех необходимых операторов можно реализовать шаблонную структуру с произвольной разрядностью , рекурсивно (делением пополам) спадающую в арифметику 64-битных встроенных чисел.

ссылка на оригинал статьи https://habr.com/ru/articles/811323/

Алгоритм деления 2W-разрядных чисел с использованием операций с числами разрядностью W

Описание алгоритма

Приложение А Итеративный алгоритм деления «широкого» числа на «узкое»

Выводы

Комментарии

Добавить комментарий Отменить ответ