rational.zig

const std = @import("std");
const builtin = @import("builtin");
const debug = std.debug;
const math = std.math;
const mem = std.mem;
const Allocator = mem.Allocator;
const ArrayList = std.ArrayList;

const TypeId = builtin.TypeId;

const bn = @import("int.zig");
const Limb = bn.Limb;
const DoubleLimb = bn.DoubleLimb;
const Int = bn.Int;

pub const Rational = struct {
    // sign of Rational is a.positive, b.positive is ignored
    p: Int,
    q: Int,

    pub fn init(a: *Allocator) !Rational {
        return Rational{
            .p = try Int.init(a),
            .q = try Int.initSet(a, 1),
        };
    }

    pub fn deinit(self: *Rational) void {
        self.p.deinit();
        self.q.deinit();
    }

    pub fn setInt(self: *Rational, a: var) !void {
        try self.p.set(a);
        try self.q.set(1);
    }

    // TODO: Accept a/b fractions and exponent form
    pub fn setFloatString(self: *Rational, str: []const u8) !void {
        if (str.len == 0) {
            return error.InvalidFloatString;
        }

        const State = enum {
            Integer,
            Fractional,
        };

        var state = State.Integer;
        var point: ?usize = null;

        var start: usize = 0;
        if (str[0] == '-') {
            start += 1;
        }

        for (str) |c, i| {
            switch (state) {
                State.Integer => {
                    switch (c) {
                        '.' => {
                            state = State.Fractional;
                            point = i;
                        },
                        '0'...'9' => {
                            // okay
                        },
                        else => {
                            return error.InvalidFloatString;
                        },
                    }
                },
                State.Fractional => {
                    switch (c) {
                        '0'...'9' => {
                            // okay
                        },
                        else => {
                            return error.InvalidFloatString;
                        },
                    }
                },
            }
        }

        // TODO: batch the multiplies by 10
        if (point) |i| {
            try self.p.setString(10, str[0..i]);

            const base = Int.initFixed(([]Limb{10})[0..]);

            var j: usize = start;
            while (j < str.len - i - 1) : (j += 1) {
                try self.p.mul(self.p, base);
            }

            try self.q.setString(10, str[i + 1 ..]);
            try self.p.add(self.p, self.q);

            try self.q.set(1);
            var k: usize = i + 1;
            while (k < str.len) : (k += 1) {
                try self.q.mul(self.q, base);
            }

            try self.reduce();
        } else {
            try self.p.setString(10, str[0..]);
            try self.q.set(1);
        }
    }

    // Translated from golang.go/src/math/big/rat.go.
    pub fn setFloat(self: *Rational, comptime T: type, f: T) !void {
        debug.assert(@typeId(T) == builtin.TypeId.Float);

        const UnsignedIntType = @IntType(false, T.bit_count);
        const f_bits = @bitCast(UnsignedIntType, f);

        const exponent_bits = math.floatExponentBits(T);
        const exponent_bias = (1 << (exponent_bits - 1)) - 1;
        const mantissa_bits = math.floatMantissaBits(T);

        const exponent_mask = (1 << exponent_bits) - 1;
        const mantissa_mask = (1 << mantissa_bits) - 1;

        var exponent = @intCast(i16, (f_bits >> mantissa_bits) & exponent_mask);
        var mantissa = f_bits & mantissa_mask;

        switch (exponent) {
            exponent_mask => {
                return error.NonFiniteFloat;
            },
            0 => {
                // denormal
                exponent -= exponent_bias - 1;
            },
            else => {
                // normal
                mantissa |= 1 << mantissa_bits;
                exponent -= exponent_bias;
            },
        }

        var shift: i16 = mantissa_bits - exponent;

        // factor out powers of two early from rational
        while (mantissa & 1 == 0 and shift > 0) {
            mantissa >>= 1;
            shift -= 1;
        }

        try self.p.set(mantissa);
        self.p.positive = f >= 0;

        try self.q.set(1);
        if (shift >= 0) {
            try self.q.shiftLeft(self.q, @intCast(usize, shift));
        } else {
            try self.p.shiftLeft(self.p, @intCast(usize, -shift));
        }

        try self.reduce();
    }

    // Translated from golang.go/src/math/big/rat.go.
    pub fn toFloat(self: Rational, comptime T: type) !T {
        debug.assert(@typeId(T) == builtin.TypeId.Float);

        const fsize = T.bit_count;
        const BitReprType = @IntType(false, T.bit_count);

        const msize = math.floatMantissaBits(T);
        const msize1 = msize + 1;
        const msize2 = msize1 + 1;

        const esize = math.floatExponentBits(T);
        const ebias = (1 << (esize - 1)) - 1;
        const emin = 1 - ebias;
        const emax = ebias;

        if (self.p.eqZero()) {
            return 0;
        }

        // 1. left-shift a or sub so that a/b is in [1 << msize1, 1 << (msize2 + 1)]
        var exp = @intCast(isize, self.p.bitCountTwosComp()) - @intCast(isize, self.q.bitCountTwosComp());

        var a2 = try self.p.clone();
        defer a2.deinit();

        var b2 = try self.q.clone();
        defer b2.deinit();

        const shift = msize2 - exp;
        if (shift >= 0) {
            try a2.shiftLeft(a2, @intCast(usize, shift));
        } else {
            try b2.shiftLeft(b2, @intCast(usize, -shift));
        }

        // 2. compute quotient and remainder
        var q = try Int.init(self.p.allocator.?);
        defer q.deinit();

        // unused
        var r = try Int.init(self.p.allocator.?);
        defer r.deinit();

        try Int.divTrunc(&q, &r, a2, b2);

        var mantissa = extractLowBits(q, BitReprType);
        var have_rem = r.len > 0;

        // 3. q didn't fit in msize2 bits, redo division b2 << 1
        if (mantissa >> msize2 == 1) {
            if (mantissa & 1 == 1) {
                have_rem = true;
            }
            mantissa >>= 1;
            exp += 1;
        }
        if (mantissa >> msize1 != 1) {
            @panic("unexpected bits in result");
        }

        // 4. Rounding
        if (emin - msize <= exp and exp <= emin) {
            // denormal
            const shift1 = @intCast(math.Log2Int(BitReprType), emin - (exp - 1));
            const lost_bits = mantissa & ((@intCast(BitReprType, 1) << shift1) - 1);
            have_rem = have_rem or lost_bits != 0;
            mantissa >>= shift1;
            exp = 2 - ebias;
        }

        // round q using round-half-to-even
        var exact = !have_rem;
        if (mantissa & 1 != 0) {
            exact = false;
            if (have_rem or (mantissa & 2 != 0)) {
                mantissa += 1;
                if (mantissa >= 1 << msize2) {
                    // 11...1 => 100...0
                    mantissa >>= 1;
                    exp += 1;
                }
            }
        }
        mantissa >>= 1;

        const f = math.scalbn(@intToFloat(T, mantissa), @intCast(i32, exp - msize1));
        if (math.isInf(f)) {
            exact = false;
        }

        return if (self.p.positive) f else -f;
    }

    pub fn setRatio(self: *Rational, p: var, q: var) !void {
        try self.p.set(p);
        try self.q.set(q);

        self.p.positive = (@boolToInt(self.p.positive) ^ @boolToInt(self.q.positive)) == 0;
        self.q.positive = true;
        try self.reduce();

        if (self.q.eqZero()) {
            @panic("cannot set rational with denominator = 0");
        }
    }

    pub fn copyInt(self: *Rational, a: Int) !void {
        try self.p.copy(a);
        try self.q.set(1);
    }

    pub fn copyRatio(self: *Rational, a: Int, b: Int) !void {
        try self.p.copy(a);
        try self.q.copy(b);

        self.p.positive = (@boolToInt(self.p.positive) ^ @boolToInt(self.q.positive)) == 0;
        self.q.positive = true;
        try self.reduce();
    }

    pub fn abs(r: *Rational) void {
        r.p.abs();
    }

    pub fn negate(r: *Rational) void {
        r.p.negate();
    }

    pub fn swap(r: *Rational, other: *Rational) void {
        r.p.swap(&other.p);
        r.q.swap(&other.q);
    }

    pub fn cmp(a: Rational, b: Rational) !i8 {
        return cmpInternal(a, b, true);
    }

    pub fn cmpAbs(a: Rational, b: Rational) !i8 {
        return cmpInternal(a, b, false);
    }

    // p/q > x/y iff p*y > x*q
    fn cmpInternal(a: Rational, b: Rational, is_abs: bool) !i8 {
        // TODO: Would a div compare algorithm of sorts be viable and quicker? Can we avoid
        // the memory allocations here?
        var q = try Int.init(a.p.allocator.?);
        defer q.deinit();

        var p = try Int.init(b.p.allocator.?);
        defer p.deinit();

        try q.mul(a.p, b.q);
        try p.mul(b.p, a.q);

        return if (is_abs) q.cmpAbs(p) else q.cmp(p);
    }

    // r/q = ap/aq + bp/bq = (ap*bq + bp*aq) / (aq*bq)
    //
    // For best performance, rma should not alias a or b.
    pub fn add(rma: *Rational, a: Rational, b: Rational) !void {
        var r = rma;
        var aliased = rma.p.limbs.ptr == a.p.limbs.ptr or rma.p.limbs.ptr == b.p.limbs.ptr;

        var sr: Rational = undefined;
        if (aliased) {
            sr = try Rational.init(rma.p.allocator.?);
            r = &sr;
            aliased = true;
        }
        defer if (aliased) {
            rma.swap(r);
            r.deinit();
        };

        try r.p.mul(a.p, b.q);
        try r.q.mul(b.p, a.q);
        try r.p.add(r.p, r.q);

        try r.q.mul(a.q, b.q);
        try r.reduce();
    }

    // r/q = ap/aq - bp/bq = (ap*bq - bp*aq) / (aq*bq)
    //
    // For best performance, rma should not alias a or b.
    pub fn sub(rma: *Rational, a: Rational, b: Rational) !void {
        var r = rma;
        var aliased = rma.p.limbs.ptr == a.p.limbs.ptr or rma.p.limbs.ptr == b.p.limbs.ptr;

        var sr: Rational = undefined;
        if (aliased) {
            sr = try Rational.init(rma.p.allocator.?);
            r = &sr;
            aliased = true;
        }
        defer if (aliased) {
            rma.swap(r);
            r.deinit();
        };

        try r.p.mul(a.p, b.q);
        try r.q.mul(b.p, a.q);
        try r.p.sub(r.p, r.q);

        try r.q.mul(a.q, b.q);
        try r.reduce();
    }

    // r/q = ap/aq * bp/bq = ap*bp / aq*bq
    pub fn mul(r: *Rational, a: Rational, b: Rational) !void {
        try r.p.mul(a.p, b.p);
        try r.q.mul(a.q, b.q);
        try r.reduce();
    }

    // r/q = (ap/aq) / (bp/bq) = ap*bq / bp*aq
    pub fn div(r: *Rational, a: Rational, b: Rational) !void {
        if (b.p.eqZero()) {
            @panic("division by zero");
        }

        try r.p.mul(a.p, b.q);
        try r.q.mul(b.p, a.q);
        try r.reduce();
    }

    // r/q = q/r
    pub fn invert(r: *Rational) void {
        Int.swap(&r.p, &r.q);
    }

    // reduce r/q such that gcd(r, q) = 1
    fn reduce(r: *Rational) !void {
        var a = try Int.init(r.p.allocator.?);
        defer a.deinit();

        const sign = r.p.positive;

        r.p.abs();
        try gcd(&a, r.p, r.q);
        r.p.positive = sign;

        const one = Int.initFixed(([]Limb{1})[0..]);
        if (a.cmp(one) != 0) {
            var unused = try Int.init(r.p.allocator.?);
            defer unused.deinit();

            // TODO: divexact would be useful here
            // TODO: don't copy r.q for div
            try Int.divTrunc(&r.p, &unused, r.p, a);
            try Int.divTrunc(&r.q, &unused, r.q, a);
        }
    }
};

