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

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

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

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

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

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

Считаем далее, что , иначе - результат деления равен нулю. Представим число в виде:

Здесь - целая часть от деления, а - остаток от деления на .

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

Пренебрежем слагаемым , чтобы получить нижнюю оценку результата деления плюс дополнительно упростить формулу:

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

Сделаем замену переменных Дело в том, что "тяжелые" случаи соответствуют параметру достаточно большому, поэтому введенный параметр "дельта" при этом будет мал и формула быстрее приведет к результату:

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

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

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

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

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

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

struct Quadrupole {
    u16 A;
    u16 B;
    u16 C;
    u16 D;
};

struct Signess {
    int s1;
    int s2;
};

static auto const seed = std::random_device{}();

/***
 * Генератор случайных чисел.
 */
auto rollu16 = [urbg = std::mt19937{seed},
                distr = std::uniform_int_distribution<u16>{}]() mutable -> u16 {
    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(u16) * CHAR_BIT) / 2;
    u16 mHigh = 0;
    u16 mLow = 0;
    int mSign = 0;
    int mOverflow = 0;

/**
 * Оператор проверки на ненулевое значение.
 */
bool operator()() {
    return (mLow != 0) || (mHigh != 0) || (mOverflow != 0);
}

U32 operator+(U32 rhs) const {
    U32 res{};
    U32 X = *this;
    if (X.mSign != 0 && rhs.mSign == 0) {
        X.mSign = 0;
        res = rhs - X;
        return res;
    }
    if (X.mSign == 0 && rhs.mSign != 0) {
        rhs.mSign = 0;
        res = X - rhs;
        return res;
    }
    res.mLow = X.mLow + rhs.mLow;
    const u16 c1 = res.mLow < std::min(X.mLow, rhs.mLow);
    res.mHigh = X.mHigh + rhs.mHigh;
    const int c2 = res.mHigh < std::min(X.mHigh, rhs.mHigh);
    auto tmp = res.mHigh;
    res.mHigh = tmp + c1;
    const int c3 = res.mHigh < std::min(tmp, c1);
    res.mOverflow = c2 || c3;
    if (X.mSign != 0 && rhs.mSign != 0) {
        res.mSign = 1;
    }
    return res;
}

U32 operator-(U32 rhs) const {
    U32 res{};
    U32 X = *this;
    if (mSign != 0 && rhs.mSign == 0) {
        rhs.mSign = 1;
        res = rhs + X;
        return res;
    }
    if (mSign == 0 && rhs.mSign != 0) {
        rhs.mSign = 0;
        res = X + rhs;
        return res;
    }
    if (mSign != 0 && rhs.mSign != 0) {
        rhs.mSign = 0;
        X.mSign = 0;
        res = rhs - X;
        return res;
    }
    res.mLow = X.mLow - rhs.mLow;
    res.mHigh = X.mHigh - rhs.mHigh;
    const int c1 = X.mLow < rhs.mLow;
    if (c1 != 0) {
        res.mLow -= 1u;
        res.mLow += 1u;
        res.mHigh -= 1u;
        const int c2 = X.mHigh < (rhs.mHigh + 1u);
        if (c2 != 0) {
            res.mHigh = -res.mHigh - 1u;
            res.mLow = -res.mLow;
            res.mSign = 1;
        }
    } else if (X.mHigh < rhs.mHigh) {
        res.mHigh = (res.mLow != 0) ? -res.mHigh - 1u : -res.mHigh;
        res.mLow = -res.mLow;
        res.mSign = 1;
    }
    return res;
}

U32 mult16(u16 x, u16 y) const {
    constexpr u16 mask = (1u << mHalfWidth) - 1;
    auto x_low = x & mask;
    auto y_low = y & mask;
    auto x_high = x >> mHalfWidth;
    auto y_high = y >> mHalfWidth;
    auto t1 = x_low * y_low;
    auto t = t1 >> mHalfWidth;
    auto t21 = x_low * y_high;
    auto q = t21 >> mHalfWidth;
    auto p = t21 & mask;
    auto t22 = x_high * y_low;
    auto s = t22 >> mHalfWidth;
    auto r = t22 & mask;
    auto t3 = x_high * y_high;
    U32 res{};
    res.mLow = t1;
    auto div = (q + s) + ((p + r + t) >> mHalfWidth);
    auto mod = (t21 << mHalfWidth) + (t22 << mHalfWidth);
    res.mLow += mod;
    res.mHigh += div;
    res.mHigh += t3;
    res.mOverflow = res.mHigh < t3 ? 1 : 0;
    return res;
}

U32 operator*(u16 rhs) {
    U32 tmpB = mult16(mLow, rhs);
    U32 tmpA = mult16(mHigh, rhs);
    tmpA.mHigh = tmpA.mLow;
    tmpA.mLow = 0;
    tmpB = tmpB + tmpA;
    if (this->mSign != 0) {
        tmpB.mSign = 1;
    }
    return tmpB;
}

