// 128-bit signed integer:
S128 :: struct {
    low:  u64;
    high: s64;
}

// 128-bit unsigned integer:
U128 :: struct {
    low:  u64;
    high: u64;
}

// These values returned for quotient and remainder of division when overflow or division by zero
U128_ERROR :: U128.{0xffffffff_ffffffff, 0xffffffff_ffffffff};
S128_ERROR :: S128.{0xffffffff_ffffffff, 0xffffffff_ffffffff};


// If you want a 128-bit result, well, we don't have any primitive types
// that hold that, so, here is some code to do it within u64.
// We might put an asm version here, but right now this is used to
// avoid asm stuff...!
multiply_64bits_with_overflow :: (a: u64, b: u64) -> (high: u64, low: u64) {
    lo_lo := (a & 0xFFFF_FFFF) * (b & 0xFFFF_FFFF);
    hi_lo := (a >> 32)         * (b & 0xFFFF_FFFF);
    lo_hi := (a & 0xFFFF_FFFF) * (b >> 32);
    hi_hi := (a >> 32)         * (b >> 32);

    cross := (lo_lo >> 32) + (hi_lo & 0xFFFF_FFFF) + lo_hi;

    high := (hi_lo >> 32) + (cross >> 32) + hi_hi;
    low  := (cross << 32) | (lo_lo & 0xFFFF_FFFF);

    return high, low;
}

to_s128 :: (a: $T) -> S128
#modify {
    return (cast(*Type_Info)T).type == .INTEGER;
} {
    return .{cast,no_check(u64) a, (cast(s64)a) >> 63};
}

to_u128 :: (a: $T) -> U128
#modify {
    return T == u64 || T == u32 || T == u16 || T == u8;
} {
    return U128.{a, 0};
}


operator == :: inline (a: S128, b: S128) -> bool {
    return (a.low == b.low) && (a.high == b.high);
}

operator == :: inline (a: U128, b: U128) -> bool {
    return (a.low == b.low) && (a.high == b.high);
}

operator ! :: inline (a: S128) -> bool {
    return a.low == 0 && a.high == 0;
}

operator ! :: inline (a: U128) -> bool {
    return a.low == 0 && a.high == 0;
}

operator + :: inline (a: S128, b: S128) -> S128 { // This can be the same code for S128 and U128.
    c: S128 = ---;

    #if CPU == .X64 {
        #asm {
            mov a_low:,  [*a + 0];
            mov a_high:, [*a + 8];

            mov b_low:,  [*b + 0];
            mov b_high:, [*b + 8];

            mov [*c + 0], a_low;
            mov [*c + 8], a_high;

            add [*c + 0], b_low;
            adc [*c + 8], b_high;
        }
    } else {
        #no_aoc {  // Let's just no_aoc this one operation for now. In theory high bits should not overflow.
            c.low  = a.low + b.low;
        }
        
        c.high = a.high + b.high;
        if c.low < a.low  c.high += 1;  // We overflowed!
    }

    return c;
}

operator - :: inline (a: S128, b: S128) -> S128 {
    return cast,force(S128) (cast,force(U128)a - cast,force(U128)b);
}

operator - :: inline (a: U128, b: U128) -> U128 {
    c: U128 = ---;

    #if CPU == .X64 {
        #asm {
            mov a_low:,  [*a + 0];
            mov a_high:, [*a + 8];

            mov b_low:,  [*b + 0];
            mov b_high:, [*b + 8];

            mov [*c + 0], a_low;
            mov [*c + 8], a_high;

            sub [*c + 0], b_low;
            sbb [*c + 8], b_high;

            // I guess these are not legal forms of 'add':
            // add [*c + 0], a_low,  b_low;
            // adc [*c + 8], a_high, b_high;
        }
    } else #no_aoc {
        c.high = a.high - b.high;
        c.low  = a.low  - b.low;
        if c.low > a.low  c.high -= 1;  // We underflowed!
    }

    return c;
}