var al = debug.global_allocator;

const SignedDoubleLimb = @IntType(true, DoubleLimb.bit_count);

fn gcd(rma: *Int, x: Int, y: Int) !void {
    var r = rma;
    var aliased = rma.limbs.ptr == x.limbs.ptr or rma.limbs.ptr == y.limbs.ptr;

    var sr: Int = undefined;
    if (aliased) {
        sr = try Int.initCapacity(rma.allocator.?, math.max(x.len, y.len));
        r = &sr;
        aliased = true;
    }
    defer if (aliased) {
        rma.swap(r);
        r.deinit();
    };

    if (x.cmp(y) > 0) {
        try gcdLehmer(r, x, y);
    } else {
        try gcdLehmer(r, y, x);
    }
}

// Storage must live for the lifetime of the returned value
fn FixedIntFromSignedDoubleLimb(A: SignedDoubleLimb, storage: []Limb) Int {
    std.debug.assert(storage.len >= 2);

    var A_is_positive = A >= 0;
    const Au = @intCast(DoubleLimb, if (A < 0) -A else A);
    storage[0] = @truncate(Limb, Au);
    storage[1] = @truncate(Limb, Au >> Limb.bit_count);
    var Ap = Int.initFixed(storage[0..2]);
    Ap.positive = A_is_positive;
    return Ap;
}

