Тестирование алгоритма деления больших чисел на С++ с использованием Python C API13.05.2024 01:30

class PythonCaller {
public:

    PythonCaller() {
        Py_Initialize();
        mMain = PyImport_AddModule("__main__");
        mGlobalDictionary = PyModule_GetDict(mMain);
        mLocalDictionary = PyDict_New();
        std::cout << "Python was initialized!\n";
    }

    ~PythonCaller() {
        Py_Finalize();
        std::cout << "~Python was finalized!\n";
    }

    /**
     * Делит два 128-битных числа.
     * @param X Делимое.
     * @param Y Делитель.
     * @return Частное от деления, объект Python.
    */
    PyObject* Divide(U128 X, U128 Y) const {
        const char* pythonScript = "quotient = nominator // denominator\n";
        PyDict_SetItemString(mLocalDictionary, "nominator", PyLong_FromString(X.value().c_str(), nullptr, 10));
        PyDict_SetItemString(mLocalDictionary, "denominator", PyLong_FromString(Y.value().c_str(), nullptr, 10));
        PyRun_String(pythonScript, Py_file_input, mGlobalDictionary, mLocalDictionary);
        return PyDict_GetItemString(mLocalDictionary, "quotient");
    }
    /**
     * Сравнивает два частных от деления.
     * @param quotient Частное от деления Python, объект Python.
     * @param reference Частное от деления С++, строка.
    */
    bool Compare(PyObject* quotient, const char* reference) const {
        PyObject* repr = PyObject_Repr(quotient);
        PyObject* str = PyUnicode_AsEncodedString(repr, "utf-8", "~E~");
        const char *bytes = PyBytes_AsString(str);
        const bool is_ok = strcmp(bytes, reference) == 0;
        if (!is_ok) {
            printf("Python: %s\n", bytes);
            printf("C++: %s\n", reference);
        }
        Py_XDECREF(repr);
        Py_XDECREF(str);
        return is_ok;
    }
private:

    PyObject* mMain = nullptr;
    PyObject* mGlobalDictionary = nullptr;
    PyObject* mLocalDictionary = nullptr;
};

#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 

using u64 = uint64_t;
using ULOW = u64; // Тип половинок: старшей и младшей частей составного числа.

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 { // Структура для задания знаков двух чисел.
    bool s1;
    bool s2;
};

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

/***
 * Генератор случайных чисел.
 */
auto roll_ulow = [urbg = std::mt19937{seed},
                distr = std::uniform_int_distribution{}]() mutable -> ULOW {
    return distr(urbg);
};

auto roll_uint = [urbg = std::mt19937{seed},
                distr = std::uniform_int_distribution{}]() mutable -> uint {
    return distr(urbg);
};

auto roll_bool = [urbg = std::mt19937{seed},
                distr = std::uniform_int_distribution{}]() mutable -> bool {
    return distr(urbg) % 2;
};

/**
 * Тест деления двух 128-битных чисел.
 * @param z1 Делимое.
 * @param z2 Делитель.
 * @return Успех/неудача.
*/
bool test_div(U128 z1, U128 z2, PythonCaller& caller) {
    const U128 z3 = z1 / z2;
    PyObject* quotient = caller.Divide(z1, z2);
    return caller.Compare(quotient, z3.value().c_str());
}

/**
 * Полуслучайный тест деления 128-битных чисел используя ограниченный набор значений вблизи угловых и граничных.
 * @param N Количество внешних итераций.
*/
void test_division_semi_randomly(long long N) {
    if (N < 1) {
        std::cout << "Skipped!\n";
        return;
    }
    PythonCaller caller;
    const std::vector choice {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 
                                    65535, 65534, 65533, 65532, 65531, 65530, 
                                    16384, 16383, 16382, 16385, 16386, 16387, 16388, 
                                    -1ull, -2ull, -3ull, -4ull, -5ull, -6ull, -7ull};
    auto make_test = [&caller](const Quadrupole q, const Signess s) -> bool {
        return test_div(U128{.mHigh = q.A, .mLow = q.B, .mSign = s.s1},
                        U128{.mHigh = q.C, .mLow = q.D, .mSign = s.s2},
                        caller);
    };
    auto get_quadrupole = [&choice]() -> Quadrupole {
        auto idx1 = roll_uint() % choice.size();
        auto idx2 = roll_uint() % choice.size();
        auto idx3 = roll_uint() % choice.size();
        auto idx4 = roll_uint() % choice.size();
        Quadrupole q {choice[idx1], choice[idx2], choice[idx3], choice[idx4]};
        return q;
    };
    long long counter = 0;
    long long external_iterations = 0;
    bool is_ok = true;
    while (external_iterations < N) {
        ++counter;
        const Quadrupole q = get_quadrupole();
        const Signess s{roll_bool(), roll_bool()};
        if (q.C == 0 && q.D == 0) {
            continue;
        }
        is_ok &= make_test(q, s);
        assert(is_ok);
        if (counter % (1ll << 20) == 0) {
            external_iterations++;
            std::cout << "... iterations: " << counter << ". External: " << 
                external_iterations << " from " << N << '\n';
        }
    }
}