operator - :: inline (a: S128) -> S128 {
    // Apollo Time is two's complement.
    // To negate it, we flip all the bits and then add 1.
    // We'll use #asm to avoid a branch.

    c: S128 = ---;
    c.low  = ~a.low;
    c.high = ~a.high;

    #if CPU == .X64 {
        #asm {
            add [*c + 0], 1;
            adc [*c + 8], 0;
        }
    } else {
        #no_aoc {
            c.low += 1;
        }
        
        if !c.low  c.high += 1;  // If c.low overflowed, carry into high.
    }

    return c;
}

// Complexity-wise, it may not be worth specializing << and >> for constant x,
// since the win is relatively low these days. But I did it anyway.
operator << :: (a: S128, $$x: u8) -> S128 {
    result : S128 = ---;

    if x >= 64 {
        new_high := cast,no_check(s64)(a.low << (x - 64));
        return .{0, new_high};
    } else {
        #if CPU == .X64 {
            result = a;

            #if !is_constant(x)  #asm { x === c; }  // Use the cl form of shift.

            #asm {
                mov low:, [*result + 0];
                shld [*result + 8], low, x;
                shl  [*result + 0], x;
            }
        } else {
            result.high = a.high << x | cast(s64)(a.low >> (64 - x));
            result.low  = a.low << x;
        }
    }

    return result;
}

operator << :: inline (a: U128, $$x: u8) -> U128 {
    result : U128 = ---;

    if x >= 64 {
        new_high := a.low << (x - 64);
        return .{0, new_high};
    } else {
        #if CPU == .X64 {
            result = a;

            #if !is_constant(x)  #asm { x === c; }  // Use the cl form of shift.

            #asm {
                mov low:, [*result + 0];
                shld [*result + 8], low, x;
                shl  [*result + 0], x;
            }
        } else {
            result.high = a.high << x | (a.low >> (64 - x));
            result.low  = a.low << x;
        }
    }

    return result;
}

operator >> :: (a: S128, $$x: u8) -> S128 {  // Arithmetic shift
    if x >= 64 {
        new_low := cast,no_check (u64) (a.high >> (x - 64));
        return .{low = new_low, high = a.high >> 63};
    } else {
        #if CPU == .X64 {
            result := a;

            #if !is_constant(x)  #asm { x === c; }  // Use the cl form of shift.

            #asm {
                mov high:, [*result + 8];
                shrd [*result + 0], high, x;
                sar  [*result + 8], x;
            }
        } else {
            result: S128 = ---;
            result.low = (a.low >> x) | cast(u64)(a.high << (64 - x));
            result.high = a.high >> x;
        }

        return result;
    }
}

operator >> :: (a: U128, $$x: u8) -> U128 {  // Logical shift.
    if !a || !x  return a;

    if x >= 64 {
        return .{low = a.high >> (x - 64), high = 0};
    } else {
        #if CPU == .X64 {
            result := a;
            #if !is_constant(x)  #asm { x === c; }  // Use the cl form of shift.

            #asm {
                mov high:, [*result + 8];
                shrd [*result + 0], high, x;
                shr  [*result + 8], x;
            }
        } else {
            result: U128 = ---;
            result.low = a.low >> x | (a.high << (64 - x));
            result.high = a.high >> x;
        }

        return result;
    }
}

// @Incomplete: Refactor logical_shift_right as >> on U128.
logical_shift_right :: (a: S128, $$x: u8) -> S128 {
    return cast,force(S128) ((cast,force(U128) a) >> x);
}