// Handbook of Applied Cryptography, 14.57
//
// r = gcd(x, y) where x, y > 0
fn gcdLehmer(r: *Int, xa: Int, ya: Int) !void {
    debug.assert(xa.positive and ya.positive);
    debug.assert(xa.cmp(ya) >= 0);

    var x = try xa.clone();
    defer x.deinit();

    var y = try ya.clone();
    defer y.deinit();

    var T = try Int.init(r.allocator.?);
    defer T.deinit();

    while (y.len > 1) {
        debug.assert(x.len >= y.len);

        // chop the leading zeros of the limbs and normalize
        const offset = @clz(x.limbs[x.len - 1]);

        var xh: SignedDoubleLimb = math.shl(Limb, x.limbs[x.len - 1], offset) |
            math.shr(Limb, x.limbs[x.len - 2], Limb.bit_count - offset);

        var yh: SignedDoubleLimb = if (y.len == x.len)
            math.shl(Limb, y.limbs[y.len - 1], offset) | math.shr(Limb, y.limbs[y.len - 2], Limb.bit_count - offset)
        else if (y.len == x.len - 1)
            math.shr(Limb, y.limbs[y.len - 2], Limb.bit_count - offset)
        else
            0;

        var A: SignedDoubleLimb = 1;
        var B: SignedDoubleLimb = 0;
        var C: SignedDoubleLimb = 0;
        var D: SignedDoubleLimb = 1;

        while (yh + C != 0 and yh + D != 0) {
            const q = @divFloor(xh + A, yh + C);
            const qp = @divFloor(xh + B, yh + D);
            if (q != qp) {
                break;
            }

            var t = A - q * C;
            A = C;
            C = t;
            t = B - q * D;
            B = D;
            D = t;

            t = xh - q * yh;
            xh = yh;
            yh = t;
        }

        if (B == 0) {
            // T = x % y, r is unused
            try Int.divTrunc(r, &T, x, y);
            debug.assert(T.positive);

            x.swap(&y);
            y.swap(&T);
        } else {
            var storage: [8]Limb = undefined;
            const Ap = FixedIntFromSignedDoubleLimb(A, storage[0..2]);
            const Bp = FixedIntFromSignedDoubleLimb(B, storage[2..4]);
            const Cp = FixedIntFromSignedDoubleLimb(C, storage[4..6]);
            const Dp = FixedIntFromSignedDoubleLimb(D, storage[6..8]);

            // T = Ax + By
            try r.mul(x, Ap);
            try T.mul(y, Bp);
            try T.add(r.*, T);

            // u = Cx + Dy, r as u
            try x.mul(x, Cp);
            try r.mul(y, Dp);
            try r.add(x, r.*);

            x.swap(&T);
            y.swap(r);
        }
    }

    // euclidean algorithm
    debug.assert(x.cmp(y) >= 0);

    while (!y.eqZero()) {
        try Int.divTrunc(&T, r, x, y);
        x.swap(&y);
        y.swap(r);
    }

    r.swap(&x);
}

