#pragma once
#include <stdint.h>
/* Algorithm due to Microsft's Coilot
function udiv_restoring(N, D, n) :
R = 0
Q = 0
negD = (~D + 1)
for i from n - 1 down to 0
{
R = (R << 1) | ((N >> i) & 1)
T = R + negD
if MSB(T) == 0:
R = T
Q = Q | (1 << i)
return (Q, R)
*/
class Arithmetic
{
public:
static bool udiv_restoring(
uint32_t numer,
uint32_t denom,
uint32_t& quo,
uint32_t& rem,
int n);
static bool umul_shift_add(
uint32_t a,
uint32_t b,
uint64_t& product,
int n);
};
#include <cstdint>
#include "Arithmetic.h"
static inline uint32_t mask_n(int bits) {
return (bits >= 32) ? 0xFFFFFFFFu : ((1u << bits) - 1u);
}
static inline uint32_t msb(uint32_t x, int bits) {
// returns top bit of a 'bits'-wide value
return (x >> (bits - 1)) & 1u;
}
bool Arithmetic::udiv_restoring(
uint32_t numer,
uint32_t denom,
uint32_t& quo,
uint32_t& rem,
int n)
{
if (denom == 0 || n <= 0 || n > 32) return false;
if (numer == 0)
{
quo = rem = 0;
return true;
}
quo = 0;
rem = 0;
if (n == 32) {
uint64_t R = 0;
uint64_t D = (uint64_t)denom;
uint64_t maskW = (1ull << 33) - 1ull; // 33-bit mask
uint64_t negD = ((~D) + 1ull) & maskW; // 33-bit two's complement
for (int i = n - 1; i >= 0; --i) {
R = ((R << 1) | ((numer >> i) & 1u)) & maskW;
uint64_t T = (R + negD) & maskW; // R - D
// Sign bit is bit 32 (the 33rd bit)
if (((T >> 32) & 1ull) == 0ull) {
R = T;
quo |= (1u << i);
}
}
rem = (uint32_t)(R & 0xFFFFFFFFu);
return true;
}
// n < 32 case: we can keep everything in uint32_t using (n+1) bits
uint32_t maskN = mask_n(n);
uint32_t maskW = mask_n(n + 1);
uint32_t N = numer & maskN;
uint32_t D = denom & maskN;
// Two's complement of D in (n+1) bits
uint32_t Dw = D; // placed in low bits of (n+1)-wide register
uint32_t negD = ((~Dw) + 1u) & maskW;
uint32_t R = 0;
for (int i = n - 1; i >= 0; --i) {
R = ((R << 1) | ((N >> i) & 1u)) & maskW;
uint32_t T = (R + negD) & maskW; // trial subtract: R - D (in w bits)
if (msb(T, n + 1) == 0) { // non-negative in (n+1) bits
R = T;
quo |= (1u << i);
}
}
rem = R & maskN; // remainder fits in n bits
return true;
}
bool Arithmetic::umul_shift_add(
uint32_t a,
uint32_t b,
uint64_t& product,
int n)
{
if (n <= 0 || n > 32) return false;
uint64_t A = a; // promote to avoid overflow
uint32_t B = b;
product = 0;
for (int i = 0; i < n; ++i) {
if (B & 1u) {
product += A;
}
A <<= 1;
B >>= 1;
}
return true;
}
#include <chrono>
#include <cstdint>
#include <iostream>
#include <limits>
#include <random>
#include <string>
#include "Arithmetic.h"
namespace {
constexpr int TESTS_PER_N = 200000;
uint32_t make_mask(int n) {
return (n == 32) ? 0xFFFFFFFFu : ((1u << n) - 1u);
}
void clear_bad_input() {
std::cin.clear();
std::cin.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
}
template <typename TrialFn>
double run_suite(const char* label, TrialFn trial, bool verbose) {
std::mt19937 rng(12345); // deterministic
auto t0 = std::chrono::high_resolution_clock::now();
for (int n = 1; n <= 32; ++n) {
const uint32_t mask = make_mask(n);
for (int i = 0; i < TESTS_PER_N; ++i) {
if (!trial(rng, mask, n)) {
std::cout << label << ": FAILED (n=" << n << ", i=" << i << ")\n";
return -1.0;
}
}
if (verbose) {
std::cout << "n=" << n << " passed\n";
}
}
auto t1 = std::chrono::high_resolution_clock::now();
double secs = std::chrono::duration<double>(t1 - t0).count();
std::cout << label << " runtime = " << secs << " sec\n";
return secs;
}
bool trial_division(std::mt19937& rng, uint32_t mask, int n) {
const uint32_t numer = rng() & mask;
const uint32_t denom = (rng() & mask) | 1u; // ensure non-zero
uint32_t q = 0, r = 0;
if (!Arithmetic::udiv_restoring(numer, denom, q, r, n)) {
std::cout << "Failure numer=" << numer << " denom=" << denom << "\n";
return false;
}
const uint32_t q2 = numer / denom;
const uint32_t r2 = numer % denom;
if (q != q2 || r != r2) {
std::cout << "Mismatch n=" << n
<< " numer=" << numer
<< " denom=" << denom
<< " got q=" << q << " r=" << r
<< " expected q=" << q2 << " r=" << r2 << "\n";
return false;
}
return true;
}
bool trial_multiplication(std::mt19937& rng, uint32_t mask, int n) {
const uint32_t a = rng() & mask;
const uint32_t b = rng() & mask;
uint64_t prod = 0;
if (!Arithmetic::umul_shift_add(a, b, prod, n)) {
std::cout << "Failure a=" << a << " b=" << b << "\n";
return false;
}
const uint64_t expected = static_cast<uint64_t>(a) * static_cast<uint64_t>(b);
if (prod != expected) {
std::cout << "Mismatch n=" << n
<< " a=" << a
<< " b=" << b
<< " got=" << prod
<< " expected=" << expected << "\n";
return false;
}
return true;
}
} // namespace
int main() {
bool verbose = true;
while (true) {
std::cout << "\nArithmetic Lab\n";
std::cout << "1. Test Division (restoring)\n";
std::cout << "2. Test Multiplication (shift-add)\n";
std::cout << "3. Run ALL tests\n";
std::cout << "4. Toggle verbose (currently " << (verbose ? "ON" : "OFF") << ")\n";
std::cout << "5. Exit\n";
std::cout << "Choice: ";
int choice = 0;
if (!(std::cin >> choice)) {
clear_bad_input();
std::cout << "Invalid input. Please enter a number.\n";
continue;
}
if (choice == 1) {
run_suite("Division test", trial_division, verbose);
}
else if (choice == 2) {
run_suite("Multiplication test", trial_multiplication, verbose);
}
else if (choice == 3) {
const double d = run_suite("Division test", trial_division, verbose);
if (d >= 0.0) run_suite("Multiplication test", trial_multiplication, verbose);
}
else if (choice == 4) {
verbose = !verbose;
}
else if (choice == 5) {
return 0;
}
else {
std::cout << "Invalid choice.\n";
}
}
}
Numerical Explorations 