// Our division procs handle divide-by-zero and overflow without aborting, returning an error flag and
// setting the result to U128_MAX.  While using operator `/` you don't get this error flag, and the expected
// result of any operator `/` divide-by-zero is boom, so here we'll assert to simulate that behaviour.
// OTOH, overflow is bit weird and not as easy to anticipate, so we'll allow it to return (and the programmer can
// test for the error return value when they notice weird behaviour with their results, but their program won't
// unexpectedly crash in production when the user puts in some bad values).
// Maybe it'd be better to halt on both conditions, with the programmer expected to validate the inputs before
// calling, but the easiest way for them to do that would be for them to call the (un)signed_divide_with_remainder
// proc themselves... at which point they don't need the `/` operator, so we're letting the `/` operator be the
// "easy use" version of divide. - IK
operator / :: (dividend: U128, divisor: U128) -> U128 {
    assert(divisor != .{0, 0}, "Cannot divide by zero!");
    return inline unsigned_divide_with_remainder(dividend, divisor);
}

operator / :: (dividend: S128, divisor: S128) -> S128 {
    assert(divisor != .{0, 0}, "Cannot divide by zero!");
    return inline signed_divide_with_remainder(dividend, divisor);
}

operator / :: inline (a: S128, factor: s64) -> S128 {
    // We'll put a faster version here later?

    b: S128;

    negate := false;
    b.low = cast,no_check(u64) factor;
    if factor < 0  b.high = -1;

    return a / b;
}


signed_divide_with_remainder :: inline (dividend: S128, divisor: S128) -> quotient: S128, remainder: S128, overflow_or_divide_by_zero: bool {
    u := dividend;
    v := divisor;

    negate := false;
    if u.high < 0 {
        negate = true;
        u = -u;
    }

    if v.high < 0 {
        negate = !negate;
        v = -v;
    }

    quotient, remainder, overflow_or_divide_by_zero := unsigned_divide_with_remainder(cast,force(U128)u, cast,force(U128)v);
    if overflow_or_divide_by_zero  return  cast,force(S128)quotient, cast,force(S128)remainder, true;
    else if negate                 return -cast,force(S128)quotient, cast,force(S128)remainder, false;
    else                           return  cast,force(S128)quotient, cast,force(S128)remainder, false;
}


unsigned_divide_with_remainder :: (u: U128, v: U128) -> quotient: U128, remainder: U128, overflow_or_divide_by_zero: bool {
    // This unsigned divide is used to help implement our operator /, which is
    // a wrapper that does sign manipulations then calls this routine.

    // Code adapted from Hacker's Delight by Henry S. Warren, Jr., Second Edition, section 9-5.
    // That listing was for 64-bit division, but here we do 128-bit division.

    if !v  return U128_ERROR, U128_ERROR, true;

    #if CPU == .X64  { // Lifted these procs out of here so they can be tested against each other
        DIVU_FF :: unsigned_divide_with_remainder_128_64_x64;
    } else {
        DIVU_FF :: unsigned_divide_with_remainder_128_64_no_asm;
    }

    if v.high == 0 {
        u1 := cast,no_check(u64) u.high;
        d  := v.low;

        if u1 < d { // If u/v cannot overflow, just do one division.
            quotient: U128 = ---;
            remainder: U128 = ---;
            quotient.high, remainder.high = 0;
            error : bool = ---;
            quotient.low, remainder.low, error = DIVU_FF(u, d);
            if error  return U128_ERROR, U128_ERROR, true;
            else      return quotient, remainder, false;
        } else {
            // If u/v would overflow: Break up u into two halves.
            u0 := u.low;
            // q1 := DIVU_FF(u1, d);   // First quotient digit.
            q1 := u1 / d;           // First quotient digit.
            k := u1 - q1*d;         // First remainder, < v.

            c: U128 = ---;
            c.high = k;
            c.low  = u0;

            q0, remainder_u64, error := DIVU_FF(c, d);
            quotient := U128.{low = q0, high = q1};
            remainder := U128.{low = remainder_u64, high = 0};

            if error  return U128_ERROR, U128_ERROR, true;
            else      return quotient, remainder, false;
        }
    }

    n := count_leading_zero_bits(v.high);

    _v1 := v << n;
    // Warren shifts right here to get the top bits, but we can just use .high
    v1 := cast(u64) _v1.high;

    // Right-shift u by 1 bit, to ensure no overflow.
    u1 := (u >> 1);

    q1, _, error := DIVU_FF(u1, v1);
    if error  return U128_ERROR, U128_ERROR, true;

    // Undo normalization and division of u by 2.
    _q1 := U128.{low = q1};
    _q1 <<= n;
    q0 := _q1 >> 63;

    // Make q0 correct or too small by 1.
    if q0.high || q0.low {  // q0 != 0
        if q0.low == 0 {
            // Borrow from the top; could use #asm.
            q0.low  -= 1;
            q0.high -= 1;
        } else {
            q0.low  -= 1;
        }
    }

    remainder := u - q0*v;
    if remainder >= v {  // Finally correct q0.
        q0.low += 1;
        if q0.low == 0  q0.high += 1;
        remainder -= v;
    }

    return q0, remainder, false;
}

