authorgravatar for marc@tiehu.isMarc Tiehuis <marc@tiehu.is> 2019-04-11 22:54:19+12:00
committergravatar for noreply@github.comGitHub <noreply@github.com> 2019-04-11 22:54:19+12:00
logb59c65e9864e20ec88ebb0c3e9e9d38d4aa2e1fc
tree462a160ce76c8d7c92766f753adef9e62841415e
parentdff201540f91bb41b44d7f73a11f4082f928cff7
parent78af62a19a87904193c59f46ac554330ba872564
signaturebadge-question-mark Signed by PGP key 4AEE18F83AFDEB23

Merge pull request #2102 from ziglang/big.int-additions

Add big.Rational type to std

5 files changed, 1303 insertions(+), 249 deletions(-)

CMakeLists.txt+1
...@@ -521,6 +521,7 @@ set(ZIG_STD_FILES...@@ -521,6 +521,7 @@ set(ZIG_STD_FILES
521 "math/atanh.zig"521 "math/atanh.zig"
522 "math/big.zig"522 "math/big.zig"
523 "math/big/int.zig"523 "math/big/int.zig"
524 "math/big/rational.zig"
524 "math/cbrt.zig"525 "math/cbrt.zig"
525 "math/ceil.zig"526 "math/ceil.zig"
526 "math/complex.zig"527 "math/complex.zig"
src-self-hosted/value.zig+4-4
...@@ -538,21 +538,21 @@ pub const Value = struct {...@@ -538,21 +538,21 @@ pub const Value = struct {
538 switch (self.base.typ.id) {538 switch (self.base.typ.id) {
539 Type.Id.Int => {539 Type.Id.Int => {
540 const type_ref = try self.base.typ.getLlvmType(ofile.arena, ofile.context);540 const type_ref = try self.base.typ.getLlvmType(ofile.arena, ofile.context);
541 if (self.big_int.len == 0) {541 if (self.big_int.len() == 0) {
542 return llvm.ConstNull(type_ref);542 return llvm.ConstNull(type_ref);
543 }543 }
544 const unsigned_val = if (self.big_int.len == 1) blk: {544 const unsigned_val = if (self.big_int.len() == 1) blk: {
545 break :blk llvm.ConstInt(type_ref, self.big_int.limbs[0], @boolToInt(false));545 break :blk llvm.ConstInt(type_ref, self.big_int.limbs[0], @boolToInt(false));
546 } else if (@sizeOf(std.math.big.Limb) == @sizeOf(u64)) blk: {546 } else if (@sizeOf(std.math.big.Limb) == @sizeOf(u64)) blk: {
547 break :blk llvm.ConstIntOfArbitraryPrecision(547 break :blk llvm.ConstIntOfArbitraryPrecision(
548 type_ref,548 type_ref,
549 @intCast(c_uint, self.big_int.len),549 @intCast(c_uint, self.big_int.len()),
550 @ptrCast([*]u64, self.big_int.limbs.ptr),550 @ptrCast([*]u64, self.big_int.limbs.ptr),
551 );551 );
552 } else {552 } else {
553 @compileError("std.math.Big.Int.Limb size does not match LLVM");553 @compileError("std.math.Big.Int.Limb size does not match LLVM");
554 };554 };
555 return if (self.big_int.positive) unsigned_val else llvm.ConstNeg(unsigned_val);555 return if (self.big_int.isPositive()) unsigned_val else llvm.ConstNeg(unsigned_val);
556 },556 },
557 Type.Id.ComptimeInt => unreachable,557 Type.Id.ComptimeInt => unreachable,
558 else => unreachable,558 else => unreachable,
std/math/big.zig+2
...@@ -1,5 +1,7 @@...@@ -1,5 +1,7 @@
1pub use @import("big/int.zig");1pub use @import("big/int.zig");
2pub use @import("big/rational.zig");
23
3test "math.big" {4test "math.big" {
4 _ = @import("big/int.zig");5 _ = @import("big/int.zig");
6 _ = @import("big/rational.zig");
5}7}
std/math/big/int.zig+400-245
...@@ -22,13 +22,18 @@ comptime {...@@ -22,13 +22,18 @@ comptime {
22}22}
2323
24pub const Int = struct {24pub const Int = struct {
25 allocator: *Allocator,25 const sign_bit: usize = 1 << (usize.bit_count - 1);
26 positive: bool,26
27 allocator: ?*Allocator,
27 // - little-endian ordered28 // - little-endian ordered
28 // - len >= 1 always29 // - len >= 1 always
29 // - zero value -> len == 1 with limbs[0] == 030 // - zero value -> len == 1 with limbs[0] == 0
30 limbs: []Limb,31 limbs: []Limb,
31 len: usize,32 // High bit is the sign bit. 1 is negative, 0 positive.
33 // Remaining bits indicate the number of used limbs.
34 //
35 // If Zig gets smarter about packing data, this can be rewritten as a u1 and usize - 1 field.
36 metadata: usize,
3237
33 const default_capacity = 4;38 const default_capacity = 4;
3439
...@@ -45,54 +50,98 @@ pub const Int = struct {...@@ -45,54 +50,98 @@ pub const Int = struct {
45 pub fn initCapacity(allocator: *Allocator, capacity: usize) !Int {50 pub fn initCapacity(allocator: *Allocator, capacity: usize) !Int {
46 return Int{51 return Int{
47 .allocator = allocator,52 .allocator = allocator,
48 .positive = true,53 .metadata = 1,
49 .limbs = block: {54 .limbs = block: {
50 var limbs = try allocator.alloc(Limb, math.max(default_capacity, capacity));55 var limbs = try allocator.alloc(Limb, math.max(default_capacity, capacity));
51 limbs[0] = 0;56 limbs[0] = 0;
52 break :block limbs;57 break :block limbs;
53 },58 },
54 .len = 1,
55 };59 };
56 }60 }
5761
62 pub fn len(self: Int) usize {
63 return self.metadata & ~sign_bit;
64 }
65
66 pub fn isPositive(self: Int) bool {
67 return self.metadata & sign_bit == 0;
68 }
69
70 pub fn setSign(self: *Int, positive: bool) void {
71 if (positive) {
72 self.metadata &= ~sign_bit;
73 } else {
74 self.metadata |= sign_bit;
75 }
76 }
77
78 pub fn setLen(self: *Int, new_len: usize) void {
79 self.metadata &= sign_bit;
80 self.metadata |= new_len;
81 }
82
83 // Initialize an Int directly from a fixed set of limb values. This is considered read-only
84 // and cannot be used as a receiver argument to any functions. If this tries to allocate
85 // at any point a panic will occur due to the null allocator.
86 pub fn initFixed(limbs: []const Limb) Int {
87 var self = Int{
88 .allocator = null,
89 .metadata = limbs.len,
90 // Cast away the const, invalid use to pass as a pointer argument.
91 .limbs = @intToPtr([*]Limb, @ptrToInt(limbs.ptr))[0..limbs.len],
92 };
93
94 self.normalize(limbs.len);
95 return self;
96 }
97
58 pub fn ensureCapacity(self: *Int, capacity: usize) !void {98 pub fn ensureCapacity(self: *Int, capacity: usize) !void {
99 self.assertWritable();
59 if (capacity <= self.limbs.len) {100 if (capacity <= self.limbs.len) {
60 return;101 return;
61 }102 }
62103
63 self.limbs = try self.allocator.realloc(self.limbs, capacity);104 self.limbs = try self.allocator.?.realloc(self.limbs, capacity);
105 }
106
107 fn assertWritable(self: Int) void {
108 if (self.allocator == null) {
109 @panic("provided Int value is read-only but must be writable");
110 }
64 }111 }
65112
66 pub fn deinit(self: *Int) void {113 pub fn deinit(self: *Int) void {
67 self.allocator.free(self.limbs);114 self.assertWritable();
115 self.allocator.?.free(self.limbs);
68 self.* = undefined;116 self.* = undefined;
69 }117 }
70118
71 pub fn clone(other: Int) !Int {119 pub fn clone(other: Int) !Int {
120 other.assertWritable();
72 return Int{121 return Int{
73 .allocator = other.allocator,122 .allocator = other.allocator,
74 .positive = other.positive,123 .metadata = other.metadata,
75 .limbs = block: {124 .limbs = block: {
76 var limbs = try other.allocator.alloc(Limb, other.len);125 var limbs = try other.allocator.?.alloc(Limb, other.len());
77 mem.copy(Limb, limbs[0..], other.limbs[0..other.len]);126 mem.copy(Limb, limbs[0..], other.limbs[0..other.len()]);
78 break :block limbs;127 break :block limbs;
79 },128 },
80 .len = other.len,
81 };129 };
82 }130 }
83131
84 pub fn copy(self: *Int, other: Int) !void {132 pub fn copy(self: *Int, other: Int) !void {
85 if (self == &other) {133 self.assertWritable();
134 if (self.limbs.ptr == other.limbs.ptr) {
86 return;135 return;
87 }136 }
88137
89 self.positive = other.positive;138 try self.ensureCapacity(other.len());
90 try self.ensureCapacity(other.len);139 mem.copy(Limb, self.limbs[0..], other.limbs[0..other.len()]);
91 mem.copy(Limb, self.limbs[0..], other.limbs[0..other.len]);140 self.metadata = other.metadata;
92 self.len = other.len;
93 }141 }
94142
95 pub fn swap(self: *Int, other: *Int) void {143 pub fn swap(self: *Int, other: *Int) void {
144 self.assertWritable();
96 mem.swap(Int, self, other);145 mem.swap(Int, self, other);
97 }146 }
98147
...@@ -103,25 +152,25 @@ pub const Int = struct {...@@ -103,25 +152,25 @@ pub const Int = struct {
103 debug.warn("\n");152 debug.warn("\n");
104 }153 }
105154
106 pub fn negate(r: *Int) void {155 pub fn negate(self: *Int) void {
107 r.positive = !r.positive;156 self.metadata ^= sign_bit;
108 }157 }
109158
110 pub fn abs(r: *Int) void {159 pub fn abs(self: *Int) void {
111 r.positive = true;160 self.metadata &= ~sign_bit;
112 }161 }
113162
114 pub fn isOdd(r: Int) bool {163 pub fn isOdd(self: Int) bool {
115 return r.limbs[0] & 1 != 0;164 return self.limbs[0] & 1 != 0;
116 }165 }
117166
118 pub fn isEven(r: Int) bool {167 pub fn isEven(self: Int) bool {
119 return !r.isOdd();168 return !self.isOdd();
120 }169 }
121170
122 // Returns the number of bits required to represent the absolute value of self.171 // Returns the number of bits required to represent the absolute value of self.
123 fn bitCountAbs(self: Int) usize {172 fn bitCountAbs(self: Int) usize {
124 return (self.len - 1) * Limb.bit_count + (Limb.bit_count - @clz(self.limbs[self.len - 1]));173 return (self.len() - 1) * Limb.bit_count + (Limb.bit_count - @clz(self.limbs[self.len() - 1]));
125 }174 }
126175
127 // Returns the number of bits required to represent the integer in twos-complement form.176 // Returns the number of bits required to represent the integer in twos-complement form.
...@@ -137,11 +186,11 @@ pub const Int = struct {...@@ -137,11 +186,11 @@ pub const Int = struct {
137186
138 // If the entire value has only one bit set (e.g. 0b100000000) then the negation in twos187 // If the entire value has only one bit set (e.g. 0b100000000) then the negation in twos
139 // complement requires one less bit.188 // complement requires one less bit.
140 if (!self.positive) block: {189 if (!self.isPositive()) block: {
141 bits += 1;190 bits += 1;
142191
143 if (@popCount(self.limbs[self.len - 1]) == 1) {192 if (@popCount(self.limbs[self.len() - 1]) == 1) {
144 for (self.limbs[0 .. self.len - 1]) |limb| {193 for (self.limbs[0 .. self.len() - 1]) |limb| {
145 if (@popCount(limb) != 0) {194 if (@popCount(limb) != 0) {
146 break :block;195 break :block;
147 }196 }
...@@ -158,11 +207,11 @@ pub const Int = struct {...@@ -158,11 +207,11 @@ pub const Int = struct {
158 if (self.eqZero()) {207 if (self.eqZero()) {
159 return true;208 return true;
160 }209 }
161 if (!is_signed and !self.positive) {210 if (!is_signed and !self.isPositive()) {
162 return false;211 return false;
163 }212 }
164213
165 const req_bits = self.bitCountTwosComp() + @boolToInt(self.positive and is_signed);214 const req_bits = self.bitCountTwosComp() + @boolToInt(self.isPositive() and is_signed);
166 return bit_count >= req_bits;215 return bit_count >= req_bits;
167 }216 }
168217
...@@ -174,11 +223,12 @@ pub const Int = struct {...@@ -174,11 +223,12 @@ pub const Int = struct {
174 // the minus sign. This is used for determining the number of characters needed to print the223 // the minus sign. This is used for determining the number of characters needed to print the
175 // value. It is inexact and will exceed the given value by 1-2 digits.224 // value. It is inexact and will exceed the given value by 1-2 digits.
176 pub fn sizeInBase(self: Int, base: usize) usize {225 pub fn sizeInBase(self: Int, base: usize) usize {
177 const bit_count = usize(@boolToInt(!self.positive)) + self.bitCountAbs();226 const bit_count = usize(@boolToInt(!self.isPositive())) + self.bitCountAbs();
178 return (bit_count / math.log2(base)) + 1;227 return (bit_count / math.log2(base)) + 1;
179 }228 }
180229
181 pub fn set(self: *Int, value: var) Allocator.Error!void {230 pub fn set(self: *Int, value: var) Allocator.Error!void {
231 self.assertWritable();
182 const T = @typeOf(value);232 const T = @typeOf(value);
183233
184 switch (@typeInfo(T)) {234 switch (@typeInfo(T)) {
...@@ -186,19 +236,19 @@ pub const Int = struct {...@@ -186,19 +236,19 @@ pub const Int = struct {
186 const UT = if (T.is_signed) @IntType(false, T.bit_count - 1) else T;236 const UT = if (T.is_signed) @IntType(false, T.bit_count - 1) else T;
187237
188 try self.ensureCapacity(@sizeOf(UT) / @sizeOf(Limb));238 try self.ensureCapacity(@sizeOf(UT) / @sizeOf(Limb));
189 self.positive = value >= 0;239 self.metadata = 0;
190 self.len = 0;240 self.setSign(value >= 0);
191241
192 var w_value: UT = if (value < 0) @intCast(UT, -value) else @intCast(UT, value);242 var w_value: UT = if (value < 0) @intCast(UT, -value) else @intCast(UT, value);
193243
194 if (info.bits <= Limb.bit_count) {244 if (info.bits <= Limb.bit_count) {
195 self.limbs[0] = Limb(w_value);245 self.limbs[0] = Limb(w_value);
196 self.len = 1;246 self.metadata += 1;
197 } else {247 } else {
198 var i: usize = 0;248 var i: usize = 0;
199 while (w_value != 0) : (i += 1) {249 while (w_value != 0) : (i += 1) {
200 self.limbs[i] = @truncate(Limb, w_value);250 self.limbs[i] = @truncate(Limb, w_value);
201 self.len += 1;251 self.metadata += 1;
202252
203 // TODO: shift == 64 at compile-time fails. Fails on u128 limbs.253 // TODO: shift == 64 at compile-time fails. Fails on u128 limbs.
204 w_value >>= Limb.bit_count / 2;254 w_value >>= Limb.bit_count / 2;
...@@ -212,8 +262,8 @@ pub const Int = struct {...@@ -212,8 +262,8 @@ pub const Int = struct {
212 const req_limbs = @divFloor(math.log2(w_value), Limb.bit_count) + 1;262 const req_limbs = @divFloor(math.log2(w_value), Limb.bit_count) + 1;
213 try self.ensureCapacity(req_limbs);263 try self.ensureCapacity(req_limbs);
214264
215 self.positive = value >= 0;265 self.metadata = req_limbs;
216 self.len = req_limbs;266 self.setSign(value >= 0);
217267
218 if (w_value <= maxInt(Limb)) {268 if (w_value <= maxInt(Limb)) {
219 self.limbs[0] = w_value;269 self.limbs[0] = w_value;
...@@ -254,17 +304,17 @@ pub const Int = struct {...@@ -254,17 +304,17 @@ pub const Int = struct {
254 if (@sizeOf(UT) <= @sizeOf(Limb)) {304 if (@sizeOf(UT) <= @sizeOf(Limb)) {
255 r = @intCast(UT, self.limbs[0]);305 r = @intCast(UT, self.limbs[0]);
256 } else {306 } else {
257 for (self.limbs[0..self.len]) |_, ri| {307 for (self.limbs[0..self.len()]) |_, ri| {
258 const limb = self.limbs[self.len - ri - 1];308 const limb = self.limbs[self.len() - ri - 1];
259 r <<= Limb.bit_count;309 r <<= Limb.bit_count;
260 r |= limb;310 r |= limb;
261 }311 }
262 }312 }
263313
264 if (!T.is_signed) {314 if (!T.is_signed) {
265 return if (self.positive) @intCast(T, r) else error.NegativeIntoUnsigned;315 return if (self.isPositive()) @intCast(T, r) else error.NegativeIntoUnsigned;
266 } else {316 } else {
267 if (self.positive) {317 if (self.isPositive()) {
268 return @intCast(T, r);318 return @intCast(T, r);
269 } else {319 } else {
270 if (math.cast(T, r)) |ok| {320 if (math.cast(T, r)) |ok| {
...@@ -304,6 +354,7 @@ pub const Int = struct {...@@ -304,6 +354,7 @@ pub const Int = struct {
304 }354 }
305355
306 pub fn setString(self: *Int, base: u8, value: []const u8) !void {356 pub fn setString(self: *Int, base: u8, value: []const u8) !void {
357 self.assertWritable();
307 if (base < 2 or base > 16) {358 if (base < 2 or base > 16) {
308 return error.InvalidBase;359 return error.InvalidBase;
309 }360 }
...@@ -315,25 +366,18 @@ pub const Int = struct {...@@ -315,25 +366,18 @@ pub const Int = struct {
315 i += 1;366 i += 1;
316 }367 }
317368
318 // TODO values less than limb size should guarantee non allocating369 const ap_base = Int.initFixed(([]Limb{base})[0..]);
319 var base_buffer: [512]u8 = undefined;
320 const base_al = &std.heap.FixedBufferAllocator.init(base_buffer[0..]).allocator;
321 const base_ap = try Int.initSet(base_al, base);
322
323 var d_buffer: [512]u8 = undefined;
324 var d_fba = std.heap.FixedBufferAllocator.init(d_buffer[0..]);
325 const d_al = &d_fba.allocator;
326
327 try self.set(0);370 try self.set(0);
371
328 for (value[i..]) |ch| {372 for (value[i..]) |ch| {
329 const d = try charToDigit(ch, base);373 const d = try charToDigit(ch, base);
330 d_fba.end_index = 0;
331 const d_ap = try Int.initSet(d_al, d);
332374
333 try self.mul(self.*, base_ap);375 const ap_d = Int.initFixed(([]Limb{d})[0..]);
334 try self.add(self.*, d_ap);376
377 try self.mul(self.*, ap_base);
378 try self.add(self.*, ap_d);
335 }379 }
336 self.positive = positive;380 self.setSign(positive);
337 }381 }
338382
339 /// TODO make this call format instead of the other way around383 /// TODO make this call format instead of the other way around
...@@ -355,7 +399,7 @@ pub const Int = struct {...@@ -355,7 +399,7 @@ pub const Int = struct {
355 if (base & (base - 1) == 0) {399 if (base & (base - 1) == 0) {
356 const base_shift = math.log2_int(Limb, base);400 const base_shift = math.log2_int(Limb, base);
357401
358 for (self.limbs[0..self.len]) |limb| {402 for (self.limbs[0..self.len()]) |limb| {
359 var shift: usize = 0;403 var shift: usize = 0;
360 while (shift < Limb.bit_count) : (shift += base_shift) {404 while (shift < Limb.bit_count) : (shift += base_shift) {
361 const r = @intCast(u8, (limb >> @intCast(Log2Limb, shift)) & Limb(base - 1));405 const r = @intCast(u8, (limb >> @intCast(Log2Limb, shift)) & Limb(base - 1));
...@@ -382,11 +426,11 @@ pub const Int = struct {...@@ -382,11 +426,11 @@ pub const Int = struct {
382 }426 }
383427
384 var q = try self.clone();428 var q = try self.clone();
385 q.positive = true;429 q.abs();
386 var r = try Int.init(allocator);430 var r = try Int.init(allocator);
387 var b = try Int.initSet(allocator, limb_base);431 var b = try Int.initSet(allocator, limb_base);
388432
389 while (q.len >= 2) {433 while (q.len() >= 2) {
390 try Int.divTrunc(&q, &r, q, b);434 try Int.divTrunc(&q, &r, q, b);
391435
392 var r_word = r.limbs[0];436 var r_word = r.limbs[0];
...@@ -399,7 +443,7 @@ pub const Int = struct {...@@ -399,7 +443,7 @@ pub const Int = struct {
399 }443 }
400444
401 {445 {
402 debug.assert(q.len == 1);446 debug.assert(q.len() == 1);
403447
404 var r_word = q.limbs[0];448 var r_word = q.limbs[0];
405 while (r_word != 0) {449 while (r_word != 0) {
...@@ -410,7 +454,7 @@ pub const Int = struct {...@@ -410,7 +454,7 @@ pub const Int = struct {
410 }454 }
411 }455 }
412456
413 if (!self.positive) {457 if (!self.isPositive()) {
414 try digits.append('-');458 try digits.append('-');
415 }459 }
416460
...@@ -428,22 +472,24 @@ pub const Int = struct {...@@ -428,22 +472,24 @@ pub const Int = struct {
428 comptime FmtError: type,472 comptime FmtError: type,
429 output: fn (@typeOf(context), []const u8) FmtError!void,473 output: fn (@typeOf(context), []const u8) FmtError!void,
430 ) FmtError!void {474 ) FmtError!void {
475 self.assertWritable();
431 // TODO look at fmt and support other bases476 // TODO look at fmt and support other bases
432 const str = self.toString(self.allocator, 10) catch @panic("TODO make this non allocating");477 // TODO support read-only fixed integers
433 defer self.allocator.free(str);478 const str = self.toString(self.allocator.?, 10) catch @panic("TODO make this non allocating");
479 defer self.allocator.?.free(str);
434 return output(context, str);480 return output(context, str);
435 }481 }
436482
437 // returns -1, 0, 1 if |a| < |b|, |a| == |b| or |a| > |b| respectively.483 // returns -1, 0, 1 if |a| < |b|, |a| == |b| or |a| > |b| respectively.
438 pub fn cmpAbs(a: Int, b: Int) i8 {484 pub fn cmpAbs(a: Int, b: Int) i8 {
439 if (a.len < b.len) {485 if (a.len() < b.len()) {
440 return -1;486 return -1;
441 }487 }
442 if (a.len > b.len) {488 if (a.len() > b.len()) {
443 return 1;489 return 1;
444 }490 }
445491
446 var i: usize = a.len - 1;492 var i: usize = a.len() - 1;
447 while (i != 0) : (i -= 1) {493 while (i != 0) : (i -= 1) {
448 if (a.limbs[i] != b.limbs[i]) {494 if (a.limbs[i] != b.limbs[i]) {
449 break;495 break;
...@@ -461,17 +507,17 @@ pub const Int = struct {...@@ -461,17 +507,17 @@ pub const Int = struct {
461507
462 // returns -1, 0, 1 if a < b, a == b or a > b respectively.508 // returns -1, 0, 1 if a < b, a == b or a > b respectively.
463 pub fn cmp(a: Int, b: Int) i8 {509 pub fn cmp(a: Int, b: Int) i8 {
464 if (a.positive != b.positive) {510 if (a.isPositive() != b.isPositive()) {
465 return if (a.positive) i8(1) else -1;511 return if (a.isPositive()) i8(1) else -1;
466 } else {512 } else {
467 const r = cmpAbs(a, b);513 const r = cmpAbs(a, b);
468 return if (a.positive) r else -r;514 return if (a.isPositive()) r else -r;
469 }515 }
470 }516 }
471517
472 // if a == 0518 // if a == 0
473 pub fn eqZero(a: Int) bool {519 pub fn eqZero(a: Int) bool {
474 return a.len == 1 and a.limbs[0] == 0;520 return a.len() == 1 and a.limbs[0] == 0;
475 }521 }
476522
477 // if |a| == |b|523 // if |a| == |b|
...@@ -484,28 +530,12 @@ pub const Int = struct {...@@ -484,28 +530,12 @@ pub const Int = struct {
484 return cmp(a, b) == 0;530 return cmp(a, b) == 0;
485 }531 }
486532
487 // Normalize for a possible single carry digit.
488 //
489 // [1, 2, 3, 4, 0] -> [1, 2, 3, 4]
490 // [1, 2, 3, 4, 5] -> [1, 2, 3, 4, 5]
491 // [0] -> [0]
492 fn norm1(r: *Int, length: usize) void {
493 debug.assert(length > 0);
494 debug.assert(length <= r.limbs.len);
495
496 if (r.limbs[length - 1] == 0) {
497 r.len = if (length > 1) length - 1 else 1;
498 } else {
499 r.len = length;
500 }
501 }
502
503 // Normalize a possible sequence of leading zeros.533 // Normalize a possible sequence of leading zeros.
504 //534 //
505 // [1, 2, 3, 4, 0] -> [1, 2, 3, 4]535 // [1, 2, 3, 4, 0] -> [1, 2, 3, 4]
506 // [1, 2, 0, 0, 0] -> [1, 2]536 // [1, 2, 0, 0, 0] -> [1, 2]
507 // [0, 0, 0, 0, 0] -> [0]537 // [0, 0, 0, 0, 0] -> [0]
508 fn normN(r: *Int, length: usize) void {538 fn normalize(r: *Int, length: usize) void {
509 debug.assert(length > 0);539 debug.assert(length > 0);
510 debug.assert(length <= r.limbs.len);540 debug.assert(length <= r.limbs.len);
511541
...@@ -517,11 +547,21 @@ pub const Int = struct {...@@ -517,11 +547,21 @@ pub const Int = struct {
517 }547 }
518548
519 // Handle zero549 // Handle zero
520 r.len = if (j != 0) j else 1;550 r.setLen(if (j != 0) j else 1);
551 }
552
553 // Cannot be used as a result argument to any function.
554 fn readOnlyPositive(a: Int) Int {
555 return Int{
556 .allocator = null,
557 .metadata = a.len(),
558 .limbs = a.limbs,
559 };
521 }560 }
522561
523 // r = a + b562 // r = a + b
524 pub fn add(r: *Int, a: Int, b: Int) Allocator.Error!void {563 pub fn add(r: *Int, a: Int, b: Int) Allocator.Error!void {
564 r.assertWritable();
525 if (a.eqZero()) {565 if (a.eqZero()) {
526 try r.copy(b);566 try r.copy(b);
527 return;567 return;
...@@ -530,38 +570,26 @@ pub const Int = struct {...@@ -530,38 +570,26 @@ pub const Int = struct {
530 return;570 return;
531 }571 }
532572
533 if (a.positive != b.positive) {573 if (a.isPositive() != b.isPositive()) {
534 if (a.positive) {574 if (a.isPositive()) {
535 // (a) + (-b) => a - b575 // (a) + (-b) => a - b
536 const bp = Int{576 try r.sub(a, readOnlyPositive(b));
537 .allocator = undefined,
538 .positive = true,
539 .limbs = b.limbs,
540 .len = b.len,
541 };
542 try r.sub(a, bp);
543 } else {577 } else {
544 // (-a) + (b) => b - a578 // (-a) + (b) => b - a
545 const ap = Int{579 try r.sub(b, readOnlyPositive(a));
546 .allocator = undefined,
547 .positive = true,
548 .limbs = a.limbs,
549 .len = a.len,
550 };
551 try r.sub(b, ap);
552 }580 }
553 } else {581 } else {
554 if (a.len >= b.len) {582 if (a.len() >= b.len()) {
555 try r.ensureCapacity(a.len + 1);583 try r.ensureCapacity(a.len() + 1);
556 lladd(r.limbs[0..], a.limbs[0..a.len], b.limbs[0..b.len]);584 lladd(r.limbs[0..], a.limbs[0..a.len()], b.limbs[0..b.len()]);
557 r.norm1(a.len + 1);585 r.normalize(a.len() + 1);
558 } else {586 } else {
559 try r.ensureCapacity(b.len + 1);587 try r.ensureCapacity(b.len() + 1);
560 lladd(r.limbs[0..], b.limbs[0..b.len], a.limbs[0..a.len]);588 lladd(r.limbs[0..], b.limbs[0..b.len()], a.limbs[0..a.len()]);
561 r.norm1(b.len + 1);589 r.normalize(b.len() + 1);
562 }590 }
563591
564 r.positive = a.positive;592 r.setSign(a.isPositive());
565 }593 }
566 }594 }
567595
...@@ -591,53 +619,42 @@ pub const Int = struct {...@@ -591,53 +619,42 @@ pub const Int = struct {
591619
592 // r = a - b620 // r = a - b
593 pub fn sub(r: *Int, a: Int, b: Int) !void {621 pub fn sub(r: *Int, a: Int, b: Int) !void {
594 if (a.positive != b.positive) {622 r.assertWritable();
595 if (a.positive) {623 if (a.isPositive() != b.isPositive()) {
624 if (a.isPositive()) {
596 // (a) - (-b) => a + b625 // (a) - (-b) => a + b
597 const bp = Int{626 try r.add(a, readOnlyPositive(b));
598 .allocator = undefined,
599 .positive = true,
600 .limbs = b.limbs,
601 .len = b.len,
602 };
603 try r.add(a, bp);
604 } else {627 } else {
605 // (-a) - (b) => -(a + b)628 // (-a) - (b) => -(a + b)
606 const ap = Int{629 try r.add(readOnlyPositive(a), b);
607 .allocator = undefined,630 r.setSign(false);
608 .positive = true,
609 .limbs = a.limbs,
610 .len = a.len,
611 };
612 try r.add(ap, b);
613 r.positive = false;
614 }631 }
615 } else {632 } else {
616 if (a.positive) {633 if (a.isPositive()) {
617 // (a) - (b) => a - b634 // (a) - (b) => a - b
618 if (a.cmp(b) >= 0) {635 if (a.cmp(b) >= 0) {
619 try r.ensureCapacity(a.len + 1);636 try r.ensureCapacity(a.len() + 1);
620 llsub(r.limbs[0..], a.limbs[0..a.len], b.limbs[0..b.len]);637 llsub(r.limbs[0..], a.limbs[0..a.len()], b.limbs[0..b.len()]);
621 r.normN(a.len);638 r.normalize(a.len());
622 r.positive = true;639 r.setSign(true);
623 } else {640 } else {
624 try r.ensureCapacity(b.len + 1);641 try r.ensureCapacity(b.len() + 1);
625 llsub(r.limbs[0..], b.limbs[0..b.len], a.limbs[0..a.len]);642 llsub(r.limbs[0..], b.limbs[0..b.len()], a.limbs[0..a.len()]);
626 r.normN(b.len);643 r.normalize(b.len());
627 r.positive = false;644 r.setSign(false);
628 }645 }
629 } else {646 } else {
630 // (-a) - (-b) => -(a - b)647 // (-a) - (-b) => -(a - b)
631 if (a.cmp(b) < 0) {648 if (a.cmp(b) < 0) {
632 try r.ensureCapacity(a.len + 1);649 try r.ensureCapacity(a.len() + 1);
633 llsub(r.limbs[0..], a.limbs[0..a.len], b.limbs[0..b.len]);650 llsub(r.limbs[0..], a.limbs[0..a.len()], b.limbs[0..b.len()]);
634 r.normN(a.len);651 r.normalize(a.len());
635 r.positive = false;652 r.setSign(false);
636 } else {653 } else {
637 try r.ensureCapacity(b.len + 1);654 try r.ensureCapacity(b.len() + 1);
638 llsub(r.limbs[0..], b.limbs[0..b.len], a.limbs[0..a.len]);655 llsub(r.limbs[0..], b.limbs[0..b.len()], a.limbs[0..a.len()]);
639 r.normN(b.len);656 r.normalize(b.len());
640 r.positive = true;657 r.setSign(true);
641 }658 }
642 }659 }
643 }660 }
...@@ -671,12 +688,14 @@ pub const Int = struct {...@@ -671,12 +688,14 @@ pub const Int = struct {
671 //688 //
672 // For greatest efficiency, ensure rma does not alias a or b.689 // For greatest efficiency, ensure rma does not alias a or b.
673 pub fn mul(rma: *Int, a: Int, b: Int) !void {690 pub fn mul(rma: *Int, a: Int, b: Int) !void {
691 rma.assertWritable();
692
674 var r = rma;693 var r = rma;
675 var aliased = rma.limbs.ptr == a.limbs.ptr or rma.limbs.ptr == b.limbs.ptr;694 var aliased = rma.limbs.ptr == a.limbs.ptr or rma.limbs.ptr == b.limbs.ptr;
676695
677 var sr: Int = undefined;696 var sr: Int = undefined;
678 if (aliased) {697 if (aliased) {
679 sr = try Int.initCapacity(rma.allocator, a.len + b.len);698 sr = try Int.initCapacity(rma.allocator.?, a.len() + b.len());
680 r = &sr;699 r = &sr;
681 aliased = true;700 aliased = true;
682 }701 }
...@@ -685,16 +704,16 @@ pub const Int = struct {...@@ -685,16 +704,16 @@ pub const Int = struct {
685 r.deinit();704 r.deinit();
686 };705 };
687706
688 try r.ensureCapacity(a.len + b.len);707 try r.ensureCapacity(a.len() + b.len());
689708
690 if (a.len >= b.len) {709 if (a.len() >= b.len()) {
691 llmul(r.limbs, a.limbs[0..a.len], b.limbs[0..b.len]);710 llmul(r.limbs, a.limbs[0..a.len()], b.limbs[0..b.len()]);
692 } else {711 } else {
693 llmul(r.limbs, b.limbs[0..b.len], a.limbs[0..a.len]);712 llmul(r.limbs, b.limbs[0..b.len()], a.limbs[0..a.len()]);
694 }713 }
695714
696 r.positive = a.positive == b.positive;715 r.normalize(a.len() + b.len());
697 r.normN(a.len + b.len);716 r.setSign(a.isPositive() == b.isPositive());
698 }717 }
699718
700 // a + b * c + *carry, sets carry to the overflow bits719 // a + b * c + *carry, sets carry to the overflow bits
...@@ -744,25 +763,24 @@ pub const Int = struct {...@@ -744,25 +763,24 @@ pub const Int = struct {
744 try div(q, r, a, b);763 try div(q, r, a, b);
745764
746 // Trunc -> Floor.765 // Trunc -> Floor.
747 if (!q.positive) {766 if (!q.isPositive()) {
748 // TODO values less than limb size should guarantee non allocating767 const one = Int.initFixed(([]Limb{1})[0..]);
749 var one_buffer: [512]u8 = undefined;768 try q.sub(q.*, one);
750 const one_al = &std.heap.FixedBufferAllocator.init(one_buffer[0..]).allocator;769 try r.add(q.*, one);
751 const one_ap = try Int.initSet(one_al, 1);
752
753 try q.sub(q.*, one_ap);
754 try r.add(q.*, one_ap);
755 }770 }
756 r.positive = b.positive;771 r.setSign(b.isPositive());
757 }772 }
758773
759 pub fn divTrunc(q: *Int, r: *Int, a: Int, b: Int) !void {774 pub fn divTrunc(q: *Int, r: *Int, a: Int, b: Int) !void {
760 try div(q, r, a, b);775 try div(q, r, a, b);
761 r.positive = a.positive;776 r.setSign(a.isPositive());
762 }777 }
763778
764 // Truncates by default.779 // Truncates by default.
765 fn div(quo: *Int, rem: *Int, a: Int, b: Int) !void {780 fn div(quo: *Int, rem: *Int, a: Int, b: Int) !void {
781 quo.assertWritable();
782 rem.assertWritable();
783
766 if (b.eqZero()) {784 if (b.eqZero()) {
767 @panic("division by zero");785 @panic("division by zero");
768 }786 }
...@@ -773,36 +791,67 @@ pub const Int = struct {...@@ -773,36 +791,67 @@ pub const Int = struct {
773 if (a.cmpAbs(b) < 0) {791 if (a.cmpAbs(b) < 0) {
774 // quo may alias a so handle rem first792 // quo may alias a so handle rem first
775 try rem.copy(a);793 try rem.copy(a);
776 rem.positive = a.positive == b.positive;794 rem.setSign(a.isPositive() == b.isPositive());
777795
778 quo.positive = true;796 quo.metadata = 1;
779 quo.len = 1;
780 quo.limbs[0] = 0;797 quo.limbs[0] = 0;
781 return;798 return;
782 }799 }
783800
784 if (b.len == 1) {801 // Handle trailing zero-words of divisor/dividend. These are not handled in the following
785 try quo.ensureCapacity(a.len);802 // algorithms.
803 const a_zero_limb_count = blk: {
804 var i: usize = 0;
805 while (i < a.len()) : (i += 1) {
806 if (a.limbs[i] != 0) break;
807 }
808 break :blk i;
809 };
810 const b_zero_limb_count = blk: {
811 var i: usize = 0;
812 while (i < b.len()) : (i += 1) {
813 if (b.limbs[i] != 0) break;
814 }
815 break :blk i;
816 };
817
818 const ab_zero_limb_count = std.math.min(a_zero_limb_count, b_zero_limb_count);
819
820 if (b.len() - ab_zero_limb_count == 1) {
821 try quo.ensureCapacity(a.len());
786822
787 lldiv1(quo.limbs[0..], &rem.limbs[0], a.limbs[0..a.len], b.limbs[0]);823 lldiv1(quo.limbs[0..], &rem.limbs[0], a.limbs[ab_zero_limb_count..a.len()], b.limbs[b.len() - 1]);
788 quo.norm1(a.len);824 quo.normalize(a.len() - ab_zero_limb_count);
789 quo.positive = a.positive == b.positive;825 quo.setSign(a.isPositive() == b.isPositive());
790826
791 rem.len = 1;827 rem.metadata = 1;
792 rem.positive = true;
793 } else {828 } else {
794 // x and y are modified during division829 // x and y are modified during division
795 var x = try a.clone();830 var x = try Int.initCapacity(quo.allocator.?, a.len());
796 defer x.deinit();831 defer x.deinit();
832 try x.copy(a);
797833
798 var y = try b.clone();834 var y = try Int.initCapacity(quo.allocator.?, b.len());
799 defer y.deinit();835 defer y.deinit();
836 try y.copy(b);
800837
801 // x may grow one limb during normalization838 // x may grow one limb during normalization
802 try quo.ensureCapacity(a.len + y.len);839 try quo.ensureCapacity(a.len() + y.len());
803 try divN(quo.allocator, quo, rem, &x, &y);840
841 // Shrink x, y such that the trailing zero limbs shared between are removed.
842 if (ab_zero_limb_count != 0) {
843 std.mem.copy(Limb, x.limbs[0..], x.limbs[ab_zero_limb_count..]);
844 std.mem.copy(Limb, y.limbs[0..], y.limbs[ab_zero_limb_count..]);
845 x.metadata -= ab_zero_limb_count;
846 y.metadata -= ab_zero_limb_count;
847 }
848
849 try divN(quo.allocator.?, quo, rem, &x, &y);
850 quo.setSign(a.isPositive() == b.isPositive());
851 }
804852
805 quo.positive = a.positive == b.positive;853 if (ab_zero_limb_count != 0) {
854 try rem.shiftLeft(rem.*, ab_zero_limb_count * Limb.bit_count);
806 }855 }
807 }856 }
808857
...@@ -837,25 +886,28 @@ pub const Int = struct {...@@ -837,25 +886,28 @@ pub const Int = struct {
837 //886 //
838 // x = qy + r where 0 <= r < y887 // x = qy + r where 0 <= r < y
839 fn divN(allocator: *Allocator, q: *Int, r: *Int, x: *Int, y: *Int) !void {888 fn divN(allocator: *Allocator, q: *Int, r: *Int, x: *Int, y: *Int) !void {
840 debug.assert(y.len >= 2);889 debug.assert(y.len() >= 2);
841 debug.assert(x.len >= y.len);890 debug.assert(x.len() >= y.len());
842 debug.assert(q.limbs.len >= x.len + y.len - 1);891 debug.assert(q.limbs.len >= x.len() + y.len() - 1);
843 debug.assert(default_capacity >= 3); // see 3.2892 debug.assert(default_capacity >= 3); // see 3.2
844893
845 var tmp = try Int.init(allocator);894 var tmp = try Int.init(allocator);
846 defer tmp.deinit();895 defer tmp.deinit();
847896
848 // Normalize so y > Limb.bit_count / 2 (i.e. leading bit is set)897 // Normalize so y > Limb.bit_count / 2 (i.e. leading bit is set) and even
849 const norm_shift = @clz(y.limbs[y.len - 1]);898 var norm_shift = @clz(y.limbs[y.len() - 1]);
899 if (norm_shift == 0 and y.isOdd()) {
900 norm_shift = Limb.bit_count;
901 }
850 try x.shiftLeft(x.*, norm_shift);902 try x.shiftLeft(x.*, norm_shift);
851 try y.shiftLeft(y.*, norm_shift);903 try y.shiftLeft(y.*, norm_shift);
852904
853 const n = x.len - 1;905 const n = x.len() - 1;
854 const t = y.len - 1;906 const t = y.len() - 1;
855907
856 // 1.908 // 1.
857 q.len = n - t + 1;909 q.metadata = n - t + 1;
858 mem.set(Limb, q.limbs[0..q.len], 0);910 mem.set(Limb, q.limbs[0..q.len()], 0);
859911
860 // 2.912 // 2.
861 try tmp.shiftLeft(y.*, Limb.bit_count * (n - t));913 try tmp.shiftLeft(y.*, Limb.bit_count * (n - t));
...@@ -880,7 +932,7 @@ pub const Int = struct {...@@ -880,7 +932,7 @@ pub const Int = struct {
880 tmp.limbs[0] = if (i >= 2) x.limbs[i - 2] else 0;932 tmp.limbs[0] = if (i >= 2) x.limbs[i - 2] else 0;
881 tmp.limbs[1] = if (i >= 1) x.limbs[i - 1] else 0;933 tmp.limbs[1] = if (i >= 1) x.limbs[i - 1] else 0;
882 tmp.limbs[2] = x.limbs[i];934 tmp.limbs[2] = x.limbs[i];
883 tmp.normN(3);935 tmp.normalize(3);
884936
885 while (true) {937 while (true) {
886 // 2x1 limb multiplication unrolled against single-limb q[i-t-1]938 // 2x1 limb multiplication unrolled against single-limb q[i-t-1]
...@@ -888,7 +940,7 @@ pub const Int = struct {...@@ -888,7 +940,7 @@ pub const Int = struct {
888 r.limbs[0] = addMulLimbWithCarry(0, if (t >= 1) y.limbs[t - 1] else 0, q.limbs[i - t - 1], &carry);940 r.limbs[0] = addMulLimbWithCarry(0, if (t >= 1) y.limbs[t - 1] else 0, q.limbs[i - t - 1], &carry);
889 r.limbs[1] = addMulLimbWithCarry(0, y.limbs[t], q.limbs[i - t - 1], &carry);941 r.limbs[1] = addMulLimbWithCarry(0, y.limbs[t], q.limbs[i - t - 1], &carry);
890 r.limbs[2] = carry;942 r.limbs[2] = carry;
891 r.normN(3);943 r.normalize(3);
892944
893 if (r.cmpAbs(tmp) <= 0) {945 if (r.cmpAbs(tmp) <= 0) {
894 break;946 break;
...@@ -903,7 +955,7 @@ pub const Int = struct {...@@ -903,7 +955,7 @@ pub const Int = struct {
903 try tmp.shiftLeft(tmp, Limb.bit_count * (i - t - 1));955 try tmp.shiftLeft(tmp, Limb.bit_count * (i - t - 1));
904 try x.sub(x.*, tmp);956 try x.sub(x.*, tmp);
905957
906 if (!x.positive) {958 if (!x.isPositive()) {
907 try tmp.shiftLeft(y.*, Limb.bit_count * (i - t - 1));959 try tmp.shiftLeft(y.*, Limb.bit_count * (i - t - 1));
908 try x.add(x.*, tmp);960 try x.add(x.*, tmp);
909 q.limbs[i - t - 1] -= 1;961 q.limbs[i - t - 1] -= 1;
...@@ -911,18 +963,20 @@ pub const Int = struct {...@@ -911,18 +963,20 @@ pub const Int = struct {
911 }963 }
912964
913 // Denormalize965 // Denormalize
914 q.normN(q.len);966 q.normalize(q.len());
915967
916 try r.shiftRight(x.*, norm_shift);968 try r.shiftRight(x.*, norm_shift);
917 r.normN(r.len);969 r.normalize(r.len());
918 }970 }
919971
920 // r = a << shift, in other words, r = a * 2^shift972 // r = a << shift, in other words, r = a * 2^shift
921 pub fn shiftLeft(r: *Int, a: Int, shift: usize) !void {973 pub fn shiftLeft(r: *Int, a: Int, shift: usize) !void {
922 try r.ensureCapacity(a.len + (shift / Limb.bit_count) + 1);974 r.assertWritable();
923 llshl(r.limbs[0..], a.limbs[0..a.len], shift);975
924 r.norm1(a.len + (shift / Limb.bit_count) + 1);976 try r.ensureCapacity(a.len() + (shift / Limb.bit_count) + 1);
925 r.positive = a.positive;977 llshl(r.limbs[0..], a.limbs[0..a.len()], shift);
978 r.normalize(a.len() + (shift / Limb.bit_count) + 1);
979 r.setSign(a.isPositive());
926 }980 }
927981
928 fn llshl(r: []Limb, a: []const Limb, shift: usize) void {982 fn llshl(r: []Limb, a: []const Limb, shift: usize) void {
...@@ -950,17 +1004,18 @@ pub const Int = struct {...@@ -950,17 +1004,18 @@ pub const Int = struct {
9501004
951 // r = a >> shift1005 // r = a >> shift
952 pub fn shiftRight(r: *Int, a: Int, shift: usize) !void {1006 pub fn shiftRight(r: *Int, a: Int, shift: usize) !void {
953 if (a.len <= shift / Limb.bit_count) {1007 r.assertWritable();
954 r.len = 1;1008
1009 if (a.len() <= shift / Limb.bit_count) {
1010 r.metadata = 1;
955 r.limbs[0] = 0;1011 r.limbs[0] = 0;
956 r.positive = true;
957 return;1012 return;
958 }1013 }
9591014
960 try r.ensureCapacity(a.len - (shift / Limb.bit_count));1015 try r.ensureCapacity(a.len() - (shift / Limb.bit_count));
961 const r_len = llshr(r.limbs[0..], a.limbs[0..a.len], shift);1016 const r_len = llshr(r.limbs[0..], a.limbs[0..a.len()], shift);
962 r.len = a.len - (shift / Limb.bit_count);1017 r.metadata = a.len() - (shift / Limb.bit_count);
963 r.positive = a.positive;1018 r.setSign(a.isPositive());
964 }1019 }
9651020
966 fn llshr(r: []Limb, a: []const Limb, shift: usize) void {1021 fn llshr(r: []Limb, a: []const Limb, shift: usize) void {
...@@ -985,14 +1040,16 @@ pub const Int = struct {...@@ -985,14 +1040,16 @@ pub const Int = struct {
9851040
986 // r = a | b1041 // r = a | b
987 pub fn bitOr(r: *Int, a: Int, b: Int) !void {1042 pub fn bitOr(r: *Int, a: Int, b: Int) !void {
988 if (a.len > b.len) {1043 r.assertWritable();
989 try r.ensureCapacity(a.len);1044
990 llor(r.limbs[0..], a.limbs[0..a.len], b.limbs[0..b.len]);1045 if (a.len() > b.len()) {
991 r.len = a.len;1046 try r.ensureCapacity(a.len());
1047 llor(r.limbs[0..], a.limbs[0..a.len()], b.limbs[0..b.len()]);
1048 r.setLen(a.len());
992 } else {1049 } else {
993 try r.ensureCapacity(b.len);1050 try r.ensureCapacity(b.len());
994 llor(r.limbs[0..], b.limbs[0..b.len], a.limbs[0..a.len]);1051 llor(r.limbs[0..], b.limbs[0..b.len()], a.limbs[0..a.len()]);
995 r.len = b.len;1052 r.setLen(b.len());
996 }1053 }
997 }1054 }
9981055
...@@ -1012,14 +1069,16 @@ pub const Int = struct {...@@ -1012,14 +1069,16 @@ pub const Int = struct {
10121069
1013 // r = a & b1070 // r = a & b
1014 pub fn bitAnd(r: *Int, a: Int, b: Int) !void {1071 pub fn bitAnd(r: *Int, a: Int, b: Int) !void {
1015 if (a.len > b.len) {1072 r.assertWritable();
1016 try r.ensureCapacity(b.len);1073
1017 lland(r.limbs[0..], a.limbs[0..a.len], b.limbs[0..b.len]);1074 if (a.len() > b.len()) {
1018 r.normN(b.len);1075 try r.ensureCapacity(b.len());
1076 lland(r.limbs[0..], a.limbs[0..a.len()], b.limbs[0..b.len()]);
1077 r.normalize(b.len());
1019 } else {1078 } else {
1020 try r.ensureCapacity(a.len);1079 try r.ensureCapacity(a.len());
1021 lland(r.limbs[0..], b.limbs[0..b.len], a.limbs[0..a.len]);1080 lland(r.limbs[0..], b.limbs[0..b.len()], a.limbs[0..a.len()]);
1022 r.normN(a.len);1081 r.normalize(a.len());
1023 }1082 }
1024 }1083 }
10251084
...@@ -1036,14 +1095,16 @@ pub const Int = struct {...@@ -1036,14 +1095,16 @@ pub const Int = struct {
10361095
1037 // r = a ^ b1096 // r = a ^ b
1038 pub fn bitXor(r: *Int, a: Int, b: Int) !void {1097 pub fn bitXor(r: *Int, a: Int, b: Int) !void {
1039 if (a.len > b.len) {1098 r.assertWritable();
1040 try r.ensureCapacity(a.len);1099
1041 llxor(r.limbs[0..], a.limbs[0..a.len], b.limbs[0..b.len]);1100 if (a.len() > b.len()) {
1042 r.normN(a.len);1101 try r.ensureCapacity(a.len());
1102 llxor(r.limbs[0..], a.limbs[0..a.len()], b.limbs[0..b.len()]);
1103 r.normalize(a.len());
1043 } else {1104 } else {
1044 try r.ensureCapacity(b.len);1105 try r.ensureCapacity(b.len());
1045 llxor(r.limbs[0..], b.limbs[0..b.len], a.limbs[0..a.len]);1106 llxor(r.limbs[0..], b.limbs[0..b.len()], a.limbs[0..a.len()]);
1046 r.normN(b.len);1107 r.normalize(b.len());
1047 }1108 }
1048 }1109 }
10491110
...@@ -1067,7 +1128,9 @@ pub const Int = struct {...@@ -1067,7 +1128,9 @@ pub const Int = struct {
1067// They will still run on larger than this and should pass, but the multi-limb code-paths1128// They will still run on larger than this and should pass, but the multi-limb code-paths
1068// may be untested in some cases.1129// may be untested in some cases.
10691130
1070const al = debug.global_allocator;1131var buffer: [64 * 8192]u8 = undefined;
1132var fixed = std.heap.FixedBufferAllocator.init(buffer[0..]);
1133const al = &fixed.allocator;
10711134
1072test "big.int comptime_int set" {1135test "big.int comptime_int set" {
1073 comptime var s = 0xefffffff00000001eeeeeeefaaaaaaab;1136 comptime var s = 0xefffffff00000001eeeeeeefaaaaaaab;
...@@ -1088,14 +1151,14 @@ test "big.int comptime_int set negative" {...@@ -1088,14 +1151,14 @@ test "big.int comptime_int set negative" {
1088 var a = try Int.initSet(al, -10);1151 var a = try Int.initSet(al, -10);
10891152
1090 testing.expect(a.limbs[0] == 10);1153 testing.expect(a.limbs[0] == 10);
1091 testing.expect(a.positive == false);1154 testing.expect(a.isPositive() == false);
1092}1155}
10931156
1094test "big.int int set unaligned small" {1157test "big.int int set unaligned small" {
1095 var a = try Int.initSet(al, u7(45));1158 var a = try Int.initSet(al, u7(45));
10961159
1097 testing.expect(a.limbs[0] == 45);1160 testing.expect(a.limbs[0] == 45);
1098 testing.expect(a.positive == true);1161 testing.expect(a.isPositive() == true);
1099}1162}
11001163
1101test "big.int comptime_int to" {1164test "big.int comptime_int to" {
...@@ -1116,7 +1179,7 @@ test "big.int to target too small error" {...@@ -1116,7 +1179,7 @@ test "big.int to target too small error" {
1116 testing.expectError(error.TargetTooSmall, a.to(u8));1179 testing.expectError(error.TargetTooSmall, a.to(u8));
1117}1180}
11181181
1119test "big.int norm1" {1182test "big.int normalize" {
1120 var a = try Int.init(al);1183 var a = try Int.init(al);
1121 try a.ensureCapacity(8);1184 try a.ensureCapacity(8);
11221185
...@@ -1124,26 +1187,26 @@ test "big.int norm1" {...@@ -1124,26 +1187,26 @@ test "big.int norm1" {
1124 a.limbs[1] = 2;1187 a.limbs[1] = 2;
1125 a.limbs[2] = 3;1188 a.limbs[2] = 3;
1126 a.limbs[3] = 0;1189 a.limbs[3] = 0;
1127 a.norm1(4);1190 a.normalize(4);
1128 testing.expect(a.len == 3);1191 testing.expect(a.len() == 3);
11291192
1130 a.limbs[0] = 1;1193 a.limbs[0] = 1;
1131 a.limbs[1] = 2;1194 a.limbs[1] = 2;
1132 a.limbs[2] = 3;1195 a.limbs[2] = 3;
1133 a.norm1(3);1196 a.normalize(3);
1134 testing.expect(a.len == 3);1197 testing.expect(a.len() == 3);
11351198
1136 a.limbs[0] = 0;1199 a.limbs[0] = 0;
1137 a.limbs[1] = 0;1200 a.limbs[1] = 0;
1138 a.norm1(2);1201 a.normalize(2);
1139 testing.expect(a.len == 1);1202 testing.expect(a.len() == 1);
11401203
1141 a.limbs[0] = 0;1204 a.limbs[0] = 0;
1142 a.norm1(1);1205 a.normalize(1);
1143 testing.expect(a.len == 1);1206 testing.expect(a.len() == 1);
1144}1207}
11451208
1146test "big.int normN" {1209test "big.int normalize multi" {
1147 var a = try Int.init(al);1210 var a = try Int.init(al);
1148 try a.ensureCapacity(8);1211 try a.ensureCapacity(8);
11491212
...@@ -1151,25 +1214,25 @@ test "big.int normN" {...@@ -1151,25 +1214,25 @@ test "big.int normN" {
1151 a.limbs[1] = 2;1214 a.limbs[1] = 2;
1152 a.limbs[2] = 0;1215 a.limbs[2] = 0;
1153 a.limbs[3] = 0;1216 a.limbs[3] = 0;
1154 a.normN(4);1217 a.normalize(4);
1155 testing.expect(a.len == 2);1218 testing.expect(a.len() == 2);
11561219
1157 a.limbs[0] = 1;1220 a.limbs[0] = 1;
1158 a.limbs[1] = 2;1221 a.limbs[1] = 2;
1159 a.limbs[2] = 3;1222 a.limbs[2] = 3;
1160 a.normN(3);1223 a.normalize(3);
1161 testing.expect(a.len == 3);1224 testing.expect(a.len() == 3);
11621225
1163 a.limbs[0] = 0;1226 a.limbs[0] = 0;
1164 a.limbs[1] = 0;1227 a.limbs[1] = 0;
1165 a.limbs[2] = 0;1228 a.limbs[2] = 0;
1166 a.limbs[3] = 0;1229 a.limbs[3] = 0;
1167 a.normN(4);1230 a.normalize(4);
1168 testing.expect(a.len == 1);1231 testing.expect(a.len() == 1);
11691232
1170 a.limbs[0] = 0;1233 a.limbs[0] = 0;
1171 a.normN(1);1234 a.normalize(1);
1172 testing.expect(a.len == 1);1235 testing.expect(a.len() == 1);
1173}1236}
11741237
1175test "big.int parity" {1238test "big.int parity" {
...@@ -1204,7 +1267,7 @@ test "big.int bitcount + sizeInBase" {...@@ -1204,7 +1267,7 @@ test "big.int bitcount + sizeInBase" {
1204 try a.shiftLeft(a, 5000);1267 try a.shiftLeft(a, 5000);
1205 testing.expect(a.bitCountAbs() == 5032);1268 testing.expect(a.bitCountAbs() == 5032);
1206 testing.expect(a.sizeInBase(2) >= 5032);1269 testing.expect(a.sizeInBase(2) >= 5032);
1207 a.positive = false;1270 a.setSign(false);
12081271
1209 testing.expect(a.bitCountAbs() == 5032);1272 testing.expect(a.bitCountAbs() == 5032);
1210 testing.expect(a.sizeInBase(2) >= 5033);1273 testing.expect(a.sizeInBase(2) >= 5033);
...@@ -1980,6 +2043,98 @@ test "big.int div multi-multi (3.1/3.3 branch)" {...@@ -1980,6 +2043,98 @@ test "big.int div multi-multi (3.1/3.3 branch)" {
1980 testing.expect((try r.to(u256)) == 0x1111111111111111111110b12222222222222222282);2043 testing.expect((try r.to(u256)) == 0x1111111111111111111110b12222222222222222282);
1981}2044}
19822045
2046test "big.int div multi-single zero-limb trailing" {
2047 var a = try Int.initSet(al, 0x60000000000000000000000000000000000000000000000000000000000000000);
2048 var b = try Int.initSet(al, 0x10000000000000000);
2049
2050 var q = try Int.init(al);
2051 var r = try Int.init(al);
2052 try Int.divTrunc(&q, &r, a, b);
2053
2054 var expected = try Int.initSet(al, 0x6000000000000000000000000000000000000000000000000);
2055 testing.expect(q.eq(expected));
2056 testing.expect(r.eqZero());
2057}
2058
2059test "big.int div multi-multi zero-limb trailing (with rem)" {
2060 var a = try Int.initSet(al, 0x86666666555555558888888777777776111111111111111100000000000000000000000000000000);
2061 var b = try Int.initSet(al, 0x8666666655555555444444443333333300000000000000000000000000000000);
2062
2063 var q = try Int.init(al);
2064 var r = try Int.init(al);
2065 try Int.divTrunc(&q, &r, a, b);
2066
2067 testing.expect((try q.to(u128)) == 0x10000000000000000);
2068
2069 const rs = try r.toString(al, 16);
2070 testing.expect(std.mem.eql(u8, rs, "4444444344444443111111111111111100000000000000000000000000000000"));
2071}
2072
2073test "big.int div multi-multi zero-limb trailing (with rem) and dividend zero-limb count > divisor zero-limb count" {
2074 var a = try Int.initSet(al, 0x8666666655555555888888877777777611111111111111110000000000000000);
2075 var b = try Int.initSet(al, 0x8666666655555555444444443333333300000000000000000000000000000000);
2076
2077 var q = try Int.init(al);
2078 var r = try Int.init(al);
2079 try Int.divTrunc(&q, &r, a, b);
2080
2081 testing.expect((try q.to(u128)) == 0x1);
2082
2083 const rs = try r.toString(al, 16);
2084 testing.expect(std.mem.eql(u8, rs, "444444434444444311111111111111110000000000000000"));
2085}
2086
2087test "big.int div multi-multi zero-limb trailing (with rem) and dividend zero-limb count < divisor zero-limb count" {
2088 var a = try Int.initSet(al, 0x86666666555555558888888777777776111111111111111100000000000000000000000000000000);
2089 var b = try Int.initSet(al, 0x866666665555555544444444333333330000000000000000);
2090
2091 var q = try Int.init(al);
2092 var r = try Int.init(al);
2093 try Int.divTrunc(&q, &r, a, b);
2094
2095 const qs = try q.toString(al, 16);
2096 testing.expect(std.mem.eql(u8, qs, "10000000000000000820820803105186f"));
2097
2098 const rs = try r.toString(al, 16);
2099 testing.expect(std.mem.eql(u8, rs, "4e11f2baa5896a321d463b543d0104e30000000000000000"));
2100}
2101
2102test "big.int div multi-multi fuzz case #1" {
2103 var a = try Int.init(al);
2104 var b = try Int.init(al);
2105
2106 try a.setString(16, "ffffffffffffffffffffffffffffc00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000003ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff80000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000003ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffc00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffe000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000");
2107 try b.setString(16, "3ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffe0000000000000000000000000000000000001ffffffffffffffffffffffffffffffffffffffffffffffffffc000000000000000000000000000000007fffffffffff");
2108
2109 var q = try Int.init(al);
2110 var r = try Int.init(al);
2111 try Int.divTrunc(&q, &r, a, b);
2112
2113 const qs = try q.toString(al, 16);
2114 testing.expect(std.mem.eql(u8, qs, "3ffffffffffffffffffffffffffff0000000000000000000000000000000000001ffffffffffffffffffffffffffff7fffffffe000000000000000000000000000180000000000000000000003fffffbfffffffdfffffffffffffeffff800000100101000000100000000020003fffffdfbfffffe3ffffffffffffeffff7fffc00800a100000017ffe000002000400007efbfff7fe9f00000037ffff3fff7fffa004006100000009ffe00000190038200bf7d2ff7fefe80400060000f7d7f8fbf9401fe38e0403ffc0bdffffa51102c300d7be5ef9df4e5060007b0127ad3fa69f97d0f820b6605ff617ddf7f32ad7a05c0d03f2e7bc78a6000e087a8bbcdc59e07a5a079128a7861f553ddebed7e8e56701756f9ead39b48cd1b0831889ea6ec1fddf643d0565b075ff07e6caea4e2854ec9227fd635ed60a2f5eef2893052ffd54718fa08604acbf6a15e78a467c4a3c53c0278af06c4416573f925491b195e8fd79302cb1aaf7caf4ecfc9aec1254cc969786363ac729f914c6ddcc26738d6b0facd54eba026580aba2eb6482a088b0d224a8852420b91ec1"));
2115
2116 const rs = try r.toString(al, 16);
2117 testing.expect(std.mem.eql(u8, rs, "310d1d4c414426b4836c2635bad1df3a424e50cbdd167ffccb4dfff57d36b4aae0d6ca0910698220171a0f3373c1060a046c2812f0027e321f72979daa5e7973214170d49e885de0c0ecc167837d44502430674a82522e5df6a0759548052420b91ec1"));
2118}
2119
2120test "big.int div multi-multi fuzz case #2" {
2121 var a = try Int.init(al);
2122 var b = try Int.init(al);
2123
2124 try a.setString(16, "3ffffffffe00000000000000000000000000fffffffffffffffffffffffffffffffffffffffffffffffffffffffffe000000000000000000000000000000000000000000000000000000000000001fffffffffffffffff800000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000ffffffffffffffffffffc000000000000000000000000000000000000000000000000000000000000000");
2125 try b.setString(16, "ffc0000000000000000000000000000000000000000000000000");
2126
2127 var q = try Int.init(al);
2128 var r = try Int.init(al);
2129 try Int.divTrunc(&q, &r, a, b);
2130
2131 const qs = try q.toString(al, 16);
2132 testing.expect(std.mem.eql(u8, qs, "40100400fe3f8fe3f8fe3f8fe3f8fe3f8fe4f93e4f93e4f93e4f93e4f93e4f93e4f93e4f93e4f93e4f93e4f93e4f91e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4992649926499264991e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4791e4792e4b92e4b92e4b92e4b92a4a92a4a92a4"));
2133
2134 const rs = try r.toString(al, 16);
2135 testing.expect(std.mem.eql(u8, rs, "a900000000000000000000000000000000000000000000000000"));
2136}
2137
1983test "big.int shift-right single" {2138test "big.int shift-right single" {
1984 var a = try Int.initSet(al, 0xffff0000);2139 var a = try Int.initSet(al, 0xffff0000);
1985 try a.shiftRight(a, 16);2140 try a.shiftRight(a, 16);
std/math/big/rational.zig created+896
...@@ -0,0 +1,896 @@
1const std = @import("../../std.zig");
2const builtin = @import("builtin");
3const debug = std.debug;
4const math = std.math;
5const mem = std.mem;
6const testing = std.testing;
7const Allocator = mem.Allocator;
8const ArrayList = std.ArrayList;
9
10const TypeId = builtin.TypeId;
11
12const bn = @import("int.zig");
13const Limb = bn.Limb;
14const DoubleLimb = bn.DoubleLimb;
15const Int = bn.Int;
16
17pub const Rational = struct {
18 // Sign of Rational is sign of p. Sign of q is ignored
19 p: Int,
20 q: Int,
21
22 pub fn init(a: *Allocator) !Rational {
23 return Rational{
24 .p = try Int.init(a),
25 .q = try Int.initSet(a, 1),
26 };
27 }
28
29 pub fn deinit(self: *Rational) void {
30 self.p.deinit();
31 self.q.deinit();
32 }
33
34 pub fn setInt(self: *Rational, a: var) !void {
35 try self.p.set(a);
36 try self.q.set(1);
37 }
38
39 // TODO: Accept a/b fractions and exponent form
40 pub fn setFloatString(self: *Rational, str: []const u8) !void {
41 if (str.len == 0) {
42 return error.InvalidFloatString;
43 }
44
45 const State = enum {
46 Integer,
47 Fractional,
48 };
49
50 var state = State.Integer;
51 var point: ?usize = null;
52
53 var start: usize = 0;
54 if (str[0] == '-') {
55 start += 1;
56 }
57
58 for (str) |c, i| {
59 switch (state) {
60 State.Integer => {
61 switch (c) {
62 '.' => {
63 state = State.Fractional;
64 point = i;
65 },
66 '0'...'9' => {
67 // okay
68 },
69 else => {
70 return error.InvalidFloatString;
71 },
72 }
73 },
74 State.Fractional => {
75 switch (c) {
76 '0'...'9' => {
77 // okay
78 },
79 else => {
80 return error.InvalidFloatString;
81 },
82 }
83 },
84 }
85 }
86
87 // TODO: batch the multiplies by 10
88 if (point) |i| {
89 try self.p.setString(10, str[0..i]);
90
91 const base = Int.initFixed(([]Limb{10})[0..]);
92
93 var j: usize = start;
94 while (j < str.len - i - 1) : (j += 1) {
95 try self.p.mul(self.p, base);
96 }
97
98 try self.q.setString(10, str[i + 1 ..]);
99 try self.p.add(self.p, self.q);
100
101 try self.q.set(1);
102 var k: usize = i + 1;
103 while (k < str.len) : (k += 1) {
104 try self.q.mul(self.q, base);
105 }
106
107 try self.reduce();
108 } else {
109 try self.p.setString(10, str[0..]);
110 try self.q.set(1);
111 }
112 }
113
114 // Translated from golang.go/src/math/big/rat.go.
115 pub fn setFloat(self: *Rational, comptime T: type, f: T) !void {
116 debug.assert(@typeId(T) == builtin.TypeId.Float);
117
118 const UnsignedIntType = @IntType(false, T.bit_count);
119 const f_bits = @bitCast(UnsignedIntType, f);
120
121 const exponent_bits = math.floatExponentBits(T);
122 const exponent_bias = (1 << (exponent_bits - 1)) - 1;
123 const mantissa_bits = math.floatMantissaBits(T);
124
125 const exponent_mask = (1 << exponent_bits) - 1;
126 const mantissa_mask = (1 << mantissa_bits) - 1;
127
128 var exponent = @intCast(i16, (f_bits >> mantissa_bits) & exponent_mask);
129 var mantissa = f_bits & mantissa_mask;
130
131 switch (exponent) {
132 exponent_mask => {
133 return error.NonFiniteFloat;
134 },
135 0 => {
136 // denormal
137 exponent -= exponent_bias - 1;
138 },
139 else => {
140 // normal
141 mantissa |= 1 << mantissa_bits;
142 exponent -= exponent_bias;
143 },
144 }
145
146 var shift: i16 = mantissa_bits - exponent;
147
148 // factor out powers of two early from rational
149 while (mantissa & 1 == 0 and shift > 0) {
150 mantissa >>= 1;
151 shift -= 1;
152 }
153
154 try self.p.set(mantissa);
155 self.p.setSign(f >= 0);
156
157 try self.q.set(1);
158 if (shift >= 0) {
159 try self.q.shiftLeft(self.q, @intCast(usize, shift));
160 } else {
161 try self.p.shiftLeft(self.p, @intCast(usize, -shift));
162 }
163
164 try self.reduce();
165 }
166
167 // Translated from golang.go/src/math/big/rat.go.
168 pub fn toFloat(self: Rational, comptime T: type) !T {
169 debug.assert(@typeId(T) == builtin.TypeId.Float);
170
171 const fsize = T.bit_count;
172 const BitReprType = @IntType(false, T.bit_count);
173
174 const msize = math.floatMantissaBits(T);
175 const msize1 = msize + 1;
176 const msize2 = msize1 + 1;
177
178 const esize = math.floatExponentBits(T);
179 const ebias = (1 << (esize - 1)) - 1;
180 const emin = 1 - ebias;
181 const emax = ebias;
182
183 if (self.p.eqZero()) {
184 return 0;
185 }
186
187 // 1. left-shift a or sub so that a/b is in [1 << msize1, 1 << (msize2 + 1)]
188 var exp = @intCast(isize, self.p.bitCountTwosComp()) - @intCast(isize, self.q.bitCountTwosComp());
189
190 var a2 = try self.p.clone();
191 defer a2.deinit();
192
193 var b2 = try self.q.clone();
194 defer b2.deinit();
195
196 const shift = msize2 - exp;
197 if (shift >= 0) {
198 try a2.shiftLeft(a2, @intCast(usize, shift));
199 } else {
200 try b2.shiftLeft(b2, @intCast(usize, -shift));
201 }
202
203 // 2. compute quotient and remainder
204 var q = try Int.init(self.p.allocator.?);
205 defer q.deinit();
206
207 // unused
208 var r = try Int.init(self.p.allocator.?);
209 defer r.deinit();
210
211 try Int.divTrunc(&q, &r, a2, b2);
212
213 var mantissa = extractLowBits(q, BitReprType);
214 var have_rem = r.len() > 0;
215
216 // 3. q didn't fit in msize2 bits, redo division b2 << 1
217 if (mantissa >> msize2 == 1) {
218 if (mantissa & 1 == 1) {
219 have_rem = true;
220 }
221 mantissa >>= 1;
222 exp += 1;
223 }
224 if (mantissa >> msize1 != 1) {
225 // NOTE: This can be hit if the limb size is small (u8/16).
226 @panic("unexpected bits in result");
227 }
228
229 // 4. Rounding
230 if (emin - msize <= exp and exp <= emin) {
231 // denormal
232 const shift1 = @intCast(math.Log2Int(BitReprType), emin - (exp - 1));
233 const lost_bits = mantissa & ((@intCast(BitReprType, 1) << shift1) - 1);
234 have_rem = have_rem or lost_bits != 0;
235 mantissa >>= shift1;
236 exp = 2 - ebias;
237 }
238
239 // round q using round-half-to-even
240 var exact = !have_rem;
241 if (mantissa & 1 != 0) {
242 exact = false;
243 if (have_rem or (mantissa & 2 != 0)) {
244 mantissa += 1;
245 if (mantissa >= 1 << msize2) {
246 // 11...1 => 100...0
247 mantissa >>= 1;
248 exp += 1;
249 }
250 }
251 }
252 mantissa >>= 1;
253
254 const f = math.scalbn(@intToFloat(T, mantissa), @intCast(i32, exp - msize1));
255 if (math.isInf(f)) {
256 exact = false;
257 }
258
259 return if (self.p.isPositive()) f else -f;
260 }
261
262 pub fn setRatio(self: *Rational, p: var, q: var) !void {
263 try self.p.set(p);
264 try self.q.set(q);
265
266 self.p.setSign(@boolToInt(self.p.isPositive()) ^ @boolToInt(self.q.isPositive()) == 0);
267 self.q.setSign(true);
268
269 try self.reduce();
270
271 if (self.q.eqZero()) {
272 @panic("cannot set rational with denominator = 0");
273 }
274 }
275
276 pub fn copyInt(self: *Rational, a: Int) !void {
277 try self.p.copy(a);
278 try self.q.set(1);
279 }
280
281 pub fn copyRatio(self: *Rational, a: Int, b: Int) !void {
282 try self.p.copy(a);
283 try self.q.copy(b);
284
285 self.p.setSign(@boolToInt(self.p.isPositive()) ^ @boolToInt(self.q.isPositive()) == 0);
286 self.q.setSign(true);
287
288 try self.reduce();
289 }
290
291 pub fn abs(r: *Rational) void {
292 r.p.abs();
293 }
294
295 pub fn negate(r: *Rational) void {
296 r.p.negate();
297 }
298
299 pub fn swap(r: *Rational, other: *Rational) void {
300 r.p.swap(&other.p);
301 r.q.swap(&other.q);
302 }
303
304 pub fn cmp(a: Rational, b: Rational) !i8 {
305 return cmpInternal(a, b, true);
306 }
307
308 pub fn cmpAbs(a: Rational, b: Rational) !i8 {
309 return cmpInternal(a, b, false);
310 }
311
312 // p/q > x/y iff p*y > x*q
313 fn cmpInternal(a: Rational, b: Rational, is_abs: bool) !i8 {
314 // TODO: Would a div compare algorithm of sorts be viable and quicker? Can we avoid
315 // the memory allocations here?
316 var q = try Int.init(a.p.allocator.?);
317 defer q.deinit();
318
319 var p = try Int.init(b.p.allocator.?);
320 defer p.deinit();
321
322 try q.mul(a.p, b.q);
323 try p.mul(b.p, a.q);
324
325 return if (is_abs) q.cmpAbs(p) else q.cmp(p);
326 }
327
328 // r/q = ap/aq + bp/bq = (ap*bq + bp*aq) / (aq*bq)
329 //
330 // For best performance, rma should not alias a or b.
331 pub fn add(rma: *Rational, a: Rational, b: Rational) !void {
332 var r = rma;
333 var aliased = rma.p.limbs.ptr == a.p.limbs.ptr or rma.p.limbs.ptr == b.p.limbs.ptr;
334
335 var sr: Rational = undefined;
336 if (aliased) {
337 sr = try Rational.init(rma.p.allocator.?);
338 r = &sr;
339 aliased = true;
340 }
341 defer if (aliased) {
342 rma.swap(r);
343 r.deinit();
344 };
345
346 try r.p.mul(a.p, b.q);
347 try r.q.mul(b.p, a.q);
348 try r.p.add(r.p, r.q);
349
350 try r.q.mul(a.q, b.q);
351 try r.reduce();
352 }
353
354 // r/q = ap/aq - bp/bq = (ap*bq - bp*aq) / (aq*bq)
355 //
356 // For best performance, rma should not alias a or b.
357 pub fn sub(rma: *Rational, a: Rational, b: Rational) !void {
358 var r = rma;
359 var aliased = rma.p.limbs.ptr == a.p.limbs.ptr or rma.p.limbs.ptr == b.p.limbs.ptr;
360
361 var sr: Rational = undefined;
362 if (aliased) {
363 sr = try Rational.init(rma.p.allocator.?);
364 r = &sr;
365 aliased = true;
366 }
367 defer if (aliased) {
368 rma.swap(r);
369 r.deinit();
370 };
371
372 try r.p.mul(a.p, b.q);
373 try r.q.mul(b.p, a.q);
374 try r.p.sub(r.p, r.q);
375
376 try r.q.mul(a.q, b.q);
377 try r.reduce();
378 }
379
380 // r/q = ap/aq * bp/bq = ap*bp / aq*bq
381 pub fn mul(r: *Rational, a: Rational, b: Rational) !void {
382 try r.p.mul(a.p, b.p);
383 try r.q.mul(a.q, b.q);
384 try r.reduce();
385 }
386
387 // r/q = (ap/aq) / (bp/bq) = ap*bq / bp*aq
388 pub fn div(r: *Rational, a: Rational, b: Rational) !void {
389 if (b.p.eqZero()) {
390 @panic("division by zero");
391 }
392
393 try r.p.mul(a.p, b.q);
394 try r.q.mul(b.p, a.q);
395 try r.reduce();
396 }
397
398 // r/q = q/r
399 pub fn invert(r: *Rational) void {
400 Int.swap(&r.p, &r.q);
401 }
402
403 // reduce r/q such that gcd(r, q) = 1
404 fn reduce(r: *Rational) !void {
405 var a = try Int.init(r.p.allocator.?);
406 defer a.deinit();
407
408 const sign = r.p.isPositive();
409 r.p.abs();
410 try gcd(&a, r.p, r.q);
411 r.p.setSign(sign);
412
413 const one = Int.initFixed(([]Limb{1})[0..]);
414 if (a.cmp(one) != 0) {
415 var unused = try Int.init(r.p.allocator.?);
416 defer unused.deinit();
417
418 // TODO: divexact would be useful here
419 // TODO: don't copy r.q for div
420 try Int.divTrunc(&r.p, &unused, r.p, a);
421 try Int.divTrunc(&r.q, &unused, r.q, a);
422 }
423 }
424};
425
426const SignedDoubleLimb = @IntType(true, DoubleLimb.bit_count);
427
428fn gcd(rma: *Int, x: Int, y: Int) !void {
429 rma.assertWritable();
430 var r = rma;
431 var aliased = rma.limbs.ptr == x.limbs.ptr or rma.limbs.ptr == y.limbs.ptr;
432
433 var sr: Int = undefined;
434 if (aliased) {
435 sr = try Int.initCapacity(rma.allocator.?, math.max(x.len(), y.len()));
436 r = &sr;
437 aliased = true;
438 }
439 defer if (aliased) {
440 rma.swap(r);
441 r.deinit();
442 };
443
444 try gcdLehmer(r, x, y);
445}
446
447// Storage must live for the lifetime of the returned value
448fn FixedIntFromSignedDoubleLimb(A: SignedDoubleLimb, storage: []Limb) Int {
449 std.debug.assert(storage.len >= 2);
450
451 var A_is_positive = A >= 0;
452 const Au = @intCast(DoubleLimb, if (A < 0) -A else A);
453 storage[0] = @truncate(Limb, Au);
454 storage[1] = @truncate(Limb, Au >> Limb.bit_count);
455 var Ap = Int.initFixed(storage[0..2]);
456 Ap.setSign(A_is_positive);
457 return Ap;
458}
459
460fn gcdLehmer(r: *Int, xa: Int, ya: Int) !void {
461 var x = try xa.clone();
462 x.abs();
463 defer x.deinit();
464
465 var y = try ya.clone();
466 y.abs();
467 defer y.deinit();
468
469 if (x.cmp(y) < 0) {
470 x.swap(&y);
471 }
472
473 var T = try Int.init(r.allocator.?);
474 defer T.deinit();
475
476 while (y.len() > 1) {
477 debug.assert(x.isPositive() and y.isPositive());
478 debug.assert(x.len() >= y.len());
479
480 var xh: SignedDoubleLimb = x.limbs[x.len() - 1];
481 var yh: SignedDoubleLimb = if (x.len() > y.len()) 0 else y.limbs[x.len() - 1];
482
483 var A: SignedDoubleLimb = 1;
484 var B: SignedDoubleLimb = 0;
485 var C: SignedDoubleLimb = 0;
486 var D: SignedDoubleLimb = 1;
487
488 while (yh + C != 0 and yh + D != 0) {
489 const q = @divFloor(xh + A, yh + C);
490 const qp = @divFloor(xh + B, yh + D);
491 if (q != qp) {
492 break;
493 }
494
495 var t = A - q * C;
496 A = C;
497 C = t;
498 t = B - q * D;
499 B = D;
500 D = t;
501
502 t = xh - q * yh;
503 xh = yh;
504 yh = t;
505 }
506
507 if (B == 0) {
508 // T = x % y, r is unused
509 try Int.divTrunc(r, &T, x, y);
510 debug.assert(T.isPositive());
511
512 x.swap(&y);
513 y.swap(&T);
514 } else {
515 var storage: [8]Limb = undefined;
516 const Ap = FixedIntFromSignedDoubleLimb(A, storage[0..2]);
517 const Bp = FixedIntFromSignedDoubleLimb(B, storage[2..4]);
518 const Cp = FixedIntFromSignedDoubleLimb(C, storage[4..6]);
519 const Dp = FixedIntFromSignedDoubleLimb(D, storage[6..8]);
520
521 // T = Ax + By
522 try r.mul(x, Ap);
523 try T.mul(y, Bp);
524 try T.add(r.*, T);
525
526 // u = Cx + Dy, r as u
527 try x.mul(x, Cp);
528 try r.mul(y, Dp);
529 try r.add(x, r.*);
530
531 x.swap(&T);
532 y.swap(r);
533 }
534 }
535
536 // euclidean algorithm
537 debug.assert(x.cmp(y) >= 0);
538
539 while (!y.eqZero()) {
540 try Int.divTrunc(&T, r, x, y);
541 x.swap(&y);
542 y.swap(r);
543 }
544
545 r.swap(&x);
546}
547
548var buffer: [64 * 8192]u8 = undefined;
549var fixed = std.heap.FixedBufferAllocator.init(buffer[0..]);
550var al = &fixed.allocator;
551
552test "big.rational gcd non-one small" {
553 var a = try Int.initSet(al, 17);
554 var b = try Int.initSet(al, 97);
555 var r = try Int.init(al);
556
557 try gcd(&r, a, b);
558
559 testing.expect((try r.to(u32)) == 1);
560}
561
562test "big.rational gcd non-one small" {
563 var a = try Int.initSet(al, 4864);
564 var b = try Int.initSet(al, 3458);
565 var r = try Int.init(al);
566
567 try gcd(&r, a, b);
568
569 testing.expect((try r.to(u32)) == 38);
570}
571
572test "big.rational gcd non-one large" {
573 var a = try Int.initSet(al, 0xffffffffffffffff);
574 var b = try Int.initSet(al, 0xffffffffffffffff7777);
575 var r = try Int.init(al);
576
577 try gcd(&r, a, b);
578
579 testing.expect((try r.to(u32)) == 4369);
580}
581
582test "big.rational gcd large multi-limb result" {
583 var a = try Int.initSet(al, 0x12345678123456781234567812345678123456781234567812345678);
584 var b = try Int.initSet(al, 0x12345671234567123456712345671234567123456712345671234567);
585 var r = try Int.init(al);
586
587 try gcd(&r, a, b);
588
589 testing.expect((try r.to(u256)) == 0xf000000ff00000fff0000ffff000fffff00ffffff1);
590}
591
592test "big.rational gcd one large" {
593 var a = try Int.initSet(al, 1897056385327307);
594 var b = try Int.initSet(al, 2251799813685248);
595 var r = try Int.init(al);
596
597 try gcd(&r, a, b);
598
599 testing.expect((try r.to(u64)) == 1);
600}
601
602fn extractLowBits(a: Int, comptime T: type) T {
603 testing.expect(@typeId(T) == builtin.TypeId.Int);
604
605 if (T.bit_count <= Limb.bit_count) {
606 return @truncate(T, a.limbs[0]);
607 } else {
608 var r: T = 0;
609 comptime var i: usize = 0;
610
611 // Remainder is always 0 since if T.bit_count >= Limb.bit_count -> Limb | T and both
612 // are powers of two.
613 inline while (i < T.bit_count / Limb.bit_count) : (i += 1) {
614 r |= math.shl(T, a.limbs[i], i * Limb.bit_count);
615 }
616
617 return r;
618 }
619}
620
621test "big.rational extractLowBits" {
622 var a = try Int.initSet(al, 0x11112222333344441234567887654321);
623
624 const a1 = extractLowBits(a, u8);
625 testing.expect(a1 == 0x21);
626
627 const a2 = extractLowBits(a, u16);
628 testing.expect(a2 == 0x4321);
629
630 const a3 = extractLowBits(a, u32);
631 testing.expect(a3 == 0x87654321);
632
633 const a4 = extractLowBits(a, u64);
634 testing.expect(a4 == 0x1234567887654321);
635
636 const a5 = extractLowBits(a, u128);
637 testing.expect(a5 == 0x11112222333344441234567887654321);
638}
639
640test "big.rational set" {
641 var a = try Rational.init(al);
642
643 try a.setInt(5);
644 testing.expect((try a.p.to(u32)) == 5);
645 testing.expect((try a.q.to(u32)) == 1);
646
647 try a.setRatio(7, 3);
648 testing.expect((try a.p.to(u32)) == 7);
649 testing.expect((try a.q.to(u32)) == 3);
650
651 try a.setRatio(9, 3);
652 testing.expect((try a.p.to(i32)) == 3);
653 testing.expect((try a.q.to(i32)) == 1);
654
655 try a.setRatio(-9, 3);
656 testing.expect((try a.p.to(i32)) == -3);
657 testing.expect((try a.q.to(i32)) == 1);
658
659 try a.setRatio(9, -3);
660 testing.expect((try a.p.to(i32)) == -3);
661 testing.expect((try a.q.to(i32)) == 1);
662
663 try a.setRatio(-9, -3);
664 testing.expect((try a.p.to(i32)) == 3);
665 testing.expect((try a.q.to(i32)) == 1);
666}
667
668test "big.rational setFloat" {
669 var a = try Rational.init(al);
670
671 try a.setFloat(f64, 2.5);
672 testing.expect((try a.p.to(i32)) == 5);
673 testing.expect((try a.q.to(i32)) == 2);
674
675 try a.setFloat(f32, -2.5);
676 testing.expect((try a.p.to(i32)) == -5);
677 testing.expect((try a.q.to(i32)) == 2);
678
679 try a.setFloat(f32, 3.141593);
680
681 // = 3.14159297943115234375
682 testing.expect((try a.p.to(u32)) == 3294199);
683 testing.expect((try a.q.to(u32)) == 1048576);
684
685 try a.setFloat(f64, 72.141593120712409172417410926841290461290467124);
686
687 // = 72.1415931207124145885245525278151035308837890625
688 testing.expect((try a.p.to(u128)) == 5076513310880537);
689 testing.expect((try a.q.to(u128)) == 70368744177664);
690}
691
692test "big.rational setFloatString" {
693 var a = try Rational.init(al);
694
695 try a.setFloatString("72.14159312071241458852455252781510353");
696
697 // = 72.1415931207124145885245525278151035308837890625
698 testing.expect((try a.p.to(u128)) == 7214159312071241458852455252781510353);
699 testing.expect((try a.q.to(u128)) == 100000000000000000000000000000000000);
700}
701
702test "big.rational toFloat" {
703 var a = try Rational.init(al);
704
705 // = 3.14159297943115234375
706 try a.setRatio(3294199, 1048576);
707 testing.expect((try a.toFloat(f64)) == 3.14159297943115234375);
708
709 // = 72.1415931207124145885245525278151035308837890625
710 try a.setRatio(5076513310880537, 70368744177664);
711 testing.expect((try a.toFloat(f64)) == 72.141593120712409172417410926841290461290467124);
712}
713
714test "big.rational set/to Float round-trip" {
715 var a = try Rational.init(al);
716 var prng = std.rand.DefaultPrng.init(0x5EED);
717 var i: usize = 0;
718 while (i < 512) : (i += 1) {
719 const r = prng.random.float(f64);
720 try a.setFloat(f64, r);
721 testing.expect((try a.toFloat(f64)) == r);
722 }
723}
724
725test "big.rational copy" {
726 var a = try Rational.init(al);
727
728 const b = try Int.initSet(al, 5);
729
730 try a.copyInt(b);
731 testing.expect((try a.p.to(u32)) == 5);
732 testing.expect((try a.q.to(u32)) == 1);
733
734 const c = try Int.initSet(al, 7);
735 const d = try Int.initSet(al, 3);
736
737 try a.copyRatio(c, d);
738 testing.expect((try a.p.to(u32)) == 7);
739 testing.expect((try a.q.to(u32)) == 3);
740
741 const e = try Int.initSet(al, 9);
742 const f = try Int.initSet(al, 3);
743
744 try a.copyRatio(e, f);
745 testing.expect((try a.p.to(u32)) == 3);
746 testing.expect((try a.q.to(u32)) == 1);
747}
748
749test "big.rational negate" {
750 var a = try Rational.init(al);
751
752 try a.setInt(-50);
753 testing.expect((try a.p.to(i32)) == -50);
754 testing.expect((try a.q.to(i32)) == 1);
755
756 a.negate();
757 testing.expect((try a.p.to(i32)) == 50);
758 testing.expect((try a.q.to(i32)) == 1);
759
760 a.negate();
761 testing.expect((try a.p.to(i32)) == -50);
762 testing.expect((try a.q.to(i32)) == 1);
763}
764
765test "big.rational abs" {
766 var a = try Rational.init(al);
767
768 try a.setInt(-50);
769 testing.expect((try a.p.to(i32)) == -50);
770 testing.expect((try a.q.to(i32)) == 1);
771
772 a.abs();
773 testing.expect((try a.p.to(i32)) == 50);
774 testing.expect((try a.q.to(i32)) == 1);
775
776 a.abs();
777 testing.expect((try a.p.to(i32)) == 50);
778 testing.expect((try a.q.to(i32)) == 1);
779}
780
781test "big.rational swap" {
782 var a = try Rational.init(al);
783 var b = try Rational.init(al);
784
785 try a.setRatio(50, 23);
786 try b.setRatio(17, 3);
787
788 testing.expect((try a.p.to(u32)) == 50);
789 testing.expect((try a.q.to(u32)) == 23);
790
791 testing.expect((try b.p.to(u32)) == 17);
792 testing.expect((try b.q.to(u32)) == 3);
793
794 a.swap(&b);
795
796 testing.expect((try a.p.to(u32)) == 17);
797 testing.expect((try a.q.to(u32)) == 3);
798
799 testing.expect((try b.p.to(u32)) == 50);
800 testing.expect((try b.q.to(u32)) == 23);
801}
802
803test "big.rational cmp" {
804 var a = try Rational.init(al);
805 var b = try Rational.init(al);
806
807 try a.setRatio(500, 231);
808 try b.setRatio(18903, 8584);
809 testing.expect((try a.cmp(b)) < 0);
810
811 try a.setRatio(890, 10);
812 try b.setRatio(89, 1);
813 testing.expect((try a.cmp(b)) == 0);
814}
815
816test "big.rational add single-limb" {
817 var a = try Rational.init(al);
818 var b = try Rational.init(al);
819
820 try a.setRatio(500, 231);
821 try b.setRatio(18903, 8584);
822 testing.expect((try a.cmp(b)) < 0);
823
824 try a.setRatio(890, 10);
825 try b.setRatio(89, 1);
826 testing.expect((try a.cmp(b)) == 0);
827}
828
829test "big.rational add" {
830 var a = try Rational.init(al);
831 var b = try Rational.init(al);
832 var r = try Rational.init(al);
833
834 try a.setRatio(78923, 23341);
835 try b.setRatio(123097, 12441414);
836 try a.add(a, b);
837
838 try r.setRatio(984786924199, 290395044174);
839 testing.expect((try a.cmp(r)) == 0);
840}
841
842test "big.rational sub" {
843 var a = try Rational.init(al);
844 var b = try Rational.init(al);
845 var r = try Rational.init(al);
846
847 try a.setRatio(78923, 23341);
848 try b.setRatio(123097, 12441414);
849 try a.sub(a, b);
850
851 try r.setRatio(979040510045, 290395044174);
852 testing.expect((try a.cmp(r)) == 0);
853}
854
855test "big.rational mul" {
856 var a = try Rational.init(al);
857 var b = try Rational.init(al);
858 var r = try Rational.init(al);
859
860 try a.setRatio(78923, 23341);
861 try b.setRatio(123097, 12441414);
862 try a.mul(a, b);
863
864 try r.setRatio(571481443, 17082061422);
865 testing.expect((try a.cmp(r)) == 0);
866}
867
868test "big.rational div" {
869 var a = try Rational.init(al);
870 var b = try Rational.init(al);
871 var r = try Rational.init(al);
872
873 try a.setRatio(78923, 23341);
874 try b.setRatio(123097, 12441414);
875 try a.div(a, b);
876
877 try r.setRatio(75531824394, 221015929);
878 testing.expect((try a.cmp(r)) == 0);
879}
880
881test "big.rational div" {
882 var a = try Rational.init(al);
883 var r = try Rational.init(al);
884
885 try a.setRatio(78923, 23341);
886 a.invert();
887
888 try r.setRatio(23341, 78923);
889 testing.expect((try a.cmp(r)) == 0);
890
891 try a.setRatio(-78923, 23341);
892 a.invert();
893
894 try r.setRatio(-23341, 78923);
895 testing.expect((try a.cmp(r)) == 0);
896}