test "big.rational gcd non-one small" {
    var a = try Int.initSet(al, 17);
    var b = try Int.initSet(al, 97);
    var r = try Int.init(al);

    try gcd(&r, a, b);

    debug.assert((try r.to(u32)) == 1);
}

test "big.rational gcd non-one small" {
    var a = try Int.initSet(al, 4864);
    var b = try Int.initSet(al, 3458);
    var r = try Int.init(al);

    try gcd(&r, a, b);

    debug.assert((try r.to(u32)) == 38);
}

test "big.rational gcd non-one large" {
    var a = try Int.initSet(al, 0xffffffffffffffff);
    var b = try Int.initSet(al, 0xffffffffffffffff7777);
    var r = try Int.init(al);

    try gcd(&r, a, b);

    debug.assert((try r.to(u32)) == 4369);
}

test "big.rational gcd large multi-limb result" {
    var a = try Int.initSet(al, 0x12345678123456781234567812345678123456781234567812345678);
    var b = try Int.initSet(al, 0x12345671234567123456712345671234567123456712345671234567);
    var r = try Int.init(al);

    try gcd(&r, a, b);

    debug.assert((try r.to(u256)) == 0xf000000ff00000fff0000ffff000fffff00ffffff1);
}

fn extractLowBits(a: Int, comptime T: type) T {
    debug.assert(@typeId(T) == builtin.TypeId.Int);

    if (T.bit_count <= Limb.bit_count) {
        return @truncate(T, a.limbs[0]);
    } else {
        var r: T = 0;
        comptime var i: usize = 0;

        // Remainder is always 0 since if T.bit_count >= Limb.bit_count -> Limb | T and both
        // are powers of two.
        inline while (i < T.bit_count / Limb.bit_count) : (i += 1) {
            r |= math.shl(T, a.limbs[i], i * Limb.bit_count);
        }

        return r;
    }
}