#scope_file

U64_MAX :: 0xFFFF_FFFF_FFFF_FFFF;

#if CPU == .X64 {
    unsigned_divide_with_remainder_128_64_x64 :: (a: U128, d: u64) -> quotient: u64, remainder: u64, overflow_or_divide_by_zero: bool {
        #assert CPU == .X64;

        if a.high >= d
            return U64_MAX, U64_MAX, true;

        quotient: u64 = ---;
        remainder: u64 = ---;
        #asm {
            high: gpr === d;
            low:  gpr === a;

            mov low,  [*a + 0];
            mov high, [*a + 8];
            div high, low, d;

            mov quotient, low;
            mov remainder, high;
        }

        return quotient, remainder, false;
    }
}

unsigned_divide_with_remainder_128_64_no_asm :: (num: U128, den: u64) -> quotient: u64, remainder: u64, overflow_or_divide_by_zero: bool #no_aoc {
    // This software version is derived from public domain/CC0 code found at:
    // https://ridiculousfish.com/blog/posts/labor-of-division-episode-v.html

    numhi := num.high;
    numlo := num.low;

    // Check for overflow and divide by 0.
    if numhi >= den
        return U64_MAX, U64_MAX, true;

    // We work in base 2**32.
    // A uint32 holds a single digit. A uint64 holds two digits.
    // Our numerator is conceptually [num3, num2, num1, num0].
    // Our denominator is [den1, den0].
    b : u64 : (1 << 32);

    // The high and low digits of our computed quotient.
    q1: u32 = ---;
    q0: u32 = ---;

    // The normalization shift factor.
    shift: int = ---;

    // The high and low digits of our denominator (after normalizing).
    // Also the low 2 digits of our numerator (after normalizing).
    den1: u32 = ---;
    den0: u32 = ---;
    num1: u32 = ---;
    num0: u32 = ---;

    // A partial remainder.
    rem: u64 = ---;

    // The estimated quotient, and its corresponding remainder (unrelated to true remainder).
    qhat: u64 = ---;
    rhat: u64 = ---;

    // Variables used to correct the estimated quotient.
    c1: u64 = ---;
    c2: u64 = ---;

    // Determine the normalization factor. We multiply den by this, so that its leading digit is at
    // least half b. In binary this means just shifting left by the number of leading zeros, so that
    // there's a 1 in the MSB.
    // We also shift numer by the same amount. This cannot overflow because numhi < den.
    // The expression (-shift & 63) is the same as (64 - shift), except it avoids the UB of shifting
    // by 64. The funny bitwise 'and' ensures that numlo does not get shifted into numhi if shift is 0.
    // clang 11 has an x86 codegen bug here: see LLVM bug 50118. The sequence below avoids it.
    shift = count_leading_zero_bits(den);
    den <<= cast(u64)shift;
    numhi <<= cast(u64)shift;
    numhi |= (numlo >> (-shift & 63)) & cast,no_check(u64)(-shift >> 63);
    numlo <<= cast(u64)shift;

    // Extract the low digits of the numerator and both digits of the denominator.
    num1 = cast(u32)(numlo >> 32);
    num0 = cast(u32)(numlo & 0xFFFF_FFFF);
    den1 = cast(u32)(den >> 32);
    den0 = cast(u32)(den & 0xFFFF_FFFF);

    // We wish to compute q1 = [n3 n2 n1] / [d1 d0].
    // Estimate q1 as [n3 n2] / [d1], and then correct it.
    // Note while qhat may be 2 digits, q1 is always 1 digit.
    qhat = numhi / den1;
    rhat = numhi % den1;
    c1 = qhat * den0;
    c2 = rhat * b + num1;
    if c1 > c2
        qhat -= ifx c1 - c2 > den then cast(u64)2 else 1;
    q1 = cast(u32)qhat;

    // Compute the true (partial) remainder.
    rem = numhi * b + num1 - q1 * den;

    // We wish to compute q0 = [rem1 rem0 n0] / [d1 d0].
    // Estimate q0 as [rem1 rem0] / [d1] and correct it.
    qhat = rem / den1;
    rhat = rem % den1;
    c1 = qhat * den0;
    c2 = rhat * b + num0;
    if c1 > c2
        qhat -= ifx c1 - c2 > den then cast(u64)2 else 1;
    q0 = cast(u32)qhat;

    quotient := ((cast(u64)q1) << 32) | q0;
    remainder := (rem * b + num0 - q0 * den) >> shift;
    return quotient, remainder, false;
}

