Meg's Mirrors

zig

1	const std = @import("../../std.zig");
2	const builtin = @import("builtin");
3	const math = std.math;
4	const Limb = std.math.big.Limb;
5	const limb_bits = @typeInfo(Limb).int.bits;
6	const HalfLimb = std.math.big.HalfLimb;
7	const half_limb_bits = @typeInfo(HalfLimb).int.bits;
8	const DoubleLimb = std.math.big.DoubleLimb;
9	const SignedDoubleLimb = std.math.big.SignedDoubleLimb;
10	const Log2Limb = std.math.big.Log2Limb;
11	const Allocator = std.mem.Allocator;
12	const mem = std.mem;
13	const maxInt = std.math.maxInt;
14	const minInt = std.math.minInt;
15	const assert = std.debug.assert;
16	const Endian = std.builtin.Endian;
17	const Signedness = std.builtin.Signedness;
18	const native_endian = builtin.cpu.arch.endian();
19
20	// Comptime-computed constants for supported bases (2 - 36)
21	// all values are set to 0 for bases 0 - 1, to make it possible to
22	// access a constant for a given base b using `constants.value[b]`
23	const Constants = struct {
24	// big_bases[b] is the biggest power of b that fit in a single Limb
25	// i.e. big_bases[b] = b^k < 2^@bitSizeOf(Limb) and b^(k+1) >= 2^@bitSizeOf(Limb)
26	big_bases: [37]Limb,
27	// digits_per_limb[b] is the value of k used in the previous field
28	digits_per_limb: [37]u8,
29	};
30	const constants: Constants = blk: {
31	@setEvalBranchQuota(2000);
32	var digits_per_limb: [37]u8 = @splat(0);
33	var bases: [37]Limb = @splat(0);
34	for (2..37) |base| {
35	digits_per_limb[base] = @intCast(math.log(Limb, base, math.maxInt(Limb)));
36	bases[base] = std.math.pow(Limb, base, digits_per_limb[base]);
37	}
38	break :blk Constants{ .big_bases = bases, .digits_per_limb = digits_per_limb };
39	};
40
41	/// Returns the number of limbs needed to store `scalar`, which must be a
42	/// primitive integer or float value.
43	/// Note: A comptime-known upper bound of this value that may be used
44	/// instead if `scalar` is not already comptime-known is
45	/// `calcTwosCompLimbCount(@typeInfo(@TypeOf(scalar)).int.bits)`
46	pub fn calcLimbLen(scalar: anytype) usize {
47	switch (@typeInfo(@TypeOf(scalar))) {
48	.int, .comptime_int => {
49	if (scalar == 0) return 1;
50	const w_value = @abs(scalar);
51	return @as(usize, @intCast(@divFloor(@as(Limb, @intCast(math.log2(w_value))), limb_bits) + 1));
52	},
53	.float => {
54	const repr: std.math.FloatRepr(@TypeOf(scalar)) = @bitCast(scalar);
55	return switch (repr.exponent) {
56	.denormal => 1,
57	else => return calcNonZeroTwosCompLimbCount(@as(usize, 2) + @max(repr.exponent.unbias(), 0)),
58	.infinite => 0,
59	};
60	},
61	.comptime_float => return calcLimbLen(@as(f128, scalar)),
62	else => @compileError("expected float or int, got " ++ @typeName(@TypeOf(scalar))),
63	}
64	}
65
66	/// Same as `calcToStringLimbsBufferLen`, without the useless base check.
67	pub fn calcLog10LimbsBufferLen(a_len: usize) usize {
68	return a_len + 2 + a_len + calcDivLimbsBufferLen(a_len, 1);
69	}
70
71	pub fn calcToStringLimbsBufferLen(a_len: usize, base: u8) usize {
72	if (math.isPowerOfTwo(base))
73	return 0;
74	return a_len + 2 + a_len + calcDivLimbsBufferLen(a_len, 1);
75	}
76
77	pub fn calcDivLimbsBufferLen(a_len: usize, b_len: usize) usize {
78	return a_len + b_len + 4;
79	}
80
81	pub fn calcMulLimbsBufferLen(a_len: usize, b_len: usize, aliases: usize) usize {
82	return aliases * @max(a_len, b_len);
83	}
84
85	pub fn calcMulWrapLimbsBufferLen(bit_count: usize, a_len: usize, b_len: usize, aliases: usize) usize {
86	const req_limbs = calcTwosCompLimbCount(bit_count);
87	return aliases * @min(req_limbs, @max(a_len, b_len));
88	}
89
90	pub fn calcSetStringLimbsBufferLen(base: u8, string_len: usize) usize {
91	const limb_count = calcSetStringLimbCount(base, string_len);
92	return calcMulLimbsBufferLen(limb_count, limb_count, 2);
93	}
94
95	/// Assumes `string_len` doesn't account for minus signs if the number is negative.
96	pub fn calcSetStringLimbCount(base: u8, string_len: usize) usize {
97	const base_f: f32 = @floatFromInt(base);
98	const string_len_f: f32 = @floatFromInt(string_len);
99	return 1 + @as(usize, @intFromFloat(@ceil(string_len_f * std.math.log2(base_f) / limb_bits)));
100	}
101
102	pub fn calcPowLimbsBufferLen(a_bit_count: usize, y: usize) usize {
103	// The 2 accounts for the minimum space requirement for llmulacc
104	return 2 + (a_bit_count * y + (limb_bits - 1)) / limb_bits;
105	}
106
107	pub fn calcSqrtLimbsBufferLen(a_bit_count: usize) usize {
108	const a_limb_count = (a_bit_count - 1) / limb_bits + 1;
109	const shift = (a_bit_count + 1) / 2;
110	const u_s_rem_limb_count = 1 + ((shift / limb_bits) + 1);
111	return a_limb_count + 3 * u_s_rem_limb_count + calcDivLimbsBufferLen(a_limb_count, u_s_rem_limb_count);
112	}
113
114	/// Compute the number of limbs required to store a 2s-complement number of `bit_count` bits.
115	pub fn calcNonZeroTwosCompLimbCount(bit_count: usize) usize {
116	assert(bit_count != 0);
117	return calcTwosCompLimbCount(bit_count);
118	}
119
120	/// Compute the number of limbs required to store a 2s-complement number of `bit_count` bits.
121	///
122	/// Special cases `bit_count == 0` to return 1. Zero-bit integers can only store the value zero
123	/// and this big integer implementation stores zero using one limb.
124	pub fn calcTwosCompLimbCount(bit_count: usize) usize {
125	return @max(@divCeil(bit_count, @bitSizeOf(Limb)), 1);
126	}
127
128	/// a + b * c + *carry, sets carry to the overflow bits
129	pub fn addMulLimbWithCarry(a: Limb, b: Limb, c: Limb, carry: *Limb) Limb {
130	// ov1[0] = a + *carry
131	const ov1 = @addWithOverflow(a, carry.*);
132
133	// r2 = b * c
134	const bc = @as(DoubleLimb, math.mulWide(Limb, b, c));
135	const r2 = @as(Limb, @truncate(bc));
136	const c2 = @as(Limb, @truncate(bc >> limb_bits));
137
138	// ov2[0] = ov1[0] + r2
139	const ov2 = @addWithOverflow(ov1[0], r2);
140
141	// This never overflows, c1, c3 are either 0 or 1 and if both are 1 then
142	// c2 is at least <= maxInt(Limb) - 2.
143	carry.* = ov1[1] + c2 + ov2[1];
144
145	return ov2[0];
146	}
147
148	/// a - b * c - *carry, sets carry to the overflow bits
149	fn subMulLimbWithBorrow(a: Limb, b: Limb, c: Limb, carry: *Limb) Limb {
150	// ov1[0] = a - *carry
151	const ov1 = @subWithOverflow(a, carry.*);
152
153	// r2 = b * c
154	const bc = @as(DoubleLimb, std.math.mulWide(Limb, b, c));
155	const r2 = @as(Limb, @truncate(bc));
156	const c2 = @as(Limb, @truncate(bc >> limb_bits));
157
158	// ov2[0] = ov1[0] - r2
159	const ov2 = @subWithOverflow(ov1[0], r2);
160	carry.* = ov1[1] + c2 + ov2[1];
161
162	return ov2[0];
163	}
164
165	/// Used to indicate either limit of a 2s-complement integer.
166	pub const TwosCompIntLimit = enum {
167	// The low limit, either 0x00 (unsigned) or (-)0x80 (signed) for an 8-bit integer.
168	min,
169
170	// The high limit, either 0xFF (unsigned) or 0x7F (signed) for an 8-bit integer.
171	max,
172	};
173
174	pub const Round = enum {
175	/// Round to the nearest representable value, with ties broken by the representation
176	/// that ends with a 0 bit.
177	nearest_even,
178	/// Round away from zero.
179	away,
180	/// Round towards zero.
181	trunc,
182	/// Round towards negative infinity.
183	floor,
184	/// Round towards positive infinity.
185	ceil,
186	};
187
188	pub const Exactness = enum { inexact, exact };
189
190	/// A arbitrary-precision big integer, with a fixed set of mutable limbs.
191	pub const Mutable = struct {
192	/// Raw digits. These are:
193	///
194	/// * Little-endian ordered
195	/// * limbs.len >= 1
196	/// * Zero is represented as limbs.len == 1 with limbs[0] == 0.
197	///
198	/// Accessing limbs directly should be avoided.
199	/// These are allocated limbs; the `len` field tells the valid range.
200	limbs: []Limb,
201	len: usize,
202	positive: bool,
203
204	pub fn toConst(self: Mutable) Const {
205	return .{
206	.limbs = self.limbs[0..self.len],
207	.positive = self.positive,
208	};
209	}
210
211	pub const ConvertError = Const.ConvertError;
212
213	/// Convert `self` to `Int`.
214	///
215	/// Returns an error if self cannot be narrowed into the requested type without truncation.
216	pub fn toInt(self: Mutable, comptime Int: type) ConvertError!Int {
217	return self.toConst().toInt(Int);
218	}
219
220	/// Convert `self` to `Float`.
221	pub fn toFloat(self: Mutable, comptime Float: type, round: Round) struct { Float, Exactness } {
222	return self.toConst().toFloat(Float, round);
223	}
224
225	/// Returns true if `a == 0`.
226	pub fn eqlZero(self: Mutable) bool {
227	return self.toConst().eqlZero();
228	}
229
230	/// Asserts that the allocator owns the limbs memory. If this is not the case,
231	/// use `toConst().toManaged()`.
232	pub fn toManaged(self: Mutable, allocator: Allocator) Managed {
233	return .{
234	.allocator = allocator,
235	.limbs = self.limbs,
236	.metadata = if (self.positive)
237	self.len & ~Managed.sign_bit
238	else
239	self.len | Managed.sign_bit,
240	};
241	}
242
243	/// `value` is a primitive integer type.
244	/// Asserts the value fits within the provided `limbs_buffer`.
245	/// Note: `calcLimbLen` can be used to figure out how big an array to allocate for `limbs_buffer`.
246	pub fn init(limbs_buffer: []Limb, value: anytype) Mutable {
247	limbs_buffer[0] = 0;
248	var self: Mutable = .{
249	.limbs = limbs_buffer,
250	.len = 1,
251	.positive = true,
252	};
253	self.set(value);
254	return self;
255	}
256
257	/// Copies the value of a Const to an existing Mutable so that they both have the same value.
258	/// Asserts the value fits in the limbs buffer.
259	pub fn copy(self: *Mutable, other: Const) void {
260	if (self.limbs.ptr != other.limbs.ptr) {
261	@memcpy(self.limbs[0..other.limbs.len], other.limbs[0..other.limbs.len]);
262	}
263	// Normalize before setting `positive` so the `eqlZero` doesn't need to iterate
264	// over the extra zero limbs.
265	self.normalize(other.limbs.len);
266	self.positive = other.positive or other.eqlZero();
267	}
268
269	/// Efficiently swap an Mutable with another. This swaps the limb pointers and a full copy is not
270	/// performed. The address of the limbs field will not be the same after this function.
271	pub fn swap(self: *Mutable, other: *Mutable) void {
272	mem.swap(Mutable, self, other);
273	}
274
275	pub fn dump(self: Mutable) void {
276	for (self.limbs[0..self.len]) |limb| {
277	std.debug.print("{x} ", .{limb});
278	}
279	std.debug.print("len={} capacity={} positive={}\n", .{ self.len, self.limbs.len, self.positive });
280	}
281
282	/// Clones an Mutable and returns a new Mutable with the same value. The new Mutable is a deep copy and
283	/// can be modified separately from the original.
284	/// Asserts that limbs is big enough to store the value.
285	pub fn clone(other: Mutable, limbs: []Limb) Mutable {
286	@memcpy(limbs[0..other.len], other.limbs[0..other.len]);
287	return .{
288	.limbs = limbs,
289	.len = other.len,
290	.positive = other.positive,
291	};
292	}
293
294	pub fn negate(self: *Mutable) void {
295	self.positive = !self.positive;
296	}
297
298	/// Modify to become the absolute value
299	pub fn abs(self: *Mutable) void {
300	self.positive = true;
301	}
302
303	/// Sets the Mutable to value. Value must be an primitive integer type.
304	/// Asserts the value fits within the limbs buffer.
305	/// Note: `calcLimbLen` can be used to figure out how big the limbs buffer
306	/// needs to be to store a specific value.
307	pub fn set(self: *Mutable, value: anytype) void {
308	const T = @TypeOf(value);
309	const needed_limbs = calcLimbLen(value);
310	assert(needed_limbs <= self.limbs.len); // value too big
311
312	self.len = needed_limbs;
313	self.positive = value >= 0;
314
315	switch (@typeInfo(T)) {
316	.int => |info| {
317	var w_value = @abs(value);
318
319	if (info.bits <= limb_bits) {
320	self.limbs[0] = w_value;
321	} else {
322	var i: usize = 0;
323	while (true) : (i += 1) {
324	self.limbs[i] = @as(Limb, @truncate(w_value));
325	w_value >>= limb_bits;
326
327	if (w_value == 0) break;
328	}
329	}
330	},
331	.comptime_int => {
332	comptime var w_value = @abs(value);
333
334	if (w_value <= maxInt(Limb)) {
335	self.limbs[0] = w_value;
336	} else {
337	const mask = (1 << limb_bits) - 1;
338
339	comptime var i = 0;
340	inline while (true) : (i += 1) {
341	self.limbs[i] = w_value & mask;
342	w_value >>= limb_bits;
343
344	if (w_value == 0) break;
345	}
346	}
347	},
348	else => @compileError("cannot set Mutable using type " ++ @typeName(T)),
349	}
350	}
351
352	/// Set self from the string representation `value`.
353	///
354	/// `value` must contain only digits <= `base` and is case insensitive. Base prefixes are
355	/// not allowed (e.g. 0x43 should simply be 43). Underscores in the input string are
356	/// ignored and can be used as digit separators.
357	///
358	/// There must be enough memory for the value in `self.limbs`. An upper bound on number of limbs can
359	/// be determined with `calcSetStringLimbCount`.
360	/// Asserts the base is in the range [2, 36].
361	///
362	/// Returns an error if the value has invalid digits for the requested base.
363	pub fn setString(
364	self: *Mutable,
365	base: u8,
366	value: []const u8,
367	) error{InvalidCharacter}!void {
368	assert(base >= 2);
369	assert(base <= 36);
370
371	var i: usize = 0;
372	var positive = true;
373	if (value.len > 0 and value[0] == '-') {
374	positive = false;
375	i += 1;
376	}
377
378	@memset(self.limbs, 0);
379	self.len = 1;
380
381	var limb: Limb = 0;
382	var j: usize = 0;
383	for (value[i..]) |ch| {
384	if (ch == '_') {
385	continue;
386	}
387	const d = try std.fmt.charToDigit(ch, base);
388	limb *= base;
389	limb += d;
390	j += 1;
391
392	if (j == constants.digits_per_limb[base]) {
393	const len = @min(self.len + 1, self.limbs.len);
394	// r = a * b = a + a * (b - 1)
395	// we assert when self.limbs is not large enough to store the number
396	assert(!llmulLimb(.add, self.limbs[0..len], self.limbs[0..len], constants.big_bases[base] - 1));
397	assert(lladdcarry(self.limbs[0..len], self.limbs[0..len], &[1]Limb{limb}) == 0);
398
399	if (self.limbs.len > self.len and self.limbs[self.len] != 0)
400	self.len += 1;
401	j = 0;
402	limb = 0;
403	}
404	}
405	if (j > 0) {
406	const len = @min(self.len + 1, self.limbs.len);
407	// we assert when self.limbs is not large enough to store the number
408	assert(!llmulLimb(.add, self.limbs[0..len], self.limbs[0..len], math.pow(Limb, base, j) - 1));
409	assert(lladdcarry(self.limbs[0..len], self.limbs[0..len], &[1]Limb{limb}) == 0);
410
411	if (self.limbs.len > self.len and self.limbs[self.len] != 0)
412	self.len += 1;
413	}
414	self.positive = positive;
415	}
416
417	/// Set self to either bound of a 2s-complement integer.
418	/// Note: The result is still sign-magnitude, not twos complement! In order to convert the
419	/// result to twos complement, it is sufficient to take the absolute value.
420	///
421	/// Asserts the result fits in `r`. An upper bound on the number of limbs needed by
422	/// r is `calcTwosCompLimbCount(bit_count)`.
423	pub fn setTwosCompIntLimit(
424	r: *Mutable,
425	limit: TwosCompIntLimit,
426	signedness: Signedness,
427	bit_count: usize,
428	) void {
429	// Handle zero-bit types.
430	if (bit_count == 0) {
431	r.set(0);
432	return;
433	}
434
435	const req_limbs = calcTwosCompLimbCount(bit_count);
436	const bit: Log2Limb = @truncate(bit_count - 1);
437	const signmask = @as(Limb, 1) << bit; // 0b0..010..0 where 1 is the sign bit.
438	const mask = (signmask << 1) -% 1; // 0b0..011..1 where the leftmost 1 is the sign bit.
439
440	r.positive = true;
441
442	switch (signedness) {
443	.signed => switch (limit) {
444	.min => {
445	// Negative bound, signed = -0x80.
446	r.len = req_limbs;
447	@memset(r.limbs[0 .. r.len - 1], 0);
448	r.limbs[r.len - 1] = signmask;
449	r.positive = false;
450	},
451	.max => {
452	// Positive bound, signed = 0x7F
453	// Note, in this branch we need to normalize because the first bit is
454	// supposed to be 0.
455
456	// Special case for 1-bit integers.
457	if (bit_count == 1) {
458	r.set(0);
459	} else {
460	const new_req_limbs = calcTwosCompLimbCount(bit_count - 1);
461	const msb = @as(Log2Limb, @truncate(bit_count - 2));
462	const new_signmask = @as(Limb, 1) << msb; // 0b0..010..0 where 1 is the sign bit.
463	const new_mask = (new_signmask << 1) -% 1; // 0b0..001..1 where the rightmost 0 is the sign bit.
464
465	r.len = new_req_limbs;
466	@memset(r.limbs[0 .. r.len - 1], maxInt(Limb));
467	r.limbs[r.len - 1] = new_mask;
468	}
469	},
470	},
471	.unsigned => switch (limit) {
472	.min => {
473	// Min bound, unsigned = 0x00
474	r.set(0);
475	},
476	.max => {
477	// Max bound, unsigned = 0xFF
478	r.len = req_limbs;
479	@memset(r.limbs[0 .. r.len - 1], maxInt(Limb));
480	r.limbs[r.len - 1] = mask;
481	},
482	},
483	}
484	}
485
486	/// Sets the Mutable to a float value rounded according to `round`.
487	/// Returns whether the conversion was exact (`round` had no effect on the result).
488	pub fn setFloat(self: *Mutable, value: anytype, round: Round) Exactness {
489	const Float = @TypeOf(value);
490	if (Float == comptime_float) return self.setFloat(@as(f128, value), round);
491	const abs_value = @abs(value);
492	if (abs_value < 1.0) {
493	if (abs_value == 0.0) {
494	self.set(0);
495	return .exact;
496	}
497	self.set(@as(i2, round: switch (round) {
498	.nearest_even => if (abs_value <= 0.5) 0 else continue :round .away,
499	.away => if (value < 0.0) -1 else 1,
500	.trunc => 0,
501	.floor => -@as(i2, @intFromBool(value < 0.0)),
502	.ceil => @intFromBool(value > 0.0),
503	}));
504	return .inexact;
505	}
506	const Repr = std.math.FloatRepr(Float);
507	const repr: Repr = @bitCast(value);
508	const exponent = repr.exponent.unbias();
509	assert(exponent >= 0);
510	const int_bit: Repr.Mantissa = 1 << (@bitSizeOf(Repr.Mantissa) - 1);
511	const mantissa = int_bit | repr.mantissa;
512	if (exponent >= @bitSizeOf(Repr.Normalized.Fraction)) {
513	self.set(mantissa);
514	self.shiftLeft(self.toConst(), @intCast(exponent - @bitSizeOf(Repr.Normalized.Fraction)));
515	self.positive = repr.sign == .positive;
516	return .exact;
517	}
518	self.set(mantissa >> @intCast(@bitSizeOf(Repr.Normalized.Fraction) - exponent));
519	const round_bits: Repr.Normalized.Fraction = @truncate(mantissa << @intCast(exponent));
520	if (round_bits == 0) {
521	self.positive = repr.sign == .positive;
522	return .exact;
523	}
524	round: switch (round) {
525	.nearest_even => {
526	const half: Repr.Normalized.Fraction = 1 << (@bitSizeOf(Repr.Normalized.Fraction) - 1);
527	if (round_bits >= half) self.addScalar(self.toConst(), 1);
528	if (round_bits == half) self.limbs[0] &= ~@as(Limb, 1);
529	},
530	.away => self.addScalar(self.toConst(), 1),
531	.trunc => {},
532	.floor => switch (repr.sign) {
533	.positive => {},
534	.negative => continue :round .away,
535	},
536	.ceil => switch (repr.sign) {
537	.positive => continue :round .away,
538	.negative => {},
539	},
540	}
541	self.positive = repr.sign == .positive;
542	return .inexact;
543	}
544
545	/// r = a + scalar
546	///
547	/// r and a may be aliases.
548	/// scalar is a primitive integer type.
549	///
550	/// Asserts the result fits in `r`. An upper bound on the number of limbs needed by
551	/// r is `@max(a.limbs.len, calcLimbLen(scalar)) + 1`.
552	pub fn addScalar(r: *Mutable, a: Const, scalar: anytype) void {
553	// Normally we could just determine the number of limbs needed with calcLimbLen,
554	// but that is not comptime-known when scalar is not a comptime_int. Instead, we
555	// use calcTwosCompLimbCount for a non-comptime_int scalar, which can be pessimistic
556	// in the case that scalar happens to be small in magnitude within its type, but it
557	// is well worth being able to use the stack and not needing an allocator passed in.
558	// Note that Mutable.init still sets len to calcLimbLen(scalar) in any case.
559	const limbs_len = comptime switch (@typeInfo(@TypeOf(scalar))) {
560	.comptime_int => calcLimbLen(scalar),
561	.int => |info| calcTwosCompLimbCount(info.bits),
562	else => @compileError("expected scalar to be an int"),
563	};
564	var limbs: [limbs_len]Limb = undefined;
565	const operand = init(&limbs, scalar).toConst();
566	return add(r, a, operand);
567	}
568
569	/// Base implementation for addition. Adds `@max(a.limbs.len, b.limbs.len)` elements from a and b,
570	/// and returns whether any overflow occurred.
571	/// r, a and b may be aliases.
572	///
573	/// Asserts r has enough elements to hold the result. The upper bound is `@max(a.limbs.len, b.limbs.len)`.
574	fn addCarry(r: *Mutable, a: Const, b: Const) bool {
575	if (a.eqlZero()) {
576	r.copy(b);
577	return false;
578	} else if (b.eqlZero()) {
579	r.copy(a);
580	return false;
581	} else if (a.positive != b.positive) {
582	if (a.positive) {
583	// (a) + (-b) => a - b
584	return r.subCarry(a, b.abs());
585	} else {
586	// (-a) + (b) => b - a
587	return r.subCarry(b, a.abs());
588	}
589	} else {
590	r.positive = a.positive;
591	if (a.limbs.len >= b.limbs.len) {
592	const c = lladdcarry(r.limbs, a.limbs, b.limbs);
593	r.normalize(a.limbs.len);
594	return c != 0;
595	} else {
596	const c = lladdcarry(r.limbs, b.limbs, a.limbs);
597	r.normalize(b.limbs.len);
598	return c != 0;
599	}
600	}
601	}
602
603	/// r = a + b
604	///
605	/// r, a and b may be aliases.
606	///
607	/// Asserts the result fits in `r`. An upper bound on the number of limbs needed by
608	/// r is `@max(a.limbs.len, b.limbs.len) + 1`.
609	pub fn add(r: *Mutable, a: Const, b: Const) void {
610	if (r.addCarry(a, b)) {
611	// Fix up the result. Note that addCarry normalizes by a.limbs.len or b.limbs.len,
612	// so we need to set the length here.
613	const msl = @max(a.limbs.len, b.limbs.len);
614	// `[add|sub]Carry` normalizes by `msl`, so we need to fix up the result manually here.
615	// Note, the fact that it normalized means that the intermediary limbs are zero here.
616	r.len = msl + 1;
617	r.limbs[msl] = 1; // If this panics, there wasn't enough space in `r`.
618	}
619	}
620
621	/// r = a + b with 2s-complement wrapping semantics. Returns whether overflow occurred.
622	/// r, a and b may be aliases
623	///
624	/// Asserts the result fits in `r`. An upper bound on the number of limbs needed by
625	/// r is `calcTwosCompLimbCount(bit_count)`.
626	pub fn addWrap(r: *Mutable, a: Const, b: Const, signedness: Signedness, bit_count: usize) bool {
627	const req_limbs = calcTwosCompLimbCount(bit_count);
628
629	// Slice of the upper bits if they exist, these will be ignored and allows us to use addCarry to determine
630	// if an overflow occurred.
631	const x: Const = .{
632	.positive = a.positive,
633	.limbs = a.limbs[0..@min(req_limbs, a.limbs.len)],
634	};
635
636	const y: Const = .{
637	.positive = b.positive,
638	.limbs = b.limbs[0..@min(req_limbs, b.limbs.len)],
639	};
640
641	var carry_truncated = false;
642	if (r.addCarry(x, y)) {
643	// There are two possibilities here:
644	// - We overflowed req_limbs. In this case, the carry is ignored, as it would be removed by
645	// truncate anyway.
646	// - a and b had less elements than req_limbs, and those were overflowed. This case needs to be handled.
647	// Note: after this we still might need to wrap.
648	const msl = @max(a.limbs.len, b.limbs.len);
649	if (msl < req_limbs) {
650	r.len = msl + 1;
651	r.limbs[msl] = 1;
652	} else {
653	carry_truncated = true;
654	}
655	}
656
657	if (!r.toConst().fitsInTwosComp(signedness, bit_count)) {
658	r.truncate(r.toConst(), signedness, bit_count);
659	return true;
660	}
661
662	return carry_truncated;
663	}
664
665	/// r = a + b with 2s-complement saturating semantics.
666	/// r, a and b may be aliases.
667	///
668	/// Assets the result fits in `r`. Upper bound on the number of limbs needed by
669	/// r is `calcTwosCompLimbCount(bit_count)`.
670	pub fn addSat(r: *Mutable, a: Const, b: Const, signedness: Signedness, bit_count: usize) void {
671	const req_limbs = calcTwosCompLimbCount(bit_count);
672
673	// Slice of the upper bits if they exist, these will be ignored and allows us to use addCarry to determine
674	// if an overflow occurred.
675	const x: Const = .{
676	.positive = a.positive,
677	.limbs = a.limbs[0..@min(req_limbs, a.limbs.len)],
678	};
679
680	const y: Const = .{
681	.positive = b.positive,
682	.limbs = b.limbs[0..@min(req_limbs, b.limbs.len)],
683	};
684
685	if (r.addCarry(x, y)) {
686	// There are two possibilities here:
687	// - We overflowed req_limbs, in which case we need to saturate.
688	// - a and b had less elements than req_limbs, and those were overflowed.
689	// Note: In this case, might _also_ need to saturate.
690	const msl = @max(a.limbs.len, b.limbs.len);
691	if (msl < req_limbs) {
692	r.len = msl + 1;
693	r.limbs[msl] = 1;
694	// Note: Saturation may still be required if msl == req_limbs - 1
695	} else {
696	// Overflowed req_limbs, definitely saturate.
697	return r.setTwosCompIntLimit(if (r.positive) .max else .min, signedness, bit_count);
698	}
699	}
700
701	// Saturate if the result didn't fit.
702	r.saturate(r.toConst(), signedness, bit_count);
703	}
704
705	/// Base implementation for subtraction. Subtracts `@max(a.limbs.len, b.limbs.len)` elements from a and b,
706	/// and returns whether any overflow occurred.
707	/// r, a and b may be aliases.
708	///
709	/// Asserts r has enough elements to hold the result. The upper bound is `@max(a.limbs.len, b.limbs.len)`.
710	fn subCarry(r: *Mutable, a: Const, b: Const) bool {
711	if (a.eqlZero()) {
712	r.copy(b);
713	r.positive = !b.positive;
714	return false;
715	} else if (b.eqlZero()) {
716	r.copy(a);
717	return false;
718	} else if (a.positive != b.positive) {
719	if (a.positive) {
720	// (a) - (-b) => a + b
721	return r.addCarry(a, b.abs());
722	} else {
723	// (-a) - (b) => -a + -b
724	return r.addCarry(a, b.negate());
725	}
726	} else if (a.positive) {
727	if (a.order(b) != .lt) {
728	// (a) - (b) => a - b
729	const c = llsubcarry(r.limbs, a.limbs, b.limbs);
730	r.normalize(a.limbs.len);
731	r.positive = true;
732	return c != 0;
733	} else {
734	// (a) - (b) => -b + a => -(b - a)
735	const c = llsubcarry(r.limbs, b.limbs, a.limbs);
736	r.normalize(b.limbs.len);
737	r.positive = false;
738	return c != 0;
739	}
740	} else {
741	if (a.order(b) == .lt) {
742	// (-a) - (-b) => -(a - b)
743	const c = llsubcarry(r.limbs, a.limbs, b.limbs);
744	r.normalize(a.limbs.len);
745	r.positive = false;
746	return c != 0;
747	} else {
748	// (-a) - (-b) => --b + -a => b - a
749	const c = llsubcarry(r.limbs, b.limbs, a.limbs);
750	r.normalize(b.limbs.len);
751	r.positive = true;
752	return c != 0;
753	}
754	}
755	}
756
757	/// r = a - b
758	///
759	/// r, a and b may be aliases.
760	///
761	/// Asserts the result fits in `r`. An upper bound on the number of limbs needed by
762	/// r is `@max(a.limbs.len, b.limbs.len) + 1`. The +1 is not needed if both operands are positive.
763	pub fn sub(r: *Mutable, a: Const, b: Const) void {
764	r.add(a, b.negate());
765	}
766
767	/// r = a - b with 2s-complement wrapping semantics. Returns whether any overflow occurred.
768	///
769	/// r, a and b may be aliases
770	/// Asserts the result fits in `r`. An upper bound on the number of limbs needed by
771	/// r is `calcTwosCompLimbCount(bit_count)`.
772	pub fn subWrap(r: *Mutable, a: Const, b: Const, signedness: Signedness, bit_count: usize) bool {
773	return r.addWrap(a, b.negate(), signedness, bit_count);
774	}
775
776	/// r = a - b with 2s-complement saturating semantics.
777	/// r, a and b may be aliases.
778	///
779	/// Assets the result fits in `r`. Upper bound on the number of limbs needed by
780	/// r is `calcTwosCompLimbCount(bit_count)`.
781	pub fn subSat(r: *Mutable, a: Const, b: Const, signedness: Signedness, bit_count: usize) void {
782	r.addSat(a, b.negate(), signedness, bit_count);
783	}
784
785	/// rma = a * b
786	///
787	/// `rma` may alias with `a` or `b`.
788	/// `a` and `b` may alias with each other.
789	///
790	/// Asserts the result fits in `rma`. An upper bound on the number of limbs needed by
791	/// rma is given by `a.limbs.len + b.limbs.len`.
792	///
793	/// `limbs_buffer` is used for temporary storage. The amount required is given by `calcMulLimbsBufferLen`.
794	pub fn mul(rma: *Mutable, a: Const, b: Const, limbs_buffer: []Limb, allocator: ?Allocator) void {
795	var buf_index: usize = 0;
796
797	const a_copy = if (rma.limbs.ptr == a.limbs.ptr) blk: {
798	const start = buf_index;
799	@memcpy(limbs_buffer[buf_index..][0..a.limbs.len], a.limbs);
800	buf_index += a.limbs.len;
801	break :blk a.toMutable(limbs_buffer[start..buf_index]).toConst();
802	} else a;
803
804	const b_copy = if (rma.limbs.ptr == b.limbs.ptr) blk: {
805	const start = buf_index;
806	@memcpy(limbs_buffer[buf_index..][0..b.limbs.len], b.limbs);
807	buf_index += b.limbs.len;
808	break :blk b.toMutable(limbs_buffer[start..buf_index]).toConst();
809	} else b;
810
811	return rma.mulNoAlias(a_copy, b_copy, allocator);
812	}
813
814	/// rma = a * b
815	///
816	/// `rma` may not alias with `a` or `b`.
817	/// `a` and `b` may alias with each other.
818	///
819	/// Asserts the result fits in `rma`. An upper bound on the number of limbs needed by
820	/// rma is given by `a.limbs.len + b.limbs.len`.
821	///
822	/// If `allocator` is provided, it will be used for temporary storage to improve
823	/// multiplication performance. `error.OutOfMemory` is handled with a fallback algorithm.
824	pub fn mulNoAlias(rma: *Mutable, a: Const, b: Const, allocator: ?Allocator) void {
825	assert(rma.limbs.ptr != a.limbs.ptr); // illegal aliasing
826	assert(rma.limbs.ptr != b.limbs.ptr); // illegal aliasing
827
828	if (a.limbs.len == 1 and b.limbs.len == 1) {
829	rma.limbs[0], const overflow_bit = @mulWithOverflow(a.limbs[0], b.limbs[0]);
830	if (overflow_bit == 0) {
831	rma.len = 1;
832	rma.positive = (a.positive == b.positive) or rma.limbs[0] == 0;
833	return;
834	}
835	}
836
837	@memset(rma.limbs[0 .. a.limbs.len + b.limbs.len], 0);
838
839	llmulacc(.add, allocator, rma.limbs, a.limbs, b.limbs);
840
841	rma.normalize(a.limbs.len + b.limbs.len);
842	rma.positive = (a.positive == b.positive);
843	}
844
845	/// rma = a * b with 2s-complement wrapping semantics.
846	///
847	/// `rma` may alias with `a` or `b`.
848	/// `a` and `b` may alias with each other.
849	///
850	/// Asserts the result fits in `rma`. An upper bound on the number of limbs needed by
851	/// rma is given by `a.limbs.len + b.limbs.len`.
852	///
853	/// `limbs_buffer` is used for temporary storage. The amount required is given by `calcMulWrapLimbsBufferLen`.
854	pub fn mulWrap(
855	rma: *Mutable,
856	a: Const,
857	b: Const,
858	signedness: Signedness,
859	bit_count: usize,
860	limbs_buffer: []Limb,
861	allocator: ?Allocator,
862	) void {
863	var buf_index: usize = 0;
864	const req_limbs = calcTwosCompLimbCount(bit_count);
865
866	const a_copy = if (rma.limbs.ptr == a.limbs.ptr) blk: {
867	const start = buf_index;
868	const a_len = @min(req_limbs, a.limbs.len);
869	@memcpy(limbs_buffer[buf_index..][0..a_len], a.limbs[0..a_len]);
870	buf_index += a_len;
871	break :blk a.toMutable(limbs_buffer[start..buf_index]).toConst();
872	} else a;
873
874	const b_copy = if (rma.limbs.ptr == b.limbs.ptr) blk: {
875	const start = buf_index;
876	const b_len = @min(req_limbs, b.limbs.len);
877	@memcpy(limbs_buffer[buf_index..][0..b_len], b.limbs[0..b_len]);
878	buf_index += b_len;
879	break :blk a.toMutable(limbs_buffer[start..buf_index]).toConst();
880	} else b;
881
882	return rma.mulWrapNoAlias(a_copy, b_copy, signedness, bit_count, allocator);
883	}
884
885	/// rma = a * b with 2s-complement wrapping semantics.
886	///
887	/// `rma` may not alias with `a` or `b`.
888	/// `a` and `b` may alias with each other.
889	///
890	/// Asserts the result fits in `rma`. An upper bound on the number of limbs needed by
891	/// rma is given by `a.limbs.len + b.limbs.len`.
892	///
893	/// If `allocator` is provided, it will be used for temporary storage to improve
894	/// multiplication performance. `error.OutOfMemory` is handled with a fallback algorithm.
895	pub fn mulWrapNoAlias(
896	rma: *Mutable,
897	a: Const,
898	b: Const,
899	signedness: Signedness,
900	bit_count: usize,
901	allocator: ?Allocator,
902	) void {
903	assert(rma.limbs.ptr != a.limbs.ptr); // illegal aliasing
904	assert(rma.limbs.ptr != b.limbs.ptr); // illegal aliasing
905
906	const req_limbs = calcTwosCompLimbCount(bit_count);
907
908	// We can ignore the upper bits here, those results will be discarded anyway.
909	const a_limbs = a.limbs[0..@min(req_limbs, a.limbs.len)];
910	const b_limbs = b.limbs[0..@min(req_limbs, b.limbs.len)];
911
912	@memset(rma.limbs[0..req_limbs], 0);
913
914	llmulacc(.add, allocator, rma.limbs, a_limbs, b_limbs);
915	rma.normalize(@min(req_limbs, a.limbs.len + b.limbs.len));
916	rma.positive = (a.positive == b.positive);
917	rma.truncate(rma.toConst(), signedness, bit_count);
918	}
919
920	/// r = @bitReverse(a) with 2s-complement semantics.
921	/// r and a may be aliases.
922	///
923	/// Asserts the result fits in `r`. Upper bound on the number of limbs needed by
924	/// r is `calcTwosCompLimbCount(bit_count)`.
925	pub fn bitReverse(r: *Mutable, a: Const, signedness: Signedness, bit_count: usize) void {
926	if (bit_count == 0) {
927	r.limbs[0] = 0;
928	r.len = 1;
929	r.positive = true;
930	return;
931	}
932
933	r.copy(a);
934
935	const limbs_required = calcTwosCompLimbCount(bit_count);
936
937	if (!a.positive) {
938	r.positive = true; // Negate.
939	r.bitNotWrap(r.toConst(), .unsigned, bit_count); // Bitwise NOT.
940	r.addScalar(r.toConst(), 1); // Add one.
941	} else if (limbs_required > a.limbs.len) {
942	// Zero-extend to our output length
943	for (r.limbs[a.limbs.len..limbs_required]) |*limb| {
944	limb.* = 0;
945	}
946	r.len = limbs_required;
947	}
948
949	// 0b0..01..1000 with @log2(@sizeOf(Limb)) consecutive ones
950	const endian_mask: usize = (@sizeOf(Limb) - 1) << 3;
951
952	const bytes = std.mem.sliceAsBytes(r.limbs);
953
954	var k: usize = 0;
955	while (k < ((bit_count + 1) / 2)) : (k += 1) {
956	var i = k;
957	var rev_i = bit_count - i - 1;
958
959	// This "endian mask" remaps a low (LE) byte to the corresponding high
960	// (BE) byte in the Limb, without changing which limbs we are indexing
961	if (native_endian == .big) {
962	i ^= endian_mask;
963	rev_i ^= endian_mask;
964	}
965
966	const bit_i = std.mem.readPackedInt(u1, bytes, i, .little);
967	const bit_rev_i = std.mem.readPackedInt(u1, bytes, rev_i, .little);
968	std.mem.writePackedInt(u1, bytes, i, bit_rev_i, .little);
969	std.mem.writePackedInt(u1, bytes, rev_i, bit_i, .little);
970	}
971
972	// Calculate signed-magnitude representation for output
973	if (signedness == .signed) {
974	const last_bit = switch (native_endian) {
975	.little => std.mem.readPackedInt(u1, bytes, bit_count - 1, .little),
976	.big => std.mem.readPackedInt(u1, bytes, (bit_count - 1) ^ endian_mask, .little),
977	};
978	if (last_bit == 1) {
979	r.bitNotWrap(r.toConst(), .unsigned, bit_count); // Bitwise NOT.
980	r.addScalar(r.toConst(), 1); // Add one.
981	r.positive = false; // Negate.
982	}
983	}
984	r.normalize(r.len);
985	}
986
987	/// r = @byteSwap(a) with 2s-complement semantics.
988	/// r and a may be aliases.
989	///
990	/// Asserts the result fits in `r`. Upper bound on the number of limbs needed by
991	/// r is `calcTwosCompLimbCount(8*byte_count)`.
992	pub fn byteSwap(r: *Mutable, a: Const, signedness: Signedness, byte_count: usize) void {
993	if (byte_count == 0) {
994	r.limbs[0] = 0;
995	r.len = 1;
996	r.positive = true;
997	return;
998	}
999
1000	r.copy(a);
1001	const limbs_required = calcTwosCompLimbCount(8 * byte_count);
1002
1003	if (!a.positive) {
1004	r.positive = true; // Negate.
1005	r.bitNotWrap(r.toConst(), .unsigned, 8 * byte_count); // Bitwise NOT.
1006	r.addScalar(r.toConst(), 1); // Add one.
1007	} else if (limbs_required > a.limbs.len) {
1008	// Zero-extend to our output length
1009	for (r.limbs[a.limbs.len..limbs_required]) |*limb| {
1010	limb.* = 0;
1011	}
1012	r.len = limbs_required;
1013	}
1014
1015	// 0b0..01..1 with @log2(@sizeOf(Limb)) trailing ones
1016	const endian_mask: usize = @sizeOf(Limb) - 1;
1017
1018	var bytes = std.mem.sliceAsBytes(r.limbs);
1019	assert(bytes.len >= byte_count);
1020
1021	var k: usize = 0;
1022	while (k < (byte_count + 1) / 2) : (k += 1) {
1023	var i = k;
1024	var rev_i = byte_count - k - 1;
1025
1026	// This "endian mask" remaps a low (LE) byte to the corresponding high
1027	// (BE) byte in the Limb, without changing which limbs we are indexing
1028	if (native_endian == .big) {
1029	i ^= endian_mask;
1030	rev_i ^= endian_mask;
1031	}
1032
1033	const byte_i = bytes[i];
1034	const byte_rev_i = bytes[rev_i];
1035	bytes[rev_i] = byte_i;
1036	bytes[i] = byte_rev_i;
1037	}
1038
1039	// Calculate signed-magnitude representation for output
1040	if (signedness == .signed) {
1041	const last_byte = switch (native_endian) {
1042	.little => bytes[byte_count - 1],
1043	.big => bytes[(byte_count - 1) ^ endian_mask],
1044	};
1045
1046	if (last_byte & (1 << 7) != 0) { // Check sign bit of last byte
1047	r.bitNotWrap(r.toConst(), .unsigned, 8 * byte_count); // Bitwise NOT.
1048	r.addScalar(r.toConst(), 1); // Add one.
1049	r.positive = false; // Negate.
1050	}
1051	}
1052	r.normalize(r.len);
1053	}
1054
1055	/// r = @popCount(a) with 2s-complement semantics.
1056	/// r and a may be aliases.
1057	///
1058	/// Assets the result fits in `r`. Upper bound on the number of limbs needed by
1059	/// r is `calcTwosCompLimbCount(bit_count)`.
1060	pub fn popCount(r: *Mutable, a: Const, bit_count: usize) void {
1061	r.copy(a);
1062
1063	if (!a.positive) {
1064	r.positive = true; // Negate.
1065	r.bitNotWrap(r.toConst(), .unsigned, bit_count); // Bitwise NOT.
1066	r.addScalar(r.toConst(), 1); // Add one.
1067	}
1068
1069	var sum: Limb = 0;
1070	for (r.limbs[0..r.len]) |limb| {
1071	sum += @popCount(limb);
1072	}
1073	r.set(sum);
1074	}
1075
1076	/// rma = a * a
1077	///
1078	/// `rma` may not alias with `a`.
1079	///
1080	/// Asserts the result fits in `rma`. An upper bound on the number of limbs needed by
1081	/// rma is given by `2 * a.limbs.len + 1`.
1082	///
1083	/// If `allocator` is provided, it will be used for temporary storage to improve
1084	/// multiplication performance. `error.OutOfMemory` is handled with a fallback algorithm.
1085	pub fn sqrNoAlias(rma: *Mutable, a: Const, opt_allocator: ?Allocator) void {
1086	_ = opt_allocator;
1087	assert(rma.limbs.ptr != a.limbs.ptr); // illegal aliasing
1088
1089	@memset(rma.limbs, 0);
1090
1091	llsquareBasecase(rma.limbs, a.limbs);
1092
1093	rma.normalize(2 * a.limbs.len + 1);
1094	rma.positive = true;
1095	}
1096
1097	/// q = a / b (rem r)
1098	///
1099	/// a / b are floored (rounded towards 0).
1100	/// q may alias with a or b.
1101	///
1102	/// Asserts there is enough memory to store q and r.
1103	/// The upper bound for r limb count is `b.limbs.len`.
1104	/// The upper bound for q limb count is given by `a.limbs`.
1105	///
1106	/// `limbs_buffer` is used for temporary storage. The amount required is given by `calcDivLimbsBufferLen`.
1107	pub fn divFloor(
1108	q: *Mutable,
1109	r: *Mutable,
1110	a: Const,
1111	b: Const,
1112	limbs_buffer: []Limb,
1113	) void {
1114	const sep = a.limbs.len + 2;
1115	var x = a.toMutable(limbs_buffer[0..sep]);
1116	var y = b.toMutable(limbs_buffer[sep..]);
1117
1118	div(q, r, &x, &y);
1119
1120	// Note, `div` performs truncating division, which satisfies
1121	// @divTrunc(a, b) * b + @rem(a, b) = a
1122	// so r = a - @divTrunc(a, b) * b
1123	// Note, @rem(a, -b) = @rem(-b, a) = -@rem(a, b) = -@rem(-a, -b)
1124	// For divTrunc, we want to perform
1125	// @divFloor(a, b) * b + @mod(a, b) = a
1126	// Note:
1127	// @divFloor(-a, b)
1128	// = @divFloor(a, -b)
1129	// = -@divCeil(a, b)
1130	// = -@divFloor(a + b - 1, b)
1131	// = -@divTrunc(a + b - 1, b)
1132
1133	// Note (1):
1134	// @divTrunc(a + b - 1, b) * b + @rem(a + b - 1, b) = a + b - 1
1135	// = @divTrunc(a + b - 1, b) * b + @rem(a - 1, b) = a + b - 1
1136	// = @divTrunc(a + b - 1, b) * b + @rem(a - 1, b) - b + 1 = a
1137
1138	if (a.positive and b.positive) {
1139	// Positive-positive case, don't need to do anything.
1140	} else if (a.positive and !b.positive) {
1141	// a/-b -> q is negative, and so we need to fix flooring.
1142	// Subtract one to make the division flooring.
1143
1144	// @divFloor(a, -b) * -b + @mod(a, -b) = a
1145	// If b divides a exactly, we have @divFloor(a, -b) * -b = a
1146	// Else, we have @divFloor(a, -b) * -b > a, so @mod(a, -b) becomes negative
1147
1148	// We have:
1149	// @divFloor(a, -b) * -b + @mod(a, -b) = a
1150	// = -@divTrunc(a + b - 1, b) * -b + @mod(a, -b) = a
1151	// = @divTrunc(a + b - 1, b) * b + @mod(a, -b) = a
1152
1153	// Substitute a for (1):
1154	// @divTrunc(a + b - 1, b) * b + @rem(a - 1, b) - b + 1 = @divTrunc(a + b - 1, b) * b + @mod(a, -b)
1155	// Yields:
1156	// @mod(a, -b) = @rem(a - 1, b) - b + 1
1157	// Note that `r` holds @rem(a, b) at this point.
1158	//
1159	// If @rem(a, b) is not 0:
1160	// @rem(a - 1, b) = @rem(a, b) - 1
1161	// => @mod(a, -b) = @rem(a, b) - 1 - b + 1 = @rem(a, b) - b
1162	// Else:
1163	// @rem(a - 1, b) = @rem(a + b - 1, b) = @rem(b - 1, b) = b - 1
1164	// => @mod(a, -b) = b - 1 - b + 1 = 0
1165	if (!r.eqlZero()) {
1166	q.addScalar(q.toConst(), -1);
1167	r.positive = true;
1168	r.sub(r.toConst(), y.toConst().abs());
1169	}
1170	} else if (!a.positive and b.positive) {
1171	// -a/b -> q is negative, and so we need to fix flooring.
1172	// Subtract one to make the division flooring.
1173
1174	// @divFloor(-a, b) * b + @mod(-a, b) = a
1175	// If b divides a exactly, we have @divFloor(-a, b) * b = -a
1176	// Else, we have @divFloor(-a, b) * b < -a, so @mod(-a, b) becomes positive
1177
1178	// We have:
1179	// @divFloor(-a, b) * b + @mod(-a, b) = -a
1180	// = -@divTrunc(a + b - 1, b) * b + @mod(-a, b) = -a
1181	// = @divTrunc(a + b - 1, b) * b - @mod(-a, b) = a
1182
1183	// Substitute a for (1):
1184	// @divTrunc(a + b - 1, b) * b + @rem(a - 1, b) - b + 1 = @divTrunc(a + b - 1, b) * b - @mod(-a, b)
1185	// Yields:
1186	// @rem(a - 1, b) - b + 1 = -@mod(-a, b)
1187	// => -@mod(-a, b) = @rem(a - 1, b) - b + 1
1188	// => @mod(-a, b) = -(@rem(a - 1, b) - b + 1) = -@rem(a - 1, b) + b - 1
1189	//
1190	// If @rem(a, b) is not 0:
1191	// @rem(a - 1, b) = @rem(a, b) - 1
1192	// => @mod(-a, b) = -(@rem(a, b) - 1) + b - 1 = -@rem(a, b) + 1 + b - 1 = -@rem(a, b) + b
1193	// Else :
1194	// @rem(a - 1, b) = b - 1
1195	// => @mod(-a, b) = -(b - 1) + b - 1 = 0
1196	if (!r.eqlZero()) {
1197	q.addScalar(q.toConst(), -1);
1198	r.positive = false;
1199	r.add(r.toConst(), y.toConst().abs());
1200	}
1201	} else if (!a.positive and !b.positive) {
1202	// a/b -> q is positive, don't need to do anything to fix flooring.
1203
1204	// @divFloor(-a, -b) * -b + @mod(-a, -b) = -a
1205	// If b divides a exactly, we have @divFloor(-a, -b) * -b = -a
1206	// Else, we have @divFloor(-a, -b) * -b > -a, so @mod(-a, -b) becomes negative
1207
1208	// We have:
1209	// @divFloor(-a, -b) * -b + @mod(-a, -b) = -a
1210	// = @divTrunc(a, b) * -b + @mod(-a, -b) = -a
1211	// = @divTrunc(a, b) * b - @mod(-a, -b) = a
1212
1213	// We also have:
1214	// @divTrunc(a, b) * b + @rem(a, b) = a
1215
1216	// Substitute a:
1217	// @divTrunc(a, b) * b + @rem(a, b) = @divTrunc(a, b) * b - @mod(-a, -b)
1218	// => @rem(a, b) = -@mod(-a, -b)
1219	// => @mod(-a, -b) = -@rem(a, b)
1220	r.positive = false;
1221	}
1222	}
1223
1224	/// q = a / b (rem r)
1225	///
1226	/// a / b are ceiled (rounded towards +inf).
1227	/// q may alias with a or b.
1228	///
1229	/// Asserts there is enough memory to store q and r.
1230	/// The upper bound for r limb count is `b.limbs.len`.
1231	/// The upper bound for q limb count is given by `a.limbs`.
1232	///
1233	/// `limbs_buffer` is used for temporary storage. The amount required is given by `calcDivLimbsBufferLen`.
1234	pub fn divCeil(
1235	q: *Mutable,
1236	r: *Mutable,
1237	a: Const,
1238	b: Const,
1239	limbs_buffer: []Limb,
1240	) void {
1241	const sep = a.limbs.len + 2;
1242	var x = a.toMutable(limbs_buffer[0..sep]);
1243	var y = b.toMutable(limbs_buffer[sep..]);
1244
1245	// div performs truncating division (@divTrunc) which rounds towards negative
1246	// infinity if the result is positive and towards positive infinity if the result is
1247	// negative.
1248	div(q, r, &x, &y);
1249
1250	// @rem gives the remainder after @divTrunc, and is defined by:
1251	// x * @divTrunc(x, y) + @rem(x, y) = x
1252	// For all integers x, y with y != 0.
1253	// In the following comments, a, b will be integers with a >= 0, b > 0, and we will take
1254	// modCeil to be the remainder after @divCeil, defined by:
1255	// x * @divCeil(x, y) + modCeil(x, y) = x
1256	// For all integers x, y with y != 0.
1257
1258	if (a.positive != b.positive or r.eqlZero()) {
1259	// In this case either the result is negative or the remainder is 0.
1260	// If the result is negative then the default truncating division already rounds
1261	// towards positive infinity, so no adjustment is needed.
1262	// If the remainder is 0 then the division is exact and no adjustment is needed.
1263	} else {
1264	// Same sign.
1265	// We have:
1266	// modCeil(a, b) != 0
1267	// => @divCeil(a, b) = @divTrunc(a, b) + 1
1268	// And:
1269	// b * @divTrunc(a, b) + @rem(a, b) = a
1270	// b * @divCeil(a, b) + modCeil(a, b) = a
1271	// => b * @divTrunc(a, b) + b + modCeil(a, b) = a
1272	// => modCeil(a, b) = @rem(a, b) - b
1273	//
1274	// This works for both positive and negative b because b keeps its sign.
1275	q.addScalar(q.toConst(), 1);
1276	r.sub(r.toConst(), y.toConst());
1277	}
1278	}
1279
1280	/// q = a / b (rem r)
1281	///
1282	/// a / b are truncated (rounded towards -inf).
1283	/// q may alias with a or b.
1284	///
1285	/// Asserts there is enough memory to store q and r.
1286	/// The upper bound for r limb count is `b.limbs.len`.
1287	/// The upper bound for q limb count is given by `a.limbs.len`.
1288	///
1289	/// `limbs_buffer` is used for temporary storage. The amount required is given by `calcDivLimbsBufferLen`.
1290	pub fn divTrunc(
1291	q: *Mutable,
1292	r: *Mutable,
1293	a: Const,
1294	b: Const,
1295	limbs_buffer: []Limb,
1296	) void {
1297	const sep = a.limbs.len + 2;
1298	var x = a.toMutable(limbs_buffer[0..sep]);
1299	var y = b.toMutable(limbs_buffer[sep..]);
1300
1301	div(q, r, &x, &y);
1302	}
1303
1304	/// r = a << shift, in other words, r = a * 2^shift
1305	///
1306	/// r and a may alias.
1307	///
1308	/// Asserts there is enough memory to fit the result. The upper bound Limb count is
1309	/// `a.limbs.len + (shift / (@sizeOf(Limb) * 8))`.
1310	pub fn shiftLeft(r: *Mutable, a: Const, shift: usize) void {
1311	const new_len = llshl(r.limbs, a.limbs, shift);
1312	r.normalize(new_len);
1313	r.positive = a.positive;
1314	}
1315
1316	/// r = a <<| shift with 2s-complement saturating semantics.
1317	///
1318	/// r and a may alias.
1319	///
1320	/// Asserts there is enough memory to fit the result. The upper bound Limb count is
1321	/// r is `calcTwosCompLimbCount(bit_count)`.
1322	pub fn shiftLeftSat(r: *Mutable, a: Const, shift: usize, signedness: Signedness, bit_count: usize) void {
1323	// Special case: When the argument is negative, but the result is supposed to be unsigned,
1324	// return 0 in all cases.
1325	if (!a.positive and signedness == .unsigned) {
1326	r.set(0);
1327	return;
1328	}
1329
1330	// Check whether the shift is going to overflow. This is the case
1331	// when (in 2s complement) any bit above `bit_count - shift` is set in the unshifted value.
1332	// Note, the sign bit is not counted here.
1333
1334	// Handle shifts larger than the target type. This also deals with
1335	// 0-bit integers.
1336	if (bit_count <= shift) {
1337	// In this case, there is only no overflow if `a` is zero.
1338	if (a.eqlZero()) {
1339	r.set(0);
1340	} else {
1341	r.setTwosCompIntLimit(if (a.positive) .max else .min, signedness, bit_count);
1342	}
1343	return;
1344	}
1345
1346	const checkbit = bit_count - shift - @intFromBool(signedness == .signed);
1347	// If `checkbit` and more significant bits are zero, no overflow will take place.
1348
1349	if (checkbit >= a.limbs.len * limb_bits) {
1350	// `checkbit` is outside the range of a, so definitely no overflow will take place. We
1351	// can defer to a normal shift.
1352	// Note that if `a` is normalized (which we assume), this checks for set bits in the upper limbs.
1353
1354	// Note, in this case r should already have enough limbs required to perform the normal shift.
1355	// In this case the shift of the most significant limb may still overflow.
1356	r.shiftLeft(a, shift);
1357	return;
1358	} else if (checkbit < (a.limbs.len - 1) * limb_bits) {
1359	// `checkbit` is not in the most significant limb. If `a` is normalized the most significant
1360	// limb will not be zero, so in this case we need to saturate. Note that `a.limbs.len` must be
1361	// at least one according to normalization rules.
1362
1363	r.setTwosCompIntLimit(if (a.positive) .max else .min, signedness, bit_count);
1364	return;
1365	}
1366
1367	// Generate a mask with the bits to check in the most significant limb. We'll need to check
1368	// all bits with equal or more significance than checkbit.
1369	// const msb = @truncate(Log2Limb, checkbit);
1370	// const checkmask = (@as(Limb, 1) << msb) -% 1;
1371
1372	if (a.limbs[a.limbs.len - 1] >> @as(Log2Limb, @truncate(checkbit)) != 0) {
1373	// Need to saturate.
1374	r.setTwosCompIntLimit(if (a.positive) .max else .min, signedness, bit_count);
1375	return;
1376	}
1377
1378	// This shift should not be able to overflow, so invoke llshl and normalize manually
1379	// to avoid the extra required limb.
1380	const new_len = llshl(r.limbs, a.limbs, shift);
1381	r.normalize(new_len);
1382	r.positive = a.positive;
1383	}
1384
1385	/// r = a >> shift
1386	/// r and a may alias.
1387	///
1388	/// Asserts there is enough memory to fit the result. The upper bound Limb count is
1389	/// `a.limbs.len - (shift / (@bitSizeOf(Limb)))`.
1390	pub fn shiftRight(r: *Mutable, a: Const, shift: usize) void {
1391	const full_limbs_shifted_out = shift / limb_bits;
1392	const remaining_bits_shifted_out = shift % limb_bits;
1393	if (a.limbs.len <= full_limbs_shifted_out) {
1394	// Shifting negative numbers converges to -1 instead of 0
1395	if (a.positive) {
1396	r.len = 1;
1397	r.positive = true;
1398	r.limbs[0] = 0;
1399	} else {
1400	r.len = 1;
1401	r.positive = false;
1402	r.limbs[0] = 1;
1403	}
1404	return;
1405	}
1406	const nonzero_negative_shiftout = if (a.positive) false else nonzero: {
1407	for (a.limbs[0..full_limbs_shifted_out]) |x| {
1408	if (x != 0)
1409	break :nonzero true;
1410	}
1411	if (remaining_bits_shifted_out == 0)
1412	break :nonzero false;
1413	const not_covered: Log2Limb = @intCast(limb_bits - remaining_bits_shifted_out);
1414	break :nonzero a.limbs[full_limbs_shifted_out] << not_covered != 0;
1415	};
1416
1417	const new_len = llshr(r.limbs, a.limbs, shift);
1418
1419	r.len = new_len;
1420	r.positive = a.positive;
1421	if (nonzero_negative_shiftout) r.addScalar(r.toConst(), -1);
1422	r.normalize(r.len);
1423	}
1424
1425	/// r = ~a under 2s complement wrapping semantics.
1426	/// r may alias with a.
1427	///
1428	/// Assets that r has enough limbs to store the result. The upper bound Limb count is
1429	/// r is `calcTwosCompLimbCount(bit_count)`.
1430	pub fn bitNotWrap(r: *Mutable, a: Const, signedness: Signedness, bit_count: usize) void {
1431	r.copy(a.negate());
1432	const negative_one: Const = .{ .limbs = &.{1}, .positive = false };
1433	_ = r.addWrap(r.toConst(), negative_one, signedness, bit_count);
1434	}
1435
1436	/// r = a | b under 2s complement semantics.
1437	/// r may alias with a or b.
1438	///
1439	/// a and b are zero-extended to the longer of a or b.
1440	///
1441	/// Asserts that r has enough limbs to store the result. Upper bound is `@max(a.limbs.len, b.limbs.len)`.
1442	pub fn bitOr(r: *Mutable, a: Const, b: Const) void {
1443	// Trivial cases, llsignedor does not support zero.
1444	if (a.eqlZero()) {
1445	r.copy(b);
1446	return;
1447	} else if (b.eqlZero()) {
1448	r.copy(a);
1449	return;
1450	}
1451
1452	if (a.limbs.len >= b.limbs.len) {
1453	r.positive = llsignedor(r.limbs, a.limbs, a.positive, b.limbs, b.positive);
1454	r.normalize(if (b.positive) a.limbs.len else b.limbs.len);
1455	} else {
1456	r.positive = llsignedor(r.limbs, b.limbs, b.positive, a.limbs, a.positive);
1457	r.normalize(if (a.positive) b.limbs.len else a.limbs.len);
1458	}
1459	}
1460
1461	/// r = a & b under 2s complement semantics.
1462	/// r may alias with a or b.
1463	///
1464	/// Asserts that r has enough limbs to store the result.
1465	/// If only a is positive, the upper bound is `a.limbs.len`.
1466	/// If only b is positive, the upper bound is `b.limbs.len`.
1467	/// If a and b are positive, the upper bound is `@min(a.limbs.len, b.limbs.len)`.
1468	/// If a and b are negative, the upper bound is `@max(a.limbs.len, b.limbs.len) + 1`.
1469	pub fn bitAnd(r: *Mutable, a: Const, b: Const) void {
1470	// Trivial cases, llsignedand does not support zero.
1471	if (a.eqlZero()) {
1472	r.copy(a);
1473	return;
1474	} else if (b.eqlZero()) {
1475	r.copy(b);
1476	return;
1477	}
1478
1479	if (a.limbs.len >= b.limbs.len) {
1480	r.positive = llsignedand(r.limbs, a.limbs, a.positive, b.limbs, b.positive);
1481	r.normalize(if (b.positive) b.limbs.len else if (a.positive) a.limbs.len else a.limbs.len + 1);
1482	} else {
1483	r.positive = llsignedand(r.limbs, b.limbs, b.positive, a.limbs, a.positive);
1484	r.normalize(if (a.positive) a.limbs.len else if (b.positive) b.limbs.len else b.limbs.len + 1);
1485	}
1486	}
1487
1488	/// r = a ^ b under 2s complement semantics.
1489	/// r may alias with a or b.
1490	///
1491	/// Asserts that r has enough limbs to store the result. If a and b share the same signedness, the
1492	/// upper bound is `@max(a.limbs.len, b.limbs.len)`. Otherwise, if either a or b is negative
1493	/// but not both, the upper bound is `@max(a.limbs.len, b.limbs.len) + 1`.
1494	pub fn bitXor(r: *Mutable, a: Const, b: Const) void {
1495	// Trivial cases, because llsignedxor does not support negative zero.
1496	if (a.eqlZero()) {
1497	r.copy(b);
1498	return;
1499	} else if (b.eqlZero()) {
1500	r.copy(a);
1501	return;
1502	}
1503
1504	if (a.limbs.len > b.limbs.len) {
1505	r.positive = llsignedxor(r.limbs, a.limbs, a.positive, b.limbs, b.positive);
1506	r.normalize(a.limbs.len + @intFromBool(a.positive != b.positive));
1507	} else {
1508	r.positive = llsignedxor(r.limbs, b.limbs, b.positive, a.limbs, a.positive);
1509	r.normalize(b.limbs.len + @intFromBool(a.positive != b.positive));
1510	}
1511	}
1512
1513	/// rma may alias x or y.
1514	/// x and y may alias each other.
1515	/// Asserts that `rma` has enough limbs to store the result. Upper bound is
1516	/// `@min(x.limbs.len, y.limbs.len)`.
1517	///
1518	/// `limbs_buffer` is used for temporary storage during the operation. When this function returns,
1519	/// it will have the same length as it had when the function was called.
1520	pub fn gcd(rma: *Mutable, x: Const, y: Const, limbs_buffer: *std.array_list.Managed(Limb)) !void {
1521	const prev_len = limbs_buffer.items.len;
1522	defer limbs_buffer.shrinkRetainingCapacity(prev_len);
1523	const x_copy = if (rma.limbs.ptr == x.limbs.ptr) blk: {
1524	const start = limbs_buffer.items.len;
1525	try limbs_buffer.appendSlice(x.limbs);
1526	break :blk x.toMutable(limbs_buffer.items[start..]).toConst();
1527	} else x;
1528	const y_copy = if (rma.limbs.ptr == y.limbs.ptr) blk: {
1529	const start = limbs_buffer.items.len;
1530	try limbs_buffer.appendSlice(y.limbs);
1531	break :blk y.toMutable(limbs_buffer.items[start..]).toConst();
1532	} else y;
1533
1534	return gcdLehmer(rma, x_copy, y_copy, limbs_buffer);
1535	}
1536
1537	/// q = a ^ b
1538	///
1539	/// r may not alias a.
1540	///
1541	/// Asserts that `r` has enough limbs to store the result. Upper bound is
1542	/// `calcPowLimbsBufferLen(a.bitCountAbs(), b)`.
1543	///
1544	/// `limbs_buffer` is used for temporary storage.
1545	/// The amount required is given by `calcPowLimbsBufferLen`.
1546	pub fn pow(r: *Mutable, a: Const, b: u32, limbs_buffer: []Limb) void {
1547	assert(r.limbs.ptr != a.limbs.ptr); // illegal aliasing
1548
1549	// Handle all the trivial cases first
1550	switch (b) {
1551	0 => {
1552	// a^0 = 1
1553	return r.set(1);
1554	},
1555	1 => {
1556	// a^1 = a
1557	return r.copy(a);
1558	},
1559	else => {},
1560	}
1561
1562	if (a.eqlZero()) {
1563	// 0^b = 0
1564	return r.set(0);
1565	} else if (a.limbs.len == 1 and a.limbs[0] == 1) {
1566	// 1^b = 1 and -1^b = ±1
1567	r.set(1);
1568	r.positive = a.positive or (b & 1) == 0;
1569	return;
1570	}
1571
1572	// Here a>1 and b>1
1573	const needed_limbs = calcPowLimbsBufferLen(a.bitCountAbs(), b);
1574	assert(r.limbs.len >= needed_limbs);
1575	assert(limbs_buffer.len >= needed_limbs);
1576
1577	llpow(r.limbs, a.limbs, b, limbs_buffer);
1578
1579	r.normalize(needed_limbs);
1580	r.positive = a.positive or (b & 1) == 0;
1581	}
1582
1583	/// r = ⌊√a⌋
1584	///
1585	/// r may alias a.
1586	///
1587	/// Asserts that `r` has enough limbs to store the result. Upper bound is
1588	/// `(a.limbs.len - 1) / 2 + 1`.
1589	///
1590	/// `limbs_buffer` is used for temporary storage.
1591	/// The amount required is given by `calcSqrtLimbsBufferLen`.
1592	pub fn sqrt(
1593	r: *Mutable,
1594	a: Const,
1595	limbs_buffer: []Limb,
1596	) void {
1597	// Brent and Zimmermann, Modern Computer Arithmetic, Algorithm 1.13 SqrtInt
1598	// https://members.loria.fr/PZimmermann/mca/pub226.html
1599	var buf_index: usize = 0;
1600	var t = b: {
1601	const start = buf_index;
1602	buf_index += a.limbs.len;
1603	break :b Mutable.init(limbs_buffer[start..buf_index], 0);
1604	};
1605	var u = b: {
1606	const start = buf_index;
1607	const shift = (a.bitCountAbs() + 1) / 2;
1608	buf_index += 1 + ((shift / limb_bits) + 1);
1609	var m = Mutable.init(limbs_buffer[start..buf_index], 1);
1610	m.shiftLeft(m.toConst(), shift); // u must be >= ⌊√a⌋, and should be as small as possible for efficiency
1611	break :b m;
1612	};
1613	var s = b: {
1614	const start = buf_index;
1615	buf_index += u.limbs.len;
1616	break :b u.toConst().toMutable(limbs_buffer[start..buf_index]);
1617	};
1618	var rem = b: {
1619	const start = buf_index;
1620	buf_index += s.limbs.len;
1621	break :b Mutable.init(limbs_buffer[start..buf_index], 0);
1622	};
1623
1624	while (true) {
1625	t.divFloor(&rem, a, s.toConst(), limbs_buffer[buf_index..]);
1626	t.add(t.toConst(), s.toConst());
1627	u.shiftRight(t.toConst(), 1);
1628
1629	if (u.toConst().order(s.toConst()).compare(.gte)) {
1630	r.copy(s.toConst());
1631	return;
1632	}
1633
1634	// Avoid copying u to s by swapping u and s
1635	const tmp_s = s;
1636	s = u;
1637	u = tmp_s;
1638	}
1639	}
1640
1641	/// rma may not alias x or y.
1642	/// x and y may alias each other.
1643	/// Asserts that `rma` has enough limbs to store the result. Upper bound is given by `calcGcdNoAliasLimbLen`.
1644	///
1645	/// `limbs_buffer` is used for temporary storage during the operation.
1646	pub fn gcdNoAlias(rma: *Mutable, x: Const, y: Const, limbs_buffer: *std.array_list.Managed(Limb)) !void {
1647	assert(rma.limbs.ptr != x.limbs.ptr); // illegal aliasing
1648	assert(rma.limbs.ptr != y.limbs.ptr); // illegal aliasing
1649	return gcdLehmer(rma, x, y, limbs_buffer);
1650	}
1651
1652	fn gcdLehmer(result: *Mutable, xa: Const, ya: Const, limbs_buffer: *std.array_list.Managed(Limb)) !void {
1653	var x = try xa.toManaged(limbs_buffer.allocator);
1654	defer x.deinit();
1655	x.abs();
1656
1657	var y = try ya.toManaged(limbs_buffer.allocator);
1658	defer y.deinit();
1659	y.abs();
1660
1661	if (x.toConst().order(y.toConst()) == .lt) {
1662	x.swap(&y);
1663	}
1664
1665	var t_big = try Managed.init(limbs_buffer.allocator);
1666	defer t_big.deinit();
1667
1668	var r = try Managed.init(limbs_buffer.allocator);
1669	defer r.deinit();
1670
1671	var tmp_x = try Managed.init(limbs_buffer.allocator);
1672	defer tmp_x.deinit();
1673
1674	while (y.len() > 1 and !y.eqlZero()) {
1675	assert(x.isPositive() and y.isPositive());
1676	assert(x.len() >= y.len());
1677
1678	var xh: SignedDoubleLimb = x.limbs[x.len() - 1];
1679	var yh: SignedDoubleLimb = if (x.len() > y.len()) 0 else y.limbs[x.len() - 1];
1680
1681	var A: SignedDoubleLimb = 1;
1682	var B: SignedDoubleLimb = 0;
1683	var C: SignedDoubleLimb = 0;
1684	var D: SignedDoubleLimb = 1;
1685
1686	while (yh + C != 0 and yh + D != 0) {
1687	const q = @divFloor(xh + A, yh + C);
1688	const qp = @divFloor(xh + B, yh + D);
1689	if (q != qp) {
1690	break;
1691	}
1692
1693	var t = A - q * C;
1694	A = C;
1695	C = t;
1696	t = B - q * D;
1697	B = D;
1698	D = t;
1699
1700	t = xh - q * yh;
1701	xh = yh;
1702	yh = t;
1703	}
1704
1705	if (B == 0) {
1706	// t_big = x % y, r is unused
1707	try r.divTrunc(&t_big, &x, &y);
1708	assert(t_big.isPositive());
1709
1710	x.swap(&y);
1711	y.swap(&t_big);
1712	} else {
1713	var storage: [8]Limb = undefined;
1714	const Ap = fixedIntFromSignedDoubleLimb(A, storage[0..2]).toManaged(limbs_buffer.allocator);
1715	const Bp = fixedIntFromSignedDoubleLimb(B, storage[2..4]).toManaged(limbs_buffer.allocator);
1716	const Cp = fixedIntFromSignedDoubleLimb(C, storage[4..6]).toManaged(limbs_buffer.allocator);
1717	const Dp = fixedIntFromSignedDoubleLimb(D, storage[6..8]).toManaged(limbs_buffer.allocator);
1718
1719	// t_big = Ax + By
1720	try r.mul(&x, &Ap);
1721	try t_big.mul(&y, &Bp);
1722	try t_big.add(&r, &t_big);
1723
1724	// u = Cx + Dy, r as u
1725	try tmp_x.copy(x.toConst());
1726	try x.mul(&tmp_x, &Cp);
1727	try r.mul(&y, &Dp);
1728	try r.add(&x, &r);
1729
1730	x.swap(&t_big);
1731	y.swap(&r);
1732	}
1733	}
1734
1735	// euclidean algorithm
1736	assert(x.toConst().order(y.toConst()) != .lt);
1737
1738	while (!y.toConst().eqlZero()) {
1739	try t_big.divTrunc(&r, &x, &y);
1740	x.swap(&y);
1741	y.swap(&r);
1742	}
1743
1744	result.copy(x.toConst());
1745	}
1746
1747	// Truncates by default.
1748	fn div(q: *Mutable, r: *Mutable, x: *Mutable, y: *Mutable) void {
1749	assert(!y.eqlZero()); // division by zero
1750	assert(q != r); // illegal aliasing
1751
1752	const q_positive = (x.positive == y.positive);
1753	const r_positive = x.positive;
1754
1755	if (x.toConst().orderAbs(y.toConst()) == .lt) {
1756	// q may alias x so handle r first.
1757	r.copy(x.toConst());
1758	r.positive = r_positive;
1759
1760	q.set(0);
1761	return;
1762	}
1763
1764	// Handle trailing zero-words of divisor/dividend. These are not handled in the following
1765	// algorithms.
1766	// Note, there must be a non-zero limb for either.
1767	// const x_trailing = std.mem.findScalar(Limb, x.limbs[0..x.len], 0).?;
1768	// const y_trailing = std.mem.findScalar(Limb, y.limbs[0..y.len], 0).?;
1769
1770	const x_trailing = for (x.limbs[0..x.len], 0..) |xi, i| {
1771	if (xi != 0) break i;
1772	} else unreachable;
1773
1774	const y_trailing = for (y.limbs[0..y.len], 0..) |yi, i| {
1775	if (yi != 0) break i;
1776	} else unreachable;
1777
1778	const xy_trailing = @min(x_trailing, y_trailing);
1779
1780	if (y.len - xy_trailing == 1) {
1781	const divisor = y.limbs[y.len - 1];
1782
1783	// Optimization for small divisor. By using a half limb we can avoid requiring DoubleLimb
1784	// divisions in the hot code path. This may often require compiler_rt software-emulation.
1785	if (divisor < maxInt(HalfLimb)) {
1786	lldiv0p5(q.limbs, &r.limbs[0], x.limbs[xy_trailing..x.len], @as(HalfLimb, @intCast(divisor)));
1787	} else {
1788	lldiv1(q.limbs, &r.limbs[0], x.limbs[xy_trailing..x.len], divisor);
1789	}
1790
1791	q.normalize(x.len - xy_trailing);
1792	q.positive = q_positive;
1793
1794	r.len = 1;
1795	r.positive = r_positive;
1796	} else {
1797	// Shrink x, y such that the trailing zero limbs shared between are removed.
1798	var x0: Mutable = .{
1799	.limbs = x.limbs[xy_trailing..],
1800	.len = x.len - xy_trailing,
1801	.positive = true,
1802	};
1803
1804	var y0: Mutable = .{
1805	.limbs = y.limbs[xy_trailing..],
1806	.len = y.len - xy_trailing,
1807	.positive = true,
1808	};
1809
1810	divmod(q, r, &x0, &y0);
1811	q.positive = q_positive;
1812
1813	r.positive = r_positive;
1814	}
1815
1816	if (xy_trailing != 0 and r.limbs[r.len - 1] != 0) {
1817	// Manually shift here since we know its limb aligned.
1818	@memmove(r.limbs[xy_trailing..][0..r.len], r.limbs[0..r.len]);
1819	@memset(r.limbs[0..xy_trailing], 0);
1820	r.len += xy_trailing;
1821	}
1822	}
1823
1824	/// Handbook of Applied Cryptography, 14.20
1825	///
1826	/// x = qy + r where 0 <= r < y
1827	/// y is modified but returned intact.
1828	fn divmod(
1829	q: *Mutable,
1830	r: *Mutable,
1831	x: *Mutable,
1832	y: *Mutable,
1833	) void {
1834	// 0.
1835	// Normalize so that y[t] > b/2
1836	const lz = @clz(y.limbs[y.len - 1]);
1837	const norm_shift = if (lz == 0 and y.toConst().isOdd())
1838	limb_bits // Force an extra limb so that y is even.
1839	else
1840	lz;
1841
1842	x.shiftLeft(x.toConst(), norm_shift);
1843	y.shiftLeft(y.toConst(), norm_shift);
1844
1845	const n = x.len - 1;
1846	const t = y.len - 1;
1847	const shift = n - t;
1848
1849	// 1.
1850	// for 0 <= j <= n - t, set q[j] to 0
1851	q.len = shift + 1;
1852	q.positive = true;
1853	@memset(q.limbs[0..q.len], 0);
1854
1855	// 2.
1856	// while x >= y * b^(n - t):
1857	// x -= y * b^(n - t)
1858	// q[n - t] += 1
1859	// Note, this algorithm is performed only once if y[t] > base/2 and y is even, which we
1860	// enforced in step 0. This means we can replace the while with an if.
1861	// Note, multiplication by b^(n - t) comes down to shifting to the right by n - t limbs.
1862	// We can also replace x >= y * b^(n - t) by x/b^(n - t) >= y, and use shifts for that.
1863	{
1864	// x >= y * b^(n - t) can be replaced by x/b^(n - t) >= y.
1865
1866	// 'divide' x by b^(n - t)
1867	var tmp: Mutable = .{
1868	.limbs = x.limbs[shift..],
1869	.len = x.len - shift,
1870	.positive = true,
1871	};
1872
1873	if (tmp.toConst().order(y.toConst()) != .lt) {
1874	// Perform x -= y * b^(n - t)
1875	// Note, we can subtract y from x[n - t..] and get the result without shifting.
1876	// We can also re-use tmp which already contains the relevant part of x. Note that
1877	// this also edits x.
1878	// Due to the check above, this cannot underflow.
1879	tmp.sub(tmp.toConst(), y.toConst());
1880
1881	// tmp.sub normalized tmp, but we need to normalize x now.
1882	x.limbs.len = tmp.limbs.len + shift;
1883
1884	q.limbs[shift] += 1;
1885	}
1886	}
1887
1888	// 3.
1889	// for i from n down to t + 1, do
1890	var i = n;
1891	while (i >= t + 1) : (i -= 1) {
1892	const k = i - t - 1;
1893	// 3.1.
1894	// if x_i == y_t:
1895	// q[i - t - 1] = b - 1
1896	// else:
1897	// q[i - t - 1] = (x[i] * b + x[i - 1]) / y[t]
1898	if (x.limbs[i] == y.limbs[t]) {
1899	q.limbs[k] = maxInt(Limb);
1900	} else {
1901	const q0 = (@as(DoubleLimb, x.limbs[i]) << limb_bits) | @as(DoubleLimb, x.limbs[i - 1]);
1902	const n0 = @as(DoubleLimb, y.limbs[t]);
1903	q.limbs[k] = @as(Limb, @intCast(q0 / n0));
1904	}
1905
1906	// 3.2
1907	// while q[i - t - 1] * (y[t] * b + y[t - 1] > x[i] * b * b + x[i - 1] + x[i - 2]:
1908	// q[i - t - 1] -= 1
1909	// Note, if y[t] > b / 2 this part is repeated no more than twice.
1910
1911	// Extract from y.
1912	const y0 = if (t > 0) y.limbs[t - 1] else 0;
1913	const y1 = y.limbs[t];
1914
1915	// Extract from x.
1916	// Note, big endian.
1917	const tmp0 = [_]Limb{
1918	x.limbs[i],
1919	if (i >= 1) x.limbs[i - 1] else 0,
1920	if (i >= 2) x.limbs[i - 2] else 0,
1921	};
1922
1923	while (true) {
1924	// Ad-hoc 2x1 multiplication with q[i - t - 1].
1925	// Note, big endian.
1926	var tmp1 = [_]Limb{ 0, undefined, undefined };
1927	tmp1[2] = addMulLimbWithCarry(0, y0, q.limbs[k], &tmp1[0]);
1928	tmp1[1] = addMulLimbWithCarry(0, y1, q.limbs[k], &tmp1[0]);
1929
1930	// Big-endian compare
1931	if (mem.order(Limb, &tmp1, &tmp0) != .gt)
1932	break;
1933
1934	q.limbs[k] -= 1;
1935	}
1936
1937	// 3.3.
1938	// x -= q[i - t - 1] * y * b^(i - t - 1)
1939	// Note, we multiply by a single limb here.
1940	// The shift doesn't need to be performed if we add the result of the first multiplication
1941	// to x[i - t - 1].
1942	const underflow = llmulLimb(.sub, x.limbs[k..x.len], y.limbs[0..y.len], q.limbs[k]);
1943
1944	// 3.4.
1945	// if x < 0:
1946	// x += y * b^(i - t - 1)
1947	// q[i - t - 1] -= 1
1948	// Note, we check for x < 0 using the underflow flag from the previous operation.
1949	if (underflow) {
1950	// While we didn't properly set the signedness of x, this operation should 'flow' it back to positive.
1951	llaccum(.add, x.limbs[k..x.len], y.limbs[0..y.len]);
1952	q.limbs[k] -= 1;
1953	}
1954	}
1955
1956	x.normalize(x.len);
1957	q.normalize(q.len);
1958
1959	// De-normalize r and y.
1960	r.shiftRight(x.toConst(), norm_shift);
1961	y.shiftRight(y.toConst(), norm_shift);
1962	}
1963
1964	/// Truncate an integer to a number of bits, following 2s-complement semantics.
1965	/// `r` may alias `a`.
1966	///
1967	/// Asserts `r` has enough storage to compute the result.
1968	/// The upper bound is `calcTwosCompLimbCount(a.len)`.
1969	pub fn truncate(r: *Mutable, a: Const, signedness: Signedness, bit_count: usize) void {
1970	// Handle 0-bit integers.
1971	if (bit_count == 0) {
1972	@branchHint(.unlikely);
1973	r.set(0);
1974	return;
1975	}
1976
1977	const max_limbs = calcTwosCompLimbCount(bit_count);
1978	const sign_bit = @as(Limb, 1) << @truncate(bit_count - 1);
1979	const mask = @as(Limb, maxInt(Limb)) >> @truncate(-%bit_count);
1980
1981	// Guess whether the result will have the same sign as `a`.
1982	// * If the result will be signed zero, the guess is `true`.
1983	// * If the result will be the minimum signed integer, the guess is `false`.
1984	// * If the result will be unsigned zero, the guess is `a.positive`.
1985	// * Otherwise the guess is correct.
1986	const same_sign_guess = switch (signedness) {
1987	.signed => max_limbs > a.limbs.len or a.limbs[max_limbs - 1] & sign_bit == 0,
1988	.unsigned => a.positive,
1989	};
1990
1991	const abs_trunc_a: Const = .{
1992	.positive = true,
1993	.limbs = a.limbs[0..llnormalize(a.limbs[0..@min(a.limbs.len, max_limbs)])],
1994	};
1995	if (same_sign_guess or abs_trunc_a.eqlZero()) {
1996	// One of the following is true:
1997	// * The result is zero.
1998	// * The result is non-zero and has the same sign as `a`.
1999	r.copy(abs_trunc_a);
2000	if (max_limbs <= r.len) r.limbs[max_limbs - 1] &= mask;
2001	r.normalize(r.len);
2002	r.positive = a.positive or r.eqlZero();
2003	} else {
2004	// One of the following is true:
2005	// * The result is the minimum signed integer.
2006	// * The result is unsigned zero.
2007	// * The result is non-zero and has the opposite sign as `a`.
2008	r.addScalar(abs_trunc_a, -1);
2009	llnot(r.limbs[0..r.len]);
2010	@memset(r.limbs[r.len..max_limbs], maxInt(Limb));
2011	r.limbs[max_limbs - 1] &= mask;
2012	r.normalize(max_limbs);
2013	r.positive = switch (signedness) {
2014	// The only value with the sign bit still set is the minimum signed integer.
2015	.signed => !a.positive and r.limbs[max_limbs - 1] & sign_bit == 0,
2016	.unsigned => !a.positive or r.eqlZero(),
2017	};
2018	}
2019	}
2020
2021	/// Saturate an integer to a number of bits, following 2s-complement semantics.
2022	/// r may alias a.
2023	///
2024	/// Asserts `r` has enough storage to store the result.
2025	/// The upper bound is `calcTwosCompLimbCount(a.len)`.
2026	pub fn saturate(r: *Mutable, a: Const, signedness: Signedness, bit_count: usize) void {
2027	if (!a.fitsInTwosComp(signedness, bit_count)) {
2028	r.setTwosCompIntLimit(if (r.positive) .max else .min, signedness, bit_count);
2029	}
2030	}
2031
2032	/// Read the value of `x` from `buffer`.
2033	/// Asserts that `buffer` is large enough to contain a value of bit-size `bit_count`.
2034	///
2035	/// The contents of `buffer` are interpreted as if they were the contents of
2036	/// @ptrCast(*[buffer.len]const u8, &x). Byte ordering is determined by `endian`
2037	/// and any required padding bits are expected on the MSB end.
2038	pub fn readTwosComplement(
2039	x: *Mutable,
2040	buffer: []const u8,
2041	bit_count: usize,
2042	endian: Endian,
2043	signedness: Signedness,
2044	) void {
2045	return readPackedTwosComplement(x, buffer, 0, bit_count, endian, signedness);
2046	}
2047
2048	/// Read the value of `x` from a packed memory `buffer`.
2049	/// Asserts that `buffer` is large enough to contain a value of bit-size `bit_count`
2050	/// at offset `bit_offset`.
2051	///
2052	/// This is equivalent to loading the value of an integer with `bit_count` bits as
2053	/// if it were a field in packed memory at the provided bit offset.
2054	pub fn readPackedTwosComplement(
2055	x: *Mutable,
2056	buffer: []const u8,
2057	bit_offset: usize,
2058	bit_count: usize,
2059	endian: Endian,
2060	signedness: Signedness,
2061	) void {
2062	if (bit_count == 0) {
2063	x.limbs[0] = 0;
2064	x.len = 1;
2065	x.positive = true;
2066	return;
2067	}
2068
2069	// Check whether the input is negative
2070	var positive = true;
2071	if (signedness == .signed) {
2072	const total_bits = bit_offset + bit_count;
2073	const last_byte = switch (endian) {
2074	.little => ((total_bits + 7) / 8) - 1,
2075	.big => buffer.len - ((total_bits + 7) / 8),
2076	};
2077
2078	const sign_bit = @as(u8, 1) << @as(u3, @intCast((total_bits - 1) % 8));
2079	positive = ((buffer[last_byte] & sign_bit) == 0);
2080	}
2081
2082	// Copy all complete limbs
2083	var carry: u1 = 1;
2084	var limb_index: usize = 0;
2085	var bit_index: usize = 0;
2086	while (limb_index < bit_count / @bitSizeOf(Limb)) : (limb_index += 1) {
2087	// Read one Limb of bits
2088	var limb = mem.readPackedInt(Limb, buffer, bit_index + bit_offset, endian);
2089	bit_index += @bitSizeOf(Limb);
2090
2091	// 2's complement (bitwise not, then add carry bit)
2092	if (!positive) {
2093	const ov = @addWithOverflow(~limb, carry);
2094	limb = ov[0];
2095	carry = ov[1];
2096	}
2097	x.limbs[limb_index] = limb;
2098	}
2099
2100	// Copy the remaining bits
2101	if (bit_count != bit_index) {
2102	// Read all remaining bits
2103	var limb = switch (signedness) {
2104	.unsigned => mem.readVarPackedInt(Limb, buffer, bit_index + bit_offset, bit_count - bit_index, endian, .unsigned),
2105	.signed => b: {
2106	const SLimb = @Int(.signed, @bitSizeOf(Limb));
2107	const limb = mem.readVarPackedInt(SLimb, buffer, bit_index + bit_offset, bit_count - bit_index, endian, .signed);
2108	break :b @as(Limb, @bitCast(limb));
2109	},
2110	};
2111
2112	// 2's complement (bitwise not, then add carry bit)
2113	if (!positive) {
2114	const ov = @addWithOverflow(~limb, carry);
2115	assert(ov[1] == 0);
2116	limb = ov[0];
2117	}
2118	x.limbs[limb_index] = limb;
2119
2120	limb_index += 1;
2121	}
2122
2123	x.positive = positive;
2124	x.len = limb_index;
2125	x.normalize(x.len);
2126	}
2127
2128	/// Normalize a possible sequence of leading zeros.
2129	///
2130	/// [1, 2, 3, 4, 0] -> [1, 2, 3, 4]
2131	/// [1, 2, 0, 0, 0] -> [1, 2]
2132	/// [0, 0, 0, 0, 0] -> [0]
2133	pub fn normalize(r: *Mutable, length: usize) void {
2134	r.len = llnormalize(r.limbs[0..length]);
2135	}
2136
2137	pub fn format(self: Mutable, w: *std.Io.Writer) std.Io.Writer.Error!void {
2138	return formatNumber(self, w, .{});
2139	}
2140
2141	/// If the absolute value of integer is greater than or equal to `pow(2, 64 * @sizeOf(usize) * 8)`,
2142	/// this function will fail to print the string, printing "(BigInt)" instead of a number.
2143	/// This is because the rendering algorithm requires reversing a string, which requires O(N) memory.
2144	/// See `Const.toString` and `Const.toStringAlloc` for a way to print big integers without failure.
2145	pub fn formatNumber(self: Mutable, w: *std.Io.Writer, n: std.fmt.Number) std.Io.Writer.Error!void {
2146	return self.toConst().formatNumber(w, n);
2147	}
2148	};
2149
2150	/// A arbitrary-precision big integer, with a fixed set of immutable limbs.
2151	pub const Const = struct {
2152	/// Raw digits. These are:
2153	///
2154	/// * Little-endian ordered
2155	/// * limbs.len >= 1
2156	/// * Zero is represented as limbs.len == 1 with limbs[0] == 0.
2157	///
2158	/// Accessing limbs directly should be avoided.
2159	limbs: []const Limb,
2160	positive: bool,
2161
2162	/// The result is an independent resource which is managed by the caller.
2163	pub fn toManaged(self: Const, allocator: Allocator) Allocator.Error!Managed {
2164	const limbs = try allocator.alloc(Limb, @max(Managed.default_capacity, self.limbs.len));
2165	@memcpy(limbs[0..self.limbs.len], self.limbs);
2166	return .{
2167	.allocator = allocator,
2168	.limbs = limbs,
2169	.metadata = if (self.positive)
2170	self.limbs.len & ~Managed.sign_bit
2171	else
2172	self.limbs.len | Managed.sign_bit,
2173	};
2174	}
2175
2176	/// Asserts `limbs` is big enough to store the value.
2177	pub fn toMutable(self: Const, limbs: []Limb) Mutable {
2178	@memcpy(limbs[0..self.limbs.len], self.limbs[0..self.limbs.len]);
2179	return .{
2180	.limbs = limbs,
2181	.positive = self.positive,
2182	.len = self.limbs.len,
2183	};
2184	}
2185
2186	pub fn dump(self: Const) void {
2187	for (self.limbs[0..self.limbs.len]) |limb| {
2188	std.debug.print("{x} ", .{limb});
2189	}
2190	std.debug.print("len={} positive={}\n", .{ self.limbs.len, self.positive });
2191	}
2192
2193	pub fn abs(self: Const) Const {
2194	return .{
2195	.limbs = self.limbs,
2196	.positive = true,
2197	};
2198	}
2199
2200	pub fn negate(self: Const) Const {
2201	return .{
2202	.limbs = self.limbs,
2203	.positive = !self.positive,
2204	};
2205	}
2206
2207	pub fn isOdd(self: Const) bool {
2208	return self.limbs[0] & 1 != 0;
2209	}
2210
2211	pub fn isEven(self: Const) bool {
2212	return !self.isOdd();
2213	}
2214
2215	/// Returns the number of bits required to represent the absolute value of an integer.
2216	pub fn bitCountAbs(self: Const) usize {
2217	return (self.limbs.len - 1) * limb_bits + (limb_bits - @clz(self.limbs[self.limbs.len - 1]));
2218	}
2219
2220	/// Returns the number of bits required to represent the integer in twos-complement form.
2221	///
2222	/// If the integer is negative the value returned is the number of bits needed by a signed
2223	/// integer to represent the value. If positive the value is the number of bits for an
2224	/// unsigned integer. Any unsigned integer will fit in the signed integer with bitcount
2225	/// one greater than the returned value.
2226	///
2227	/// e.g. -127 returns 8 as it will fit in an i8. 127 returns 7 since it fits in a u7.
2228	pub fn bitCountTwosComp(self: Const) usize {
2229	var bits = self.bitCountAbs();
2230
2231	// If the entire value has only one bit set (e.g. 0b100000000) then the negation in twos
2232	// complement requires one less bit.
2233	if (!self.positive) block: {
2234	bits += 1;
2235
2236	if (@popCount(self.limbs[self.limbs.len - 1]) == 1) {
2237	for (self.limbs[0 .. self.limbs.len - 1]) |limb| {
2238	if (@popCount(limb) != 0) {
2239	break :block;
2240	}
2241	}
2242
2243	bits -= 1;
2244	}
2245	}
2246
2247	return bits;
2248	}
2249
2250	/// Returns the number of bits required to represent the integer in twos-complement form
2251	/// with the given signedness.
2252	pub fn bitCountTwosCompForSignedness(self: Const, signedness: std.builtin.Signedness) usize {
2253	return self.bitCountTwosComp() + @intFromBool(self.positive and signedness == .signed);
2254	}
2255
2256	/// @popCount with two's complement semantics.
2257	///
2258	/// This returns the number of 1 bits set when the value would be represented in
2259	/// two's complement with the given integer width (bit_count).
2260	/// This includes the leading sign bit, which will be set for negative values.
2261	///
2262	/// Asserts that bit_count is enough to represent value in two's compliment
2263	/// and that the final result fits in a usize.
2264	/// Asserts that there are no trailing empty limbs on the most significant end,
2265	/// i.e. that limb count matches `calcLimbLen()` and zero is not negative.
2266	pub fn popCount(self: Const, bit_count: usize) usize {
2267	var sum: usize = 0;
2268	if (self.positive) {
2269	for (self.limbs) |limb| {
2270	sum += @popCount(limb);
2271	}
2272	} else {
2273	assert(self.fitsInTwosComp(.signed, bit_count));
2274	assert(self.limbs[self.limbs.len - 1] != 0);
2275
2276	var remaining_bits = bit_count;
2277	var carry: u1 = 1;
2278	var add_res: Limb = undefined;
2279
2280	// All but the most significant limb.
2281	for (self.limbs[0 .. self.limbs.len - 1]) |limb| {
2282	const ov = @addWithOverflow(~limb, carry);
2283	add_res = ov[0];
2284	carry = ov[1];
2285	sum += @popCount(add_res);
2286	remaining_bits -= limb_bits; // Asserted not to underflow by fitsInTwosComp
2287	}
2288
2289	// The most significant limb may have fewer than @bitSizeOf(Limb) meaningful bits,
2290	// which we can detect with @clz().
2291	// There may also be fewer limbs than needed to fill bit_count.
2292	const limb = self.limbs[self.limbs.len - 1];
2293	const leading_zeroes = @clz(limb);
2294	// The most significant limb is asserted not to be all 0s (above),
2295	// so ~limb cannot be all 1s, and ~limb + 1 cannot overflow.
2296	sum += @popCount(~limb + carry);
2297	sum -= leading_zeroes; // All leading zeroes were flipped and added to sum, so undo those
2298	const remaining_ones = remaining_bits - (limb_bits - leading_zeroes); // All bits not covered by limbs
2299	sum += remaining_ones;
2300	}
2301	return sum;
2302	}
2303
2304	pub fn fitsInTwosComp(self: Const, signedness: Signedness, bit_count: usize) bool {
2305	if (self.eqlZero()) {
2306	return true;
2307	}
2308	if (signedness == .unsigned and !self.positive) {
2309	return false;
2310	}
2311	return bit_count >= self.bitCountTwosCompForSignedness(signedness);
2312	}
2313
2314	/// Returns whether self can fit into an integer of the requested type.
2315	pub fn fits(self: Const, comptime T: type) bool {
2316	const info = @typeInfo(T).int;
2317	return self.fitsInTwosComp(info.signedness, info.bits);
2318	}
2319
2320	/// Returns the approximate size of the integer in the given base. Negative values accommodate for
2321	/// the minus sign. This is used for determining the number of characters needed to print the
2322	/// value. It is inexact and may exceed the given value by ~1-2 bytes.
2323	/// TODO See if we can make this exact.
2324	pub fn sizeInBaseUpperBound(self: Const, base: usize) usize {
2325	const bit_count = @as(usize, @intFromBool(!self.positive)) + self.bitCountAbs();
2326	return (bit_count / math.log2(base)) + 2;
2327	}
2328
2329	pub const ConvertError = error{
2330	NegativeIntoUnsigned,
2331	TargetTooSmall,
2332	};
2333
2334	/// Convert `self` to `Int`.
2335	///
2336	/// Returns an error if self cannot be narrowed into the requested type without truncation.
2337	pub fn toInt(self: Const, comptime Int: type) ConvertError!Int {
2338	switch (@typeInfo(Int)) {
2339	.int => |info| {
2340	// Make sure -0 is handled correctly.
2341	if (self.eqlZero()) return 0;
2342
2343	const Unsigned = @Int(.unsigned, info.bits);
2344
2345	if (!self.fitsInTwosComp(info.signedness, info.bits)) {
2346	return error.TargetTooSmall;
2347	}
2348
2349	var r: Unsigned = 0;
2350
2351	if (@sizeOf(Unsigned) <= @sizeOf(Limb)) {
2352	r = @intCast(self.limbs[0]);
2353	} else {
2354	for (self.limbs[0..self.limbs.len], 0..) |_, ri| {
2355	const limb = self.limbs[self.limbs.len - ri - 1];
2356	r <<= limb_bits;
2357	r |= limb;
2358	}
2359	}
2360
2361	if (info.signedness == .unsigned) {
2362	return if (self.positive) @intCast(r) else error.NegativeIntoUnsigned;
2363	} else {
2364	if (self.positive) {
2365	return @intCast(r);
2366	} else {
2367	if (math.cast(Int, r)) |ok| {
2368	return -ok;
2369	} else {
2370	return minInt(Int);
2371	}
2372	}
2373	}
2374	},
2375	else => @compileError("expected int type, found '" ++ @typeName(Int) ++ "'"),
2376	}
2377	}
2378
2379	/// Convert self to `Float`.
2380	pub fn toFloat(self: Const, comptime Float: type, round: Round) struct { Float, Exactness } {
2381	if (Float == comptime_float) return self.toFloat(f128, round);
2382	const normalized_abs: Const = .{
2383	.limbs = self.limbs[0..llnormalize(self.limbs)],
2384	.positive = true,
2385	};
2386	if (normalized_abs.eqlZero()) return .{ if (self.positive) 0.0 else -0.0, .exact };
2387
2388	const Repr = std.math.FloatRepr(Float);
2389	var mantissa_limbs: [calcNonZeroTwosCompLimbCount(1 + @bitSizeOf(Repr.Mantissa))]Limb = undefined;
2390	var mantissa: Mutable = .{
2391	.limbs = &mantissa_limbs,
2392	.positive = undefined,
2393	.len = undefined,
2394	};
2395	var exponent = normalized_abs.bitCountAbs() - 1;
2396	const exactness: Exactness = exactness: {
2397	if (exponent <= @bitSizeOf(Repr.Normalized.Fraction)) {
2398	mantissa.shiftLeft(normalized_abs, @intCast(@bitSizeOf(Repr.Normalized.Fraction) - exponent));
2399	break :exactness .exact;
2400	}
2401	const shift: usize = @intCast(exponent - @bitSizeOf(Repr.Normalized.Fraction));
2402	mantissa.shiftRight(normalized_abs, shift);
2403	const final_limb_index = (shift - 1) / limb_bits;
2404	const round_bits = normalized_abs.limbs[final_limb_index] << @truncate(-%shift) |
2405	@intFromBool(!std.mem.allEqual(Limb, normalized_abs.limbs[0..final_limb_index], 0));
2406	if (round_bits == 0) break :exactness .exact;
2407	round: switch (round) {
2408	.nearest_even => {
2409	const half: Limb = 1 << (limb_bits - 1);
2410	if (round_bits >= half) mantissa.addScalar(mantissa.toConst(), 1);
2411	if (round_bits == half) mantissa.limbs[0] &= ~@as(Limb, 1);
2412	},
2413	.away => mantissa.addScalar(mantissa.toConst(), 1),
2414	.trunc => {},
2415	.floor => if (!self.positive) continue :round .away,
2416	.ceil => if (self.positive) continue :round .away,
2417	}
2418	break :exactness .inexact;
2419	};
2420	const normalized_res: Repr.Normalized = .{
2421	.fraction = @truncate(mantissa.toInt(Repr.Mantissa) catch |err| switch (err) {
2422	error.NegativeIntoUnsigned => unreachable,
2423	error.TargetTooSmall => fraction: {
2424	assert(mantissa.toConst().orderAgainstScalar(1 << @bitSizeOf(Repr.Mantissa)).compare(.eq));
2425	exponent += 1;
2426	break :fraction 1 << (@bitSizeOf(Repr.Mantissa) - 1);
2427	},
2428	}),
2429	.exponent = std.math.lossyCast(Repr.Normalized.Exponent, exponent),
2430	};
2431	return .{ normalized_res.reconstruct(if (self.positive) .positive else .negative), exactness };
2432	}
2433
2434	pub fn format(self: Const, w: *std.Io.Writer) std.Io.Writer.Error!void {
2435	return self.formatNumber(w, .{});
2436	}
2437
2438	/// If the absolute value of integer is greater than or equal to `pow(2, 64 * @sizeOf(usize) * 8)`,
2439	/// this function will fail to print the string, printing "(BigInt)" instead of a number.
2440	/// This is because the rendering algorithm requires reversing a string, which requires O(N) memory.
2441	/// See `toString` and `toStringAlloc` for a way to print big integers without failure.
2442	pub fn formatNumber(self: Const, w: *std.Io.Writer, number: std.fmt.Number) std.Io.Writer.Error!void {
2443	const available_len = 64;
2444	if (self.limbs.len > available_len)
2445	return w.writeAll("(BigInt)");
2446
2447	var limbs: [calcToStringLimbsBufferLen(available_len, 10)]Limb = undefined;
2448
2449	const biggest: Const = .{
2450	.limbs = &@as([available_len]Limb, @splat(comptime math.maxInt(Limb))),
2451	.positive = false,
2452	};
2453	var buf: [biggest.sizeInBaseUpperBound(2)]u8 = undefined;
2454	const base: u8 = number.mode.base() orelse @panic("TODO print big int in scientific form");
2455	const len = self.toString(&buf, base, number.case, &limbs);
2456	return w.writeAll(buf[0..len]);
2457	}
2458
2459	/// Converts self to a string in the requested base.
2460	/// Caller owns returned memory.
2461	/// Asserts that `base` is in the range [2, 36].
2462	/// See also `toString`, a lower level function than this.
2463	pub fn toStringAlloc(self: Const, allocator: Allocator, base: u8, case: std.fmt.Case) Allocator.Error![]u8 {
2464	assert(base >= 2);
2465	assert(base <= 36);
2466
2467	if (self.eqlZero()) {
2468	return allocator.dupe(u8, "0");
2469	}
2470	const string = try allocator.alloc(u8, self.sizeInBaseUpperBound(base));
2471	errdefer allocator.free(string);
2472
2473	const limbs = try allocator.alloc(Limb, calcToStringLimbsBufferLen(self.limbs.len, base));
2474	defer allocator.free(limbs);
2475
2476	return allocator.realloc(string, self.toString(string, base, case, limbs));
2477	}
2478
2479	/// Converts self to a string in the requested base.
2480	/// Asserts that `base` is in the range [2, 36].
2481	/// `string` is a caller-provided slice of at least `sizeInBaseUpperBound` bytes,
2482	/// where the result is written to.
2483	/// Returns the length of the string.
2484	/// `limbs_buffer` is caller-provided memory for `toString` to use as a working area. It must have
2485	/// length of at least `calcToStringLimbsBufferLen`.
2486	/// In the case of power-of-two base, `limbs_buffer` is ignored.
2487	/// See also `toStringAlloc`, a higher level function than this.
2488	pub fn toString(self: Const, string: []u8, base: u8, case: std.fmt.Case, limbs_buffer: []Limb) usize {
2489	assert(base >= 2);
2490	assert(base <= 36);
2491
2492	if (self.eqlZero()) {
2493	string[0] = '0';
2494	return 1;
2495	}
2496
2497	var digits_len: usize = 0;
2498
2499	// Power of two: can do a single pass and use masks to extract digits.
2500	if (math.isPowerOfTwo(base)) {
2501	const base_shift = math.log2_int(Limb, base);
2502
2503	outer: for (self.limbs[0..self.limbs.len]) |limb| {
2504	var shift: usize = 0;
2505	while (shift < limb_bits) : (shift += base_shift) {
2506	const r = @as(u8, @intCast((limb >> @as(Log2Limb, @intCast(shift))) & @as(Limb, base - 1)));
2507	const ch = std.fmt.digitToChar(r, case);
2508	string[digits_len] = ch;
2509	digits_len += 1;
2510	// If we hit the end, it must be all zeroes from here.
2511	if (digits_len == string.len) break :outer;
2512	}
2513	}
2514
2515	// Always will have a non-zero digit somewhere.
2516	while (string[digits_len - 1] == '0') {
2517	digits_len -= 1;
2518	}
2519	} else {
2520	// Non power-of-two: batch divisions per word size.
2521	// We use a HalfLimb here so the division uses the faster lldiv0p5 over lldiv1 codepath.
2522	const digits_per_limb = math.log(HalfLimb, base, maxInt(HalfLimb));
2523	var limb_base: Limb = 1;
2524	var j: usize = 0;
2525	while (j < digits_per_limb) : (j += 1) {
2526	limb_base *= base;
2527	}
2528	const b: Const = .{ .limbs = &[_]Limb{limb_base}, .positive = true };
2529
2530	var q: Mutable = .{
2531	.limbs = limbs_buffer[0 .. self.limbs.len + 2],
2532	.positive = true, // Make absolute by ignoring self.positive.
2533	.len = self.limbs.len,
2534	};
2535	@memcpy(q.limbs[0..self.limbs.len], self.limbs);
2536
2537	var r: Mutable = .{
2538	.limbs = limbs_buffer[q.limbs.len..][0..self.limbs.len],
2539	.positive = true,
2540	.len = 1,
2541	};
2542	r.limbs[0] = 0;
2543
2544	const rest_of_the_limbs_buf = limbs_buffer[q.limbs.len + r.limbs.len ..];
2545
2546	while (q.len >= 2) {
2547	// Passing an allocator here would not be helpful since this division is destroying
2548	// information, not creating it. [TODO citation needed]
2549	q.divTrunc(&r, q.toConst(), b, rest_of_the_limbs_buf);
2550
2551	var r_word = r.limbs[0];
2552	var i: usize = 0;
2553	while (i < digits_per_limb) : (i += 1) {
2554	const ch = std.fmt.digitToChar(@as(u8, @intCast(r_word % base)), case);
2555	r_word /= base;
2556	string[digits_len] = ch;
2557	digits_len += 1;
2558	}
2559	}
2560
2561	{
2562	assert(q.len == 1);
2563
2564	var r_word = q.limbs[0];
2565	while (r_word != 0) {
2566	const ch = std.fmt.digitToChar(@as(u8, @intCast(r_word % base)), case);
2567	r_word /= base;
2568	string[digits_len] = ch;
2569	digits_len += 1;
2570	}
2571	}
2572	}
2573
2574	if (!self.positive) {
2575	string[digits_len] = '-';
2576	digits_len += 1;
2577	}
2578
2579	const s = string[0..digits_len];
2580	mem.reverse(u8, s);
2581	return s.len;
2582	}
2583
2584	/// Write the value of `x` into `buffer`
2585	/// Asserts that `buffer` is large enough to store the value.
2586	///
2587	/// `buffer` is filled so that its contents match what would be observed via
2588	/// @ptrCast(*[buffer.len]const u8, &x). Byte ordering is determined by `endian`,
2589	/// and any required padding bits are added on the MSB end.
2590	pub fn writeTwosComplement(x: Const, buffer: []u8, endian: Endian) void {
2591	return writePackedTwosComplement(x, buffer, 0, 8 * buffer.len, endian);
2592	}
2593
2594	/// Write the value of `x` to a packed memory `buffer`.
2595	/// Asserts that `buffer` is large enough to contain a value of bit-size `bit_count`
2596	/// at offset `bit_offset`.
2597	///
2598	/// This is equivalent to storing the value of an integer with `bit_count` bits as
2599	/// if it were a field in packed memory at the provided bit offset.
2600	pub fn writePackedTwosComplement(x: Const, buffer: []u8, bit_offset: usize, bit_count: usize, endian: Endian) void {
2601	assert(x.fitsInTwosComp(if (x.positive) .unsigned else .signed, bit_count));
2602
2603	// Copy all complete limbs
2604	var carry: u1 = 1;
2605	var limb_index: usize = 0;
2606	var bit_index: usize = 0;
2607	while (limb_index < bit_count / @bitSizeOf(Limb)) : (limb_index += 1) {
2608	var limb: Limb = if (limb_index < x.limbs.len) x.limbs[limb_index] else 0;
2609
2610	// 2's complement (bitwise not, then add carry bit)
2611	if (!x.positive) {
2612	const ov = @addWithOverflow(~limb, carry);
2613	limb = ov[0];
2614	carry = ov[1];
2615	}
2616
2617	// Write one Limb of bits
2618	mem.writePackedInt(Limb, buffer, bit_index + bit_offset, limb, endian);
2619	bit_index += @bitSizeOf(Limb);
2620	}
2621
2622	// Copy the remaining bits
2623	if (bit_count != bit_index) {
2624	var limb: Limb = if (limb_index < x.limbs.len) x.limbs[limb_index] else 0;
2625
2626	// 2's complement (bitwise not, then add carry bit)
2627	if (!x.positive) limb = ~limb +% carry;
2628
2629	// Write all remaining bits
2630	mem.writeVarPackedInt(buffer, bit_index + bit_offset, bit_count - bit_index, limb, endian);
2631	}
2632	}
2633
2634	/// Returns `math.Order.lt`, `math.Order.eq`, `math.Order.gt` if
2635	/// `|a| < |b|`, `|a| == |b|`, or `|a| > |b|` respectively.
2636	pub fn orderAbs(a: Const, b: Const) math.Order {
2637	if (a.limbs.len < b.limbs.len) {
2638	return .lt;
2639	}
2640	if (a.limbs.len > b.limbs.len) {
2641	return .gt;
2642	}
2643
2644	var i: usize = a.limbs.len - 1;
2645	while (i != 0) : (i -= 1) {
2646	if (a.limbs[i] != b.limbs[i]) {
2647	break;
2648	}
2649	}
2650
2651	if (a.limbs[i] < b.limbs[i]) {
2652	return .lt;
2653	} else if (a.limbs[i] > b.limbs[i]) {
2654	return .gt;
2655	} else {
2656	return .eq;
2657	}
2658	}
2659
2660	/// Returns `math.Order.lt`, `math.Order.eq`, `math.Order.gt` if `a < b`, `a == b` or `a > b` respectively.
2661	pub fn order(a: Const, b: Const) math.Order {
2662	if (a.positive != b.positive) {
2663	if (eqlZero(a) and eqlZero(b)) {
2664	return .eq;
2665	} else {
2666	return if (a.positive) .gt else .lt;
2667	}
2668	} else {
2669	const r = orderAbs(a, b);
2670	return if (a.positive) r else switch (r) {
2671	.lt => math.Order.gt,
2672	.eq => math.Order.eq,
2673	.gt => math.Order.lt,
2674	};
2675	}
2676	}
2677
2678	/// Same as `order` but the right-hand operand is a primitive integer.
2679	pub fn orderAgainstScalar(lhs: Const, scalar: anytype) math.Order {
2680	// Normally we could just determine the number of limbs needed with calcLimbLen,
2681	// but that is not comptime-known when scalar is not a comptime_int. Instead, we
2682	// use calcTwosCompLimbCount for a non-comptime_int scalar, which can be pessimistic
2683	// in the case that scalar happens to be small in magnitude within its type, but it
2684	// is well worth being able to use the stack and not needing an allocator passed in.
2685	// Note that Mutable.init still sets len to calcLimbLen(scalar) in any case.
2686	const limbs_len = comptime switch (@typeInfo(@TypeOf(scalar))) {
2687	.comptime_int => calcLimbLen(scalar),
2688	.int => |info| calcTwosCompLimbCount(info.bits),
2689	else => @compileError("expected scalar to be an int"),
2690	};
2691	var limbs: [limbs_len]Limb = undefined;
2692	const rhs = Mutable.init(&limbs, scalar);
2693	return order(lhs, rhs.toConst());
2694	}
2695
2696	/// Returns true if `a == 0`.
2697	pub fn eqlZero(a: Const) bool {
2698	var d: Limb = 0;
2699	for (a.limbs) |limb| d |= limb;
2700	return d == 0;
2701	}
2702
2703	/// Returns true if `|a| == |b|`.
2704	pub fn eqlAbs(a: Const, b: Const) bool {
2705	return orderAbs(a, b) == .eq;
2706	}
2707
2708	/// Returns true if `a == b`.
2709	pub fn eql(a: Const, b: Const) bool {
2710	return order(a, b) == .eq;
2711	}
2712
2713	/// Returns the number of leading zeros in twos-complement form.
2714	pub fn clz(a: Const, bits: Limb) Limb {
2715	// Limbs are stored in little-endian order but we need to iterate big-endian.
2716	if (!a.positive and !a.eqlZero()) return 0;
2717	var total_limb_lz: Limb = 0;
2718	var i: usize = a.limbs.len;
2719	const bits_per_limb = @bitSizeOf(Limb);
2720	while (i != 0) {
2721	i -= 1;
2722	const this_limb_lz = @clz(a.limbs[i]);
2723	total_limb_lz += this_limb_lz;
2724	if (this_limb_lz != bits_per_limb) break;
2725	}
2726	const total_limb_bits = a.limbs.len * bits_per_limb;
2727	return total_limb_lz + bits - total_limb_bits;
2728	}
2729
2730	/// Returns the number of trailing zeros in twos-complement form.
2731	pub fn ctz(a: Const, bits: Limb) Limb {
2732	// Limbs are stored in little-endian order. Converting a negative number to twos-complement
2733	// flips all bits above the lowest set bit, which does not affect the trailing zero count.
2734	if (a.eqlZero()) return bits;
2735	var result: Limb = 0;
2736	for (a.limbs) |limb| {
2737	const limb_tz = @ctz(limb);
2738	result += limb_tz;
2739	if (limb_tz != @bitSizeOf(Limb)) break;
2740	}
2741	return @min(result, bits);
2742	}
2743
2744	/// Calculate the base 2 logarithm, rounded down.
2745	pub fn log2(a: Const) Limb {
2746	assert(a.positive);
2747	assert(!a.eqlZero());
2748	return a.bitCountAbs() - 1;
2749	}
2750
2751	/// Calculate the base 10 logarithm, rounded down.
2752	///
2753	/// The allocator is used to allocate a temporary buffer.
2754	pub fn log10Alloc(a: Const, allocator: Allocator) Allocator.Error!Limb {
2755	const limbs_buffer = try allocator.alloc(Limb, calcLog10LimbsBufferLen(a.limbs.len));
2756	defer allocator.free(limbs_buffer);
2757
2758	return a.log10(limbs_buffer);
2759	}
2760
2761	/// Calculate the base 10 logarithm, rounded down.
2762	///
2763	/// `limbs_buffer` is used for temporary storage. The amount required is given by `calcLog10LimbsBufferLen`.
2764	pub fn log10(a: Const, limbs_buffer: []Limb) Limb {
2765	assert(a.positive);
2766	assert(!a.eqlZero());
2767	const limb_base_as_bigint: Const = .{ .limbs = &.{constants.big_bases[10]}, .positive = true };
2768
2769	var q: Mutable = .{
2770	.limbs = limbs_buffer[0 .. a.limbs.len + 2],
2771	.positive = true,
2772	.len = a.limbs.len,
2773	};
2774	@memcpy(q.limbs[0..a.limbs.len], a.limbs);
2775
2776	var remainder: Mutable = .{
2777	.limbs = limbs_buffer[q.limbs.len..][0..a.limbs.len],
2778	.positive = true,
2779	.len = 1,
2780	};
2781
2782	const division_buf = limbs_buffer[q.limbs.len + remainder.limbs.len ..];
2783
2784	var num_digits: Limb = 0;
2785	while (q.len >= 2) {
2786	q.divTrunc(&remainder, q.toConst(), limb_base_as_bigint, division_buf);
2787	num_digits += constants.digits_per_limb[10];
2788	}
2789	var remaining_limb = q.limbs[0];
2790	while (remaining_limb != 0) {
2791	remaining_limb /= 10;
2792	num_digits += 1;
2793	}
2794
2795	return num_digits - 1;
2796	}
2797	};
2798
2799	/// An arbitrary-precision big integer along with an allocator which manages the memory.
2800	///
2801	/// Memory is allocated as needed to ensure operations never overflow. The range
2802	/// is bounded only by available memory.
2803	pub const Managed = struct {
2804	pub const sign_bit: usize = 1 << (@typeInfo(usize).int.bits - 1);
2805
2806	/// Default number of limbs to allocate on creation of a `Managed`.
2807	pub const default_capacity = 4;
2808
2809	/// Allocator used by the Managed when requesting memory.
2810	allocator: Allocator,
2811
2812	/// Raw digits. These are:
2813	///
2814	/// * Little-endian ordered
2815	/// * limbs.len >= 1
2816	/// * Zero is represent as Managed.len() == 1 with limbs[0] == 0.
2817	///
2818	/// Accessing limbs directly should be avoided.
2819	limbs: []Limb,
2820
2821	/// High bit is the sign bit. If set, Managed is negative, else Managed is positive.
2822	/// The remaining bits represent the number of limbs used by Managed.
2823	metadata: usize,
2824
2825	/// Creates a new `Managed`. `default_capacity` limbs will be allocated immediately.
2826	/// The integer value after initializing is `0`.
2827	pub fn init(allocator: Allocator) !Managed {
2828	return initCapacity(allocator, default_capacity);
2829	}
2830
2831	pub fn toMutable(self: Managed) Mutable {
2832	return .{
2833	.limbs = self.limbs,
2834	.positive = self.isPositive(),
2835	.len = self.len(),
2836	};
2837	}
2838
2839	pub fn toConst(self: Managed) Const {
2840	return .{
2841	.limbs = self.limbs[0..self.len()],
2842	.positive = self.isPositive(),
2843	};
2844	}
2845
2846	/// Creates a new `Managed` with value `value`.
2847	///
2848	/// This is identical to an `init`, followed by a `set`.
2849	pub fn initSet(allocator: Allocator, value: anytype) !Managed {
2850	var s = try Managed.init(allocator);
2851	errdefer s.deinit();
2852	try s.set(value);
2853	return s;
2854	}
2855
2856	/// Creates a new Managed with a specific capacity. If capacity < default_capacity then the
2857	/// default capacity will be used instead.
2858	/// The integer value after initializing is `0`.
2859	pub fn initCapacity(allocator: Allocator, capacity: usize) !Managed {
2860	return .{
2861	.allocator = allocator,
2862	.metadata = 1,
2863	.limbs = block: {
2864	const limbs = try allocator.alloc(Limb, @max(default_capacity, capacity));
2865	limbs[0] = 0;
2866	break :block limbs;
2867	},
2868	};
2869	}
2870
2871	/// Returns the number of limbs currently in use.
2872	pub fn len(self: Managed) usize {
2873	return self.metadata & ~sign_bit;
2874	}
2875
2876	/// Returns whether an Managed is positive.
2877	pub fn isPositive(self: Managed) bool {
2878	return self.metadata & sign_bit == 0;
2879	}
2880
2881	/// Sets the sign of an Managed.
2882	pub fn setSign(self: *Managed, positive: bool) void {
2883	if (positive) {
2884	self.metadata &= ~sign_bit;
2885	} else {
2886	self.metadata |= sign_bit;
2887	}
2888	}
2889
2890	/// Sets the length of an Managed.
2891	///
2892	/// If setLen is used, then the Managed must be normalized to suit.
2893	pub fn setLen(self: *Managed, new_len: usize) void {
2894	self.metadata &= sign_bit;
2895	self.metadata |= new_len;
2896	}
2897
2898	pub fn setMetadata(self: *Managed, positive: bool, length: usize) void {
2899	self.metadata = if (positive) length & ~sign_bit else length | sign_bit;
2900	}
2901
2902	/// Ensures an Managed has enough space allocated for capacity limbs. If the Managed does not have
2903	/// sufficient capacity, the exact amount will be allocated. This occurs even if the requested
2904	/// capacity is only greater than the current capacity by one limb.
2905	pub fn ensureCapacity(self: *Managed, capacity: usize) !void {
2906	if (capacity <= self.limbs.len) {
2907	return;
2908	}
2909	self.limbs = try self.allocator.realloc(self.limbs, capacity);
2910	}
2911
2912	/// Frees all associated memory.
2913	pub fn deinit(self: *Managed) void {
2914	self.allocator.free(self.limbs);
2915	self.* = undefined;
2916	}
2917
2918	/// Returns a `Managed` with the same value. The returned `Managed` is a deep copy and
2919	/// can be modified separately from the original, and its resources are managed
2920	/// separately from the original.
2921	pub fn clone(other: Managed) !Managed {
2922	return other.cloneWithDifferentAllocator(other.allocator);
2923	}
2924
2925	pub fn cloneWithDifferentAllocator(other: Managed, allocator: Allocator) !Managed {
2926	return .{
2927	.allocator = allocator,
2928	.metadata = other.metadata,
2929	.limbs = block: {
2930	const limbs = try allocator.alloc(Limb, other.len());
2931	@memcpy(limbs, other.limbs[0..other.len()]);
2932	break :block limbs;
2933	},
2934	};
2935	}
2936
2937	/// Copies the value of the integer to an existing `Managed` so that they both have the same value.
2938	/// Extra memory will be allocated if the receiver does not have enough capacity.
2939	pub fn copy(self: *Managed, other: Const) !void {
2940	if (self.limbs.ptr == other.limbs.ptr) return;
2941
2942	try self.ensureCapacity(other.limbs.len);
2943	@memcpy(self.limbs[0..other.limbs.len], other.limbs[0..other.limbs.len]);
2944	self.setMetadata(other.positive, other.limbs.len);
2945	}
2946
2947	/// Efficiently swap a `Managed` with another. This swaps the limb pointers and a full copy is not
2948	/// performed. The address of the limbs field will not be the same after this function.
2949	pub fn swap(self: *Managed, other: *Managed) void {
2950	mem.swap(Managed, self, other);
2951	}
2952
2953	/// Debugging tool: prints the state to stderr.
2954	pub fn dump(self: Managed) void {
2955	for (self.limbs[0..self.len()]) |limb| {
2956	std.debug.print("{x} ", .{limb});
2957	}
2958	std.debug.print("len={} capacity={} positive={}\n", .{ self.len(), self.limbs.len, self.isPositive() });
2959	}
2960
2961	/// Negate the sign.
2962	pub fn negate(self: *Managed) void {
2963	self.metadata ^= sign_bit;
2964	}
2965
2966	/// Make positive.
2967	pub fn abs(self: *Managed) void {
2968	self.metadata &= ~sign_bit;
2969	}
2970
2971	pub fn isOdd(self: Managed) bool {
2972	return self.limbs[0] & 1 != 0;
2973	}
2974
2975	pub fn isEven(self: Managed) bool {
2976	return !self.isOdd();
2977	}
2978
2979	/// Returns the number of bits required to represent the absolute value of an integer.
2980	pub fn bitCountAbs(self: Managed) usize {
2981	return self.toConst().bitCountAbs();
2982	}
2983
2984	/// Returns the number of bits required to represent the integer in twos-complement form.
2985	///
2986	/// If the integer is negative the value returned is the number of bits needed by a signed
2987	/// integer to represent the value. If positive the value is the number of bits for an
2988	/// unsigned integer. Any unsigned integer will fit in the signed integer with bitcount
2989	/// one greater than the returned value.
2990	///
2991	/// e.g. -127 returns 8 as it will fit in an i8. 127 returns 7 since it fits in a u7.
2992	pub fn bitCountTwosComp(self: Managed) usize {
2993	return self.toConst().bitCountTwosComp();
2994	}
2995
2996	pub fn fitsInTwosComp(self: Managed, signedness: Signedness, bit_count: usize) bool {
2997	return self.toConst().fitsInTwosComp(signedness, bit_count);
2998	}
2999
3000	/// Returns whether self can fit into an integer of the requested type.
3001	pub fn fits(self: Managed, comptime T: type) bool {
3002	return self.toConst().fits(T);
3003	}
3004
3005	/// Returns the approximate size of the integer in the given base. Negative values accommodate for
3006	/// the minus sign. This is used for determining the number of characters needed to print the
3007	/// value. It is inexact and may exceed the given value by ~1-2 bytes.
3008	pub fn sizeInBaseUpperBound(self: Managed, base: usize) usize {
3009	return self.toConst().sizeInBaseUpperBound(base);
3010	}
3011
3012	/// Sets an Managed to value. Value must be an primitive integer type.
3013	pub fn set(self: *Managed, value: anytype) Allocator.Error!void {
3014	try self.ensureCapacity(calcLimbLen(value));
3015	var m = self.toMutable();
3016	m.set(value);
3017	self.setMetadata(m.positive, m.len);
3018	}
3019
3020	pub const ConvertError = Const.ConvertError;
3021
3022	/// Convert `self` to `Int`.
3023	///
3024	/// Returns an error if self cannot be narrowed into the requested type without truncation.
3025	pub fn toInt(self: Managed, comptime Int: type) ConvertError!Int {
3026	return self.toConst().toInt(Int);
3027	}
3028
3029	/// Convert `self` to `Float`.
3030	pub fn toFloat(self: Managed, comptime Float: type, round: Round) struct { Float, Exactness } {
3031	return self.toConst().toFloat(Float, round);
3032	}
3033
3034	/// Set self from the string representation `value`.
3035	///
3036	/// `value` must contain only digits <= `base` and is case insensitive. Base prefixes are
3037	/// not allowed (e.g. 0x43 should simply be 43). Underscores in the input string are
3038	/// ignored and can be used as digit separators.
3039	///
3040	/// Returns an error if memory could not be allocated or `value` has invalid digits for the
3041	/// requested base.
3042	///
3043	/// self's allocator is used for temporary storage to boost multiplication performance.
3044	pub fn setString(self: *Managed, base: u8, value: []const u8) !void {
3045	if (base < 2 or base > 36) return error.InvalidBase;
3046	try self.ensureCapacity(calcSetStringLimbCount(base, value.len));
3047	var m = self.toMutable();
3048	try m.setString(base, value);
3049	self.setMetadata(m.positive, m.len);
3050	}
3051
3052	/// Set self to either bound of a 2s-complement integer.
3053	/// Note: The result is still sign-magnitude, not twos complement! In order to convert the
3054	/// result to twos complement, it is sufficient to take the absolute value.
3055	pub fn setTwosCompIntLimit(
3056	r: *Managed,
3057	limit: TwosCompIntLimit,
3058	signedness: Signedness,
3059	bit_count: usize,
3060	) !void {
3061	try r.ensureCapacity(calcTwosCompLimbCount(bit_count));
3062	var m = r.toMutable();
3063	m.setTwosCompIntLimit(limit, signedness, bit_count);
3064	r.setMetadata(m.positive, m.len);
3065	}
3066
3067	/// Converts self to a string in the requested base. Memory is allocated from the provided
3068	/// allocator and not the one present in self.
3069	pub fn toString(self: Managed, allocator: Allocator, base: u8, case: std.fmt.Case) ![]u8 {
3070	if (base < 2 or base > 36) return error.InvalidBase;
3071	return self.toConst().toStringAlloc(allocator, base, case);
3072	}
3073
3074	/// To allow `std.fmt.format` to work with `Managed`.
3075	pub fn format(self: Managed, w: *std.Io.Writer) std.Io.Writer.Error!void {
3076	return formatNumber(self, w, .{});
3077	}
3078
3079	/// If the absolute value of integer is greater than or equal to `pow(2, 64 * @sizeOf(usize) * 8)`,
3080	/// this function will fail to print the string, printing "(BigInt)" instead of a number.
3081	/// This is because the rendering algorithm requires reversing a string, which requires O(N) memory.
3082	/// See `toString` and `toStringAlloc` for a way to print big integers without failure.
3083	pub fn formatNumber(self: Managed, w: *std.Io.Writer, n: std.fmt.Number) std.Io.Writer.Error!void {
3084	return self.toConst().formatNumber(w, n);
3085	}
3086
3087	/// Returns math.Order.lt, math.Order.eq, math.Order.gt if |a| < |b|, |a| ==
3088	/// |b| or |a| > |b| respectively.
3089	pub fn orderAbs(a: Managed, b: Managed) math.Order {
3090	return a.toConst().orderAbs(b.toConst());
3091	}
3092
3093	/// Returns math.Order.lt, math.Order.eq, math.Order.gt if a < b, a == b or a > b
3094	/// respectively.
3095	pub fn order(a: Managed, b: Managed) math.Order {
3096	return a.toConst().order(b.toConst());
3097	}
3098
3099	/// Returns true if a == 0.
3100	pub fn eqlZero(a: Managed) bool {
3101	return a.toConst().eqlZero();
3102	}
3103
3104	/// Returns true if |a| == |b|.
3105	pub fn eqlAbs(a: Managed, b: Managed) bool {
3106	return a.toConst().eqlAbs(b.toConst());
3107	}
3108
3109	/// Returns true if a == b.
3110	pub fn eql(a: Managed, b: Managed) bool {
3111	return a.toConst().eql(b.toConst());
3112	}
3113
3114	/// Normalize a possible sequence of leading zeros.
3115	///
3116	/// [1, 2, 3, 4, 0] -> [1, 2, 3, 4]
3117	/// [1, 2, 0, 0, 0] -> [1, 2]
3118	/// [0, 0, 0, 0, 0] -> [0]
3119	pub fn normalize(r: *Managed, length: usize) void {
3120	assert(length > 0);
3121	assert(length <= r.limbs.len);
3122
3123	var j = length;
3124	while (j > 0) : (j -= 1) {
3125	if (r.limbs[j - 1] != 0) {
3126	break;
3127	}
3128	}
3129
3130	// Handle zero
3131	r.setLen(if (j != 0) j else 1);
3132	}
3133
3134	/// r = a + scalar
3135	///
3136	/// r and a may be aliases.
3137	///
3138	/// Returns an error if memory could not be allocated.
3139	pub fn addScalar(r: *Managed, a: *const Managed, scalar: anytype) Allocator.Error!void {
3140	const needed = @max(a.len(), calcLimbLen(scalar)) + 1;
3141	const aliased = limbsAliasDistinct(r, a);
3142	try r.ensureAliasAwareCapacity(needed, aliased);
3143	var m = r.toMutable();
3144	m.addScalar(a.toConst(), scalar);
3145	r.setMetadata(m.positive, m.len);
3146	}
3147
3148	/// r = a + b
3149	///
3150	/// r, a and b may be aliases.
3151	///
3152	/// Returns an error if memory could not be allocated.
3153	pub fn add(r: *Managed, a: *const Managed, b: *const Managed) Allocator.Error!void {
3154	const needed = @max(a.len(), b.len()) + 1;
3155	const aliased = limbsAliasDistinct(r, a) or limbsAliasDistinct(r, b);
3156	try r.ensureAliasAwareCapacity(needed, aliased);
3157	var m = r.toMutable();
3158	m.add(a.toConst(), b.toConst());
3159	r.setMetadata(m.positive, m.len);
3160	}
3161
3162	/// r = a + b with 2s-complement wrapping semantics. Returns whether any overflow occurred.
3163	///
3164	/// r, a and b may be aliases.
3165	///
3166	/// Returns an error if memory could not be allocated.
3167	pub fn addWrap(
3168	r: *Managed,
3169	a: *const Managed,
3170	b: *const Managed,
3171	signedness: Signedness,
3172	bit_count: usize,
3173	) Allocator.Error!bool {
3174	const aliased = limbsAliasDistinct(r, a) or limbsAliasDistinct(r, b);
3175	const needed = calcTwosCompLimbCount(bit_count);
3176	try r.ensureAliasAwareCapacity(needed, aliased);
3177	var m = r.toMutable();
3178	const wrapped = m.addWrap(a.toConst(), b.toConst(), signedness, bit_count);
3179	r.setMetadata(m.positive, m.len);
3180	return wrapped;
3181	}
3182
3183	/// r = a + b with 2s-complement saturating semantics.
3184	///
3185	/// r, a and b may be aliases.
3186	///
3187	/// Returns an error if memory could not be allocated.
3188	pub fn addSat(r: *Managed, a: *const Managed, b: *const Managed, signedness: Signedness, bit_count: usize) Allocator.Error!void {
3189	const aliased = limbsAliasDistinct(r, a) or limbsAliasDistinct(r, b);
3190	const needed = calcTwosCompLimbCount(bit_count);
3191	try r.ensureAliasAwareCapacity(needed, aliased);
3192	var m = r.toMutable();
3193	m.addSat(a.toConst(), b.toConst(), signedness, bit_count);
3194	r.setMetadata(m.positive, m.len);
3195	}
3196
3197	/// r = a - b
3198	///
3199	/// r, a and b may be aliases.
3200	///
3201	/// Returns an error if memory could not be allocated.
3202	pub fn sub(r: *Managed, a: *const Managed, b: *const Managed) !void {
3203	const aliased = limbsAliasDistinct(r, a) or limbsAliasDistinct(r, b);
3204	const needed = @max(a.len(), b.len()) + 1;
3205	try r.ensureAliasAwareCapacity(needed, aliased);
3206	var m = r.toMutable();
3207	m.sub(a.toConst(), b.toConst());
3208	r.setMetadata(m.positive, m.len);
3209	}
3210
3211	/// r = a - b with 2s-complement wrapping semantics. Returns whether any overflow occurred.
3212	///
3213	/// r, a and b may be aliases.
3214	///
3215	/// Returns an error if memory could not be allocated.
3216	pub fn subWrap(
3217	r: *Managed,
3218	a: *const Managed,
3219	b: *const Managed,
3220	signedness: Signedness,
3221	bit_count: usize,
3222	) Allocator.Error!bool {
3223	const aliased = limbsAliasDistinct(r, a) or limbsAliasDistinct(r, b);
3224	const needed = calcTwosCompLimbCount(bit_count);
3225	try r.ensureAliasAwareCapacity(needed, aliased);
3226	var m = r.toMutable();
3227	const wrapped = m.subWrap(a.toConst(), b.toConst(), signedness, bit_count);
3228	r.setMetadata(m.positive, m.len);
3229	return wrapped;
3230	}
3231
3232	/// r = a - b with 2s-complement saturating semantics.
3233	///
3234	/// r, a and b may be aliases.
3235	///
3236	/// Returns an error if memory could not be allocated.
3237	pub fn subSat(
3238	r: *Managed,
3239	a: *const Managed,
3240	b: *const Managed,
3241	signedness: Signedness,
3242	bit_count: usize,
3243	) Allocator.Error!void {
3244	const aliased = limbsAliasDistinct(r, a) or limbsAliasDistinct(r, b);
3245	const needed = calcTwosCompLimbCount(bit_count);
3246	try r.ensureAliasAwareCapacity(needed, aliased);
3247	var m = r.toMutable();
3248	m.subSat(a.toConst(), b.toConst(), signedness, bit_count);
3249	r.setMetadata(m.positive, m.len);
3250	}
3251
3252	/// rma = a * b
3253	///
3254	/// rma, a and b may be aliases. However, it is more efficient if rma does not alias a or b.
3255	///
3256	/// Returns an error if memory could not be allocated.
3257	///
3258	/// rma's allocator is used for temporary storage to speed up the multiplication.
3259	pub fn mul(rma: *Managed, a: *const Managed, b: *const Managed) !void {
3260	var alias_count: usize = 0;
3261	if (rma.limbs.ptr == a.limbs.ptr)
3262	alias_count += 1;
3263	if (rma.limbs.ptr == b.limbs.ptr)
3264	alias_count += 1;
3265	const needed = a.len() + b.len() + 1;
3266	const capacity_alias = limbsAliasDistinct(rma, a) or limbsAliasDistinct(rma, b);
3267	try rma.ensureAliasAwareCapacity(needed, capacity_alias);
3268	var m = rma.toMutable();
3269	if (alias_count == 0) {
3270	m.mulNoAlias(a.toConst(), b.toConst(), rma.allocator);
3271	} else {
3272	const limb_count = calcMulLimbsBufferLen(a.len(), b.len(), alias_count);
3273	const limbs_buffer = try rma.allocator.alloc(Limb, limb_count);
3274	defer rma.allocator.free(limbs_buffer);
3275	m.mul(a.toConst(), b.toConst(), limbs_buffer, rma.allocator);
3276	}
3277	rma.setMetadata(m.positive, m.len);
3278	}
3279
3280	/// rma = a * b with 2s-complement wrapping semantics.
3281	///
3282	/// rma, a and b may be aliases. However, it is more efficient if rma does not alias a or b.
3283	///
3284	/// Returns an error if memory could not be allocated.
3285	///
3286	/// rma's allocator is used for temporary storage to speed up the multiplication.
3287	pub fn mulWrap(
3288	rma: *Managed,
3289	a: *const Managed,
3290	b: *const Managed,
3291	signedness: Signedness,
3292	bit_count: usize,
3293	) !void {
3294	var alias_count: usize = 0;
3295	if (rma.limbs.ptr == a.limbs.ptr)
3296	alias_count += 1;
3297	if (rma.limbs.ptr == b.limbs.ptr)
3298	alias_count += 1;
3299	const needed = calcTwosCompLimbCount(bit_count);
3300	const capacity_alias = limbsAliasDistinct(rma, a) or limbsAliasDistinct(rma, b);
3301	try rma.ensureAliasAwareCapacity(needed, capacity_alias);
3302	var m = rma.toMutable();
3303	if (alias_count == 0) {
3304	m.mulWrapNoAlias(a.toConst(), b.toConst(), signedness, bit_count, rma.allocator);
3305	} else {
3306	const limb_count = calcMulWrapLimbsBufferLen(bit_count, a.len(), b.len(), alias_count);
3307	const limbs_buffer = try rma.allocator.alloc(Limb, limb_count);
3308	defer rma.allocator.free(limbs_buffer);
3309	m.mulWrap(a.toConst(), b.toConst(), signedness, bit_count, limbs_buffer, rma.allocator);
3310	}
3311	rma.setMetadata(m.positive, m.len);
3312	}
3313
3314	pub fn ensureTwosCompCapacity(r: *Managed, bit_count: usize) !void {
3315	try r.ensureCapacity(calcTwosCompLimbCount(bit_count));
3316	}
3317
3318	/// True if two distinct `Managed` parameters share the same limbs buffer.
3319	///
3320	/// We specifically exclude the case where `@intFromPtr(a) == @intFromPtr(b)` (same object).
3321	/// When both pointers refer to the same `Managed` instance, `ensureCapacity` can reallocate
3322	/// the buffer (if needed) without creating dangling pointers for that object.
3323	fn limbsAliasDistinct(a: *const Managed, b: *const Managed) bool {
3324	return @intFromPtr(a) != @intFromPtr(b) and a.limbs.ptr == b.limbs.ptr;
3325	}
3326
3327	/// When `aliased` is false (including when both pointers refer to the same object),
3328	/// `ensureCapacity` may reallocate; callers who rely on distinct `Managed` instances
3329	/// aliasing must ensure capacity before aliasing.
3330	/// See https://github.com/ziglang/zig/issues/6167
3331	fn ensureAliasAwareCapacity(r: *Managed, needed: usize, aliased: bool) !void {
3332	if (aliased) {
3333	assert(needed <= r.limbs.len);
3334	} else {
3335	try r.ensureCapacity(needed);
3336	}
3337	}
3338
3339	/// Use this function before doing `addScalar` if some of your parameters alias each other
3340	pub fn ensureAddScalarCapacity(r: *Managed, a: *const Managed, scalar: anytype) !void {
3341	try r.ensureCapacity(@max(a.len(), calcLimbLen(scalar)) + 1);
3342	}
3343
3344	/// Use this function before doing `add` if some of your parameters alias each other
3345	pub fn ensureAddCapacity(r: *Managed, a: *const Managed, b: *const Managed) !void {
3346	try r.ensureCapacity(@max(a.len(), b.len()) + 1);
3347	}
3348
3349	/// Use this function before doing `mul` if some of your parameters alias each other
3350	pub fn ensureMulCapacity(rma: *Managed, a: *const Managed, b: *const Managed) !void {
3351	try rma.ensureCapacity(a.len() + b.len() + 1);
3352	}
3353
3354	/// q = a / b (rem r)
3355	///
3356	/// a / b are floored (rounded towards 0).
3357	///
3358	/// Returns an error if memory could not be allocated.
3359	pub fn divFloor(q: *Managed, r: *Managed, a: *const Managed, b: *const Managed) !void {
3360	const q_alias = limbsAliasDistinct(q, a) or limbsAliasDistinct(q, b);
3361	const r_alias = limbsAliasDistinct(r, a) or limbsAliasDistinct(r, b);
3362	try q.ensureAliasAwareCapacity(a.len(), q_alias);
3363	try r.ensureAliasAwareCapacity(b.len(), r_alias);
3364	var mq = q.toMutable();
3365	var mr = r.toMutable();
3366	const limbs_buffer = try q.allocator.alloc(Limb, calcDivLimbsBufferLen(a.len(), b.len()));
3367	defer q.allocator.free(limbs_buffer);
3368	mq.divFloor(&mr, a.toConst(), b.toConst(), limbs_buffer);
3369	q.setMetadata(mq.positive, mq.len);
3370	r.setMetadata(mr.positive, mr.len);
3371	}
3372
3373	/// q = a / b (rem r)
3374	///
3375	/// a / b are ceiled (rounded towards positive infinity).
3376	///
3377	/// Returns an error if memory could not be allocated.
3378	pub fn divCeil(q: *Managed, r: *Managed, a: *const Managed, b: *const Managed) !void {
3379	const q_alias = limbsAliasDistinct(q, a) or limbsAliasDistinct(q, b);
3380	const r_alias = limbsAliasDistinct(r, a) or limbsAliasDistinct(r, b);
3381	try q.ensureAliasAwareCapacity(a.len(), q_alias);
3382	try r.ensureAliasAwareCapacity(b.len(), r_alias);
3383	var mq = q.toMutable();
3384	var mr = r.toMutable();
3385	const limbs_buffer = try q.allocator.alloc(Limb, calcDivLimbsBufferLen(a.len(), b.len()));
3386	defer q.allocator.free(limbs_buffer);
3387	mq.divCeil(&mr, a.toConst(), b.toConst(), limbs_buffer);
3388	q.setMetadata(mq.positive, mq.len);
3389	r.setMetadata(mr.positive, mr.len);
3390	}
3391
3392	/// q = a / b (rem r)
3393	///
3394	/// a / b are truncated (rounded towards -inf).
3395	///
3396	/// Returns an error if memory could not be allocated.
3397	pub fn divTrunc(q: *Managed, r: *Managed, a: *const Managed, b: *const Managed) !void {
3398	const q_alias = limbsAliasDistinct(q, a) or limbsAliasDistinct(q, b);
3399	const r_alias = limbsAliasDistinct(r, a) or limbsAliasDistinct(r, b);
3400	try q.ensureAliasAwareCapacity(a.len(), q_alias);
3401	try r.ensureAliasAwareCapacity(b.len(), r_alias);
3402	var mq = q.toMutable();
3403	var mr = r.toMutable();
3404	const limbs_buffer = try q.allocator.alloc(Limb, calcDivLimbsBufferLen(a.len(), b.len()));
3405	defer q.allocator.free(limbs_buffer);
3406	mq.divTrunc(&mr, a.toConst(), b.toConst(), limbs_buffer);
3407	q.setMetadata(mq.positive, mq.len);
3408	r.setMetadata(mr.positive, mr.len);
3409	}
3410
3411	/// r = a << shift, in other words, r = a * 2^shift
3412	/// r and a may alias.
3413	pub fn shiftLeft(r: *Managed, a: *const Managed, shift: usize) !void {
3414	const aliased = limbsAliasDistinct(r, a);
3415	const needed = a.len() + (shift / limb_bits) + 1;
3416	try r.ensureAliasAwareCapacity(needed, aliased);
3417	var m = r.toMutable();
3418	m.shiftLeft(a.toConst(), shift);
3419	r.setMetadata(m.positive, m.len);
3420	}
3421
3422	/// r = a <<| shift with 2s-complement saturating semantics.
3423	/// r and a may alias.
3424	pub fn shiftLeftSat(r: *Managed, a: *const Managed, shift: usize, signedness: Signedness, bit_count: usize) !void {
3425	const aliased = limbsAliasDistinct(r, a);
3426	const needed = calcTwosCompLimbCount(bit_count);
3427	try r.ensureAliasAwareCapacity(needed, aliased);
3428	var m = r.toMutable();
3429	m.shiftLeftSat(a.toConst(), shift, signedness, bit_count);
3430	r.setMetadata(m.positive, m.len);
3431	}
3432
3433	/// r = a >> shift
3434	/// r and a may alias.
3435	pub fn shiftRight(r: *Managed, a: *const Managed, shift: usize) !void {
3436	if (a.len() <= shift / limb_bits) {
3437	// Shifting negative numbers converges to -1 instead of 0
3438	if (a.isPositive()) {
3439	r.metadata = 1;
3440	r.limbs[0] = 0;
3441	} else {
3442	r.metadata = 1;
3443	r.setSign(false);
3444	r.limbs[0] = 1;
3445	}
3446	return;
3447	}
3448
3449	const aliased = limbsAliasDistinct(r, a);
3450	const needed = a.len() - (shift / limb_bits);
3451	try r.ensureAliasAwareCapacity(needed, aliased);
3452	var m = r.toMutable();
3453	m.shiftRight(a.toConst(), shift);
3454	r.setMetadata(m.positive, m.len);
3455	}
3456
3457	/// r = ~a under 2s-complement wrapping semantics.
3458	/// r and a may alias.
3459	pub fn bitNotWrap(r: *Managed, a: *const Managed, signedness: Signedness, bit_count: usize) !void {
3460	const aliased = limbsAliasDistinct(r, a);
3461	const needed = calcTwosCompLimbCount(bit_count);
3462	try r.ensureAliasAwareCapacity(needed, aliased);
3463	var m = r.toMutable();
3464	m.bitNotWrap(a.toConst(), signedness, bit_count);
3465	r.setMetadata(m.positive, m.len);
3466	}
3467
3468	/// r = a | b
3469	///
3470	/// a and b are zero-extended to the longer of a or b.
3471	pub fn bitOr(r: *Managed, a: *const Managed, b: *const Managed) !void {
3472	const aliased = limbsAliasDistinct(r, a) or limbsAliasDistinct(r, b);
3473	const needed = @max(a.len(), b.len());
3474	try r.ensureAliasAwareCapacity(needed, aliased);
3475	var m = r.toMutable();
3476	m.bitOr(a.toConst(), b.toConst());
3477	r.setMetadata(m.positive, m.len);
3478	}
3479
3480	/// r = a & b
3481	pub fn bitAnd(r: *Managed, a: *const Managed, b: *const Managed) !void {
3482	const cap = if (a.len() >= b.len())
3483	if (b.isPositive()) b.len() else if (a.isPositive()) a.len() else a.len() + 1
3484	else if (a.isPositive()) a.len() else if (b.isPositive()) b.len() else b.len() + 1;
3485
3486	const aliased = limbsAliasDistinct(r, a) or limbsAliasDistinct(r, b);
3487	try r.ensureAliasAwareCapacity(cap, aliased);
3488	var m = r.toMutable();
3489	m.bitAnd(a.toConst(), b.toConst());
3490	r.setMetadata(m.positive, m.len);
3491	}
3492
3493	/// r = a ^ b
3494	pub fn bitXor(r: *Managed, a: *const Managed, b: *const Managed) !void {
3495	const cap = @max(a.len(), b.len()) + @intFromBool(a.isPositive() != b.isPositive());
3496	const aliased = limbsAliasDistinct(r, a) or limbsAliasDistinct(r, b);
3497	try r.ensureAliasAwareCapacity(cap, aliased);
3498
3499	var m = r.toMutable();
3500	m.bitXor(a.toConst(), b.toConst());
3501	r.setMetadata(m.positive, m.len);
3502	}
3503
3504	/// rma may alias x or y.
3505	/// x and y may alias each other.
3506	///
3507	/// rma's allocator is used for temporary storage to boost multiplication performance.
3508	pub fn gcd(rma: *Managed, x: *const Managed, y: *const Managed) !void {
3509	const aliased = limbsAliasDistinct(rma, x) or limbsAliasDistinct(rma, y);
3510	const needed = @min(x.len(), y.len());
3511	try rma.ensureAliasAwareCapacity(needed, aliased);
3512	var m = rma.toMutable();
3513	var limbs_buffer = std.array_list.Managed(Limb).init(rma.allocator);
3514	defer limbs_buffer.deinit();
3515	try m.gcd(x.toConst(), y.toConst(), &limbs_buffer);
3516	rma.setMetadata(m.positive, m.len);
3517	}
3518
3519	/// r = a * a
3520	pub fn sqr(rma: *Managed, a: *const Managed) !void {
3521	const needed_limbs = 2 * a.len() + 1;
3522	const capacity_alias = limbsAliasDistinct(rma, a);
3523	const same_buffer = rma.limbs.ptr == a.limbs.ptr;
3524	try rma.ensureAliasAwareCapacity(needed_limbs, capacity_alias);
3525
3526	if (same_buffer) {
3527	const a_len = a.len();
3528	const tmp = try rma.allocator.alloc(Limb, a_len);
3529	defer rma.allocator.free(tmp);
3530	@memcpy(tmp[0..a_len], a.limbs[0..a_len]);
3531	const a_const: Const = .{ .limbs = tmp[0..a_len], .positive = a.isPositive() };
3532	var rma_mut = rma.toMutable();
3533	rma_mut.sqrNoAlias(a_const, rma.allocator);
3534	rma.setMetadata(rma_mut.positive, rma_mut.len);
3535	} else {
3536	var rma_mut = rma.toMutable();
3537	rma_mut.sqrNoAlias(a.toConst(), rma.allocator);
3538	rma.setMetadata(rma_mut.positive, rma_mut.len);
3539	}
3540	}
3541
3542	pub fn pow(rma: *Managed, a: *const Managed, b: u32) !void {
3543	const needed_limbs = calcPowLimbsBufferLen(a.bitCountAbs(), b);
3544	const capacity_alias = limbsAliasDistinct(rma, a);
3545	const same_buffer = rma.limbs.ptr == a.limbs.ptr;
3546
3547	try rma.ensureAliasAwareCapacity(needed_limbs, capacity_alias);
3548	const limbs_buffer = try rma.allocator.alloc(Limb, needed_limbs);
3549	defer rma.allocator.free(limbs_buffer);
3550
3551	if (same_buffer) {
3552	const a_len = a.len();
3553	const tmp = try rma.allocator.alloc(Limb, a_len);
3554	defer rma.allocator.free(tmp);
3555	@memcpy(tmp[0..a_len], a.limbs[0..a_len]);
3556	const a_const: Const = .{ .limbs = tmp[0..a_len], .positive = a.isPositive() };
3557	var rma_mut = rma.toMutable();
3558	rma_mut.pow(a_const, b, limbs_buffer);
3559	rma.setMetadata(rma_mut.positive, rma_mut.len);
3560	} else {
3561	var rma_mut = rma.toMutable();
3562	rma_mut.pow(a.toConst(), b, limbs_buffer);
3563	rma.setMetadata(rma_mut.positive, rma_mut.len);
3564	}
3565	}
3566
3567	/// r = ⌊√a⌋
3568	pub fn sqrt(rma: *Managed, a: *const Managed) !void {
3569	const bit_count = a.bitCountAbs();
3570	const aliased = limbsAliasDistinct(rma, a);
3571
3572	if (bit_count == 0) {
3573	try rma.set(0);
3574	rma.setMetadata(a.isPositive(), rma.len());
3575	return;
3576	}
3577
3578	if (!a.isPositive()) {
3579	return error.SqrtOfNegativeNumber;
3580	}
3581
3582	const needed_limbs = calcSqrtLimbsBufferLen(bit_count);
3583	const limbs_buffer = try rma.allocator.alloc(Limb, needed_limbs);
3584	defer rma.allocator.free(limbs_buffer);
3585
3586	const needed = (a.len() - 1) / 2 + 1;
3587	try rma.ensureAliasAwareCapacity(needed, aliased);
3588	var m = rma.toMutable();
3589	m.sqrt(a.toConst(), limbs_buffer);
3590	rma.setMetadata(m.positive, m.len);
3591	}
3592
3593	/// r = truncate(Int(signedness, bit_count), a)
3594	pub fn truncate(r: *Managed, a: *const Managed, signedness: Signedness, bit_count: usize) !void {
3595	const aliased = limbsAliasDistinct(r, a);
3596	const needed = calcTwosCompLimbCount(bit_count);
3597	try r.ensureAliasAwareCapacity(needed, aliased);
3598	var m = r.toMutable();
3599	m.truncate(a.toConst(), signedness, bit_count);
3600	r.setMetadata(m.positive, m.len);
3601	}
3602
3603	/// r = saturate(Int(signedness, bit_count), a)
3604	pub fn saturate(r: *Managed, a: *const Managed, signedness: Signedness, bit_count: usize) !void {
3605	const aliased = limbsAliasDistinct(r, a);
3606	const needed = calcTwosCompLimbCount(bit_count);
3607	try r.ensureAliasAwareCapacity(needed, aliased);
3608	var m = r.toMutable();
3609	m.saturate(a.toConst(), signedness, bit_count);
3610	r.setMetadata(m.positive, m.len);
3611	}
3612
3613	/// r = @popCount(a) with 2s-complement semantics.
3614	/// r and a may be aliases.
3615	pub fn popCount(r: *Managed, a: *const Managed, bit_count: usize) !void {
3616	const aliased = limbsAliasDistinct(r, a);
3617	const needed = calcTwosCompLimbCount(bit_count);
3618	try r.ensureAliasAwareCapacity(needed, aliased);
3619	var m = r.toMutable();
3620	m.popCount(a.toConst(), bit_count);
3621	r.setMetadata(m.positive, m.len);
3622	}
3623	};
3624
3625	/// Different operators which can be used in accumulation style functions
3626	/// (llmulacc, llmulaccKaratsuba, llmulaccLong, llmulLimb). In all these functions,
3627	/// a computed value is accumulated with an existing result.
3628	const AccOp = enum {
3629	/// The computed value is added to the result.
3630	add,
3631
3632	/// The computed value is subtracted from the result.
3633	sub,
3634	};
3635
3636	/// Knuth 4.3.1, Algorithm M.
3637	///
3638	/// r = r (op) a * b
3639	/// r MUST NOT alias any of a or b.
3640	///
3641	/// The result is computed modulo `r.len`. When `r.len >= a.len + b.len`, no overflow occurs.
3642	fn llmulacc(comptime op: AccOp, opt_allocator: ?Allocator, r: []Limb, a: []const Limb, b: []const Limb) void {
3643	assert(r.len >= a.len);
3644	assert(r.len >= b.len);
3645	assert(!slicesOverlap(r, a));
3646	assert(!slicesOverlap(r, b));
3647
3648	// Order greatest first.
3649	var x = a;
3650	var y = b;
3651	if (a.len < b.len) {
3652	x = b;
3653	y = a;
3654	}
3655
3656	k_mul: {
3657	if (y.len > 48) {
3658	if (opt_allocator) |allocator| {
3659	llmulaccKaratsuba(op, allocator, r, x, y) catch |err| switch (err) {
3660	error.OutOfMemory => break :k_mul, // handled below
3661	};
3662	return;
3663	}
3664	}
3665	}
3666
3667	llmulaccLong(op, r, x, y);
3668	}
3669
3670	/// Knuth 4.3.1, Algorithm M.
3671	///
3672	/// r = r (op) a * b
3673	/// r MUST NOT alias any of a or b.
3674	///
3675	/// The result is computed modulo `r.len`. When `r.len >= a.len + b.len`, no overflow occurs.
3676	fn llmulaccKaratsuba(
3677	comptime op: AccOp,
3678	allocator: Allocator,
3679	r: []Limb,
3680	a: []const Limb,
3681	b: []const Limb,
3682	) error{OutOfMemory}!void {
3683	assert(r.len >= a.len);
3684	assert(a.len >= b.len);
3685	assert(!slicesOverlap(r, a));
3686	assert(!slicesOverlap(r, b));
3687
3688	// Classical karatsuba algorithm:
3689	// a = a1 * B + a0
3690	// b = b1 * B + b0
3691	// Where a0, b0 < B
3692	//
3693	// We then have:
3694	// ab = a * b
3695	// = (a1 * B + a0) * (b1 * B + b0)
3696	// = a1 * b1 * B * B + a1 * B * b0 + a0 * b1 * B + a0 * b0
3697	// = a1 * b1 * B * B + (a1 * b0 + a0 * b1) * B + a0 * b0
3698	//
3699	// Note that:
3700	// a1 * b0 + a0 * b1
3701	// = (a1 + a0)(b1 + b0) - a1 * b1 - a0 * b0
3702	// = (a0 - a1)(b1 - b0) + a1 * b1 + a0 * b0
3703	//
3704	// This yields:
3705	// ab = p2 * B^2 + (p0 + p1 + p2) * B + p0
3706	//
3707	// Where:
3708	// p0 = a0 * b0
3709	// p1 = (a0 - a1)(b1 - b0)
3710	// p2 = a1 * b1
3711	//
3712	// Note, (a0 - a1) and (b1 - b0) produce values -B < x < B, and so we need to mind the sign here.
3713	// We also have:
3714	// 0 <= p0 <= 2B
3715	// -2B <= p1 <= 2B
3716	//
3717	// Note, when B is a multiple of the limb size, multiplies by B amount to shifts or
3718	// slices of a limbs array.
3719	//
3720	// This function computes the result of the multiplication modulo r.len. This means:
3721	// - p2 and p1 only need to be computed modulo r.len - B.
3722	// - In the case of p2, p2 * B^2 needs to be added modulo r.len - 2 * B.
3723
3724	const split = b.len / 2; // B
3725
3726	const limbs_after_split = r.len - split; // Limbs to compute for p1 and p2.
3727	const limbs_after_split2 = r.len - split * 2; // Limbs to add for p2 * B^2.
3728
3729	// For a0 and b0 we need the full range.
3730	const a0 = a[0..llnormalize(a[0..split])];
3731	const b0 = b[0..llnormalize(b[0..split])];
3732
3733	// For a1 and b1 we only need `limbs_after_split` limbs.
3734	const a1 = blk: {
3735	var a1 = a[split..];
3736	a1.len = @min(llnormalize(a1), limbs_after_split);
3737	break :blk a1;
3738	};
3739
3740	const b1 = blk: {
3741	var b1 = b[split..];
3742	b1.len = @min(llnormalize(b1), limbs_after_split);
3743	break :blk b1;
3744	};
3745
3746	// Note that the above slices relative to `split` work because we have a.len > b.len.
3747
3748	// We need some temporary memory to store intermediate results.
3749	// Note, we can reduce the amount of temporaries we need by reordering the computation here:
3750	// ab = p2 * B^2 + (p0 + p1 + p2) * B + p0
3751	// = p2 * B^2 + (p0 * B + p1 * B + p2 * B) + p0
3752	// = (p2 * B^2 + p2 * B) + (p0 * B + p0) + p1 * B
3753
3754	// Allocate at least enough memory to be able to multiply the upper two segments of a and b, assuming
3755	// no overflow.
3756	const tmp = try allocator.alloc(Limb, a.len - split + b.len - split);
3757	defer allocator.free(tmp);
3758
3759	// Compute p2.
3760	// Note, we don't need to compute all of p2, just enough limbs to satisfy r.
3761	const p2_limbs = @min(limbs_after_split, a1.len + b1.len);
3762
3763	@memset(tmp[0..p2_limbs], 0);
3764	llmulacc(.add, allocator, tmp[0..p2_limbs], a1[0..@min(a1.len, p2_limbs)], b1[0..@min(b1.len, p2_limbs)]);
3765	const p2 = tmp[0..llnormalize(tmp[0..p2_limbs])];
3766
3767	// Add p2 * B to the result.
3768	llaccum(op, r[split..], p2);
3769
3770	// Add p2 * B^2 to the result if required.
3771	if (limbs_after_split2 > 0) {
3772	llaccum(op, r[split * 2 ..], p2[0..@min(p2.len, limbs_after_split2)]);
3773	}
3774
3775	// Compute p0.
3776	// Since a0.len, b0.len <= split and r.len >= split * 2, the full width of p0 needs to be computed.
3777	const p0_limbs = a0.len + b0.len;
3778	@memset(tmp[0..p0_limbs], 0);
3779	llmulacc(.add, allocator, tmp[0..p0_limbs], a0, b0);
3780	const p0 = tmp[0..llnormalize(tmp[0..p0_limbs])];
3781
3782	// Add p0 to the result.
3783	llaccum(op, r, p0);
3784
3785	// Add p0 * B to the result. In this case, we may not need all of it.
3786	llaccum(op, r[split..], p0[0..@min(limbs_after_split, p0.len)]);
3787
3788	// Finally, compute and add p1.
3789	// From now on we only need `limbs_after_split` limbs for a0 and b0, since the result of the
3790	// following computation will be added * B.
3791	const a0x = a0[0..@min(a0.len, limbs_after_split)];
3792	const b0x = b0[0..@min(b0.len, limbs_after_split)];
3793
3794	const j0_sign = llcmp(a0x, a1);
3795	const j1_sign = llcmp(b1, b0x);
3796
3797	if (j0_sign * j1_sign == 0) {
3798	// p1 is zero, we don't need to do any computation at all.
3799	return;
3800	}
3801
3802	@memset(tmp, 0);
3803
3804	// p1 is nonzero, so compute the intermediary terms j0 = a0 - a1 and j1 = b1 - b0.
3805	// Note that in this case, we again need some storage for intermediary results
3806	// j0 and j1. Since we have tmp.len >= 2B, we can store both
3807	// intermediaries in the already allocated array.
3808	const j0 = tmp[0 .. a.len - split];
3809	const j1 = tmp[a.len - split ..];
3810
3811	// Ensure that no subtraction overflows.
3812	if (j0_sign == 1) {
3813	// a0 > a1.
3814	_ = llsubcarry(j0, a0x, a1);
3815	} else {
3816	// a0 < a1.
3817	_ = llsubcarry(j0, a1, a0x);
3818	}
3819
3820	if (j1_sign == 1) {
3821	// b1 > b0.
3822	_ = llsubcarry(j1, b1, b0x);
3823	} else {
3824	// b1 > b0.
3825	_ = llsubcarry(j1, b0x, b1);
3826	}
3827
3828	if (j0_sign * j1_sign == 1) {
3829	// If j0 and j1 are both positive, we now have:
3830	// p1 = j0 * j1
3831	// If j0 and j1 are both negative, we now have:
3832	// p1 = -j0 * -j1 = j0 * j1
3833	// In this case we can add p1 to the result using llmulacc.
3834	llmulacc(op, allocator, r[split..], j0[0..llnormalize(j0)], j1[0..llnormalize(j1)]);
3835	} else {
3836	// In this case either j0 or j1 is negative, an we have:
3837	// p1 = -(j0 * j1)
3838	// Now we need to subtract instead of accumulate.
3839	const inverted_op = if (op == .add) .sub else .add;
3840	llmulacc(inverted_op, allocator, r[split..], j0[0..llnormalize(j0)], j1[0..llnormalize(j1)]);
3841	}
3842	}
3843
3844	/// r = r (op) a.
3845	/// The result is computed modulo `r.len`.
3846	fn llaccum(comptime op: AccOp, r: []Limb, a: []const Limb) void {
3847	assert(!slicesOverlap(r, a) or @intFromPtr(r.ptr) <= @intFromPtr(a.ptr));
3848	if (op == .sub) {
3849	_ = llsubcarry(r, r, a);
3850	return;
3851	}
3852
3853	assert(r.len != 0 and a.len != 0);
3854	assert(r.len >= a.len);
3855
3856	var i: usize = 0;
3857	var carry: Limb = 0;
3858
3859	while (i < a.len) : (i += 1) {
3860	const ov1 = @addWithOverflow(r[i], a[i]);
3861	r[i] = ov1[0];
3862	const ov2 = @addWithOverflow(r[i], carry);
3863	r[i] = ov2[0];
3864	carry = @as(Limb, ov1[1]) + ov2[1];
3865	}
3866
3867	while ((carry != 0) and i < r.len) : (i += 1) {
3868	const ov = @addWithOverflow(r[i], carry);
3869	r[i] = ov[0];
3870	carry = ov[1];
3871	}
3872	}
3873
3874	/// Returns -1, 0, 1 if |a| < |b|, |a| == |b| or |a| > |b| respectively for limbs.
3875	pub fn llcmp(a: []const Limb, b: []const Limb) i8 {
3876	const a_len = llnormalize(a);
3877	const b_len = llnormalize(b);
3878	if (a_len < b_len) {
3879	return -1;
3880	}
3881	if (a_len > b_len) {
3882	return 1;
3883	}
3884
3885	var i: usize = a_len - 1;
3886	while (i != 0) : (i -= 1) {
3887	if (a[i] != b[i]) {
3888	break;
3889	}
3890	}
3891
3892	if (a[i] < b[i]) {
3893	return -1;
3894	} else if (a[i] > b[i]) {
3895	return 1;
3896	} else {
3897	return 0;
3898	}
3899	}
3900
3901	/// r = r (op) y * xi
3902	/// The result is computed modulo `r.len`. When `r.len >= a.len + b.len`, no overflow occurs.
3903	fn llmulaccLong(comptime op: AccOp, r: []Limb, a: []const Limb, b: []const Limb) void {
3904	assert(r.len >= a.len);
3905	assert(a.len >= b.len);
3906
3907	var i: usize = 0;
3908	while (i < b.len) : (i += 1) {
3909	_ = llmulLimb(op, r[i..], a, b[i]);
3910	}
3911	}
3912
3913	/// r = r (op) y * xi
3914	/// The result is computed modulo `r.len`.
3915	/// Returns whether the operation overflowed.
3916	fn llmulLimb(comptime op: AccOp, acc: []Limb, y: []const Limb, xi: Limb) bool {
3917	assert(!slicesOverlap(acc, y) or @intFromPtr(acc.ptr) <= @intFromPtr(y.ptr));
3918
3919	if (xi == 0) {
3920	return false;
3921	}
3922
3923	const split = @min(y.len, acc.len);
3924	var a_lo = acc[0..split];
3925	var a_hi = acc[split..];
3926
3927	switch (op) {
3928	.add => {
3929	var carry: Limb = 0;
3930	var j: usize = 0;
3931	while (j < a_lo.len) : (j += 1) {
3932	a_lo[j] = addMulLimbWithCarry(a_lo[j], y[j], xi, &carry);
3933	}
3934
3935	j = 0;
3936	while ((carry != 0) and (j < a_hi.len)) : (j += 1) {
3937	const ov = @addWithOverflow(a_hi[j], carry);
3938	a_hi[j] = ov[0];
3939	carry = ov[1];
3940	}
3941
3942	return carry != 0;
3943	},
3944	.sub => {
3945	var borrow: Limb = 0;
3946	var j: usize = 0;
3947	while (j < a_lo.len) : (j += 1) {
3948	a_lo[j] = subMulLimbWithBorrow(a_lo[j], y[j], xi, &borrow);
3949	}
3950
3951	j = 0;
3952	while ((borrow != 0) and (j < a_hi.len)) : (j += 1) {
3953	const ov = @subWithOverflow(a_hi[j], borrow);
3954	a_hi[j] = ov[0];
3955	borrow = ov[1];
3956	}
3957
3958	return borrow != 0;
3959	},
3960	}
3961	}
3962
3963	/// returns the min length the limb could be.
3964	fn llnormalize(a: []const Limb) usize {
3965	var j = a.len;
3966	while (j > 0) : (j -= 1) {
3967	if (a[j - 1] != 0) {
3968	break;
3969	}
3970	}
3971
3972	// Handle zero
3973	return if (j != 0) j else 1;
3974	}
3975
3976	/// Knuth 4.3.1, Algorithm S.
3977	fn llsubcarry(r: []Limb, a: []const Limb, b: []const Limb) Limb {
3978	assert(a.len != 0 and b.len != 0);
3979	assert(a.len >= b.len);
3980	assert(r.len >= a.len);
3981	assert(!slicesOverlap(r, a) or @intFromPtr(r.ptr) <= @intFromPtr(a.ptr));
3982	assert(!slicesOverlap(r, b) or @intFromPtr(r.ptr) <= @intFromPtr(b.ptr));
3983
3984	var i: usize = 0;
3985	var borrow: Limb = 0;
3986
3987	while (i < b.len) : (i += 1) {
3988	const ov1 = @subWithOverflow(a[i], b[i]);
3989	r[i] = ov1[0];
3990	const ov2 = @subWithOverflow(r[i], borrow);
3991	r[i] = ov2[0];
3992	borrow = @as(Limb, ov1[1]) + ov2[1];
3993	}
3994
3995	while (i < a.len) : (i += 1) {
3996	const ov = @subWithOverflow(a[i], borrow);
3997	r[i] = ov[0];
3998	borrow = ov[1];
3999	}
4000
4001	return borrow;
4002	}
4003
4004	fn llsub(r: []Limb, a: []const Limb, b: []const Limb) void {
4005	assert(a.len > b.len or (a.len == b.len and a[a.len - 1] >= b[b.len - 1]));
4006	assert(llsubcarry(r, a, b) == 0);
4007	}
4008
4009	/// Knuth 4.3.1, Algorithm A.
4010	fn lladdcarry(r: []Limb, a: []const Limb, b: []const Limb) Limb {
4011	assert(a.len != 0 and b.len != 0);
4012	assert(a.len >= b.len);
4013	assert(r.len >= a.len);
4014	assert(!slicesOverlap(r, a) or @intFromPtr(r.ptr) <= @intFromPtr(a.ptr));
4015	assert(!slicesOverlap(r, b) or @intFromPtr(r.ptr) <= @intFromPtr(b.ptr));
4016
4017	var i: usize = 0;
4018	var carry: Limb = 0;
4019
4020	while (i < b.len) : (i += 1) {
4021	const ov1 = @addWithOverflow(a[i], b[i]);
4022	r[i] = ov1[0];
4023	const ov2 = @addWithOverflow(r[i], carry);
4024	r[i] = ov2[0];
4025	carry = @as(Limb, ov1[1]) + ov2[1];
4026	}
4027
4028	while (i < a.len) : (i += 1) {
4029	const ov = @addWithOverflow(a[i], carry);
4030	r[i] = ov[0];
4031	carry = ov[1];
4032	}
4033
4034	return carry;
4035	}
4036
4037	fn lladd(r: []Limb, a: []const Limb, b: []const Limb) void {
4038	assert(r.len >= a.len + 1);
4039	r[a.len] = lladdcarry(r, a, b);
4040	}
4041
4042	/// Knuth 4.3.1, Exercise 16.
4043	fn lldiv1(quo: []Limb, rem: *Limb, a: []const Limb, b: Limb) void {
4044	assert(a.len > 1 or a[0] >= b);
4045	assert(quo.len >= a.len);
4046
4047	rem.* = 0;
4048	for (a, 0..) |_, ri| {
4049	const i = a.len - ri - 1;
4050	const pdiv = ((@as(DoubleLimb, rem.*) << limb_bits) | a[i]);
4051
4052	if (pdiv == 0) {
4053	quo[i] = 0;
4054	rem.* = 0;
4055	} else if (pdiv < b) {
4056	quo[i] = 0;
4057	rem.* = @as(Limb, @truncate(pdiv));
4058	} else if (pdiv == b) {
4059	quo[i] = 1;
4060	rem.* = 0;
4061	} else {
4062	quo[i] = @as(Limb, @truncate(@divTrunc(pdiv, b)));
4063	rem.* = @as(Limb, @truncate(pdiv - (quo[i] *% b)));
4064	}
4065	}
4066	}
4067
4068	fn lldiv0p5(quo: []Limb, rem: *Limb, a: []const Limb, b: HalfLimb) void {
4069	assert(a.len > 1 or a[0] >= b);
4070	assert(quo.len >= a.len);
4071
4072	rem.* = 0;
4073	for (a, 0..) |_, ri| {
4074	const i = a.len - ri - 1;
4075	const ai_high = a[i] >> half_limb_bits;
4076	const ai_low = a[i] & ((1 << half_limb_bits) - 1);
4077
4078	// Split the division into two divisions acting on half a limb each. Carry remainder.
4079	const ai_high_with_carry = (rem.* << half_limb_bits) | ai_high;
4080	const ai_high_quo = ai_high_with_carry / b;
4081	rem.* = ai_high_with_carry % b;
4082
4083	const ai_low_with_carry = (rem.* << half_limb_bits) | ai_low;
4084	const ai_low_quo = ai_low_with_carry / b;
4085	rem.* = ai_low_with_carry % b;
4086
4087	quo[i] = (ai_high_quo << half_limb_bits) | ai_low_quo;
4088	}
4089	}
4090
4091	/// Performs r = a << shift and returns the amount of limbs affected
4092	///
4093	/// if a and r overlaps, then r.ptr >= a.ptr is asserted
4094	/// r must have the capacity to store a << shift
4095	fn llshl(r: []Limb, a: []const Limb, shift: usize) usize {
4096	std.debug.assert(a.len >= 1);
4097	if (slicesOverlap(a, r))
4098	std.debug.assert(@intFromPtr(r.ptr) >= @intFromPtr(a.ptr));
4099
4100	if (shift == 0) {
4101	if (a.ptr != r.ptr) @memmove(r[0..a.len], a);
4102	return a.len;
4103	}
4104	if (shift >= limb_bits) {
4105	const limb_shift = shift / limb_bits;
4106
4107	const affected = llshl(r[limb_shift..], a, shift % limb_bits);
4108	@memset(r[0..limb_shift], 0);
4109
4110	return limb_shift + affected;
4111	}
4112
4113	// shift is guaranteed to be < limb_bits
4114	const bit_shift: Log2Limb = @truncate(shift);
4115	const opposite_bit_shift: Log2Limb = @truncate(limb_bits - bit_shift);
4116
4117	// We only need the extra limb if the shift of the last element overflows.
4118	// This is useful for the implementation of `shiftLeftSat`.
4119	const overflows = a[a.len - 1] >> opposite_bit_shift != 0;
4120	if (overflows) {
4121	std.debug.assert(r.len >= a.len + 1);
4122	} else {
4123	std.debug.assert(r.len >= a.len);
4124	}
4125
4126	var i: usize = a.len;
4127	if (overflows) {
4128	// r is asserted to be large enough above
4129	r[a.len] = a[a.len - 1] >> opposite_bit_shift;
4130	}
4131	while (i > 1) {
4132	i -= 1;
4133	r[i] = (a[i - 1] >> opposite_bit_shift) | (a[i] << bit_shift);
4134	}
4135	r[0] = a[0] << bit_shift;
4136
4137	return a.len + @intFromBool(overflows);
4138	}
4139
4140	/// Performs r = a >> shift and returns the amount of limbs affected
4141	///
4142	/// if a and r overlaps, then r.ptr <= a.ptr is asserted
4143	/// r must have the capacity to store a >> shift
4144	///
4145	/// See tests below for examples of behaviour
4146	fn llshr(r: []Limb, a: []const Limb, shift: usize) usize {
4147	if (slicesOverlap(a, r))
4148	std.debug.assert(@intFromPtr(r.ptr) <= @intFromPtr(a.ptr));
4149
4150	if (a.len == 0) return 0;
4151
4152	if (shift == 0) {
4153	std.debug.assert(r.len >= a.len);
4154
4155	if (a.ptr != r.ptr) @memmove(r[0..a.len], a);
4156	return a.len;
4157	}
4158	if (shift >= limb_bits) {
4159	if (shift / limb_bits >= a.len) {
4160	r[0] = 0;
4161	return 1;
4162	}
4163	return llshr(r, a[shift / limb_bits ..], shift % limb_bits);
4164	}
4165
4166	// shift is guaranteed to be < limb_bits
4167	const bit_shift: Log2Limb = @truncate(shift);
4168	const opposite_bit_shift: Log2Limb = @truncate(limb_bits - bit_shift);
4169
4170	// special case, where there is a risk to set r to 0
4171	if (a.len == 1) {
4172	r[0] = a[0] >> bit_shift;
4173	return 1;
4174	}
4175	if (a.len == 0) {
4176	r[0] = 0;
4177	return 1;
4178	}
4179
4180	// if the most significant limb becomes 0 after the shift
4181	const shrink = a[a.len - 1] >> bit_shift == 0;
4182	std.debug.assert(r.len >= a.len - @intFromBool(shrink));
4183
4184	var i: usize = 0;
4185	while (i < a.len - 1) : (i += 1) {
4186	r[i] = (a[i] >> bit_shift) | (a[i + 1] << opposite_bit_shift);
4187	}
4188
4189	if (!shrink)
4190	r[i] = a[i] >> bit_shift;
4191
4192	return a.len - @intFromBool(shrink);
4193	}
4194
4195	// r = ~r
4196	fn llnot(r: []Limb) void {
4197	for (r) |*elem| {
4198	elem.* = ~elem.*;
4199	}
4200	}
4201
4202	// r = a | b with 2s complement semantics.
4203	// r may alias.
4204	// a and b must not be 0.
4205	// Returns `true` when the result is positive.
4206	// When b is positive, r requires at least `a.len` limbs of storage.
4207	// When b is negative, r requires at least `b.len` limbs of storage.
4208	fn llsignedor(r: []Limb, a: []const Limb, a_positive: bool, b: []const Limb, b_positive: bool) bool {
4209	assert(r.len >= a.len);
4210	assert(a.len >= b.len);
4211
4212	if (a_positive and b_positive) {
4213	// Trivial case, result is positive.
4214	var i: usize = 0;
4215	while (i < b.len) : (i += 1) {
4216	r[i] = a[i] | b[i];
4217	}
4218	while (i < a.len) : (i += 1) {
4219	r[i] = a[i];
4220	}
4221
4222	return true;
4223	} else if (!a_positive and b_positive) {
4224	// Result is negative.
4225	// r = (--a) | b
4226	// = ~(-a - 1) | b
4227	// = ~(-a - 1) | ~~b
4228	// = ~((-a - 1) & ~b)
4229	// = -(((-a - 1) & ~b) + 1)
4230
4231	var i: usize = 0;
4232	var a_borrow: u1 = 1;
4233	var r_carry: u1 = 1;
4234
4235	while (i < b.len) : (i += 1) {
4236	const ov1 = @subWithOverflow(a[i], a_borrow);
4237	a_borrow = ov1[1];
4238	const ov2 = @addWithOverflow(ov1[0] & ~b[i], r_carry);
4239	r[i] = ov2[0];
4240	r_carry = ov2[1];
4241	}
4242
4243	// In order for r_carry to be nonzero at this point, ~b[i] would need to be
4244	// all ones, which would require b[i] to be zero. This cannot be when
4245	// b is normalized, so there cannot be a carry here.
4246	// Also, x & ~b can only clear bits, so (x & ~b) <= x, meaning (-a - 1) + 1 never overflows.
4247	assert(r_carry == 0);
4248
4249	// With b = 0, we get (-a - 1) & ~0 = -a - 1.
4250	// Note, if a_borrow is zero we do not need to compute anything for
4251	// the higher limbs so we can early return here.
4252	while (i < a.len and a_borrow == 1) : (i += 1) {
4253	const ov = @subWithOverflow(a[i], a_borrow);
4254	r[i] = ov[0];
4255	a_borrow = ov[1];
4256	}
4257
4258	assert(a_borrow == 0); // a was 0.
4259
4260	return false;
4261	} else if (a_positive and !b_positive) {
4262	// Result is negative.
4263	// r = a | (--b)
4264	// = a | ~(-b - 1)
4265	// = ~~a | ~(-b - 1)
4266	// = ~(~a & (-b - 1))
4267	// = -((~a & (-b - 1)) + 1)
4268
4269	var i: usize = 0;
4270	var b_borrow: u1 = 1;
4271	var r_carry: u1 = 1;
4272
4273	while (i < b.len) : (i += 1) {
4274	const ov1 = @subWithOverflow(b[i], b_borrow);
4275	b_borrow = ov1[1];
4276	const ov2 = @addWithOverflow(~a[i] & ov1[0], r_carry);
4277	r[i] = ov2[0];
4278	r_carry = ov2[1];
4279	}
4280
4281	// b is at least 1, so this should never underflow.
4282	assert(b_borrow == 0); // b was 0
4283
4284	// x & ~a can only clear bits, so (x & ~a) <= x, meaning (-b - 1) + 1 never overflows.
4285	assert(r_carry == 0);
4286
4287	// With b = 0 and b_borrow = 0, we get ~a & (0 - 0) = ~a & 0 = 0.
4288	// Omit setting the upper bytes, just deal with those when calling llsignedor.
4289
4290	return false;
4291	} else {
4292	// Result is negative.
4293	// r = (--a) | (--b)
4294	// = ~(-a - 1) | ~(-b - 1)
4295	// = ~((-a - 1) & (-b - 1))
4296	// = -(~(~((-a - 1) & (-b - 1))) + 1)
4297	// = -((-a - 1) & (-b - 1) + 1)
4298
4299	var i: usize = 0;
4300	var a_borrow: u1 = 1;
4301	var b_borrow: u1 = 1;
4302	var r_carry: u1 = 1;
4303
4304	while (i < b.len) : (i += 1) {
4305	const ov1 = @subWithOverflow(a[i], a_borrow);
4306	a_borrow = ov1[1];
4307	const ov2 = @subWithOverflow(b[i], b_borrow);
4308	b_borrow = ov2[1];
4309	const ov3 = @addWithOverflow(ov1[0] & ov2[0], r_carry);
4310	r[i] = ov3[0];
4311	r_carry = ov3[1];
4312	}
4313
4314	// b is at least 1, so this should never underflow.
4315	assert(b_borrow == 0); // b was 0
4316
4317	// Can never overflow because in order for b_limb to be maxInt(Limb),
4318	// b_borrow would need to equal 1.
4319
4320	// x & y can only clear bits, meaning x & y <= x and x & y <= y. This implies that
4321	// for x = a - 1 and y = b - 1, the +1 term would never cause an overflow.
4322	assert(r_carry == 0);
4323
4324	// With b = 0 and b_borrow = 0 we get (-a - 1) & (0 - 0) = (-a - 1) & 0 = 0.
4325	// Omit setting the upper bytes, just deal with those when calling llsignedor.
4326	return false;
4327	}
4328	}
4329
4330	// r = a & b with 2s complement semantics.
4331	// r may alias.
4332	// a and b must not be 0.
4333	// Returns `true` when the result is positive.
4334	// We assume `a.len >= b.len` here, so:
4335	// 1. when b is positive, r requires at least `b.len` limbs of storage,
4336	// 2. when b is negative but a is positive, r requires at least `a.len` limbs of storage,
4337	// 3. when both a and b are negative, r requires at least `a.len + 1` limbs of storage.
4338	fn llsignedand(r: []Limb, a: []const Limb, a_positive: bool, b: []const Limb, b_positive: bool) bool {
4339	assert(a.len != 0 and b.len != 0);
4340	assert(a.len >= b.len);
4341	assert(r.len >= if (b_positive) b.len else if (a_positive) a.len else a.len + 1);
4342
4343	if (a_positive and b_positive) {
4344	// Trivial case, result is positive.
4345	var i: usize = 0;
4346	while (i < b.len) : (i += 1) {
4347	r[i] = a[i] & b[i];
4348	}
4349
4350	// With b = 0 we have a & 0 = 0, so the upper bytes are zero.
4351	// Omit setting them here and simply discard them whenever
4352	// llsignedand is called.
4353
4354	return true;
4355	} else if (!a_positive and b_positive) {
4356	// Result is positive.
4357	// r = (--a) & b
4358	// = ~(-a - 1) & b
4359
4360	var i: usize = 0;
4361	var a_borrow: u1 = 1;
4362
4363	while (i < b.len) : (i += 1) {
4364	const ov = @subWithOverflow(a[i], a_borrow);
4365	a_borrow = ov[1];
4366	r[i] = ~ov[0] & b[i];
4367	}
4368
4369	// With b = 0 we have ~(a - 1) & 0 = 0, so the upper bytes are zero.
4370	// Omit setting them here and simply discard them whenever
4371	// llsignedand is called.
4372
4373	return true;
4374	} else if (a_positive and !b_positive) {
4375	// Result is positive.
4376	// r = a & (--b)
4377	// = a & ~(-b - 1)
4378
4379	var i: usize = 0;
4380	var b_borrow: u1 = 1;
4381
4382	while (i < b.len) : (i += 1) {
4383	const ov = @subWithOverflow(b[i], b_borrow);
4384	b_borrow = ov[1];
4385	r[i] = a[i] & ~ov[0];
4386	}
4387
4388	assert(b_borrow == 0); // b was 0
4389
4390	// With b = 0 and b_borrow = 0 we have a & ~(0 - 0) = a & ~0 = a, so
4391	// the upper bytes are the same as those of a.
4392
4393	while (i < a.len) : (i += 1) {
4394	r[i] = a[i];
4395	}
4396
4397	return true;
4398	} else {
4399	// Result is negative.
4400	// r = (--a) & (--b)
4401	// = ~(-a - 1) & ~(-b - 1)
4402	// = ~((-a - 1) | (-b - 1))
4403	// = -(((-a - 1) | (-b - 1)) + 1)
4404
4405	var i: usize = 0;
4406	var a_borrow: u1 = 1;
4407	var b_borrow: u1 = 1;
4408	var r_carry: u1 = 1;
4409
4410	while (i < b.len) : (i += 1) {
4411	const ov1 = @subWithOverflow(a[i], a_borrow);
4412	a_borrow = ov1[1];
4413	const ov2 = @subWithOverflow(b[i], b_borrow);
4414	b_borrow = ov2[1];
4415	const ov3 = @addWithOverflow(ov1[0] | ov2[0], r_carry);
4416	r[i] = ov3[0];
4417	r_carry = ov3[1];
4418	}
4419
4420	// b is at least 1, so this should never underflow.
4421	assert(b_borrow == 0); // b was 0
4422
4423	// With b = 0 and b_borrow = 0 we get (-a - 1) | (0 - 0) = (-a - 1) | 0 = -a - 1.
4424	while (i < a.len) : (i += 1) {
4425	const ov1 = @subWithOverflow(a[i], a_borrow);
4426	a_borrow = ov1[1];
4427	const ov2 = @addWithOverflow(ov1[0], r_carry);
4428	r[i] = ov2[0];
4429	r_carry = ov2[1];
4430	}
4431
4432	assert(a_borrow == 0); // a was 0.
4433
4434	// The final addition can overflow here, so we need to keep that in mind.
4435	r[i] = r_carry;
4436
4437	return false;
4438	}
4439	}
4440
4441	// r = a ^ b with 2s complement semantics.
4442	// r may alias.
4443	// a and b must not be -0.
4444	// Returns `true` when the result is positive.
4445	// If the sign of a and b is equal, then r requires at least `@max(a.len, b.len)` limbs are required.
4446	// Otherwise, r requires at least `@max(a.len, b.len) + 1` limbs.
4447	fn llsignedxor(r: []Limb, a: []const Limb, a_positive: bool, b: []const Limb, b_positive: bool) bool {
4448	assert(a.len != 0 and b.len != 0);
4449	assert(r.len >= a.len);
4450	assert(a.len >= b.len);
4451
4452	// If a and b are positive, the result is positive and r = a ^ b.
4453	// If a negative, b positive, result is negative and we have
4454	// r = --(--a ^ b)
4455	// = --(~(-a - 1) ^ b)
4456	// = -(~(~(-a - 1) ^ b) + 1)
4457	// = -(((-a - 1) ^ b) + 1)
4458	// Same if a is positive and b is negative, sides switched.
4459	// If both a and b are negative, the result is positive and we have
4460	// r = (--a) ^ (--b)
4461	// = ~(-a - 1) ^ ~(-b - 1)
4462	// = (-a - 1) ^ (-b - 1)
4463	// These operations can be made more generic as follows:
4464	// - If a is negative, subtract 1 from |a| before the xor.
4465	// - If b is negative, subtract 1 from |b| before the xor.
4466	// - if the result is supposed to be negative, add 1.
4467
4468	var i: usize = 0;
4469	var a_borrow = @intFromBool(!a_positive);
4470	var b_borrow = @intFromBool(!b_positive);
4471	var r_carry = @intFromBool(a_positive != b_positive);
4472
4473	while (i < b.len) : (i += 1) {
4474	const ov1 = @subWithOverflow(a[i], a_borrow);
4475	a_borrow = ov1[1];
4476	const ov2 = @subWithOverflow(b[i], b_borrow);
4477	b_borrow = ov2[1];
4478	const ov3 = @addWithOverflow(ov1[0] ^ ov2[0], r_carry);
4479	r[i] = ov3[0];
4480	r_carry = ov3[1];
4481	}
4482
4483	while (i < a.len) : (i += 1) {
4484	const ov1 = @subWithOverflow(a[i], a_borrow);
4485	a_borrow = ov1[1];
4486	const ov2 = @addWithOverflow(ov1[0], r_carry);
4487	r[i] = ov2[0];
4488	r_carry = ov2[1];
4489	}
4490
4491	// If both inputs don't share the same sign, an extra limb is required.
4492	if (a_positive != b_positive) {
4493	r[i] = r_carry;
4494	} else {
4495	assert(r_carry == 0);
4496	}
4497
4498	assert(a_borrow == 0);
4499	assert(b_borrow == 0);
4500
4501	return a_positive == b_positive;
4502	}
4503
4504	/// r MUST NOT alias x.
4505	fn llsquareBasecase(r: []Limb, x: []const Limb) void {
4506	const x_norm = x;
4507	assert(r.len >= 2 * x_norm.len + 1);
4508	assert(!slicesOverlap(r, x));
4509
4510	// Compute the square of a N-limb bigint with only (N^2 + N)/2
4511	// multiplications by exploiting the symmetry of the coefficients around the
4512	// diagonal:
4513	//
4514	// a b c *
4515	// a b c =
4516	// -------------------
4517	// ca cb cc +
4518	// ba bb bc +
4519	// aa ab ac
4520	//
4521	// Note that:
4522	// - Each mixed-product term appears twice for each column,
4523	// - Squares are always in the 2k (0 <= k < N) column
4524
4525	for (x_norm, 0..) |v, i| {
4526	// Accumulate all the x[i]*x[j] (with x!=j) products
4527	const overflow = llmulLimb(.add, r[2 * i + 1 ..], x_norm[i + 1 ..], v);
4528	assert(!overflow);
4529	}
4530
4531	// Each product appears twice, multiply by 2
4532	_ = llshl(r, r[0 .. 2 * x_norm.len], 1);
4533
4534	for (x_norm, 0..) |v, i| {
4535	// Compute and add the squares
4536	const overflow = llmulLimb(.add, r[2 * i ..], x[i..][0..1], v);
4537	assert(!overflow);
4538	}
4539	}
4540
4541	/// Knuth 4.6.3
4542	fn llpow(r: []Limb, a: []const Limb, b: u32, tmp_limbs: []Limb) void {
4543	var tmp1: []Limb = undefined;
4544	var tmp2: []Limb = undefined;
4545
4546	// Multiplication requires no aliasing between the operand and the result
4547	// variable, use the output limbs and another temporary set to overcome this
4548	// limitation.
4549	// The initial assignment makes the result end in `r` so an extra memory
4550	// copy is saved, each 1 flips the index twice so it's only the zeros that
4551	// matter.
4552	const b_leading_zeros = @clz(b);
4553	const exp_zeros = @popCount(~b) - b_leading_zeros;
4554	if (exp_zeros & 1 != 0) {
4555	tmp1 = tmp_limbs;
4556	tmp2 = r;
4557	} else {
4558	tmp1 = r;
4559	tmp2 = tmp_limbs;
4560	}
4561
4562	@memcpy(tmp1[0..a.len], a);
4563	@memset(tmp1[a.len..], 0);
4564
4565	// Scan the exponent as a binary number, from left to right, dropping the
4566	// most significant bit set.
4567	// Square the result if the current bit is zero, square and multiply by a if
4568	// it is one.
4569	const exp_bits = 32 - 1 - b_leading_zeros;
4570	var exp = b << @as(u5, @intCast(1 + b_leading_zeros));
4571
4572	var i: usize = 0;
4573	while (i < exp_bits) : (i += 1) {
4574	// Square
4575	@memset(tmp2, 0);
4576	llsquareBasecase(tmp2, tmp1[0..llnormalize(tmp1)]);
4577	mem.swap([]Limb, &tmp1, &tmp2);
4578	// Multiply by a
4579	const ov = @shlWithOverflow(exp, 1);
4580	exp = ov[0];
4581	if (ov[1] != 0) {
4582	@memset(tmp2, 0);
4583	llmulacc(.add, null, tmp2, tmp1[0..llnormalize(tmp1)], a);
4584	mem.swap([]Limb, &tmp1, &tmp2);
4585	}
4586	}
4587	}
4588
4589	// Storage must live for the lifetime of the returned value
4590	fn fixedIntFromSignedDoubleLimb(A: SignedDoubleLimb, storage: []Limb) Mutable {
4591	assert(storage.len >= 2);
4592
4593	const A_is_positive = A >= 0;
4594	const Au = @as(DoubleLimb, @intCast(if (A < 0) -A else A));
4595	storage[0] = @as(Limb, @truncate(Au));
4596	storage[1] = @as(Limb, @truncate(Au >> limb_bits));
4597	return .{
4598	.limbs = storage[0..2],
4599	.positive = A_is_positive,
4600	.len = 2,
4601	};
4602	}
4603
4604	fn slicesOverlap(a: []const Limb, b: []const Limb) bool {
4605	// there is no overlap if a.ptr + a.len <= b.ptr or b.ptr + b.len <= a.ptr
4606	return @intFromPtr(a.ptr + a.len) > @intFromPtr(b.ptr) and @intFromPtr(b.ptr + b.len) > @intFromPtr(a.ptr);
4607	}
4608
4609	test {
4610	_ = @import("int_test.zig");
4611	}
4612
4613	const testing_allocator = std.testing.allocator;
4614	test "llshl shift by whole number of limb" {
4615	const padding = maxInt(Limb);
4616
4617	var r: [10]Limb = @splat(padding);
4618
4619	const A: Limb = @truncate(0xCCCCCCCCCCCCCCCCCCCCCCC);
4620	const B: Limb = @truncate(0x22222222222222222222222);
4621
4622	const data = [2]Limb{ A, B };
4623	for (0..9) |i| {
4624	@memset(&r, padding);
4625	const len = llshl(&r, &data, i * @bitSizeOf(Limb));
4626
4627	try std.testing.expectEqual(i + 2, len);
4628	try std.testing.expectEqualSlices(Limb, &data, r[i .. i + 2]);
4629	for (r[0..i]) |x|
4630	try std.testing.expectEqual(0, x);
4631	for (r[i + 2 ..]) |x|
4632	try std.testing.expectEqual(padding, x);
4633	}
4634	}
4635
4636	test llshl {
4637	if (limb_bits != 64) return error.SkipZigTest;
4638
4639	// 1 << 63
4640	const left_one = 0x8000000000000000;
4641	const maxint: Limb = 0xFFFFFFFFFFFFFFFF;
4642
4643	// zig fmt: off
4644	try testOneShiftCase(.llshl, .{0, &.{0}, &.{0}});
4645	try testOneShiftCase(.llshl, .{0, &.{1}, &.{1}});
4646	try testOneShiftCase(.llshl, .{0, &.{125484842448}, &.{125484842448}});
4647	try testOneShiftCase(.llshl, .{0, &.{0xdeadbeef}, &.{0xdeadbeef}});
4648	try testOneShiftCase(.llshl, .{0, &.{maxint}, &.{maxint}});
4649	try testOneShiftCase(.llshl, .{0, &.{left_one}, &.{left_one}});
4650	try testOneShiftCase(.llshl, .{0, &.{0, 1}, &.{0, 1}});
4651	try testOneShiftCase(.llshl, .{0, &.{1, 2}, &.{1, 2}});
4652	try testOneShiftCase(.llshl, .{0, &.{left_one, 1}, &.{left_one, 1}});
4653	try testOneShiftCase(.llshl, .{1, &.{0}, &.{0}});
4654	try testOneShiftCase(.llshl, .{1, &.{2}, &.{1}});
4655	try testOneShiftCase(.llshl, .{1, &.{250969684896}, &.{125484842448}});
4656	try testOneShiftCase(.llshl, .{1, &.{0x1bd5b7dde}, &.{0xdeadbeef}});
4657	try testOneShiftCase(.llshl, .{1, &.{0xfffffffffffffffe, 1}, &.{maxint}});
4658	try testOneShiftCase(.llshl, .{1, &.{0, 1}, &.{left_one}});
4659	try testOneShiftCase(.llshl, .{1, &.{0, 2}, &.{0, 1}});
4660	try testOneShiftCase(.llshl, .{1, &.{2, 4}, &.{1, 2}});
4661	try testOneShiftCase(.llshl, .{1, &.{0, 3}, &.{left_one, 1}});
4662	try testOneShiftCase(.llshl, .{5, &.{32}, &.{1}});
4663	try testOneShiftCase(.llshl, .{5, &.{4015514958336}, &.{125484842448}});
4664	try testOneShiftCase(.llshl, .{5, &.{0x1bd5b7dde0}, &.{0xdeadbeef}});
4665	try testOneShiftCase(.llshl, .{5, &.{0xffffffffffffffe0, 0x1f}, &.{maxint}});
4666	try testOneShiftCase(.llshl, .{5, &.{0, 16}, &.{left_one}});
4667	try testOneShiftCase(.llshl, .{5, &.{0, 32}, &.{0, 1}});
4668	try testOneShiftCase(.llshl, .{5, &.{32, 64}, &.{1, 2}});
4669	try testOneShiftCase(.llshl, .{5, &.{0, 48}, &.{left_one, 1}});
4670	try testOneShiftCase(.llshl, .{64, &.{0, 1}, &.{1}});
4671	try testOneShiftCase(.llshl, .{64, &.{0, 125484842448}, &.{125484842448}});
4672	try testOneShiftCase(.llshl, .{64, &.{0, 0xdeadbeef}, &.{0xdeadbeef}});
4673	try testOneShiftCase(.llshl, .{64, &.{0, maxint}, &.{maxint}});
4674	try testOneShiftCase(.llshl, .{64, &.{0, left_one}, &.{left_one}});
4675	try testOneShiftCase(.llshl, .{64, &.{0, 0, 1}, &.{0, 1}});
4676	try testOneShiftCase(.llshl, .{64, &.{0, 1, 2}, &.{1, 2}});
4677	try testOneShiftCase(.llshl, .{64, &.{0, left_one, 1}, &.{left_one, 1}});
4678	try testOneShiftCase(.llshl, .{35, &.{0x800000000}, &.{1}});
4679	try testOneShiftCase(.llshl, .{35, &.{13534986488655118336, 233}, &.{125484842448}});
4680	try testOneShiftCase(.llshl, .{35, &.{0xf56df77800000000, 6}, &.{0xdeadbeef}});
4681	try testOneShiftCase(.llshl, .{35, &.{0xfffffff800000000, 0x7ffffffff}, &.{maxint}});
4682	try testOneShiftCase(.llshl, .{35, &.{0, 17179869184}, &.{left_one}});
4683	try testOneShiftCase(.llshl, .{35, &.{0, 0x800000000}, &.{0, 1}});
4684	try testOneShiftCase(.llshl, .{35, &.{0x800000000, 0x1000000000}, &.{1, 2}});
4685	try testOneShiftCase(.llshl, .{35, &.{0, 0xc00000000}, &.{left_one, 1}});
4686	try testOneShiftCase(.llshl, .{70, &.{0, 64}, &.{1}});
4687	try testOneShiftCase(.llshl, .{70, &.{0, 8031029916672}, &.{125484842448}});
4688	try testOneShiftCase(.llshl, .{70, &.{0, 0x37ab6fbbc0}, &.{0xdeadbeef}});
4689	try testOneShiftCase(.llshl, .{70, &.{0, 0xffffffffffffffc0, 63}, &.{maxint}});
4690	try testOneShiftCase(.llshl, .{70, &.{0, 0, 32}, &.{left_one}});
4691	try testOneShiftCase(.llshl, .{70, &.{0, 0, 64}, &.{0, 1}});
4692	try testOneShiftCase(.llshl, .{70, &.{0, 64, 128}, &.{1, 2}});
4693	try testOneShiftCase(.llshl, .{70, &.{0, 0, 0x60}, &.{left_one, 1}});
4694	// zig fmt: on
4695	}
4696
4697	test "llshl shift 0" {
4698	const n = @bitSizeOf(Limb);
4699	if (n <= 20) return error.SkipZigTest;
4700
4701	// zig fmt: off
4702	try testOneShiftCase(.llshl, .{0, &.{0}, &.{0}});
4703	try testOneShiftCase(.llshl, .{1, &.{0}, &.{0}});
4704	try testOneShiftCase(.llshl, .{5, &.{0}, &.{0}});
4705	try testOneShiftCase(.llshl, .{13, &.{0}, &.{0}});
4706	try testOneShiftCase(.llshl, .{20, &.{0}, &.{0}});
4707	try testOneShiftCase(.llshl, .{0, &.{0, 0}, &.{0, 0}});
4708	try testOneShiftCase(.llshl, .{2, &.{0, 0}, &.{0, 0}});
4709	try testOneShiftCase(.llshl, .{7, &.{0, 0}, &.{0, 0}});
4710	try testOneShiftCase(.llshl, .{11, &.{0, 0}, &.{0, 0}});
4711	try testOneShiftCase(.llshl, .{19, &.{0, 0}, &.{0, 0}});
4712
4713	try testOneShiftCase(.llshl, .{0, &.{0}, &.{0}});
4714	try testOneShiftCase(.llshl, .{n, &.{0, 0}, &.{0}});
4715	try testOneShiftCase(.llshl, .{2*n, &.{0, 0, 0}, &.{0}});
4716	try testOneShiftCase(.llshl, .{3*n, &.{0, 0, 0, 0}, &.{0}});
4717	try testOneShiftCase(.llshl, .{4*n, &.{0, 0, 0, 0, 0}, &.{0}});
4718	try testOneShiftCase(.llshl, .{0, &.{0, 0}, &.{0, 0}});
4719	try testOneShiftCase(.llshl, .{n, &.{0, 0, 0}, &.{0, 0}});
4720	try testOneShiftCase(.llshl, .{2*n, &.{0, 0, 0, 0}, &.{0, 0}});
4721	try testOneShiftCase(.llshl, .{3*n, &.{0, 0, 0, 0, 0}, &.{0, 0}});
4722	try testOneShiftCase(.llshl, .{4*n, &.{0, 0, 0, 0, 0, 0}, &.{0, 0}});
4723	// zig fmt: on
4724	}
4725
4726	test "llshr shift 0" {
4727	const n = @bitSizeOf(Limb);
4728
4729	// zig fmt: off
4730	try testOneShiftCase(.llshr, .{0, &.{0}, &.{0}});
4731	try testOneShiftCase(.llshr, .{1, &.{0}, &.{0}});
4732	try testOneShiftCase(.llshr, .{5, &.{0}, &.{0}});
4733	try testOneShiftCase(.llshr, .{13, &.{0}, &.{0}});
4734	try testOneShiftCase(.llshr, .{20, &.{0}, &.{0}});
4735	try testOneShiftCase(.llshr, .{0, &.{0, 0}, &.{0, 0}});
4736	try testOneShiftCase(.llshr, .{2, &.{0}, &.{0, 0}});
4737	try testOneShiftCase(.llshr, .{7, &.{0}, &.{0, 0}});
4738	try testOneShiftCase(.llshr, .{11, &.{0}, &.{0, 0}});
4739	try testOneShiftCase(.llshr, .{19, &.{0}, &.{0, 0}});
4740
4741	try testOneShiftCase(.llshr, .{n, &.{0}, &.{0}});
4742	try testOneShiftCase(.llshr, .{2*n, &.{0}, &.{0}});
4743	try testOneShiftCase(.llshr, .{3*n, &.{0}, &.{0}});
4744	try testOneShiftCase(.llshr, .{4*n, &.{0}, &.{0}});
4745	try testOneShiftCase(.llshr, .{n, &.{0}, &.{0, 0}});
4746	try testOneShiftCase(.llshr, .{2*n, &.{0}, &.{0, 0}});
4747	try testOneShiftCase(.llshr, .{3*n, &.{0}, &.{0, 0}});
4748	try testOneShiftCase(.llshr, .{4*n, &.{0}, &.{0, 0}});
4749
4750	try testOneShiftCase(.llshr, .{1, &.{}, &.{}});
4751	try testOneShiftCase(.llshr, .{2, &.{}, &.{}});
4752	try testOneShiftCase(.llshr, .{64, &.{}, &.{}});
4753	// zig fmt: on
4754	}
4755
4756	test "llshr to 0" {
4757	const n = @bitSizeOf(Limb);
4758	if (n != 64 and n != 32) return error.SkipZigTest;
4759
4760	// zig fmt: off
4761	try testOneShiftCase(.llshr, .{1, &.{0}, &.{0}});
4762	try testOneShiftCase(.llshr, .{1, &.{0}, &.{1}});
4763	try testOneShiftCase(.llshr, .{5, &.{0}, &.{1}});
4764	try testOneShiftCase(.llshr, .{65, &.{0}, &.{0, 1}});
4765	try testOneShiftCase(.llshr, .{193, &.{0}, &.{0, 0, maxInt(Limb)}});
4766	try testOneShiftCase(.llshr, .{193, &.{0}, &.{maxInt(Limb), 1, maxInt(Limb)}});
4767	try testOneShiftCase(.llshr, .{193, &.{0}, &.{0xdeadbeef, 0xabcdefab, 0x1234}});
4768	// zig fmt: on
4769	}
4770
4771	test "llshr single" {
4772	if (limb_bits != 64) return error.SkipZigTest;
4773
4774	// 1 << 63
4775	const left_one = 0x8000000000000000;
4776	const maxint: Limb = 0xFFFFFFFFFFFFFFFF;
4777
4778	// zig fmt: off
4779	try testOneShiftCase(.llshr, .{0, &.{0}, &.{0}});
4780	try testOneShiftCase(.llshr, .{0, &.{1}, &.{1}});
4781	try testOneShiftCase(.llshr, .{0, &.{125484842448}, &.{125484842448}});
4782	try testOneShiftCase(.llshr, .{0, &.{0xdeadbeef}, &.{0xdeadbeef}});
4783	try testOneShiftCase(.llshr, .{0, &.{maxint}, &.{maxint}});
4784	try testOneShiftCase(.llshr, .{0, &.{left_one}, &.{left_one}});
4785	try testOneShiftCase(.llshr, .{1, &.{0}, &.{0}});
4786	try testOneShiftCase(.llshr, .{1, &.{1}, &.{2}});
4787	try testOneShiftCase(.llshr, .{1, &.{62742421224}, &.{125484842448}});
4788	try testOneShiftCase(.llshr, .{1, &.{62742421223}, &.{125484842447}});
4789	try testOneShiftCase(.llshr, .{1, &.{0x6f56df77}, &.{0xdeadbeef}});
4790	try testOneShiftCase(.llshr, .{1, &.{0x7fffffffffffffff}, &.{maxint}});
4791	try testOneShiftCase(.llshr, .{1, &.{0x4000000000000000}, &.{left_one}});
4792	try testOneShiftCase(.llshr, .{8, &.{1}, &.{256}});
4793	try testOneShiftCase(.llshr, .{8, &.{490175165}, &.{125484842448}});
4794	try testOneShiftCase(.llshr, .{8, &.{0xdeadbe}, &.{0xdeadbeef}});
4795	try testOneShiftCase(.llshr, .{8, &.{0xffffffffffffff}, &.{maxint}});
4796	try testOneShiftCase(.llshr, .{8, &.{0x80000000000000}, &.{left_one}});
4797	// zig fmt: on
4798	}
4799
4800	test llshr {
4801	if (limb_bits != 64) return error.SkipZigTest;
4802
4803	// 1 << 63
4804	const left_one = 0x8000000000000000;
4805	const maxint: Limb = 0xFFFFFFFFFFFFFFFF;
4806
4807	// zig fmt: off
4808	try testOneShiftCase(.llshr, .{0, &.{0, 0}, &.{0, 0}});
4809	try testOneShiftCase(.llshr, .{0, &.{0, 1}, &.{0, 1}});
4810	try testOneShiftCase(.llshr, .{0, &.{15, 1}, &.{15, 1}});
4811	try testOneShiftCase(.llshr, .{0, &.{987656565, 123456789456}, &.{987656565, 123456789456}});
4812	try testOneShiftCase(.llshr, .{0, &.{0xfeebdaed, 0xdeadbeef}, &.{0xfeebdaed, 0xdeadbeef}});
4813	try testOneShiftCase(.llshr, .{0, &.{1, maxint}, &.{1, maxint}});
4814	try testOneShiftCase(.llshr, .{0, &.{0, left_one}, &.{0, left_one}});
4815	try testOneShiftCase(.llshr, .{1, &.{0}, &.{0, 0}});
4816	try testOneShiftCase(.llshr, .{1, &.{left_one}, &.{0, 1}});
4817	try testOneShiftCase(.llshr, .{1, &.{0x8000000000000007}, &.{15, 1}});
4818	try testOneShiftCase(.llshr, .{1, &.{493828282, 61728394728}, &.{987656565, 123456789456}});
4819	try testOneShiftCase(.llshr, .{1, &.{0x800000007f75ed76, 0x6f56df77}, &.{0xfeebdaed, 0xdeadbeef}});
4820	try testOneShiftCase(.llshr, .{1, &.{left_one, 0x7fffffffffffffff}, &.{1, maxint}});
4821	try testOneShiftCase(.llshr, .{1, &.{0, 0x4000000000000000}, &.{0, left_one}});
4822	try testOneShiftCase(.llshr, .{64, &.{0}, &.{0, 0}});
4823	try testOneShiftCase(.llshr, .{64, &.{1}, &.{0, 1}});
4824	try testOneShiftCase(.llshr, .{64, &.{1}, &.{15, 1}});
4825	try testOneShiftCase(.llshr, .{64, &.{123456789456}, &.{987656565, 123456789456}});
4826	try testOneShiftCase(.llshr, .{64, &.{0xdeadbeef}, &.{0xfeebdaed, 0xdeadbeef}});
4827	try testOneShiftCase(.llshr, .{64, &.{maxint}, &.{1, maxint}});
4828	try testOneShiftCase(.llshr, .{64, &.{left_one}, &.{0, left_one}});
4829	try testOneShiftCase(.llshr, .{72, &.{0}, &.{0, 0}});
4830	try testOneShiftCase(.llshr, .{72, &.{0}, &.{0, 1}});
4831	try testOneShiftCase(.llshr, .{72, &.{0}, &.{15, 1}});
4832	try testOneShiftCase(.llshr, .{72, &.{482253083}, &.{987656565, 123456789456}});
4833	try testOneShiftCase(.llshr, .{72, &.{0xdeadbe}, &.{0xfeebdaed, 0xdeadbeef}});
4834	try testOneShiftCase(.llshr, .{72, &.{0xffffffffffffff}, &.{1, maxint}});
4835	try testOneShiftCase(.llshr, .{72, &.{0x80000000000000}, &.{0, left_one}});
4836	// zig fmt: on
4837	}
4838
4839	const Case = struct { usize, []const Limb, []const Limb };
4840
4841	fn testOneShiftCase(comptime function: enum { llshr, llshl }, case: Case) !void {
4842	const func = if (function == .llshl) llshl else llshr;
4843	const shift_direction = if (function == .llshl) -1 else 1;
4844
4845	try testOneShiftCaseNoAliasing(func, case);
4846	try testOneShiftCaseAliasing(func, case, shift_direction);
4847	}
4848
4849	fn testOneShiftCaseNoAliasing(func: fn ([]Limb, []const Limb, usize) usize, case: Case) !void {
4850	const padding = maxInt(Limb);
4851	var r: [20]Limb = @splat(padding);
4852
4853	const shift = case[0];
4854	const expected = case[1];
4855	const data = case[2];
4856
4857	std.debug.assert(expected.len <= 20);
4858
4859	const len = func(&r, data, shift);
4860
4861	try std.testing.expectEqual(expected.len, len);
4862	try std.testing.expectEqualSlices(Limb, expected, r[0..len]);
4863	try std.testing.expect(mem.allEqual(Limb, r[len..], padding));
4864	}
4865
4866	fn testOneShiftCaseAliasing(func: fn ([]Limb, []const Limb, usize) usize, case: Case, shift_direction: isize) !void {
4867	const padding = maxInt(Limb);
4868	var r: [60]Limb = @splat(padding);
4869	const base = 20;
4870
4871	assert(shift_direction == 1 or shift_direction == -1);
4872
4873	for (0..10) |limb_shift| {
4874	const shift = case[0];
4875	const expected = case[1];
4876	const data = case[2];
4877
4878	std.debug.assert(expected.len <= 20);
4879
4880	@memset(&r, padding);
4881	const final_limb_base: usize = @intCast(base + shift_direction * @as(isize, @intCast(limb_shift)));
4882	const written_data = r[final_limb_base..][0..data.len];
4883	@memcpy(written_data, data);
4884
4885	const len = func(r[base..], written_data, shift);
4886
4887	try std.testing.expectEqual(expected.len, len);
4888	try std.testing.expectEqualSlices(Limb, expected, r[base .. base + len]);
4889	}
4890	}
4891
4892	test "format" {
4893	var a: Managed = try .init(std.testing.allocator);
4894	defer a.deinit();
4895
4896	try a.set(123);
4897	try testFormat(a, "123");
4898
4899	try a.set(-123);
4900	try testFormat(a, "-123");
4901
4902	try a.set(20000000000000000000); // > maxInt(u64)
4903	try testFormat(a, "20000000000000000000");
4904
4905	try a.set(1 << 64 * @sizeOf(usize) * 8);
4906	try testFormat(a, "(BigInt)");
4907
4908	try a.set(-(1 << 64 * @sizeOf(usize) * 8));
4909	try testFormat(a, "(BigInt)");
4910	}
4911
4912	fn testFormat(a: Managed, expected: []const u8) !void {
4913	try std.testing.expectFmt(expected, "{f}", .{a});
4914	try std.testing.expectFmt(expected, "{f}", .{a.toMutable()});
4915	try std.testing.expectFmt(expected, "{f}", .{a.toConst()});
4916	}