test "big.rational extractLowBits" {
    var a = try Int.initSet(al, 0x11112222333344441234567887654321);

    const a1 = extractLowBits(a, u8);
    debug.assert(a1 == 0x21);

    const a2 = extractLowBits(a, u16);
    debug.assert(a2 == 0x4321);

    const a3 = extractLowBits(a, u32);
    debug.assert(a3 == 0x87654321);

    const a4 = extractLowBits(a, u64);
    debug.assert(a4 == 0x1234567887654321);

    const a5 = extractLowBits(a, u128);
    debug.assert(a5 == 0x11112222333344441234567887654321);
}

test "big.rational set" {
    var a = try Rational.init(al);

    try a.setInt(5);
    debug.assert((try a.p.to(u32)) == 5);
    debug.assert((try a.q.to(u32)) == 1);

    try a.setRatio(7, 3);
    debug.assert((try a.p.to(u32)) == 7);
    debug.assert((try a.q.to(u32)) == 3);

    try a.setRatio(9, 3);
    debug.assert((try a.p.to(i32)) == 3);
    debug.assert((try a.q.to(i32)) == 1);

    try a.setRatio(-9, 3);
    debug.assert((try a.p.to(i32)) == -3);
    debug.assert((try a.q.to(i32)) == 1);

    try a.setRatio(9, -3);
    debug.assert((try a.p.to(i32)) == -3);
    debug.assert((try a.q.to(i32)) == 1);

    try a.setRatio(-9, -3);
    debug.assert((try a.p.to(i32)) == 3);
    debug.assert((try a.q.to(i32)) == 1);
}