count_leading_zero_bits :: inline (x: u64) -> u8 {
    // Maybe Basic should have a Simple_Math component with stuff like this in it? - IK
    n: u8 = ---;
    #if CPU == .X64 {
        #asm {
            lzcnt n, x;
        }
    } else {
        if x == 0 return(64);
        n = 0;
        if x <= 0x00000000FFFFFFFF { n = n + 32; x = x << 32; }
        if x <= 0x0000FFFFFFFFFFFF { n = n + 16; x = x << 16; }
        if x <= 0x00FFFFFFFFFFFFFF { n = n +  8; x = x <<  8; }
        if x <= 0x0FFFFFFFFFFFFFFF { n = n +  4; x = x <<  4; }
        if x <= 0x3FFFFFFFFFFFFFFF { n = n +  2; x = x <<  2; }
        if x <= 0x7FFFFFFFFFFFFFFF { n = n +  1; }
    }
    return n;
}


#scope_export

//
// The sign handling in these operator * below is a little uncertain given
// that I don't really understand the difference between imul and mul.
// @Incomplete: Consult Hacker's Delight.
//
operator * :: (_a: S128, _b: S128) -> S128 {
    a := _a; // @Speed @Temporary
    b := _b;

    negate := false;
    if a.high < 0 {
        a = -a;
        negate = true;
    }

    if b.high < 0 {
        b = -b;
        negate = !negate;
    }

    u := cast,force(U128)a * cast,force(U128) b;
    result := cast,force(S128) u;

    if negate return -result;
    else      return result;
}

operator * :: (_a: U128, _b: U128) -> U128 {
    a := _a; // @Speed @Temporary
    b := _b;

    result: U128 = ---;

    // The code below is really convoluted because I am not sure if we can
    // bind multiple names to the same physical register? Would really like to.
    #if CPU == .X64 {
        #asm {
            result_high: gpr === d;
            result_low:  gpr === a;

            factor:   gpr;
            high_sum: gpr;

            // Low times low:
            mov result_low, [*a + 0];  // This is input as well as result.
            mov factor,     [*b + 0];

            mul result_high, result_low, factor;

            mov [*result + 0],  result_low;   // Store the low part.
            mov high_sum, result_high;  // Store the high part.

            // High a times low b. 'factor' still holds low b.
            mov result_low, [*a + 8];
            imul result_high, result_low, factor;
            add high_sum, result_low;

            // Low a times high b.
            mov result_low, [*a + 0];
            mov factor,     [*b + 8];
            imul result_high, result_low, factor;
            add high_sum, result_low;

            // High a times high b is discarded.

            mov [*result + 8], high_sum;
        }
    } else {
        // Intentionally leaving aoc enabled here, since these should not in principle overflow.
        result.high, result.low = multiply_64bits_with_overflow(a.low, b.low);
        result.high += a.high * b.low + a.low * b.high;
    }

    return result;
}

