From ea1d2a240956b2b25d945d361fd274e67ec7849c Mon Sep 17 00:00:00 2001 From: Marc Tiehuis Date: Tue, 26 Mar 2019 19:47:26 +1300 Subject: [PATCH 1/7] Add read-only, non-allocating Int for internal constants A constant Int is one which has a value of null for its allocator field. It cannot be resized or have its limbs written. Any attempt made to write to it will be caught with a runtime panic. --- std/math/big/int.zig | 146 ++++++++++++++++++++++++------------------- 1 file changed, 83 insertions(+), 63 deletions(-) diff --git a/std/math/big/int.zig b/std/math/big/int.zig index 8800c2c7a959..7812d1ee0664 100644 --- a/std/math/big/int.zig +++ b/std/math/big/int.zig @@ -22,7 +22,7 @@ comptime { } pub const Int = struct { - allocator: *Allocator, + allocator: ?*Allocator, positive: bool, // - little-endian ordered // - len >= 1 always @@ -55,16 +55,40 @@ pub const Int = struct { }; } + // Initialize an Int directly from a fixed set of limb values. This is considered read-only + // and cannot be used as a receiver argument to any functions. If this tries to allocate + // at any point a panic will occur due to the null allocator. + pub fn initFixed(limbs: []const Limb) Int { + var self = Int{ + .allocator = null, + .positive = true, + // Cast away the const, invalid use to pass as a pointer argument. + .limbs = @intToPtr([*]Limb, @ptrToInt(limbs.ptr))[0..limbs.len], + .len = limbs.len, + }; + + self.normN(limbs.len); + return self; + } + pub fn ensureCapacity(self: *Int, capacity: usize) !void { + self.assertWritable(); if (capacity <= self.limbs.len) { return; } - self.limbs = try self.allocator.realloc(self.limbs, capacity); + self.limbs = try self.allocator.?.realloc(self.limbs, capacity); + } + + fn assertWritable(self: Int) void { + if (self.allocator == null) { + @panic("provided Int value is read-only but must be writable"); + } } pub fn deinit(self: *Int) void { - self.allocator.free(self.limbs); + self.assertWritable(); + self.allocator.?.free(self.limbs); self.* = undefined; } @@ -73,7 +97,7 @@ pub const Int = struct { .allocator = other.allocator, .positive = other.positive, .limbs = block: { - var limbs = try other.allocator.alloc(Limb, other.len); + var limbs = try other.allocator.?.alloc(Limb, other.len); mem.copy(Limb, limbs[0..], other.limbs[0..other.len]); break :block limbs; }, @@ -82,7 +106,8 @@ pub const Int = struct { } pub fn copy(self: *Int, other: Int) !void { - if (self == &other) { + self.assertWritable(); + if (self.limbs.ptr == other.limbs.ptr) { return; } @@ -93,6 +118,7 @@ pub const Int = struct { } pub fn swap(self: *Int, other: *Int) void { + self.assertWritable(); mem.swap(Int, self, other); } @@ -103,20 +129,20 @@ pub const Int = struct { debug.warn("\n"); } - pub fn negate(r: *Int) void { - r.positive = !r.positive; + pub fn negate(self: *Int) void { + self.positive = !self.positive; } - pub fn abs(r: *Int) void { - r.positive = true; + pub fn abs(self: *Int) void { + self.positive = true; } - pub fn isOdd(r: Int) bool { - return r.limbs[0] & 1 != 0; + pub fn isOdd(self: Int) bool { + return self.limbs[0] & 1 != 0; } - pub fn isEven(r: Int) bool { - return !r.isOdd(); + pub fn isEven(self: Int) bool { + return !self.isOdd(); } // Returns the number of bits required to represent the absolute value of self. @@ -179,6 +205,7 @@ pub const Int = struct { } pub fn set(self: *Int, value: var) Allocator.Error!void { + self.assertWritable(); const T = @typeOf(value); switch (@typeInfo(T)) { @@ -304,6 +331,7 @@ pub const Int = struct { } pub fn setString(self: *Int, base: u8, value: []const u8) !void { + self.assertWritable(); if (base < 2 or base > 16) { return error.InvalidBase; } @@ -315,23 +343,16 @@ pub const Int = struct { i += 1; } - // TODO values less than limb size should guarantee non allocating - var base_buffer: [512]u8 = undefined; - const base_al = &std.heap.FixedBufferAllocator.init(base_buffer[0..]).allocator; - const base_ap = try Int.initSet(base_al, base); - - var d_buffer: [512]u8 = undefined; - var d_fba = std.heap.FixedBufferAllocator.init(d_buffer[0..]); - const d_al = &d_fba.allocator; - + const ap_base = Int.initFixed(([]Limb{base})[0..]); try self.set(0); + for (value[i..]) |ch| { const d = try charToDigit(ch, base); - d_fba.end_index = 0; - const d_ap = try Int.initSet(d_al, d); - try self.mul(self.*, base_ap); - try self.add(self.*, d_ap); + const ap_d = Int.initFixed(([]Limb{d})[0..]); + + try self.mul(self.*, ap_base); + try self.add(self.*, ap_d); } self.positive = positive; } @@ -520,8 +541,19 @@ pub const Int = struct { r.len = if (j != 0) j else 1; } + // Cannot be used as a result argument to any function. + fn readOnlyPositive(a: Int) Int { + return Int{ + .allocator = null, + .positive = true, + .limbs = a.limbs, + .len = a.len, + }; + } + // r = a + b pub fn add(r: *Int, a: Int, b: Int) Allocator.Error!void { + r.assertWritable(); if (a.eqZero()) { try r.copy(b); return; @@ -533,22 +565,10 @@ pub const Int = struct { if (a.positive != b.positive) { if (a.positive) { // (a) + (-b) => a - b - const bp = Int{ - .allocator = undefined, - .positive = true, - .limbs = b.limbs, - .len = b.len, - }; - try r.sub(a, bp); + try r.sub(a, readOnlyPositive(b)); } else { // (-a) + (b) => b - a - const ap = Int{ - .allocator = undefined, - .positive = true, - .limbs = a.limbs, - .len = a.len, - }; - try r.sub(b, ap); + try r.sub(b, readOnlyPositive(a)); } } else { if (a.len >= b.len) { @@ -591,25 +611,14 @@ pub const Int = struct { // r = a - b pub fn sub(r: *Int, a: Int, b: Int) !void { + r.assertWritable(); if (a.positive != b.positive) { if (a.positive) { // (a) - (-b) => a + b - const bp = Int{ - .allocator = undefined, - .positive = true, - .limbs = b.limbs, - .len = b.len, - }; - try r.add(a, bp); + try r.add(a, readOnlyPositive(b)); } else { // (-a) - (b) => -(a + b) - const ap = Int{ - .allocator = undefined, - .positive = true, - .limbs = a.limbs, - .len = a.len, - }; - try r.add(ap, b); + try r.add(readOnlyPositive(a), b); r.positive = false; } } else { @@ -671,12 +680,14 @@ pub const Int = struct { // // For greatest efficiency, ensure rma does not alias a or b. pub fn mul(rma: *Int, a: Int, b: Int) !void { + rma.assertWritable(); + var r = rma; var aliased = rma.limbs.ptr == a.limbs.ptr or rma.limbs.ptr == b.limbs.ptr; var sr: Int = undefined; if (aliased) { - sr = try Int.initCapacity(rma.allocator, a.len + b.len); + sr = try Int.initCapacity(rma.allocator.?, a.len + b.len); r = &sr; aliased = true; } @@ -745,13 +756,9 @@ pub const Int = struct { // Trunc -> Floor. if (!q.positive) { - // TODO values less than limb size should guarantee non allocating - var one_buffer: [512]u8 = undefined; - const one_al = &std.heap.FixedBufferAllocator.init(one_buffer[0..]).allocator; - const one_ap = try Int.initSet(one_al, 1); - - try q.sub(q.*, one_ap); - try r.add(q.*, one_ap); + const one = Int.initFixed(([]Limb{1})[0..]); + try q.sub(q.*, one); + try r.add(q.*, one); } r.positive = b.positive; } @@ -763,6 +770,9 @@ pub const Int = struct { // Truncates by default. fn div(quo: *Int, rem: *Int, a: Int, b: Int) !void { + quo.assertWritable(); + rem.assertWritable(); + if (b.eqZero()) { @panic("division by zero"); } @@ -800,7 +810,7 @@ pub const Int = struct { // x may grow one limb during normalization try quo.ensureCapacity(a.len + y.len); - try divN(quo.allocator, quo, rem, &x, &y); + try divN(quo.allocator.?, quo, rem, &x, &y); quo.positive = a.positive == b.positive; } @@ -919,6 +929,8 @@ pub const Int = struct { // r = a << shift, in other words, r = a * 2^shift pub fn shiftLeft(r: *Int, a: Int, shift: usize) !void { + r.assertWritable(); + try r.ensureCapacity(a.len + (shift / Limb.bit_count) + 1); llshl(r.limbs[0..], a.limbs[0..a.len], shift); r.norm1(a.len + (shift / Limb.bit_count) + 1); @@ -950,6 +962,8 @@ pub const Int = struct { // r = a >> shift pub fn shiftRight(r: *Int, a: Int, shift: usize) !void { + r.assertWritable(); + if (a.len <= shift / Limb.bit_count) { r.len = 1; r.limbs[0] = 0; @@ -985,6 +999,8 @@ pub const Int = struct { // r = a | b pub fn bitOr(r: *Int, a: Int, b: Int) !void { + r.assertWritable(); + if (a.len > b.len) { try r.ensureCapacity(a.len); llor(r.limbs[0..], a.limbs[0..a.len], b.limbs[0..b.len]); @@ -1012,6 +1028,8 @@ pub const Int = struct { // r = a & b pub fn bitAnd(r: *Int, a: Int, b: Int) !void { + r.assertWritable(); + if (a.len > b.len) { try r.ensureCapacity(b.len); lland(r.limbs[0..], a.limbs[0..a.len], b.limbs[0..b.len]); @@ -1036,6 +1054,8 @@ pub const Int = struct { // r = a ^ b pub fn bitXor(r: *Int, a: Int, b: Int) !void { + r.assertWritable(); + if (a.len > b.len) { try r.ensureCapacity(a.len); llxor(r.limbs[0..], a.limbs[0..a.len], b.limbs[0..b.len]); From 5f4fcd5030fa93cf2518f98ffc666e4e4dc1052b Mon Sep 17 00:00:00 2001 From: Marc Tiehuis Date: Tue, 26 Mar 2019 19:53:02 +1300 Subject: [PATCH 2/7] Add initial big.Rational type --- CMakeLists.txt | 1 + std/math/big.zig | 2 + std/math/big/rational.zig | 898 ++++++++++++++++++++++++++++++++++++++ 3 files changed, 901 insertions(+) create mode 100644 std/math/big/rational.zig diff --git a/CMakeLists.txt b/CMakeLists.txt index 8588c3ae9ef7..2c34f334f97c 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 94b6d864e783..44b5ce675fad 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 000000000000..188d36d909ba --- /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); +} From b3ecdfd7bfb9c6b3bd085b489fe20dd0ca1e12d8 Mon Sep 17 00:00:00 2001 From: Marc Tiehuis Date: Tue, 26 Mar 2019 20:31:16 +1300 Subject: [PATCH 3/7] Fix big.Int toString maybe-null allocator --- std/math/big/int.zig | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/std/math/big/int.zig b/std/math/big/int.zig index 7812d1ee0664..89be663a4a5c 100644 --- a/std/math/big/int.zig +++ b/std/math/big/int.zig @@ -449,9 +449,11 @@ pub const Int = struct { comptime FmtError: type, output: fn (@typeOf(context), []const u8) FmtError!void, ) FmtError!void { + self.assertWritable(); // TODO look at fmt and support other bases - const str = self.toString(self.allocator, 10) catch @panic("TODO make this non allocating"); - defer self.allocator.free(str); + // TODO support read-only fixed integers + const str = self.toString(self.allocator.?, 10) catch @panic("TODO make this non allocating"); + defer self.allocator.?.free(str); return output(context, str); } From d3e1f323626ba5955aa98628c9c817a19baa35fb Mon Sep 17 00:00:00 2001 From: Marc Tiehuis Date: Wed, 27 Mar 2019 22:28:38 +1300 Subject: [PATCH 4/7] Small fixes for big.Rational and corrections for gcdLehmer The int div code still causes some edge cases to fail right now. --- std/math/big/int.zig | 7 +- std/math/big/rational.zig | 160 ++++++++++++++++++++------------------ 2 files changed, 91 insertions(+), 76 deletions(-) diff --git a/std/math/big/int.zig b/std/math/big/int.zig index 89be663a4a5c..6f4ce1b54ed9 100644 --- a/std/math/big/int.zig +++ b/std/math/big/int.zig @@ -93,6 +93,7 @@ pub const Int = struct { } pub fn clone(other: Int) !Int { + other.assertWritable(); return Int{ .allocator = other.allocator, .positive = other.positive, @@ -804,11 +805,13 @@ pub const Int = struct { rem.positive = true; } else { // x and y are modified during division - var x = try a.clone(); + var x = try Int.initCapacity(quo.allocator.?, a.len); defer x.deinit(); + try x.copy(a); - var y = try b.clone(); + var y = try Int.initCapacity(quo.allocator.?, b.len); defer y.deinit(); + try y.copy(b); // x may grow one limb during normalization try quo.ensureCapacity(a.len + y.len); diff --git a/std/math/big/rational.zig b/std/math/big/rational.zig index 188d36d909ba..ca3b3a5926fb 100644 --- a/std/math/big/rational.zig +++ b/std/math/big/rational.zig @@ -3,6 +3,7 @@ const builtin = @import("builtin"); const debug = std.debug; const math = std.math; const mem = std.mem; +const testing = std.testing; const Allocator = mem.Allocator; const ArrayList = std.ArrayList; @@ -425,6 +426,7 @@ var al = debug.global_allocator; const SignedDoubleLimb = @IntType(true, DoubleLimb.bit_count); fn gcd(rma: *Int, x: Int, y: Int) !void { + rma.assertWritable(); var r = rma; var aliased = rma.limbs.ptr == x.limbs.ptr or rma.limbs.ptr == y.limbs.ptr; @@ -463,19 +465,23 @@ fn FixedIntFromSignedDoubleLimb(A: SignedDoubleLimb, storage: []Limb) Int { // // 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(); + x.abs(); defer x.deinit(); var y = try ya.clone(); + y.abs(); defer y.deinit(); + if (x.cmp(y) < 0) { + x.swap(&y); + } + var T = try Int.init(r.allocator.?); defer T.deinit(); while (y.len > 1) { + debug.assert(x.positive and y.positive); debug.assert(x.len >= y.len); // chop the leading zeros of the limbs and normalize @@ -540,7 +546,13 @@ fn gcdLehmer(r: *Int, xa: Int, ya: Int) !void { try r.add(x, r.*); x.swap(&T); + x.abs(); y.swap(r); + y.abs(); + + if (x.cmp(y) < 0) { + x.swap(&y); + } } } @@ -563,7 +575,7 @@ test "big.rational gcd non-one small" { try gcd(&r, a, b); - debug.assert((try r.to(u32)) == 1); + testing.expect((try r.to(u32)) == 1); } test "big.rational gcd non-one small" { @@ -573,7 +585,7 @@ test "big.rational gcd non-one small" { try gcd(&r, a, b); - debug.assert((try r.to(u32)) == 38); + testing.expect((try r.to(u32)) == 38); } test "big.rational gcd non-one large" { @@ -583,7 +595,7 @@ test "big.rational gcd non-one large" { try gcd(&r, a, b); - debug.assert((try r.to(u32)) == 4369); + testing.expect((try r.to(u32)) == 4369); } test "big.rational gcd large multi-limb result" { @@ -593,11 +605,11 @@ test "big.rational gcd large multi-limb result" { try gcd(&r, a, b); - debug.assert((try r.to(u256)) == 0xf000000ff00000fff0000ffff000fffff00ffffff1); + testing.expect((try r.to(u256)) == 0xf000000ff00000fff0000ffff000fffff00ffffff1); } fn extractLowBits(a: Int, comptime T: type) T { - debug.assert(@typeId(T) == builtin.TypeId.Int); + testing.expect(@typeId(T) == builtin.TypeId.Int); if (T.bit_count <= Limb.bit_count) { return @truncate(T, a.limbs[0]); @@ -619,71 +631,71 @@ test "big.rational extractLowBits" { var a = try Int.initSet(al, 0x11112222333344441234567887654321); const a1 = extractLowBits(a, u8); - debug.assert(a1 == 0x21); + testing.expect(a1 == 0x21); const a2 = extractLowBits(a, u16); - debug.assert(a2 == 0x4321); + testing.expect(a2 == 0x4321); const a3 = extractLowBits(a, u32); - debug.assert(a3 == 0x87654321); + testing.expect(a3 == 0x87654321); const a4 = extractLowBits(a, u64); - debug.assert(a4 == 0x1234567887654321); + testing.expect(a4 == 0x1234567887654321); const a5 = extractLowBits(a, u128); - debug.assert(a5 == 0x11112222333344441234567887654321); + testing.expect(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); + testing.expect((try a.p.to(u32)) == 5); + testing.expect((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); + testing.expect((try a.p.to(u32)) == 7); + testing.expect((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); + testing.expect((try a.p.to(i32)) == 3); + testing.expect((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); + testing.expect((try a.p.to(i32)) == -3); + testing.expect((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); + testing.expect((try a.p.to(i32)) == -3); + testing.expect((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); + testing.expect((try a.p.to(i32)) == 3); + testing.expect((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); + testing.expect((try a.p.to(i32)) == 5); + testing.expect((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); + testing.expect((try a.p.to(i32)) == -5); + testing.expect((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); + testing.expect((try a.p.to(u32)) == 3294199); + testing.expect((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); + testing.expect((try a.p.to(u128)) == 5076513310880537); + testing.expect((try a.q.to(u128)) == 70368744177664); } test "big.rational setFloatString" { @@ -692,8 +704,8 @@ test "big.rational setFloatString" { try a.setFloatString("72.14159312071241458852455252781510353"); // = 72.1415931207124145885245525278151035308837890625 - debug.assert((try a.p.to(u128)) == 7214159312071241458852455252781510353); - debug.assert((try a.q.to(u128)) == 100000000000000000000000000000000000); + testing.expect((try a.p.to(u128)) == 7214159312071241458852455252781510353); + testing.expect((try a.q.to(u128)) == 100000000000000000000000000000000000); } test "big.rational toFloat" { @@ -701,11 +713,11 @@ test "big.rational toFloat" { // = 3.14159297943115234375 try a.setRatio(3294199, 1048576); - debug.assert((try a.toFloat(f64)) == 3.14159297943115234375); + testing.expect((try a.toFloat(f64)) == 3.14159297943115234375); // = 72.1415931207124145885245525278151035308837890625 try a.setRatio(5076513310880537, 70368744177664); - debug.assert((try a.toFloat(f64)) == 72.141593120712409172417410926841290461290467124); + testing.expect((try a.toFloat(f64)) == 72.141593120712409172417410926841290461290467124); } test "big.rational set/to Float round-trip" { @@ -720,7 +732,7 @@ test "big.rational set/to Float round-trip" { while (i < 512) : (i += 1) { const r = prng.random.float(f64); try a.setFloat(f64, r); - debug.assert((try a.toFloat(f64)) == r); + testing.expect((try a.toFloat(f64)) == r); } } @@ -730,54 +742,54 @@ test "big.rational copy" { 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); + testing.expect((try a.p.to(u32)) == 5); + testing.expect((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); + testing.expect((try a.p.to(u32)) == 7); + testing.expect((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); + testing.expect((try a.p.to(u32)) == 3); + testing.expect((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); + testing.expect((try a.p.to(i32)) == -50); + testing.expect((try a.q.to(i32)) == 1); a.negate(); - debug.assert((try a.p.to(i32)) == 50); - debug.assert((try a.q.to(i32)) == 1); + testing.expect((try a.p.to(i32)) == 50); + testing.expect((try a.q.to(i32)) == 1); a.negate(); - debug.assert((try a.p.to(i32)) == -50); - debug.assert((try a.q.to(i32)) == 1); + testing.expect((try a.p.to(i32)) == -50); + testing.expect((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); + testing.expect((try a.p.to(i32)) == -50); + testing.expect((try a.q.to(i32)) == 1); a.abs(); - debug.assert((try a.p.to(i32)) == 50); - debug.assert((try a.q.to(i32)) == 1); + testing.expect((try a.p.to(i32)) == 50); + testing.expect((try a.q.to(i32)) == 1); a.abs(); - debug.assert((try a.p.to(i32)) == 50); - debug.assert((try a.q.to(i32)) == 1); + testing.expect((try a.p.to(i32)) == 50); + testing.expect((try a.q.to(i32)) == 1); } test "big.rational swap" { @@ -787,19 +799,19 @@ test "big.rational swap" { 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); + testing.expect((try a.p.to(u32)) == 50); + testing.expect((try a.q.to(u32)) == 23); - debug.assert((try b.p.to(u32)) == 17); - debug.assert((try b.q.to(u32)) == 3); + testing.expect((try b.p.to(u32)) == 17); + testing.expect((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); + testing.expect((try a.p.to(u32)) == 17); + testing.expect((try a.q.to(u32)) == 3); - debug.assert((try b.p.to(u32)) == 50); - debug.assert((try b.q.to(u32)) == 23); + testing.expect((try b.p.to(u32)) == 50); + testing.expect((try b.q.to(u32)) == 23); } test "big.rational cmp" { @@ -808,11 +820,11 @@ test "big.rational cmp" { try a.setRatio(500, 231); try b.setRatio(18903, 8584); - debug.assert((try a.cmp(b)) < 0); + testing.expect((try a.cmp(b)) < 0); try a.setRatio(890, 10); try b.setRatio(89, 1); - debug.assert((try a.cmp(b)) == 0); + testing.expect((try a.cmp(b)) == 0); } test "big.rational add single-limb" { @@ -821,11 +833,11 @@ test "big.rational add single-limb" { try a.setRatio(500, 231); try b.setRatio(18903, 8584); - debug.assert((try a.cmp(b)) < 0); + testing.expect((try a.cmp(b)) < 0); try a.setRatio(890, 10); try b.setRatio(89, 1); - debug.assert((try a.cmp(b)) == 0); + testing.expect((try a.cmp(b)) == 0); } test "big.rational add" { @@ -838,7 +850,7 @@ test "big.rational add" { try a.add(a, b); try r.setRatio(984786924199, 290395044174); - debug.assert((try a.cmp(r)) == 0); + testing.expect((try a.cmp(r)) == 0); } test "big.rational sub" { @@ -851,7 +863,7 @@ test "big.rational sub" { try a.sub(a, b); try r.setRatio(979040510045, 290395044174); - debug.assert((try a.cmp(r)) == 0); + testing.expect((try a.cmp(r)) == 0); } test "big.rational mul" { @@ -864,7 +876,7 @@ test "big.rational mul" { try a.mul(a, b); try r.setRatio(571481443, 17082061422); - debug.assert((try a.cmp(r)) == 0); + testing.expect((try a.cmp(r)) == 0); } test "big.rational div" { @@ -877,7 +889,7 @@ test "big.rational div" { try a.div(a, b); try r.setRatio(75531824394, 221015929); - debug.assert((try a.cmp(r)) == 0); + testing.expect((try a.cmp(r)) == 0); } test "big.rational div" { @@ -888,11 +900,11 @@ test "big.rational div" { a.invert(); try r.setRatio(23341, 78923); - debug.assert((try a.cmp(r)) == 0); + testing.expect((try a.cmp(r)) == 0); try a.setRatio(-78923, 23341); a.invert(); try r.setRatio(-23341, 78923); - debug.assert((try a.cmp(r)) == 0); + testing.expect((try a.cmp(r)) == 0); } From 30788a98b1bb58886d409464baf7eb01e428d939 Mon Sep 17 00:00:00 2001 From: Marc Tiehuis Date: Thu, 28 Mar 2019 20:39:57 +1300 Subject: [PATCH 5/7] Handle zero-limb trailing div case in big.Int --- std/math/big/int.zig | 100 ++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 95 insertions(+), 5 deletions(-) diff --git a/std/math/big/int.zig b/std/math/big/int.zig index 6f4ce1b54ed9..5b84c460da86 100644 --- a/std/math/big/int.zig +++ b/std/math/big/int.zig @@ -794,11 +794,30 @@ pub const Int = struct { return; } - if (b.len == 1) { + // Handle trailing zero-words of divisor/dividend. These are not handled in the following + // algorithms. + const a_zero_limb_count = blk: { + var i: usize = 0; + while (i < a.len) : (i += 1) { + if (a.limbs[i] != 0) break; + } + break :blk i; + }; + const b_zero_limb_count = blk: { + var i: usize = 0; + while (i < b.len) : (i += 1) { + if (b.limbs[i] != 0) break; + } + break :blk i; + }; + + const ab_zero_limb_count = std.math.min(a_zero_limb_count, b_zero_limb_count); + + if (b.len - ab_zero_limb_count == 1) { try quo.ensureCapacity(a.len); - lldiv1(quo.limbs[0..], &rem.limbs[0], a.limbs[0..a.len], b.limbs[0]); - quo.norm1(a.len); + lldiv1(quo.limbs[0..], &rem.limbs[0], a.limbs[ab_zero_limb_count..a.len], b.limbs[b.len - 1]); + quo.norm1(a.len - ab_zero_limb_count); quo.positive = a.positive == b.positive; rem.len = 1; @@ -815,9 +834,25 @@ pub const Int = struct { // x may grow one limb during normalization try quo.ensureCapacity(a.len + y.len); + + // Shrink x, y such that the trailing zero limbs shared between are removed. + if (ab_zero_limb_count != 0) { + std.mem.copy(Limb, x.limbs[0..], x.limbs[ab_zero_limb_count..]); + std.mem.copy(Limb, y.limbs[0..], y.limbs[ab_zero_limb_count..]); + x.len -= ab_zero_limb_count; + y.len -= ab_zero_limb_count; + } + try divN(quo.allocator.?, quo, rem, &x, &y); quo.positive = a.positive == b.positive; + + // If dividend had trailing zeros beyond divisor, add extra trailing limbs. + // Single-limb division never has multi-limb remainder so nothing to add. + if (a_zero_limb_count > b_zero_limb_count) { + const shift = a_zero_limb_count - b_zero_limb_count; + try rem.shiftLeft(rem.*, shift * Limb.bit_count); + } } } @@ -860,8 +895,11 @@ pub const Int = struct { var tmp = try Int.init(allocator); defer tmp.deinit(); - // Normalize so y > Limb.bit_count / 2 (i.e. leading bit is set) - const norm_shift = @clz(y.limbs[y.len - 1]); + // Normalize so y > Limb.bit_count / 2 (i.e. leading bit is set) and even + var norm_shift = @clz(y.limbs[y.len - 1]); + if (norm_shift == 0 and y.isOdd()) { + norm_shift = Limb.bit_count; + } try x.shiftLeft(x.*, norm_shift); try y.shiftLeft(y.*, norm_shift); @@ -2005,6 +2043,58 @@ test "big.int div multi-multi (3.1/3.3 branch)" { testing.expect((try r.to(u256)) == 0x1111111111111111111110b12222222222222222282); } +test "big.int div multi-single zero-limb trailing" { + var a = try Int.initSet(al, 0x60000000000000000000000000000000000000000000000000000000000000000); + var b = try Int.initSet(al, 0x10000000000000000); + + var q = try Int.init(al); + var r = try Int.init(al); + try Int.divTrunc(&q, &r, a, b); + + var expected = try Int.initSet(al, 0x6000000000000000000000000000000000000000000000000); + testing.expect(q.eq(expected)); + testing.expect(r.eqZero()); +} + +test "big.int div multi-multi zero-limb trailing (with rem)" { + var a = try Int.initSet(al, 0x86666666555555558888888777777776111111111111111100000000000000000000000000000000); + var b = try Int.initSet(al, 0x8666666655555555444444443333333300000000000000000000000000000000); + + var q = try Int.init(al); + var r = try Int.init(al); + try Int.divTrunc(&q, &r, a, b); + + testing.expect((try q.to(u128)) == 0x10000000000000000); + testing.expect((try r.to(u128)) == 0x44444443444444431111111111111111); +} + +test "big.int div multi-multi zero-limb trailing (with rem) and dividend zero-limb count > divisor zero-limb count" { + var a = try Int.initSet(al, 0x8666666655555555888888877777777611111111111111110000000000000000); + var b = try Int.initSet(al, 0x8666666655555555444444443333333300000000000000000000000000000000); + + var q = try Int.init(al); + var r = try Int.init(al); + try Int.divTrunc(&q, &r, a, b); + + testing.expect((try q.to(u128)) == 0x1); + testing.expect((try r.to(u128)) == 0x44444443444444431111111111111111); +} + +test "big.int div multi-multi zero-limb trailing (with rem) and dividend zero-limb count < divisor zero-limb count" { + var a = try Int.initSet(al, 0x86666666555555558888888777777776111111111111111100000000000000000000000000000000); + var b = try Int.initSet(al, 0x866666665555555544444444333333330000000000000000); + + var q = try Int.init(al); + var r = try Int.init(al); + try Int.divTrunc(&q, &r, a, b); + + const qs = try q.toString(al, 16); + testing.expect(std.mem.eql(u8, qs, "10000000000000000820820803105186f")); + + const rs = try r.toString(al, 16); + testing.expect(std.mem.eql(u8, rs, "4e11f2baa5896a321d463b543d0104e30000000000000000")); +} + test "big.int shift-right single" { var a = try Int.initSet(al, 0xffff0000); try a.shiftRight(a, 16); From 87d8ecda462688c597c726b1da5dbd5f8478e0fc Mon Sep 17 00:00:00 2001 From: Marc Tiehuis Date: Wed, 3 Apr 2019 17:20:23 +1300 Subject: [PATCH 6/7] Fix math.big.Int divN/gcdLehmer and fuzz-test failures --- std/math/big/int.zig | 131 ++++++++++++++++++++++---------------- std/math/big/rational.zig | 53 ++++++--------- 2 files changed, 96 insertions(+), 88 deletions(-) diff --git a/std/math/big/int.zig b/std/math/big/int.zig index 5b84c460da86..bd17ed16f727 100644 --- a/std/math/big/int.zig +++ b/std/math/big/int.zig @@ -67,7 +67,7 @@ pub const Int = struct { .len = limbs.len, }; - self.normN(limbs.len); + self.normalize(limbs.len); return self; } @@ -508,28 +508,12 @@ pub const Int = struct { return cmp(a, b) == 0; } - // Normalize for a possible single carry digit. - // - // [1, 2, 3, 4, 0] -> [1, 2, 3, 4] - // [1, 2, 3, 4, 5] -> [1, 2, 3, 4, 5] - // [0] -> [0] - fn norm1(r: *Int, length: usize) void { - debug.assert(length > 0); - debug.assert(length <= r.limbs.len); - - if (r.limbs[length - 1] == 0) { - r.len = if (length > 1) length - 1 else 1; - } else { - r.len = length; - } - } - // Normalize a possible sequence of leading zeros. // // [1, 2, 3, 4, 0] -> [1, 2, 3, 4] // [1, 2, 0, 0, 0] -> [1, 2] // [0, 0, 0, 0, 0] -> [0] - fn normN(r: *Int, length: usize) void { + fn normalize(r: *Int, length: usize) void { debug.assert(length > 0); debug.assert(length <= r.limbs.len); @@ -577,11 +561,11 @@ pub const Int = struct { if (a.len >= b.len) { try r.ensureCapacity(a.len + 1); lladd(r.limbs[0..], a.limbs[0..a.len], b.limbs[0..b.len]); - r.norm1(a.len + 1); + r.normalize(a.len + 1); } else { try r.ensureCapacity(b.len + 1); lladd(r.limbs[0..], b.limbs[0..b.len], a.limbs[0..a.len]); - r.norm1(b.len + 1); + r.normalize(b.len + 1); } r.positive = a.positive; @@ -630,12 +614,12 @@ pub const Int = struct { if (a.cmp(b) >= 0) { try r.ensureCapacity(a.len + 1); llsub(r.limbs[0..], a.limbs[0..a.len], b.limbs[0..b.len]); - r.normN(a.len); + r.normalize(a.len); r.positive = true; } else { try r.ensureCapacity(b.len + 1); llsub(r.limbs[0..], b.limbs[0..b.len], a.limbs[0..a.len]); - r.normN(b.len); + r.normalize(b.len); r.positive = false; } } else { @@ -643,12 +627,12 @@ pub const Int = struct { if (a.cmp(b) < 0) { try r.ensureCapacity(a.len + 1); llsub(r.limbs[0..], a.limbs[0..a.len], b.limbs[0..b.len]); - r.normN(a.len); + r.normalize(a.len); r.positive = false; } else { try r.ensureCapacity(b.len + 1); llsub(r.limbs[0..], b.limbs[0..b.len], a.limbs[0..a.len]); - r.normN(b.len); + r.normalize(b.len); r.positive = true; } } @@ -708,7 +692,7 @@ pub const Int = struct { } r.positive = a.positive == b.positive; - r.normN(a.len + b.len); + r.normalize(a.len + b.len); } // a + b * c + *carry, sets carry to the overflow bits @@ -817,7 +801,7 @@ pub const Int = struct { try quo.ensureCapacity(a.len); lldiv1(quo.limbs[0..], &rem.limbs[0], a.limbs[ab_zero_limb_count..a.len], b.limbs[b.len - 1]); - quo.norm1(a.len - ab_zero_limb_count); + quo.normalize(a.len - ab_zero_limb_count); quo.positive = a.positive == b.positive; rem.len = 1; @@ -846,13 +830,10 @@ pub const Int = struct { try divN(quo.allocator.?, quo, rem, &x, &y); quo.positive = a.positive == b.positive; + } - // If dividend had trailing zeros beyond divisor, add extra trailing limbs. - // Single-limb division never has multi-limb remainder so nothing to add. - if (a_zero_limb_count > b_zero_limb_count) { - const shift = a_zero_limb_count - b_zero_limb_count; - try rem.shiftLeft(rem.*, shift * Limb.bit_count); - } + if (ab_zero_limb_count != 0) { + try rem.shiftLeft(rem.*, ab_zero_limb_count * Limb.bit_count); } } @@ -933,7 +914,7 @@ pub const Int = struct { tmp.limbs[0] = if (i >= 2) x.limbs[i - 2] else 0; tmp.limbs[1] = if (i >= 1) x.limbs[i - 1] else 0; tmp.limbs[2] = x.limbs[i]; - tmp.normN(3); + tmp.normalize(3); while (true) { // 2x1 limb multiplication unrolled against single-limb q[i-t-1] @@ -941,7 +922,7 @@ pub const Int = struct { r.limbs[0] = addMulLimbWithCarry(0, if (t >= 1) y.limbs[t - 1] else 0, q.limbs[i - t - 1], &carry); r.limbs[1] = addMulLimbWithCarry(0, y.limbs[t], q.limbs[i - t - 1], &carry); r.limbs[2] = carry; - r.normN(3); + r.normalize(3); if (r.cmpAbs(tmp) <= 0) { break; @@ -964,10 +945,10 @@ pub const Int = struct { } // Denormalize - q.normN(q.len); + q.normalize(q.len); try r.shiftRight(x.*, norm_shift); - r.normN(r.len); + r.normalize(r.len); } // r = a << shift, in other words, r = a * 2^shift @@ -976,7 +957,7 @@ pub const Int = struct { try r.ensureCapacity(a.len + (shift / Limb.bit_count) + 1); llshl(r.limbs[0..], a.limbs[0..a.len], shift); - r.norm1(a.len + (shift / Limb.bit_count) + 1); + r.normalize(a.len + (shift / Limb.bit_count) + 1); r.positive = a.positive; } @@ -1076,11 +1057,11 @@ pub const Int = struct { if (a.len > b.len) { try r.ensureCapacity(b.len); lland(r.limbs[0..], a.limbs[0..a.len], b.limbs[0..b.len]); - r.normN(b.len); + r.normalize(b.len); } else { try r.ensureCapacity(a.len); lland(r.limbs[0..], b.limbs[0..b.len], a.limbs[0..a.len]); - r.normN(a.len); + r.normalize(a.len); } } @@ -1102,11 +1083,11 @@ pub const Int = struct { if (a.len > b.len) { try r.ensureCapacity(a.len); llxor(r.limbs[0..], a.limbs[0..a.len], b.limbs[0..b.len]); - r.normN(a.len); + r.normalize(a.len); } else { try r.ensureCapacity(b.len); llxor(r.limbs[0..], b.limbs[0..b.len], a.limbs[0..a.len]); - r.normN(b.len); + r.normalize(b.len); } } @@ -1130,7 +1111,9 @@ pub const Int = struct { // They will still run on larger than this and should pass, but the multi-limb code-paths // may be untested in some cases. -const al = debug.global_allocator; +var buffer: [64 * 8192]u8 = undefined; +var fixed = std.heap.FixedBufferAllocator.init(buffer[0..]); +const al = &fixed.allocator; test "big.int comptime_int set" { comptime var s = 0xefffffff00000001eeeeeeefaaaaaaab; @@ -1179,7 +1162,7 @@ test "big.int to target too small error" { testing.expectError(error.TargetTooSmall, a.to(u8)); } -test "big.int norm1" { +test "big.int normalize" { var a = try Int.init(al); try a.ensureCapacity(8); @@ -1187,26 +1170,26 @@ test "big.int norm1" { a.limbs[1] = 2; a.limbs[2] = 3; a.limbs[3] = 0; - a.norm1(4); + a.normalize(4); testing.expect(a.len == 3); a.limbs[0] = 1; a.limbs[1] = 2; a.limbs[2] = 3; - a.norm1(3); + a.normalize(3); testing.expect(a.len == 3); a.limbs[0] = 0; a.limbs[1] = 0; - a.norm1(2); + a.normalize(2); testing.expect(a.len == 1); a.limbs[0] = 0; - a.norm1(1); + a.normalize(1); testing.expect(a.len == 1); } -test "big.int normN" { +test "big.int normalize multi" { var a = try Int.init(al); try a.ensureCapacity(8); @@ -1214,24 +1197,24 @@ test "big.int normN" { a.limbs[1] = 2; a.limbs[2] = 0; a.limbs[3] = 0; - a.normN(4); + a.normalize(4); testing.expect(a.len == 2); a.limbs[0] = 1; a.limbs[1] = 2; a.limbs[2] = 3; - a.normN(3); + a.normalize(3); testing.expect(a.len == 3); a.limbs[0] = 0; a.limbs[1] = 0; a.limbs[2] = 0; a.limbs[3] = 0; - a.normN(4); + a.normalize(4); testing.expect(a.len == 1); a.limbs[0] = 0; - a.normN(1); + a.normalize(1); testing.expect(a.len == 1); } @@ -2065,7 +2048,9 @@ test "big.int div multi-multi zero-limb trailing (with rem)" { try Int.divTrunc(&q, &r, a, b); testing.expect((try q.to(u128)) == 0x10000000000000000); - testing.expect((try r.to(u128)) == 0x44444443444444431111111111111111); + + const rs = try r.toString(al, 16); + testing.expect(std.mem.eql(u8, rs, "4444444344444443111111111111111100000000000000000000000000000000")); } test "big.int div multi-multi zero-limb trailing (with rem) and dividend zero-limb count > divisor zero-limb count" { @@ -2077,7 +2062,9 @@ test "big.int div multi-multi zero-limb trailing (with rem) and dividend zero-li try Int.divTrunc(&q, &r, a, b); testing.expect((try q.to(u128)) == 0x1); - testing.expect((try r.to(u128)) == 0x44444443444444431111111111111111); + + const rs = try r.toString(al, 16); + testing.expect(std.mem.eql(u8, rs, "444444434444444311111111111111110000000000000000")); } test "big.int div multi-multi zero-limb trailing (with rem) and dividend zero-limb count < divisor zero-limb count" { @@ -2095,6 +2082,42 @@ test "big.int div multi-multi zero-limb trailing (with rem) and dividend zero-li testing.expect(std.mem.eql(u8, rs, "4e11f2baa5896a321d463b543d0104e30000000000000000")); } +test "big.int div multi-multi fuzz case #1" { + var a = try Int.init(al); + var b = try Int.init(al); + + try a.setString(16, "ffffffffffffffffffffffffffffc00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000003ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff80000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000003ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffc00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffe000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000"); + try b.setString(16, "3ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffe0000000000000000000000000000000000001ffffffffffffffffffffffffffffffffffffffffffffffffffc000000000000000000000000000000007fffffffffff"); + + var q = try Int.init(al); + var r = try Int.init(al); + try Int.divTrunc(&q, &r, a, b); + + const qs = try q.toString(al, 16); + testing.expect(std.mem.eql(u8, qs, "3ffffffffffffffffffffffffffff0000000000000000000000000000000000001ffffffffffffffffffffffffffff7fffffffe000000000000000000000000000180000000000000000000003fffffbfffffffdfffffffffffffeffff800000100101000000100000000020003fffffdfbfffffe3ffffffffffffeffff7fffc00800a100000017ffe000002000400007efbfff7fe9f00000037ffff3fff7fffa004006100000009ffe00000190038200bf7d2ff7fefe80400060000f7d7f8fbf9401fe38e0403ffc0bdffffa51102c300d7be5ef9df4e5060007b0127ad3fa69f97d0f820b6605ff617ddf7f32ad7a05c0d03f2e7bc78a6000e087a8bbcdc59e07a5a079128a7861f553ddebed7e8e56701756f9ead39b48cd1b0831889ea6ec1fddf643d0565b075ff07e6caea4e2854ec9227fd635ed60a2f5eef2893052ffd54718fa08604acbf6a15e78a467c4a3c53c0278af06c4416573f925491b195e8fd79302cb1aaf7caf4ecfc9aec1254cc969786363ac729f914c6ddcc26738d6b0facd54eba026580aba2eb6482a088b0d224a8852420b91ec1")); + + const rs = try r.toString(al, 16); + testing.expect(std.mem.eql(u8, rs, "310d1d4c414426b4836c2635bad1df3a424e50cbdd167ffccb4dfff57d36b4aae0d6ca0910698220171a0f3373c1060a046c2812f0027e321f72979daa5e7973214170d49e885de0c0ecc167837d44502430674a82522e5df6a0759548052420b91ec1")); +} + +test "big.int div multi-multi fuzz case #2" { + var a = try Int.init(al); + var b = try Int.init(al); + + try a.setString(16, "3ffffffffe00000000000000000000000000fffffffffffffffffffffffffffffffffffffffffffffffffffffffffe000000000000000000000000000000000000000000000000000000000000001fffffffffffffffff800000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000ffffffffffffffffffffc000000000000000000000000000000000000000000000000000000000000000"); + try b.setString(16, "ffc0000000000000000000000000000000000000000000000000"); + + var q = try Int.init(al); + var r = try Int.init(al); + try Int.divTrunc(&q, &r, a, b); + + const qs = try q.toString(al, 16); + testing.expect(std.mem.eql(u8, qs, "40100400fe3f8fe3f8fe3f8fe3f8fe3f8fe4f93e4f93e4f93e4f93e4f93e4f93e4f93e4f93e4f93e4f93e4f93e4f91e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4992649926499264991e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4792e4b92e4b92e4b92e4b92a4a92a4a92a4")); + + const rs = try r.toString(al, 16); + testing.expect(std.mem.eql(u8, rs, "a900000000000000000000000000000000000000000000000000")); +} + test "big.int shift-right single" { var a = try Int.initSet(al, 0xffff0000); try a.shiftRight(a, 16); diff --git a/std/math/big/rational.zig b/std/math/big/rational.zig index ca3b3a5926fb..ab356d980fd2 100644 --- a/std/math/big/rational.zig +++ b/std/math/big/rational.zig @@ -222,6 +222,7 @@ pub const Rational = struct { exp += 1; } if (mantissa >> msize1 != 1) { + // NOTE: This can be hit if the limb size is small (u8/16). @panic("unexpected bits in result"); } @@ -421,8 +422,6 @@ pub const Rational = struct { } }; -var al = debug.global_allocator; - const SignedDoubleLimb = @IntType(true, DoubleLimb.bit_count); fn gcd(rma: *Int, x: Int, y: Int) !void { @@ -441,11 +440,7 @@ fn gcd(rma: *Int, x: Int, y: Int) !void { r.deinit(); }; - if (x.cmp(y) > 0) { - try gcdLehmer(r, x, y); - } else { - try gcdLehmer(r, y, x); - } + try gcdLehmer(r, x, y); } // Storage must live for the lifetime of the returned value @@ -461,9 +456,6 @@ fn FixedIntFromSignedDoubleLimb(A: SignedDoubleLimb, storage: []Limb) Int { 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 { var x = try xa.clone(); x.abs(); @@ -484,18 +476,8 @@ fn gcdLehmer(r: *Int, xa: Int, ya: Int) !void { debug.assert(x.positive and y.positive); 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 xh: SignedDoubleLimb = x.limbs[x.len - 1]; + var yh: SignedDoubleLimb = if (x.len > y.len) 0 else y.limbs[x.len - 1]; var A: SignedDoubleLimb = 1; var B: SignedDoubleLimb = 0; @@ -546,13 +528,7 @@ fn gcdLehmer(r: *Int, xa: Int, ya: Int) !void { try r.add(x, r.*); x.swap(&T); - x.abs(); y.swap(r); - y.abs(); - - if (x.cmp(y) < 0) { - x.swap(&y); - } } } @@ -568,6 +544,10 @@ fn gcdLehmer(r: *Int, xa: Int, ya: Int) !void { r.swap(&x); } +var buffer: [64 * 8192]u8 = undefined; +var fixed = std.heap.FixedBufferAllocator.init(buffer[0..]); +var al = &fixed.allocator; + test "big.rational gcd non-one small" { var a = try Int.initSet(al, 17); var b = try Int.initSet(al, 97); @@ -608,6 +588,16 @@ test "big.rational gcd large multi-limb result" { testing.expect((try r.to(u256)) == 0xf000000ff00000fff0000ffff000fffff00ffffff1); } +test "big.rational gcd one large" { + var a = try Int.initSet(al, 1897056385327307); + var b = try Int.initSet(al, 2251799813685248); + var r = try Int.init(al); + + try gcd(&r, a, b); + + testing.expect((try r.to(u64)) == 1); +} + fn extractLowBits(a: Int, comptime T: type) T { testing.expect(@typeId(T) == builtin.TypeId.Int); @@ -721,12 +711,7 @@ test "big.rational toFloat" { } 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 a = try Rational.init(al); var prng = std.rand.DefaultPrng.init(0x5EED); var i: usize = 0; while (i < 512) : (i += 1) { From 78af62a19a87904193c59f46ac554330ba872564 Mon Sep 17 00:00:00 2001 From: Marc Tiehuis Date: Tue, 9 Apr 2019 17:44:49 +1200 Subject: [PATCH 7/7] Pack big.Int sign and length fields This effectively takes one-bit from the length field and uses it as the sign bit. It reduces the size of an Int from 40 bits to 32 bits on a 64-bit arch. This also reduces std.Rational from 80 bits to 64 bits. --- src-self-hosted/value.zig | 8 +- std/math/big/int.zig | 341 ++++++++++++++++++++------------------ std/math/big/rational.zig | 39 ++--- 3 files changed, 203 insertions(+), 185 deletions(-) diff --git a/src-self-hosted/value.zig b/src-self-hosted/value.zig index d8c0f7b5c87c..f1c7fc0b5056 100644 --- a/src-self-hosted/value.zig +++ b/src-self-hosted/value.zig @@ -538,21 +538,21 @@ pub const Value = struct { switch (self.base.typ.id) { Type.Id.Int => { const type_ref = try self.base.typ.getLlvmType(ofile.arena, ofile.context); - if (self.big_int.len == 0) { + if (self.big_int.len() == 0) { return llvm.ConstNull(type_ref); } - const unsigned_val = if (self.big_int.len == 1) blk: { + const unsigned_val = if (self.big_int.len() == 1) blk: { break :blk llvm.ConstInt(type_ref, self.big_int.limbs[0], @boolToInt(false)); } else if (@sizeOf(std.math.big.Limb) == @sizeOf(u64)) blk: { break :blk llvm.ConstIntOfArbitraryPrecision( type_ref, - @intCast(c_uint, self.big_int.len), + @intCast(c_uint, self.big_int.len()), @ptrCast([*]u64, self.big_int.limbs.ptr), ); } else { @compileError("std.math.Big.Int.Limb size does not match LLVM"); }; - return if (self.big_int.positive) unsigned_val else llvm.ConstNeg(unsigned_val); + return if (self.big_int.isPositive()) unsigned_val else llvm.ConstNeg(unsigned_val); }, Type.Id.ComptimeInt => unreachable, else => unreachable, diff --git a/std/math/big/int.zig b/std/math/big/int.zig index bd17ed16f727..aced892e186b 100644 --- a/std/math/big/int.zig +++ b/std/math/big/int.zig @@ -22,13 +22,18 @@ comptime { } pub const Int = struct { + const sign_bit: usize = 1 << (usize.bit_count - 1); + allocator: ?*Allocator, - positive: bool, // - little-endian ordered // - len >= 1 always // - zero value -> len == 1 with limbs[0] == 0 limbs: []Limb, - len: usize, + // High bit is the sign bit. 1 is negative, 0 positive. + // Remaining bits indicate the number of used limbs. + // + // If Zig gets smarter about packing data, this can be rewritten as a u1 and usize - 1 field. + metadata: usize, const default_capacity = 4; @@ -45,26 +50,45 @@ pub const Int = struct { pub fn initCapacity(allocator: *Allocator, capacity: usize) !Int { return Int{ .allocator = allocator, - .positive = true, + .metadata = 1, .limbs = block: { var limbs = try allocator.alloc(Limb, math.max(default_capacity, capacity)); limbs[0] = 0; break :block limbs; }, - .len = 1, }; } + pub fn len(self: Int) usize { + return self.metadata & ~sign_bit; + } + + pub fn isPositive(self: Int) bool { + return self.metadata & sign_bit == 0; + } + + pub fn setSign(self: *Int, positive: bool) void { + if (positive) { + self.metadata &= ~sign_bit; + } else { + self.metadata |= sign_bit; + } + } + + pub fn setLen(self: *Int, new_len: usize) void { + self.metadata &= sign_bit; + self.metadata |= new_len; + } + // Initialize an Int directly from a fixed set of limb values. This is considered read-only // and cannot be used as a receiver argument to any functions. If this tries to allocate // at any point a panic will occur due to the null allocator. pub fn initFixed(limbs: []const Limb) Int { var self = Int{ .allocator = null, - .positive = true, + .metadata = limbs.len, // Cast away the const, invalid use to pass as a pointer argument. .limbs = @intToPtr([*]Limb, @ptrToInt(limbs.ptr))[0..limbs.len], - .len = limbs.len, }; self.normalize(limbs.len); @@ -96,13 +120,12 @@ pub const Int = struct { other.assertWritable(); return Int{ .allocator = other.allocator, - .positive = other.positive, + .metadata = other.metadata, .limbs = block: { - var limbs = try other.allocator.?.alloc(Limb, other.len); - mem.copy(Limb, limbs[0..], other.limbs[0..other.len]); + var limbs = try other.allocator.?.alloc(Limb, other.len()); + mem.copy(Limb, limbs[0..], other.limbs[0..other.len()]); break :block limbs; }, - .len = other.len, }; } @@ -112,10 +135,9 @@ pub const Int = struct { return; } - self.positive = other.positive; - try self.ensureCapacity(other.len); - mem.copy(Limb, self.limbs[0..], other.limbs[0..other.len]); - self.len = other.len; + try self.ensureCapacity(other.len()); + mem.copy(Limb, self.limbs[0..], other.limbs[0..other.len()]); + self.metadata = other.metadata; } pub fn swap(self: *Int, other: *Int) void { @@ -131,11 +153,11 @@ pub const Int = struct { } pub fn negate(self: *Int) void { - self.positive = !self.positive; + self.metadata ^= sign_bit; } pub fn abs(self: *Int) void { - self.positive = true; + self.metadata &= ~sign_bit; } pub fn isOdd(self: Int) bool { @@ -148,7 +170,7 @@ pub const Int = struct { // Returns the number of bits required to represent the absolute value of self. fn bitCountAbs(self: Int) usize { - return (self.len - 1) * Limb.bit_count + (Limb.bit_count - @clz(self.limbs[self.len - 1])); + return (self.len() - 1) * Limb.bit_count + (Limb.bit_count - @clz(self.limbs[self.len() - 1])); } // Returns the number of bits required to represent the integer in twos-complement form. @@ -164,11 +186,11 @@ pub const Int = struct { // If the entire value has only one bit set (e.g. 0b100000000) then the negation in twos // complement requires one less bit. - if (!self.positive) block: { + if (!self.isPositive()) block: { bits += 1; - if (@popCount(self.limbs[self.len - 1]) == 1) { - for (self.limbs[0 .. self.len - 1]) |limb| { + if (@popCount(self.limbs[self.len() - 1]) == 1) { + for (self.limbs[0 .. self.len() - 1]) |limb| { if (@popCount(limb) != 0) { break :block; } @@ -185,11 +207,11 @@ pub const Int = struct { if (self.eqZero()) { return true; } - if (!is_signed and !self.positive) { + if (!is_signed and !self.isPositive()) { return false; } - const req_bits = self.bitCountTwosComp() + @boolToInt(self.positive and is_signed); + const req_bits = self.bitCountTwosComp() + @boolToInt(self.isPositive() and is_signed); return bit_count >= req_bits; } @@ -201,7 +223,7 @@ pub const Int = struct { // the minus sign. This is used for determining the number of characters needed to print the // value. It is inexact and will exceed the given value by 1-2 digits. pub fn sizeInBase(self: Int, base: usize) usize { - const bit_count = usize(@boolToInt(!self.positive)) + self.bitCountAbs(); + const bit_count = usize(@boolToInt(!self.isPositive())) + self.bitCountAbs(); return (bit_count / math.log2(base)) + 1; } @@ -214,19 +236,19 @@ pub const Int = struct { const UT = if (T.is_signed) @IntType(false, T.bit_count - 1) else T; try self.ensureCapacity(@sizeOf(UT) / @sizeOf(Limb)); - self.positive = value >= 0; - self.len = 0; + self.metadata = 0; + self.setSign(value >= 0); var w_value: UT = if (value < 0) @intCast(UT, -value) else @intCast(UT, value); if (info.bits <= Limb.bit_count) { self.limbs[0] = Limb(w_value); - self.len = 1; + self.metadata += 1; } else { var i: usize = 0; while (w_value != 0) : (i += 1) { self.limbs[i] = @truncate(Limb, w_value); - self.len += 1; + self.metadata += 1; // TODO: shift == 64 at compile-time fails. Fails on u128 limbs. w_value >>= Limb.bit_count / 2; @@ -240,8 +262,8 @@ pub const Int = struct { const req_limbs = @divFloor(math.log2(w_value), Limb.bit_count) + 1; try self.ensureCapacity(req_limbs); - self.positive = value >= 0; - self.len = req_limbs; + self.metadata = req_limbs; + self.setSign(value >= 0); if (w_value <= maxInt(Limb)) { self.limbs[0] = w_value; @@ -282,17 +304,17 @@ pub const Int = struct { if (@sizeOf(UT) <= @sizeOf(Limb)) { r = @intCast(UT, self.limbs[0]); } else { - for (self.limbs[0..self.len]) |_, ri| { - const limb = self.limbs[self.len - ri - 1]; + for (self.limbs[0..self.len()]) |_, ri| { + const limb = self.limbs[self.len() - ri - 1]; r <<= Limb.bit_count; r |= limb; } } if (!T.is_signed) { - return if (self.positive) @intCast(T, r) else error.NegativeIntoUnsigned; + return if (self.isPositive()) @intCast(T, r) else error.NegativeIntoUnsigned; } else { - if (self.positive) { + if (self.isPositive()) { return @intCast(T, r); } else { if (math.cast(T, r)) |ok| { @@ -355,7 +377,7 @@ pub const Int = struct { try self.mul(self.*, ap_base); try self.add(self.*, ap_d); } - self.positive = positive; + self.setSign(positive); } /// TODO make this call format instead of the other way around @@ -377,7 +399,7 @@ pub const Int = struct { if (base & (base - 1) == 0) { const base_shift = math.log2_int(Limb, base); - for (self.limbs[0..self.len]) |limb| { + for (self.limbs[0..self.len()]) |limb| { var shift: usize = 0; while (shift < Limb.bit_count) : (shift += base_shift) { const r = @intCast(u8, (limb >> @intCast(Log2Limb, shift)) & Limb(base - 1)); @@ -404,11 +426,11 @@ pub const Int = struct { } var q = try self.clone(); - q.positive = true; + q.abs(); var r = try Int.init(allocator); var b = try Int.initSet(allocator, limb_base); - while (q.len >= 2) { + while (q.len() >= 2) { try Int.divTrunc(&q, &r, q, b); var r_word = r.limbs[0]; @@ -421,7 +443,7 @@ pub const Int = struct { } { - debug.assert(q.len == 1); + debug.assert(q.len() == 1); var r_word = q.limbs[0]; while (r_word != 0) { @@ -432,7 +454,7 @@ pub const Int = struct { } } - if (!self.positive) { + if (!self.isPositive()) { try digits.append('-'); } @@ -460,14 +482,14 @@ pub const Int = struct { // returns -1, 0, 1 if |a| < |b|, |a| == |b| or |a| > |b| respectively. pub fn cmpAbs(a: Int, b: Int) i8 { - if (a.len < b.len) { + if (a.len() < b.len()) { return -1; } - if (a.len > b.len) { + if (a.len() > b.len()) { return 1; } - var i: usize = a.len - 1; + var i: usize = a.len() - 1; while (i != 0) : (i -= 1) { if (a.limbs[i] != b.limbs[i]) { break; @@ -485,17 +507,17 @@ pub const Int = struct { // returns -1, 0, 1 if a < b, a == b or a > b respectively. pub fn cmp(a: Int, b: Int) i8 { - if (a.positive != b.positive) { - return if (a.positive) i8(1) else -1; + if (a.isPositive() != b.isPositive()) { + return if (a.isPositive()) i8(1) else -1; } else { const r = cmpAbs(a, b); - return if (a.positive) r else -r; + return if (a.isPositive()) r else -r; } } // if a == 0 pub fn eqZero(a: Int) bool { - return a.len == 1 and a.limbs[0] == 0; + return a.len() == 1 and a.limbs[0] == 0; } // if |a| == |b| @@ -525,16 +547,15 @@ pub const Int = struct { } // Handle zero - r.len = if (j != 0) j else 1; + r.setLen(if (j != 0) j else 1); } // Cannot be used as a result argument to any function. fn readOnlyPositive(a: Int) Int { return Int{ .allocator = null, - .positive = true, + .metadata = a.len(), .limbs = a.limbs, - .len = a.len, }; } @@ -549,8 +570,8 @@ pub const Int = struct { return; } - if (a.positive != b.positive) { - if (a.positive) { + if (a.isPositive() != b.isPositive()) { + if (a.isPositive()) { // (a) + (-b) => a - b try r.sub(a, readOnlyPositive(b)); } else { @@ -558,17 +579,17 @@ pub const Int = struct { try r.sub(b, readOnlyPositive(a)); } } else { - if (a.len >= b.len) { - try r.ensureCapacity(a.len + 1); - lladd(r.limbs[0..], a.limbs[0..a.len], b.limbs[0..b.len]); - r.normalize(a.len + 1); + if (a.len() >= b.len()) { + try r.ensureCapacity(a.len() + 1); + lladd(r.limbs[0..], a.limbs[0..a.len()], b.limbs[0..b.len()]); + r.normalize(a.len() + 1); } else { - try r.ensureCapacity(b.len + 1); - lladd(r.limbs[0..], b.limbs[0..b.len], a.limbs[0..a.len]); - r.normalize(b.len + 1); + try r.ensureCapacity(b.len() + 1); + lladd(r.limbs[0..], b.limbs[0..b.len()], a.limbs[0..a.len()]); + r.normalize(b.len() + 1); } - r.positive = a.positive; + r.setSign(a.isPositive()); } } @@ -599,41 +620,41 @@ pub const Int = struct { // r = a - b pub fn sub(r: *Int, a: Int, b: Int) !void { r.assertWritable(); - if (a.positive != b.positive) { - if (a.positive) { + if (a.isPositive() != b.isPositive()) { + if (a.isPositive()) { // (a) - (-b) => a + b try r.add(a, readOnlyPositive(b)); } else { // (-a) - (b) => -(a + b) try r.add(readOnlyPositive(a), b); - r.positive = false; + r.setSign(false); } } else { - if (a.positive) { + if (a.isPositive()) { // (a) - (b) => a - b if (a.cmp(b) >= 0) { - try r.ensureCapacity(a.len + 1); - llsub(r.limbs[0..], a.limbs[0..a.len], b.limbs[0..b.len]); - r.normalize(a.len); - r.positive = true; + try r.ensureCapacity(a.len() + 1); + llsub(r.limbs[0..], a.limbs[0..a.len()], b.limbs[0..b.len()]); + r.normalize(a.len()); + r.setSign(true); } else { - try r.ensureCapacity(b.len + 1); - llsub(r.limbs[0..], b.limbs[0..b.len], a.limbs[0..a.len]); - r.normalize(b.len); - r.positive = false; + try r.ensureCapacity(b.len() + 1); + llsub(r.limbs[0..], b.limbs[0..b.len()], a.limbs[0..a.len()]); + r.normalize(b.len()); + r.setSign(false); } } else { // (-a) - (-b) => -(a - b) if (a.cmp(b) < 0) { - try r.ensureCapacity(a.len + 1); - llsub(r.limbs[0..], a.limbs[0..a.len], b.limbs[0..b.len]); - r.normalize(a.len); - r.positive = false; + try r.ensureCapacity(a.len() + 1); + llsub(r.limbs[0..], a.limbs[0..a.len()], b.limbs[0..b.len()]); + r.normalize(a.len()); + r.setSign(false); } else { - try r.ensureCapacity(b.len + 1); - llsub(r.limbs[0..], b.limbs[0..b.len], a.limbs[0..a.len]); - r.normalize(b.len); - r.positive = true; + try r.ensureCapacity(b.len() + 1); + llsub(r.limbs[0..], b.limbs[0..b.len()], a.limbs[0..a.len()]); + r.normalize(b.len()); + r.setSign(true); } } } @@ -674,7 +695,7 @@ pub const Int = struct { var sr: Int = undefined; if (aliased) { - sr = try Int.initCapacity(rma.allocator.?, a.len + b.len); + sr = try Int.initCapacity(rma.allocator.?, a.len() + b.len()); r = &sr; aliased = true; } @@ -683,16 +704,16 @@ pub const Int = struct { r.deinit(); }; - try r.ensureCapacity(a.len + b.len); + try r.ensureCapacity(a.len() + b.len()); - if (a.len >= b.len) { - llmul(r.limbs, a.limbs[0..a.len], b.limbs[0..b.len]); + if (a.len() >= b.len()) { + llmul(r.limbs, a.limbs[0..a.len()], b.limbs[0..b.len()]); } else { - llmul(r.limbs, b.limbs[0..b.len], a.limbs[0..a.len]); + llmul(r.limbs, b.limbs[0..b.len()], a.limbs[0..a.len()]); } - r.positive = a.positive == b.positive; - r.normalize(a.len + b.len); + r.normalize(a.len() + b.len()); + r.setSign(a.isPositive() == b.isPositive()); } // a + b * c + *carry, sets carry to the overflow bits @@ -742,17 +763,17 @@ pub const Int = struct { try div(q, r, a, b); // Trunc -> Floor. - if (!q.positive) { + if (!q.isPositive()) { const one = Int.initFixed(([]Limb{1})[0..]); try q.sub(q.*, one); try r.add(q.*, one); } - r.positive = b.positive; + r.setSign(b.isPositive()); } pub fn divTrunc(q: *Int, r: *Int, a: Int, b: Int) !void { try div(q, r, a, b); - r.positive = a.positive; + r.setSign(a.isPositive()); } // Truncates by default. @@ -770,10 +791,9 @@ pub const Int = struct { if (a.cmpAbs(b) < 0) { // quo may alias a so handle rem first try rem.copy(a); - rem.positive = a.positive == b.positive; + rem.setSign(a.isPositive() == b.isPositive()); - quo.positive = true; - quo.len = 1; + quo.metadata = 1; quo.limbs[0] = 0; return; } @@ -782,14 +802,14 @@ pub const Int = struct { // algorithms. const a_zero_limb_count = blk: { var i: usize = 0; - while (i < a.len) : (i += 1) { + while (i < a.len()) : (i += 1) { if (a.limbs[i] != 0) break; } break :blk i; }; const b_zero_limb_count = blk: { var i: usize = 0; - while (i < b.len) : (i += 1) { + while (i < b.len()) : (i += 1) { if (b.limbs[i] != 0) break; } break :blk i; @@ -797,39 +817,37 @@ pub const Int = struct { const ab_zero_limb_count = std.math.min(a_zero_limb_count, b_zero_limb_count); - if (b.len - ab_zero_limb_count == 1) { - try quo.ensureCapacity(a.len); + if (b.len() - ab_zero_limb_count == 1) { + try quo.ensureCapacity(a.len()); - lldiv1(quo.limbs[0..], &rem.limbs[0], a.limbs[ab_zero_limb_count..a.len], b.limbs[b.len - 1]); - quo.normalize(a.len - ab_zero_limb_count); - quo.positive = a.positive == b.positive; + lldiv1(quo.limbs[0..], &rem.limbs[0], a.limbs[ab_zero_limb_count..a.len()], b.limbs[b.len() - 1]); + quo.normalize(a.len() - ab_zero_limb_count); + quo.setSign(a.isPositive() == b.isPositive()); - rem.len = 1; - rem.positive = true; + rem.metadata = 1; } else { // x and y are modified during division - var x = try Int.initCapacity(quo.allocator.?, a.len); + var x = try Int.initCapacity(quo.allocator.?, a.len()); defer x.deinit(); try x.copy(a); - var y = try Int.initCapacity(quo.allocator.?, b.len); + var y = try Int.initCapacity(quo.allocator.?, b.len()); defer y.deinit(); try y.copy(b); // x may grow one limb during normalization - try quo.ensureCapacity(a.len + y.len); + try quo.ensureCapacity(a.len() + y.len()); // Shrink x, y such that the trailing zero limbs shared between are removed. if (ab_zero_limb_count != 0) { std.mem.copy(Limb, x.limbs[0..], x.limbs[ab_zero_limb_count..]); std.mem.copy(Limb, y.limbs[0..], y.limbs[ab_zero_limb_count..]); - x.len -= ab_zero_limb_count; - y.len -= ab_zero_limb_count; + x.metadata -= ab_zero_limb_count; + y.metadata -= ab_zero_limb_count; } try divN(quo.allocator.?, quo, rem, &x, &y); - - quo.positive = a.positive == b.positive; + quo.setSign(a.isPositive() == b.isPositive()); } if (ab_zero_limb_count != 0) { @@ -868,28 +886,28 @@ pub const Int = struct { // // x = qy + r where 0 <= r < y fn divN(allocator: *Allocator, q: *Int, r: *Int, x: *Int, y: *Int) !void { - debug.assert(y.len >= 2); - debug.assert(x.len >= y.len); - debug.assert(q.limbs.len >= x.len + y.len - 1); + debug.assert(y.len() >= 2); + debug.assert(x.len() >= y.len()); + debug.assert(q.limbs.len >= x.len() + y.len() - 1); debug.assert(default_capacity >= 3); // see 3.2 var tmp = try Int.init(allocator); defer tmp.deinit(); // Normalize so y > Limb.bit_count / 2 (i.e. leading bit is set) and even - var norm_shift = @clz(y.limbs[y.len - 1]); + var norm_shift = @clz(y.limbs[y.len() - 1]); if (norm_shift == 0 and y.isOdd()) { norm_shift = Limb.bit_count; } try x.shiftLeft(x.*, norm_shift); try y.shiftLeft(y.*, norm_shift); - const n = x.len - 1; - const t = y.len - 1; + const n = x.len() - 1; + const t = y.len() - 1; // 1. - q.len = n - t + 1; - mem.set(Limb, q.limbs[0..q.len], 0); + q.metadata = n - t + 1; + mem.set(Limb, q.limbs[0..q.len()], 0); // 2. try tmp.shiftLeft(y.*, Limb.bit_count * (n - t)); @@ -937,7 +955,7 @@ pub const Int = struct { try tmp.shiftLeft(tmp, Limb.bit_count * (i - t - 1)); try x.sub(x.*, tmp); - if (!x.positive) { + if (!x.isPositive()) { try tmp.shiftLeft(y.*, Limb.bit_count * (i - t - 1)); try x.add(x.*, tmp); q.limbs[i - t - 1] -= 1; @@ -945,20 +963,20 @@ pub const Int = struct { } // Denormalize - q.normalize(q.len); + q.normalize(q.len()); try r.shiftRight(x.*, norm_shift); - r.normalize(r.len); + r.normalize(r.len()); } // r = a << shift, in other words, r = a * 2^shift pub fn shiftLeft(r: *Int, a: Int, shift: usize) !void { r.assertWritable(); - try r.ensureCapacity(a.len + (shift / Limb.bit_count) + 1); - llshl(r.limbs[0..], a.limbs[0..a.len], shift); - r.normalize(a.len + (shift / Limb.bit_count) + 1); - r.positive = a.positive; + try r.ensureCapacity(a.len() + (shift / Limb.bit_count) + 1); + llshl(r.limbs[0..], a.limbs[0..a.len()], shift); + r.normalize(a.len() + (shift / Limb.bit_count) + 1); + r.setSign(a.isPositive()); } fn llshl(r: []Limb, a: []const Limb, shift: usize) void { @@ -988,17 +1006,16 @@ pub const Int = struct { pub fn shiftRight(r: *Int, a: Int, shift: usize) !void { r.assertWritable(); - if (a.len <= shift / Limb.bit_count) { - r.len = 1; + if (a.len() <= shift / Limb.bit_count) { + r.metadata = 1; r.limbs[0] = 0; - r.positive = true; return; } - try r.ensureCapacity(a.len - (shift / Limb.bit_count)); - const r_len = llshr(r.limbs[0..], a.limbs[0..a.len], shift); - r.len = a.len - (shift / Limb.bit_count); - r.positive = a.positive; + try r.ensureCapacity(a.len() - (shift / Limb.bit_count)); + const r_len = llshr(r.limbs[0..], a.limbs[0..a.len()], shift); + r.metadata = a.len() - (shift / Limb.bit_count); + r.setSign(a.isPositive()); } fn llshr(r: []Limb, a: []const Limb, shift: usize) void { @@ -1025,14 +1042,14 @@ pub const Int = struct { pub fn bitOr(r: *Int, a: Int, b: Int) !void { r.assertWritable(); - if (a.len > b.len) { - try r.ensureCapacity(a.len); - llor(r.limbs[0..], a.limbs[0..a.len], b.limbs[0..b.len]); - r.len = a.len; + if (a.len() > b.len()) { + try r.ensureCapacity(a.len()); + llor(r.limbs[0..], a.limbs[0..a.len()], b.limbs[0..b.len()]); + r.setLen(a.len()); } else { - try r.ensureCapacity(b.len); - llor(r.limbs[0..], b.limbs[0..b.len], a.limbs[0..a.len]); - r.len = b.len; + try r.ensureCapacity(b.len()); + llor(r.limbs[0..], b.limbs[0..b.len()], a.limbs[0..a.len()]); + r.setLen(b.len()); } } @@ -1054,14 +1071,14 @@ pub const Int = struct { pub fn bitAnd(r: *Int, a: Int, b: Int) !void { r.assertWritable(); - if (a.len > b.len) { - try r.ensureCapacity(b.len); - lland(r.limbs[0..], a.limbs[0..a.len], b.limbs[0..b.len]); - r.normalize(b.len); + if (a.len() > b.len()) { + try r.ensureCapacity(b.len()); + lland(r.limbs[0..], a.limbs[0..a.len()], b.limbs[0..b.len()]); + r.normalize(b.len()); } else { - try r.ensureCapacity(a.len); - lland(r.limbs[0..], b.limbs[0..b.len], a.limbs[0..a.len]); - r.normalize(a.len); + try r.ensureCapacity(a.len()); + lland(r.limbs[0..], b.limbs[0..b.len()], a.limbs[0..a.len()]); + r.normalize(a.len()); } } @@ -1080,14 +1097,14 @@ pub const Int = struct { pub fn bitXor(r: *Int, a: Int, b: Int) !void { r.assertWritable(); - if (a.len > b.len) { - try r.ensureCapacity(a.len); - llxor(r.limbs[0..], a.limbs[0..a.len], b.limbs[0..b.len]); - r.normalize(a.len); + if (a.len() > b.len()) { + try r.ensureCapacity(a.len()); + llxor(r.limbs[0..], a.limbs[0..a.len()], b.limbs[0..b.len()]); + r.normalize(a.len()); } else { - try r.ensureCapacity(b.len); - llxor(r.limbs[0..], b.limbs[0..b.len], a.limbs[0..a.len]); - r.normalize(b.len); + try r.ensureCapacity(b.len()); + llxor(r.limbs[0..], b.limbs[0..b.len()], a.limbs[0..a.len()]); + r.normalize(b.len()); } } @@ -1134,14 +1151,14 @@ test "big.int comptime_int set negative" { var a = try Int.initSet(al, -10); testing.expect(a.limbs[0] == 10); - testing.expect(a.positive == false); + testing.expect(a.isPositive() == false); } test "big.int int set unaligned small" { var a = try Int.initSet(al, u7(45)); testing.expect(a.limbs[0] == 45); - testing.expect(a.positive == true); + testing.expect(a.isPositive() == true); } test "big.int comptime_int to" { @@ -1171,22 +1188,22 @@ test "big.int normalize" { a.limbs[2] = 3; a.limbs[3] = 0; a.normalize(4); - testing.expect(a.len == 3); + testing.expect(a.len() == 3); a.limbs[0] = 1; a.limbs[1] = 2; a.limbs[2] = 3; a.normalize(3); - testing.expect(a.len == 3); + testing.expect(a.len() == 3); a.limbs[0] = 0; a.limbs[1] = 0; a.normalize(2); - testing.expect(a.len == 1); + testing.expect(a.len() == 1); a.limbs[0] = 0; a.normalize(1); - testing.expect(a.len == 1); + testing.expect(a.len() == 1); } test "big.int normalize multi" { @@ -1198,24 +1215,24 @@ test "big.int normalize multi" { a.limbs[2] = 0; a.limbs[3] = 0; a.normalize(4); - testing.expect(a.len == 2); + testing.expect(a.len() == 2); a.limbs[0] = 1; a.limbs[1] = 2; a.limbs[2] = 3; a.normalize(3); - testing.expect(a.len == 3); + testing.expect(a.len() == 3); a.limbs[0] = 0; a.limbs[1] = 0; a.limbs[2] = 0; a.limbs[3] = 0; a.normalize(4); - testing.expect(a.len == 1); + testing.expect(a.len() == 1); a.limbs[0] = 0; a.normalize(1); - testing.expect(a.len == 1); + testing.expect(a.len() == 1); } test "big.int parity" { @@ -1250,7 +1267,7 @@ test "big.int bitcount + sizeInBase" { try a.shiftLeft(a, 5000); testing.expect(a.bitCountAbs() == 5032); testing.expect(a.sizeInBase(2) >= 5032); - a.positive = false; + a.setSign(false); testing.expect(a.bitCountAbs() == 5032); testing.expect(a.sizeInBase(2) >= 5033); diff --git a/std/math/big/rational.zig b/std/math/big/rational.zig index ab356d980fd2..34cf3bf7644b 100644 --- a/std/math/big/rational.zig +++ b/std/math/big/rational.zig @@ -15,7 +15,7 @@ const DoubleLimb = bn.DoubleLimb; const Int = bn.Int; pub const Rational = struct { - // sign of Rational is a.positive, b.positive is ignored + // Sign of Rational is sign of p. Sign of q is ignored p: Int, q: Int, @@ -152,7 +152,7 @@ pub const Rational = struct { } try self.p.set(mantissa); - self.p.positive = f >= 0; + self.p.setSign(f >= 0); try self.q.set(1); if (shift >= 0) { @@ -211,7 +211,7 @@ pub const Rational = struct { try Int.divTrunc(&q, &r, a2, b2); var mantissa = extractLowBits(q, BitReprType); - var have_rem = r.len > 0; + var have_rem = r.len() > 0; // 3. q didn't fit in msize2 bits, redo division b2 << 1 if (mantissa >> msize2 == 1) { @@ -256,15 +256,16 @@ pub const Rational = struct { exact = false; } - return if (self.p.positive) f else -f; + return if (self.p.isPositive()) 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; + self.p.setSign(@boolToInt(self.p.isPositive()) ^ @boolToInt(self.q.isPositive()) == 0); + self.q.setSign(true); + try self.reduce(); if (self.q.eqZero()) { @@ -281,8 +282,9 @@ pub const Rational = struct { 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; + self.p.setSign(@boolToInt(self.p.isPositive()) ^ @boolToInt(self.q.isPositive()) == 0); + self.q.setSign(true); + try self.reduce(); } @@ -403,11 +405,10 @@ pub const Rational = struct { var a = try Int.init(r.p.allocator.?); defer a.deinit(); - const sign = r.p.positive; - + const sign = r.p.isPositive(); r.p.abs(); try gcd(&a, r.p, r.q); - r.p.positive = sign; + r.p.setSign(sign); const one = Int.initFixed(([]Limb{1})[0..]); if (a.cmp(one) != 0) { @@ -431,7 +432,7 @@ fn gcd(rma: *Int, x: Int, y: Int) !void { var sr: Int = undefined; if (aliased) { - sr = try Int.initCapacity(rma.allocator.?, math.max(x.len, y.len)); + sr = try Int.initCapacity(rma.allocator.?, math.max(x.len(), y.len())); r = &sr; aliased = true; } @@ -452,7 +453,7 @@ fn FixedIntFromSignedDoubleLimb(A: SignedDoubleLimb, storage: []Limb) Int { 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; + Ap.setSign(A_is_positive); return Ap; } @@ -472,12 +473,12 @@ fn gcdLehmer(r: *Int, xa: Int, ya: Int) !void { var T = try Int.init(r.allocator.?); defer T.deinit(); - while (y.len > 1) { - debug.assert(x.positive and y.positive); - debug.assert(x.len >= y.len); + while (y.len() > 1) { + debug.assert(x.isPositive() and y.isPositive()); + debug.assert(x.len() >= y.len()); - var xh: SignedDoubleLimb = x.limbs[x.len - 1]; - var yh: SignedDoubleLimb = if (x.len > y.len) 0 else y.limbs[x.len - 1]; + var xh: SignedDoubleLimb = x.limbs[x.len() - 1]; + var yh: SignedDoubleLimb = if (x.len() > y.len()) 0 else y.limbs[x.len() - 1]; var A: SignedDoubleLimb = 1; var B: SignedDoubleLimb = 0; @@ -506,7 +507,7 @@ fn gcdLehmer(r: *Int, xa: Int, ya: Int) !void { if (B == 0) { // T = x % y, r is unused try Int.divTrunc(r, &T, x, y); - debug.assert(T.positive); + debug.assert(T.isPositive()); x.swap(&y); y.swap(&T);