test "big.rational setFloat" {
    var a = try Rational.init(al);

    try a.setFloat(f64, 2.5);
    debug.assert((try a.p.to(i32)) == 5);
    debug.assert((try a.q.to(i32)) == 2);

    try a.setFloat(f32, -2.5);
    debug.assert((try a.p.to(i32)) == -5);
    debug.assert((try a.q.to(i32)) == 2);

    try a.setFloat(f32, 3.141593);

    //                = 3.14159297943115234375
    debug.assert((try a.p.to(u32)) == 3294199);
    debug.assert((try a.q.to(u32)) == 1048576);

    try a.setFloat(f64, 72.141593120712409172417410926841290461290467124);

    //                = 72.1415931207124145885245525278151035308837890625
    debug.assert((try a.p.to(u128)) == 5076513310880537);
    debug.assert((try a.q.to(u128)) == 70368744177664);
}

test "big.rational setFloatString" {
    var a = try Rational.init(al);

    try a.setFloatString("72.14159312071241458852455252781510353");

    //                  = 72.1415931207124145885245525278151035308837890625
    debug.assert((try a.p.to(u128)) == 7214159312071241458852455252781510353);
    debug.assert((try a.q.to(u128)) == 100000000000000000000000000000000000);
}

test "big.rational toFloat" {
    var a = try Rational.init(al);

    // = 3.14159297943115234375
    try a.setRatio(3294199, 1048576);
    debug.assert((try a.toFloat(f64)) == 3.14159297943115234375);

    // = 72.1415931207124145885245525278151035308837890625
    try a.setRatio(5076513310880537, 70368744177664);
    debug.assert((try a.toFloat(f64)) == 72.141593120712409172417410926841290461290467124);
}

test "big.rational set/to Float round-trip" {
    // toFloat allocates memory in a loop so we need to free it
    var buf: [512 * 1024]u8 = undefined;
    var fixed = std.heap.FixedBufferAllocator.init(buf[0..]);

    var a = try Rational.init(&fixed.allocator);

    var prng = std.rand.DefaultPrng.init(0x5EED);
    var i: usize = 0;
    while (i < 512) : (i += 1) {
        const r = prng.random.float(f64);
        try a.setFloat(f64, r);
        debug.assert((try a.toFloat(f64)) == r);
    }
}

test "big.rational copy" {
    var a = try Rational.init(al);

    const b = try Int.initSet(al, 5);

    try a.copyInt(b);
    debug.assert((try a.p.to(u32)) == 5);
    debug.assert((try a.q.to(u32)) == 1);

    const c = try Int.initSet(al, 7);
    const d = try Int.initSet(al, 3);

    try a.copyRatio(c, d);
    debug.assert((try a.p.to(u32)) == 7);
    debug.assert((try a.q.to(u32)) == 3);

    const e = try Int.initSet(al, 9);
    const f = try Int.initSet(al, 3);

    try a.copyRatio(e, f);
    debug.assert((try a.p.to(u32)) == 3);
    debug.assert((try a.q.to(u32)) == 1);
}

test "big.rational negate" {
    var a = try Rational.init(al);

    try a.setInt(-50);
    debug.assert((try a.p.to(i32)) == -50);
    debug.assert((try a.q.to(i32)) == 1);

    a.negate();
    debug.assert((try a.p.to(i32)) == 50);
    debug.assert((try a.q.to(i32)) == 1);

    a.negate();
    debug.assert((try a.p.to(i32)) == -50);
    debug.assert((try a.q.to(i32)) == 1);
}