operator * :: (a: S128, factor: s64) -> S128 #symmetric {
    at := a; // @Speed @Temporary

    negate := false;
    if at.high < 0 {
        at = -at;
        negate = true;
    }

    if factor < 0 {
        factor = -factor;
        negate = !negate;
    }

    result: S128 = ---;

    // We would like to be able to bind multiple names to the same machine register
    // to make code like the following less confusing.
    #if CPU == .X64 {
        #asm {
            result_high: gpr === d;
            result_low:  gpr === a;

            mov result_low, [*at + 0];  // This is input as well as result.
            mul result_high, result_low, factor;

            mov [*result + 0], result_low;

            mov carry:, result_high;

            mov result_low, [*at + 8];  // This is input as well as result.
            imul result_high, result_low, factor;

            add result_low, carry;

            mov [*result + 8], result_low;  // This is actually the high.
        }
    } else #no_aoc {
        u128 : U128 = ---;
        u128.high, u128.low = multiply_64bits_with_overflow(at.low, cast(u64)factor);
        u128.high += cast(u64)at.high * cast(u64)factor;
        result = cast,force(S128)u128;
    }

    if negate return -result;
    else      return result;
}

operator < :: (a: S128, b: S128) -> bool {
    result: bool = ---;

    #if CPU == .X64 {
        #asm {
            mov a_low:,  [*a + 0];
            mov a_high:, [*a + 8];

            mov b_low:,  [*b + 0];
            mov b_high:, [*b + 8];

            sub a_low,  b_low;
            sbb a_high, b_high;

            setl  result;
        }
    } else {
        if a.high < b.high return true;
        if a.high > b.high return false;
        return a.low < b.low;
    }

    return result;
}

operator < :: (a: U128, b: U128) -> bool {  // @Speed: Do #asm.
    if (a.high < b.high) return true;
    if (a.high > b.high) return false;

    return (a.low < b.low);
}

operator <= :: (a: S128, b: S128) -> bool {
    // We would like the compiler to generate this for us, probably, like it does with +=.
    return (a == b) || (a < b);  // Probably not very efficient.
}

operator > :: (a: S128, b: S128) -> bool {
    // We would like the compiler to generate this for us, probably, like it does with +=.
    return !(a <= b);
}

operator >= :: (a: S128, b: S128) -> bool {
    // We would like the compiler to generate this for us, probably, like it does with +=.
    return !(a < b);
}

operator <= :: (a: U128, b: U128) -> bool {
    // We would like the compiler to generate this for us, probably, like it does with +=.
    return (a == b) || (a < b);  // Probably not very efficient.
}

operator > :: (a: U128, b: U128) -> bool {
    // We would like the compiler to generate this for us, probably, like it does with +=.
    return !(a <= b);
}

operator >= :: (a: U128, b: U128) -> bool {
    // We would like the compiler to generate this for us, probably, like it does with +=.
    return !(a < b);
}


operator * :: (a: U128, b: s64) -> U128 {
    return a * cast(u64)b;
}

operator * :: (a: U128, b: u64) -> U128 {
    b128 : U128 = ---;
    b128.high = 0;
    b128.low = b;
    return a * b128;
}

operator + :: (a: U128, b: U128) -> U128 {
    return cast,force(U128)(cast,force(S128)a + cast,force(S128)b);
}

