authorgravatar for robin@voetter.nlRobin Voetter <robin@voetter.nl> 2021-10-03 16:17:46+02:00
committergravatar for robin@voetter.nlRobin Voetter <robin@voetter.nl> 2021-10-04 11:25:29+02:00
logfdf13fb81940ada1673914cab4ff46a1ddae42dd
tree4f7e11fbf3416c3ca139be4e6c29c9a9c1e05ef2
parent41e9c1bac1c447fe42a191bf16ee25ddb3bba97a

big ints: Wrapping multiplication


1 files changed, 90 insertions(+), 24 deletions(-)

lib/std/math/big/int.zig+90-24
...@@ -44,6 +44,11 @@ pub fn calcMulLimbsBufferLen(a_len: usize, b_len: usize, aliases: usize) usize {...@@ -44,6 +44,11 @@ pub fn calcMulLimbsBufferLen(a_len: usize, b_len: usize, aliases: usize) usize {
44 return aliases * math.max(a_len, b_len);44 return aliases * math.max(a_len, b_len);
45}45}
4646
47pub fn calcMulWrapLimbsBufferLen(bit_count: usize, a_len: usize, b_len: usize, aliases: usize) usize {
48 const req_limbs = calcTwosCompLimbCount(bit_count);
49 return aliases * math.min(req_limbs, math.max(a_len, b_len));
50}
51
47pub fn calcSetStringLimbsBufferLen(base: u8, string_len: usize) usize {52pub fn calcSetStringLimbsBufferLen(base: u8, string_len: usize) usize {
48 const limb_count = calcSetStringLimbCount(base, string_len);53 const limb_count = calcSetStringLimbCount(base, string_len);
49 return calcMulLimbsBufferLen(limb_count, limb_count, 2);54 return calcMulLimbsBufferLen(limb_count, limb_count, 2);
...@@ -611,7 +616,7 @@ pub const Mutable = struct {...@@ -611,7 +616,7 @@ pub const Mutable = struct {
611 /// `a` and `b` may alias with each other.616 /// `a` and `b` may alias with each other.
612 ///617 ///
613 /// Asserts the result fits in `rma`. An upper bound on the number of limbs needed by618 /// Asserts the result fits in `rma`. An upper bound on the number of limbs needed by
614 /// rma is given by `a.limbs.len + b.limbs.len + 1`.619 /// rma is given by `a.limbs.len + b.limbs.len`.
615 ///620 ///
616 /// `limbs_buffer` is used for temporary storage. The amount required is given by `calcMulLimbsBufferLen`.621 /// `limbs_buffer` is used for temporary storage. The amount required is given by `calcMulLimbsBufferLen`.
617 pub fn mul(rma: *Mutable, a: Const, b: Const, limbs_buffer: []Limb, allocator: ?*Allocator) void {622 pub fn mul(rma: *Mutable, a: Const, b: Const, limbs_buffer: []Limb, allocator: ?*Allocator) void {
...@@ -640,7 +645,7 @@ pub const Mutable = struct {...@@ -640,7 +645,7 @@ pub const Mutable = struct {
640 /// `a` and `b` may alias with each other.645 /// `a` and `b` may alias with each other.
641 ///646 ///
642 /// Asserts the result fits in `rma`. An upper bound on the number of limbs needed by647 /// Asserts the result fits in `rma`. An upper bound on the number of limbs needed by
643 /// rma is given by `a.limbs.len + b.limbs.len + 1`.648 /// rma is given by `a.limbs.len + b.limbs.len`.
644 ///649 ///
645 /// If `allocator` is provided, it will be used for temporary storage to improve650 /// If `allocator` is provided, it will be used for temporary storage to improve
646 /// multiplication performance. `error.OutOfMemory` is handled with a fallback algorithm.651 /// multiplication performance. `error.OutOfMemory` is handled with a fallback algorithm.
...@@ -658,18 +663,69 @@ pub const Mutable = struct {...@@ -658,18 +663,69 @@ pub const Mutable = struct {
658663
659 mem.set(Limb, rma.limbs[0 .. a.limbs.len + b.limbs.len + 1], 0);664 mem.set(Limb, rma.limbs[0 .. a.limbs.len + b.limbs.len + 1], 0);
660665
661 _ = llmulacc(.add, allocator, rma.limbs, a.limbs, b.limbs);666 llmulacc(.add, allocator, rma.limbs, a.limbs, b.limbs);
662667
663 rma.normalize(a.limbs.len + b.limbs.len);668 rma.normalize(a.limbs.len + b.limbs.len);
664 rma.positive = (a.positive == b.positive);669 rma.positive = (a.positive == b.positive);
665 }670 }
666671
667 pub fn mulNoAliasWrap(672 /// rma = a * b with 2s-complement wrapping semantics.
673 ///
674 /// `rma` may alias with `a` or `b`.
675 /// `a` and `b` may alias with each other.
676 ///
677 /// Asserts the result fits in `rma`. An upper bound on the number of limbs needed by
678 /// rma is given by `a.limbs.len + b.limbs.len`.
679 ///
680 /// `limbs_buffer` is used for temporary storage. The amount required is given by `calcMulWrapLimbsBufferLen`.
681 pub fn mulWrap(
682 rma: *Mutable,
683 a: Const,
684 b: Const,
685 signedness: std.builtin.Signedness,
686 bit_count: usize,
687 limbs_buffer: []Limb,
688 allocator: ?*Allocator,
689 ) void {
690 var buf_index: usize = 0;
691 const req_limbs = calcTwosCompLimbCount(bit_count);
692
693 const a_copy = if (rma.limbs.ptr == a.limbs.ptr) blk: {
694 const start = buf_index;
695 const a_len = math.min(req_limbs, a.limbs.len);
696 mem.copy(Limb, limbs_buffer[buf_index..], a.limbs[0..a_len]);
697 buf_index += a_len;
698 break :blk a.toMutable(limbs_buffer[start..buf_index]).toConst();
699 } else a;
700
701 const b_copy = if (rma.limbs.ptr == b.limbs.ptr) blk: {
702 const start = buf_index;
703 const b_len = math.min(req_limbs, b.limbs.len);
704 mem.copy(Limb, limbs_buffer[buf_index..], b.limbs[0..b_len]);
705 buf_index += b_len;
706 break :blk a.toMutable(limbs_buffer[start..buf_index]).toConst();
707 } else b;
708
709 return rma.mulWrapNoAlias(a_copy, b_copy, signedness, bit_count, allocator);
710 }
711
712 /// rma = a * b with 2s-complement wrapping semantics.
713 ///
714 /// `rma` may not alias with `a` or `b`.
715 /// `a` and `b` may alias with each other.
716 ///
717 /// Asserts the result fits in `rma`. An upper bound on the number of limbs needed by
718 /// rma is given by `a.limbs.len + b.limbs.len`.
719 ///
720 /// If `allocator` is provided, it will be used for temporary storage to improve
721 /// multiplication performance. `error.OutOfMemory` is handled with a fallback algorithm.
722 pub fn mulWrapNoAlias(
668 rma: *Mutable,723 rma: *Mutable,
669 a: Const,724 a: Const,
670 b: Const,725 b: Const,
671 signedness: std.builtin.Signedness,726 signedness: std.builtin.Signedness,
672 bit_count: usize,727 bit_count: usize,
728 allocator: ?*Allocator,
673 ) void {729 ) void {
674 assert(rma.limbs.ptr != a.limbs.ptr); // illegal aliasing730 assert(rma.limbs.ptr != a.limbs.ptr); // illegal aliasing
675 assert(rma.limbs.ptr != b.limbs.ptr); // illegal aliasing731 assert(rma.limbs.ptr != b.limbs.ptr); // illegal aliasing
...@@ -682,15 +738,9 @@ pub const Mutable = struct {...@@ -682,15 +738,9 @@ pub const Mutable = struct {
682738
683 mem.set(Limb, rma.limbs[0..req_limbs], 0);739 mem.set(Limb, rma.limbs[0..req_limbs], 0);
684740
685 if (a_limbs.len >= b_limbs.len) {741 llmulacc(.add, allocator, rma.limbs, a_limbs, b_limbs);
686 llmulaccLow(rma.limbs, a_limbs, b_limbs);
687 } else {
688 llmulaccLow(rma.limbs, b_limbs, a_limbs);
689 }
690
691 rma.normalize(math.min(req_limbs, a.limbs.len + b.limbs.len));742 rma.normalize(math.min(req_limbs, a.limbs.len + b.limbs.len));
692 rma.positive = (a.positive == b.positive);743 rma.positive = (a.positive == b.positive);
693
694 rma.truncate(rma.toConst(), signedness, bit_count);744 rma.truncate(rma.toConst(), signedness, bit_count);
695 }745 }
696746
...@@ -2167,6 +2217,35 @@ pub const Managed = struct {...@@ -2167,6 +2217,35 @@ pub const Managed = struct {
2167 rma.setMetadata(m.positive, m.len);2217 rma.setMetadata(m.positive, m.len);
2168 }2218 }
21692219
2220 /// rma = a * b with 2s-complement wrapping semantics.
2221 ///
2222 /// rma, a and b may be aliases. However, it is more efficient if rma does not alias a or b.
2223 /// If rma aliases a or b, then caller must call `rma.ensureCapacity(calcTwosCompLimbCount(bit_count))`
2224 /// prior to calling `mul`.
2225 ///
2226 /// Returns an error if memory could not be allocated.
2227 ///
2228 /// rma's allocator is used for temporary storage to speed up the multiplication.
2229 pub fn mulWrap(rma: *Managed, a: Const, b: Const, signedness: std.builtin.Signedness, bit_count: usize) !void {
2230 var alias_count: usize = 0;
2231 if (rma.limbs.ptr == a.limbs.ptr)
2232 alias_count += 1;
2233 if (rma.limbs.ptr == b.limbs.ptr)
2234 alias_count += 1;
2235
2236 try rma.ensureCapacity(calcTwosCompLimbCount(bit_count));
2237 var m = rma.toMutable();
2238 if (alias_count == 0) {
2239 m.mulWrapNoAlias(a, b, signedness, bit_count, rma.allocator);
2240 } else {
2241 const limb_count = calcMulWrapLimbsBufferLen(bit_count, a.limbs.len, b.limbs.len, alias_count);
2242 const limbs_buffer = try rma.allocator.alloc(Limb, limb_count);
2243 defer rma.allocator.free(limbs_buffer);
2244 m.mulWrap(a, b, signedness, bit_count, limbs_buffer, rma.allocator);
2245 }
2246 rma.setMetadata(m.positive, m.len);
2247 }
2248
2170 pub fn ensureAddScalarCapacity(r: *Managed, a: Const, scalar: anytype) !void {2249 pub fn ensureAddScalarCapacity(r: *Managed, a: Const, scalar: anytype) !void {
2171 try r.ensureCapacity(math.max(a.limbs.len, calcLimbLen(scalar)) + 1);2250 try r.ensureCapacity(math.max(a.limbs.len, calcLimbLen(scalar)) + 1);
2172 }2251 }
...@@ -2337,19 +2416,6 @@ pub const Managed = struct {...@@ -2337,19 +2416,6 @@ pub const Managed = struct {
2337 }2416 }
2338};2417};
23392418
2340/// r = r + a * b, ignoring overflow
2341fn llmulaccLow(r: []Limb, a: []const Limb, b: []const Limb) void {
2342 assert(r.len >= a.len);
2343 assert(a.len >= b.len);
2344
2345 // TODO: Improve performance.
2346
2347 var i: usize = 0;
2348 while (i < b.len) : (i += 1) {
2349 llmulLimb(.add, r[i..], a, b[i]);
2350 }
2351}
2352
2353/// Different operators which can be used in accumulation style functions2419/// Different operators which can be used in accumulation style functions
2354/// (llmulacc, llmulaccKaratsuba, llmulaccLong, llmulLimb). In all these functions,2420/// (llmulacc, llmulaccKaratsuba, llmulaccLong, llmulLimb). In all these functions,
2355/// a computed value is accumulated with an existing result.2421/// a computed value is accumulated with an existing result.