test "big.rational abs" {
    var a = try Rational.init(al);

    try a.setInt(-50);
    debug.assert((try a.p.to(i32)) == -50);
    debug.assert((try a.q.to(i32)) == 1);

    a.abs();
    debug.assert((try a.p.to(i32)) == 50);
    debug.assert((try a.q.to(i32)) == 1);

    a.abs();
    debug.assert((try a.p.to(i32)) == 50);
    debug.assert((try a.q.to(i32)) == 1);
}

test "big.rational swap" {
    var a = try Rational.init(al);
    var b = try Rational.init(al);

    try a.setRatio(50, 23);
    try b.setRatio(17, 3);

    debug.assert((try a.p.to(u32)) == 50);
    debug.assert((try a.q.to(u32)) == 23);

    debug.assert((try b.p.to(u32)) == 17);
    debug.assert((try b.q.to(u32)) == 3);

    a.swap(&b);

    debug.assert((try a.p.to(u32)) == 17);
    debug.assert((try a.q.to(u32)) == 3);

    debug.assert((try b.p.to(u32)) == 50);
    debug.assert((try b.q.to(u32)) == 23);
}

test "big.rational cmp" {
    var a = try Rational.init(al);
    var b = try Rational.init(al);

    try a.setRatio(500, 231);
    try b.setRatio(18903, 8584);
    debug.assert((try a.cmp(b)) < 0);

    try a.setRatio(890, 10);
    try b.setRatio(89, 1);
    debug.assert((try a.cmp(b)) == 0);
}

test "big.rational add single-limb" {
    var a = try Rational.init(al);
    var b = try Rational.init(al);

    try a.setRatio(500, 231);
    try b.setRatio(18903, 8584);
    debug.assert((try a.cmp(b)) < 0);

    try a.setRatio(890, 10);
    try b.setRatio(89, 1);
    debug.assert((try a.cmp(b)) == 0);
}

test "big.rational add" {
    var a = try Rational.init(al);
    var b = try Rational.init(al);
    var r = try Rational.init(al);

    try a.setRatio(78923, 23341);
    try b.setRatio(123097, 12441414);
    try a.add(a, b);

    try r.setRatio(984786924199, 290395044174);
    debug.assert((try a.cmp(r)) == 0);
}

test "big.rational sub" {
    var a = try Rational.init(al);
    var b = try Rational.init(al);
    var r = try Rational.init(al);

    try a.setRatio(78923, 23341);
    try b.setRatio(123097, 12441414);
    try a.sub(a, b);

    try r.setRatio(979040510045, 290395044174);
    debug.assert((try a.cmp(r)) == 0);
}

test "big.rational mul" {
    var a = try Rational.init(al);
    var b = try Rational.init(al);
    var r = try Rational.init(al);

    try a.setRatio(78923, 23341);
    try b.setRatio(123097, 12441414);
    try a.mul(a, b);

    try r.setRatio(571481443, 17082061422);
    debug.assert((try a.cmp(r)) == 0);
}

test "big.rational div" {
    var a = try Rational.init(al);
    var b = try Rational.init(al);
    var r = try Rational.init(al);

    try a.setRatio(78923, 23341);
    try b.setRatio(123097, 12441414);
    try a.div(a, b);

    try r.setRatio(75531824394, 221015929);
    debug.assert((try a.cmp(r)) == 0);
}

test "big.rational div" {
    var a = try Rational.init(al);
    var r = try Rational.init(al);

    try a.setRatio(78923, 23341);
    a.invert();

    try r.setRatio(23341, 78923);
    debug.assert((try a.cmp(r)) == 0);

    try a.setRatio(-78923, 23341);
    a.invert();

    try r.setRatio(-23341, 78923);
    debug.assert((try a.cmp(r)) == 0);
}