operator + :: (a: U128, b: u64) -> U128 {
    return a + U128.{b, 0};
}

#if CPU == .X64 {
    test_unsigned_division :: (random_iterations := 0, random_get : () -> u64 = null) {
        to_string :: (x: U128) -> string {
            if x.high
                return tprint("0x%0_%0", formatInt(x.high, base=16), formatInt(x.low, base=16, minimum_digits=16));
            else
                return tprint("%0", x.low);
        }

        test_pattern : u64 : 0xAAAAAAAA_AAAAAAAA;
        dividend := U128.{test_pattern, test_pattern};

        // Test x64 version against no_asm version - both should give identical results!
        // Also test the result played backwards through multiply gives the start point.
        while true  {
            divisor := test_pattern;
            while true {
                x64_quotient, x64_remainder, x64_error := unsigned_divide_with_remainder_128_64_x64(dividend, divisor);
                jai_quotient, jai_remainder, jai_error := unsigned_divide_with_remainder_128_64_no_asm(dividend, divisor);
                assert(x64_quotient == jai_quotient);
                assert(x64_remainder == jai_remainder);
                assert(x64_error == jai_error);
                if !x64_error {
                    print("[OK] % / % = % remainder %\n", to_string(dividend), divisor, x64_quotient, x64_remainder);
                    assert(U128.{divisor, 0} * x64_quotient + x64_remainder == dividend);
                }
                if !divisor  break;
                divisor >>= 1;
            }
            if !dividend  break;
            dividend >>= 1;
        }

        // Test `/` operator
        dividend = U128.{test_pattern, test_pattern};
        while true {
            divisor := U128.{test_pattern, test_pattern};
            while true {
                quotient, remainder, overflow_or_div_zero := unsigned_divide_with_remainder(dividend, divisor);
                operator_quotient := dividend / divisor;
                if overflow_or_div_zero {
                    print("[OVR] % / %\n", to_string(dividend), to_string(divisor));
                    assert(operator_quotient.high == U64_MAX && operator_quotient.low == U64_MAX);
                } else {
                    print("[OK] % / % = % remainder %\n", to_string(dividend), to_string(divisor), to_string(quotient), to_string(remainder));
                    assert(quotient * divisor + remainder == dividend);
                    assert(quotient == operator_quotient);
                    operator_remainder := dividend - (operator_quotient * divisor);
                    assert(remainder == operator_remainder);
                    if quotient == .{}
                        assert(remainder == dividend);
                }
                divisor >>= 1; // Currently operator `/` halts on divide-by-zero - if that behaviour changes then swap this line and the one below it around
                if !divisor  break;
            }
            if !dividend  break;
            dividend >>= 1;
        }

        // Test randomly
        for 1 .. random_iterations {
            assert(random_get != null);
            x := U128.{low = random_get(), high = 0};
            y := U128.{low = random_get(), high = 0};
            r := U128.{low = random_get() % x.low, high = 0};
            dividend = x * y + r;
            divisor := x;

            quotient, remainder, overflow_or_div_zero := unsigned_divide_with_remainder(dividend, divisor);
            operator_quotient := dividend / divisor;
            if overflow_or_div_zero {
                print("[OVR] % / %\n", to_string(dividend), to_string(divisor));
                assert(operator_quotient.high == U64_MAX && operator_quotient.low == U64_MAX);
            } else {
                print("[OK] % / % = % remainder %\n", to_string(dividend), to_string(divisor), to_string(quotient), to_string(remainder));
                assert(quotient * divisor + remainder == dividend);
                assert(quotient == operator_quotient);
                operator_remainder := dividend - (operator_quotient * divisor);
                assert(remainder == operator_remainder);
                assert(divisor == x);
                assert(quotient == y);
                assert(remainder == r);
                if quotient == .{}
                    assert(remainder == dividend);
            }
        }
    }
}