Алгоритм деления 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.$

Численный эксперимент показал, что достаточно одной вышеописанной итерации. Физически это объясняется тем, что алгоритм разработан специально для , что приводит к достаточно большому числу в знаменателе и сравнительно небольшому частному. Однако, для получения окончательного результата требуется "fine tuning" - последовательное прибавление или вычитание единицы в зависимости от текущей ошибки деления, что можно сделать в цикле и за одно получить остаток от деления. Для этого необходима реализация оператора умножения "широкого" числа на "узкое", при этом дополнительно следует отслеживать переполнение, в результате которого придется уменьшить частное на единицу.

Заметим, что предложенный алгоритм не привязан к типу чисел: целое, с плавающей точкой, в альтернативной системе счисления и т. п. Основа алгоритма предполагает возможность представления произвольного числа в виде двух половинок: старшей и младшей. Знаковые биты обрабатываются отдельной логикой, так, чтобы фактически обрабатывались числа без знака.

Пример реализации алгоритма деления на С++

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

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

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

$Z = {{X} \over {C}} = {{ A \cdot {2}^{w} + B } \over {С}}, ~~C>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-битных встроенных чисел.