/**
 * Случайный тест деления 128-битных чисел.
 * @param N Количество внешних итераций.
*/
void test_division_randomly(long long N) {
    if (N < 1) {
        std::cout << "Skipped!\n";
        return;
    }
    PythonCaller caller;
    auto make_test = [&caller](const Quadrupole q, const Signess s) -> bool {
        return test_div(U128{.mHigh = q.A, .mLow = q.B, .mSign = s.s1},
                        U128{.mHigh = q.C, .mLow = q.D, .mSign = s.s2},
                        caller);
    };
    auto get_quadrupole = []() -> Quadrupole {
        Quadrupole q {roll_ulow(), roll_ulow(), roll_ulow(), roll_ulow()};
        return q;
    };
    long long counter = 0;
    long long external_iterations = 0;
    bool is_ok = true;
    while (external_iterations < N) {
        ++counter;
        const Quadrupole q = get_quadrupole();
        const Signess s{roll_bool(), roll_bool()};
        if (q.C == 0 && q.D == 0) {
            continue;
        }
        is_ok &= make_test(q, s);
        assert(is_ok);
        if (counter % (1ll << 20) == 0) {
            external_iterations++;
            std::cout << "... iterations: " << counter << ". External: " << 
                external_iterations << " from " << N << '\n';
        }
    }
}


int main(int argc, char* argv[]) {
    long long N = 10;
    if (argc > 1) {
        N = atoi(argv[1]);
        std::cout << "You set the number of external iterations N: " << N << '\n';
    }
    std::cout << "Run semi-random test...\n";
    test_division_semi_randomly(N);
    std::cout << "Ok\n";

    std::cout << "Run random test...\n";
    test_division_randomly(N);
    std::cout << "Ok\n";
    return 0;
}

static constexpr char DIGITS[10]{'0', '1', '2', '3', '4',
                                 '5', '6', '7', '8', '9'};

