diff --git a/CMakeLists.txt b/CMakeLists.txt index 8588c3ae9..2c34f334f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -521,6 +521,7 @@ set(ZIG_STD_FILES "math/atanh.zig" "math/big.zig" "math/big/int.zig" + "math/big/rational.zig" "math/cbrt.zig" "math/ceil.zig" "math/complex.zig" diff --git a/std/math/big.zig b/std/math/big.zig index 94b6d864e..44b5ce675 100644 --- a/std/math/big.zig +++ b/std/math/big.zig @@ -1,5 +1,7 @@ pub use @import("big/int.zig"); +pub use @import("big/rational.zig"); test "math.big" { _ = @import("big/int.zig"); + _ = @import("big/rational.zig"); } diff --git a/std/math/big/rational.zig b/std/math/big/rational.zig new file mode 100644 index 000000000..188d36d90 --- /dev/null +++ b/std/math/big/rational.zig @@ -0,0 +1,898 @@ +const std = @import("../../std.zig"); +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); +}