U32 operator/(u16 y) const {
    U32 rem = *this;
    auto iter1 = [&rem](u16 v) {
        U32 res{};
        auto q0 = rem.mLow / v;
        auto r0 = rem.mLow % v;
        auto q1 = rem.mHigh / v;
        auto r1 = rem.mHigh % v;
        res.mLow = q0;
        res.mHigh = q1;
        rem.mLow = r0;
        rem.mHigh = r1;
        return res;
    };
    auto iter2 = [&rem](u16 v) {
        U32 res{};
        auto d = (1 u << mHalfWidth) / v;
        auto t = (1 u << mHalfWidth) % v;
        res.mLow = d * rem.mHigh << mHalfWidth;
        res.mLow += ((t * rem.mHigh << mHalfWidth) + rem.mLow) / v;
        rem.mLow = 0;
        rem.mHigh = 0;
        return res;
    };
    U32 result{};
    result = result + iter1(y);
    result = result + iter2(y);
    if (this->mSign != 0) {
        result.mSign = 1;
    }
    return result;
}

U32& operator/=(u16 y) {
    *this = *this / y;
    return *this;
}

U32 operator/(const U32 other) const {
    U32 X = *this;
    U32 Y = other;
    constexpr U32 ZERO{.mHigh = 0, .mLow = 0, .mSign = 0, .mOverflow = 0};
    constexpr U32 UNIT{.mHigh = 0, .mLow = 1, .mSign = 0, .mOverflow = 0};
    constexpr U32 UNIT_NEG{
        .mHigh = 0, .mLow = 1, .mSign = 1, .mOverflow = 0};
    if (mHigh < Y.mHigh) {
        return ZERO;
    }
    if (Y.mHigh == 0) {
        U32 res = X / Y.mLow;
        if (Y.mSign != 0) {
            res.mSign = res.mSign != 0 ? 0 : 1;
        }
        return res;
    }
    if (X.mSign != 0 && Y.mSign != 0) {
        X.mSign = 0;
        Y.mSign = 0;
    }
    bool make_sign_inverse = false;
    if ((X.mSign != 0 && Y.mSign == 0) || (X.mSign == 0 && Y.mSign != 0)) {
        X.mSign = 0;
        Y.mSign = 0;
        make_sign_inverse = true;
    }
    assert((Y.mHigh != 0) || (Y.mLow != 0));
    const u16 Q = mHigh / Y.mHigh;
    const u16 R = mHigh % Y.mHigh;
    const u16 Delta = (1u << 2 * mHalfWidth) - 1 - Y.mLow;
    const U32 DeltaQ = mult16(Delta, Q);
    U32 W1{};
    W1.mHigh = R >= Q ? R - Q : Q - R;
    W1.mLow = 0;
    W1.mSign = R >= Q ? 0 : 1;
    W1 = W1 + DeltaQ;
    u16 C1 = Y.mHigh + 1u;
    C1 = C1 != 0 ? C1 : (1 u << 2 * mHalfWidth) - 1;
    u16 W2 = (1 u << 2 * mHalfWidth) - 1 - 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.mSign == 0) ? X - N : X + N;
    U32 More = Error - Y;
    bool do_inc = More.mSign == 0;
    bool do_dec = Error.mSign == 1;
    while (do_dec || do_inc) {
        result = result + (do_inc ? UNIT : (do_dec ? UNIT_NEG : ZERO));
        if (do_dec) {
            Error = Error + Y;
        }
        if (do_inc) {
            Error = Error - Y;
        }
        More = Error - Y;
        do_inc = More.mSign == 0;
        do_dec = Error.mSign == 1;
    }
    if (make_sign_inverse) {
        result.mSign = result.mSign != 0 ? 0 : 1;
    }
    return result;
}

/**
 * Возвращает строковое представление числа.
 */
std::string value() const {
    std::string result{};
    if (this->mOverflow != 0) {
        result = "Overflow";
        return result;
    }
    U32 tmp = *this;
    constexpr int multiplier_mod10 =
        (1u << 2 * mHalfWidth) % 10;  // = 65536 mod 10.
    while (tmp()) {
        const int d =
            ((tmp.mLow % 10u) + multiplier_mod10 * (tmp.mHigh % 10u)) %
            10;
        result.push_back(DIGITS[d]);
        tmp /= 10u;
    }
    if (mSign != 0 && (mHigh != 0 || mLow != 0)) {
        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.mSign != 0 ? -x : x;
    y = z2.mSign != 0 ? -y : y;
    const i64 z_built_in = x / y;
    i64 z = ((u64)z3.mHigh << 16) + (u64)z3.mLow;
    z = z3.mSign != 0 ? -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> args = {
        {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}};
    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& arg : args) {
        is_ok &= make_test(arg);
    }
    assert(is_ok);
}

void test_division_randomly(long long N) {
    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 (true) {
        ++counter;
        const Quadrupole q{rollu16(), rollu16(), rollu16(), rollu16()};
        const Signess s{rollu16() % 2, rollu16() % 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 << '\n';
        }
        if (ext >= N) {
            break;
        }
    }
}

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";
    }
    {
        U32 x{.mHigh = 1, .mLow = 0};
        std::cout << x.value() << '\n';
    }
    {
        U32 x{.mHigh = 1, .mLow = 65535};
        std::cout << x.value() << '\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;
}

Выводы

Предложен и протестирован алгоритм деления чисел, состоящих из старшей и младшей половинок, масштабирующийся на произвольную разрядность кратно . Данный алгоритм в некотором смысле является вариантом "умного" деления в столбик: сначала вычисляется первое приближение, равное делению старших половинок числа, а затем - второе, равное скорректированному первому. Корректор равен последовательному делению некоторого широкого числа на два узких.

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