Алгоритм деления 2W-разрядных чисел с использованием операций с числами разрядностью W
На примере 32-битных целых чисел рассматривается масштабируемый алгоритм деления, использующий числа с двукратно меньшей (16 бит) разрядностью. Для иллюстрации работоспособности алгоритма приведен код тестового приложения на языке С++.
Описание алгоритма
Представим некоторое
где
Тогда деление двух чисел,
Заметим, что если
Считаем далее, что
Здесь
Перепишем дробь:
Пренебрежем слагаемым
Выделим слагаемое
Сделаем замену переменных
В числителе и знаменателе дроби стоят "широкие" числа (разрядностью
Таким образом, "широкий" числитель последовательно делим на "узкие" знаменатели. Первый знаменатель иногда может равняться множителю
Численный эксперимент показал, что достаточно одной вышеописанной итерации. Физически это объясняется тем, что алгоритм разработан специально для
Заметим, что предложенный алгоритм не привязан к типу чисел: целое, с плавающей точкой, в альтернативной системе счисления и т. п. Основа алгоритма предполагает возможность представления произвольного числа в виде двух половинок: старшей и младшей. Знаковые биты обрабатываются отдельной логикой, так, чтобы фактически обрабатывались числа без знака.
Пример реализации алгоритма деления на С++
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;
}
Приложение А Итеративный алгоритм деления "широкого" числа на "узкое"
Пусть имеем дробь
Представим число
где
Заменим
Далее вычисляем ошибку приближения:
которую подаем на следующий шаг итерации если она больше или равна знаменателю
Выводы
Предложен и протестирован алгоритм деления чисел, состоящих из старшей и младшей половинок, масштабирующийся на произвольную разрядность кратно
Предложенный алгоритм легко масштабируется на 128-битный вариант с использованием встроенной 64-битной арифметики. Однако, вариант с масштабированием, например, на 256-бит, требует реализации в структуре U128 полноценного умножения, что можно сделать, масштабируя реализованный оператор "половинчатого" умножения: U32 на u16. Также потребуется реализация оператора побитового сдвига. В конечном итоге, при реализации всех необходимых операторов можно реализовать шаблонную структуру с произвольной разрядностью