// High/Low структура 128-битного числа со знаком и флагом переполнения.
// Для иллюстрации алгоритма деления двух U128 чисел реализованы основные
// арифметические операторы, кроме умножения двух U128 чисел.
struct U128 {
    // Битовая полуширина половинок.
    static constexpr int mHalfWidth = (sizeof(ULOW) * CHAR_BIT) / 2;
    // Наибольшее значение половинок, M-1.
    static constexpr ULOW mMaxULOW = ULOW(-1);
    ULOW mHigh = 0;
    ULOW mLow = 0;
    bool mSign = 0;
    bool mOverflow = 0;

bool is_zero() const {
    return (mLow | mHigh | mOverflow) == 0;
}

bool is_negative() const {
    return mSign;
}

bool is_non_negative() const {
    return !mSign && !mOverflow;
}

bool is_nonzero_negative() const {
    return (mLow | mHigh) && mSign && !mOverflow;
}

bool is_overflow() const {
    return mOverflow;
}

constexpr static U128 get_zero() {
    return U128 {.mHigh = 0, .mLow = 0, .mSign = 0, .mOverflow = 0};
}

constexpr static U128 get_unit() {
    return U128 {.mHigh = 0, .mLow = 1, .mSign = 0, .mOverflow = 0};
}

constexpr static U128 get_unit_neg() {
    return U128 {.mHigh = 0, .mLow = 1, .mSign = 1, .mOverflow = 0};
}

U128 operator+(U128 rhs) const {
    U128 result{};
    U128 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;
}

U128& operator+=(U128 other) {
    *this = *this + other;
    return *this;
}

U128 operator-(U128 rhs) const {
    U128 result{};
    U128 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;
}

U128& operator-=(U128 other) {
    *this = *this - other;
    return *this;
}

U128 mult64(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;
    U128 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;
}

U128 operator*(ULOW rhs) const {
    U128 result = mult64(mLow, rhs);
    U128 tmp = mult64(mHigh, rhs);
    tmp.mHigh = tmp.mLow;
    tmp.mLow = 0;
    result += tmp;
    result.mSign = !result.is_zero() ? this->mSign : 0;
    return result;
}

U128 div10() const { // Специальный метод деления на 10 для строкового представления числа.
    U128 X = *this;
    const int sign = X.mSign;
    X.mSign = 0;
    ULOW Q = X.mHigh / 10;
    ULOW R = X.mHigh % 10;
    ULOW N = R * (mMaxULOW / 10) + (X.mLow / 10);
    U128 result{.mHigh = Q, .mLow = N};
    U128 E = X - result * 10;
    while (E.mHigh != 0 || E.mLow >= 10) {
        Q = E.mHigh / 10;
        R = E.mHigh % 10;
        N = R * (mMaxULOW / 10) + (E.mLow / 10);
        const U128 tmp {.mHigh = Q, .mLow = N};
        result += tmp;
        E -= tmp * 10;
    }
    result.mSign = sign;
    return result;
}

U128 operator/(ULOW y) const {
    assert(y != 0);
    const U128 X = *this;
    ULOW Q = X.mHigh / y;
    ULOW R = X.mHigh % y;
    ULOW N = R * (mMaxULOW / y) + (X.mLow / y);
    U128 result{.mHigh = Q, .mLow = N, .mSign = X.mSign};
    U128 E = X - result * y; // Ошибка от деления: остаток от деления.
    while (1) {
        Q = E.mHigh / y;
        R = E.mHigh % y;
        N = R * (mMaxULOW / y) + (E.mLow / y);
        const U128 tmp {.mHigh = Q, .mLow = N, .mSign = E.mSign};
        if (tmp.is_zero()) {
            break;
        }
        result += tmp;
        E -= tmp * y;
    }
    if (E.is_nonzero_negative()) {
        result -= get_unit();
        E += U128{.mHigh = 0, .mLow = y};
    }
    return result;
}

U128& operator/=(ULOW y) {
    *this = *this / y;
    return *this;
}

U128 operator/(const U128 other) const {
    U128 X = *this;
    U128 Y = other;
    constexpr U128 ZERO = get_zero();
    constexpr U128 UNIT = get_unit();
    constexpr U128 UNIT_NEG = get_unit_neg();
    if (Y.mHigh == 0) {
        X.mSign ^= Y.mSign;
        U128 result = X / Y.mLow;
        return result;
    }
    const bool make_sign_inverse = X.mSign != Y.mSign;
    X.mSign = make_sign_inverse;
    Y.mSign = 0;
    const ULOW Q = X.mHigh / Y.mHigh;
    const ULOW R = X.mHigh % Y.mHigh;
    const ULOW Delta = mMaxULOW - Y.mLow;
    const U128 DeltaQ = mult64(Delta, Q);
    U128 W1 = U128{.mHigh = R, .mLow = 0} - U128{.mHigh = Q, .mLow = 0};
    W1 = W1 + DeltaQ;
    const ULOW C1 = (Y.mHigh < mMaxULOW) ? Y.mHigh + ULOW(1) : mMaxULOW;
    const ULOW W2 = mMaxULOW - Delta / C1;
    U128 Quotient = W1 / W2;
    Quotient = Quotient / C1;
    U128 result = U128{.mHigh = 0, .mLow = Q} + Quotient;
    assert(result.mHigh == 0);
    result.mSign ^= make_sign_inverse;
    U128 N = Y * result.mLow;
    N.mSign ^= make_sign_inverse;
    assert(!N.is_overflow());
    U128 Error = X - N;
    U128 More = Error - Y;
    bool do_inc = More.is_non_negative();
    bool do_dec = Error.is_nonzero_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_non_negative();
        do_dec = Error.is_nonzero_negative();
    }
    return result;
}

/**
 * Возвращает строковое представление числа.
 */
std::string value() const {
    std::string result{};
    if (this->is_overflow()) {
        result = "Overflow";
        return result;
    }
    U128 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 = X.div10();
    }
    if (this->is_negative() && !this->is_zero()) {
        result.push_back('-');
    }
    std::reverse(result.begin(), result.end());
    return result.length() != 0 ? result : "0";
}
};