A arbitrary-precision big integer, with a fixed set of immutable limbs.
pub const Const = structpub const Const = struct {
/// Raw digits. These are:
///
/// * Little-endian ordered
/// * limbs.len >= 1
/// * Zero is represented as limbs.len == 1 with limbs[0] == 0.
///
/// Accessing limbs directly should be avoided.
limbs: []const Limb,
positive: bool,
/// The result is an independent resource which is managed by the caller.
pub fn toManaged(self: Const, allocator: Allocator) Allocator.Error!Managed {
const limbs = try allocator.alloc(Limb, @max(Managed.default_capacity, self.limbs.len));
@memcpy(limbs[0..self.limbs.len], self.limbs);
return .{
.allocator = allocator,
.limbs = limbs,
.metadata = if (self.positive)
self.limbs.len & ~Managed.sign_bit
else
self.limbs.len | Managed.sign_bit,
};
}
/// Asserts `limbs` is big enough to store the value.
pub fn toMutable(self: Const, limbs: []Limb) Mutable {
@memcpy(limbs[0..self.limbs.len], self.limbs[0..self.limbs.len]);
return .{
.limbs = limbs,
.positive = self.positive,
.len = self.limbs.len,
};
}
pub fn dump(self: Const) void {
for (self.limbs[0..self.limbs.len]) |limb| {
std.debug.print("{x} ", .{limb});
}
std.debug.print("len={} positive={}\n", .{ self.limbs.len, self.positive });
}
pub fn abs(self: Const) Const {
return .{
.limbs = self.limbs,
.positive = true,
};
}
pub fn negate(self: Const) Const {
return .{
.limbs = self.limbs,
.positive = !self.positive,
};
}
pub fn isOdd(self: Const) bool {
return self.limbs[0] & 1 != 0;
}
pub fn isEven(self: Const) bool {
return !self.isOdd();
}
/// Returns the number of bits required to represent the absolute value of an integer.
pub fn bitCountAbs(self: Const) usize {
return (self.limbs.len - 1) * limb_bits + (limb_bits - @clz(self.limbs[self.limbs.len - 1]));
}
/// Returns the number of bits required to represent the integer in twos-complement form.
///
/// If the integer is negative the value returned is the number of bits needed by a signed
/// integer to represent the value. If positive the value is the number of bits for an
/// unsigned integer. Any unsigned integer will fit in the signed integer with bitcount
/// one greater than the returned value.
///
/// e.g. -127 returns 8 as it will fit in an i8. 127 returns 7 since it fits in a u7.
pub fn bitCountTwosComp(self: Const) usize {
var bits = self.bitCountAbs();
// If the entire value has only one bit set (e.g. 0b100000000) then the negation in twos
// complement requires one less bit.
if (!self.positive) block: {
bits += 1;
if (@popCount(self.limbs[self.limbs.len - 1]) == 1) {
for (self.limbs[0 .. self.limbs.len - 1]) |limb| {
if (@popCount(limb) != 0) {
break :block;
}
}
bits -= 1;
}
}
return bits;
}
/// Returns the number of bits required to represent the integer in twos-complement form
/// with the given signedness.
pub fn bitCountTwosCompForSignedness(self: Const, signedness: std.builtin.Signedness) usize {
return self.bitCountTwosComp() + @intFromBool(self.positive and signedness == .signed);
}
/// @popCount with two's complement semantics.
///
/// This returns the number of 1 bits set when the value would be represented in
/// two's complement with the given integer width (bit_count).
/// This includes the leading sign bit, which will be set for negative values.
///
/// Asserts that bit_count is enough to represent value in two's compliment
/// and that the final result fits in a usize.
/// Asserts that there are no trailing empty limbs on the most significant end,
/// i.e. that limb count matches `calcLimbLen()` and zero is not negative.
pub fn popCount(self: Const, bit_count: usize) usize {
var sum: usize = 0;
if (self.positive) {
for (self.limbs) |limb| {
sum += @popCount(limb);
}
} else {
assert(self.fitsInTwosComp(.signed, bit_count));
assert(self.limbs[self.limbs.len - 1] != 0);
var remaining_bits = bit_count;
var carry: u1 = 1;
var add_res: Limb = undefined;
// All but the most significant limb.
for (self.limbs[0 .. self.limbs.len - 1]) |limb| {
const ov = @addWithOverflow(~limb, carry);
add_res = ov[0];
carry = ov[1];
sum += @popCount(add_res);
remaining_bits -= limb_bits; // Asserted not to underflow by fitsInTwosComp
}
// The most significant limb may have fewer than @bitSizeOf(Limb) meaningful bits,
// which we can detect with @clz().
// There may also be fewer limbs than needed to fill bit_count.
const limb = self.limbs[self.limbs.len - 1];
const leading_zeroes = @clz(limb);
// The most significant limb is asserted not to be all 0s (above),
// so ~limb cannot be all 1s, and ~limb + 1 cannot overflow.
sum += @popCount(~limb + carry);
sum -= leading_zeroes; // All leading zeroes were flipped and added to sum, so undo those
const remaining_ones = remaining_bits - (limb_bits - leading_zeroes); // All bits not covered by limbs
sum += remaining_ones;
}
return sum;
}
pub fn fitsInTwosComp(self: Const, signedness: Signedness, bit_count: usize) bool {
if (self.eqlZero()) {
return true;
}
if (signedness == .unsigned and !self.positive) {
return false;
}
return bit_count >= self.bitCountTwosCompForSignedness(signedness);
}
/// Returns whether self can fit into an integer of the requested type.
pub fn fits(self: Const, comptime T: type) bool {
const info = @typeInfo(T).int;
return self.fitsInTwosComp(info.signedness, info.bits);
}
/// Returns the approximate size of the integer in the given base. Negative values accommodate for
/// the minus sign. This is used for determining the number of characters needed to print the
/// value. It is inexact and may exceed the given value by ~1-2 bytes.
/// TODO See if we can make this exact.
pub fn sizeInBaseUpperBound(self: Const, base: usize) usize {
const bit_count = @as(usize, @intFromBool(!self.positive)) + self.bitCountAbs();
return (bit_count / math.log2(base)) + 2;
}
pub const ConvertError = error{
NegativeIntoUnsigned,
TargetTooSmall,
};
/// Convert `self` to `Int`.
///
/// Returns an error if self cannot be narrowed into the requested type without truncation.
pub fn toInt(self: Const, comptime Int: type) ConvertError!Int {
switch (@typeInfo(Int)) {
.int => |info| {
// Make sure -0 is handled correctly.
if (self.eqlZero()) return 0;
const Unsigned = @Int(.unsigned, info.bits);
if (!self.fitsInTwosComp(info.signedness, info.bits)) {
return error.TargetTooSmall;
}
var r: Unsigned = 0;
if (@sizeOf(Unsigned) <= @sizeOf(Limb)) {
r = @intCast(self.limbs[0]);
} else {
for (self.limbs[0..self.limbs.len], 0..) |_, ri| {
const limb = self.limbs[self.limbs.len - ri - 1];
r <<= limb_bits;
r |= limb;
}
}
if (info.signedness == .unsigned) {
return if (self.positive) @intCast(r) else error.NegativeIntoUnsigned;
} else {
if (self.positive) {
return @intCast(r);
} else {
if (math.cast(Int, r)) |ok| {
return -ok;
} else {
return minInt(Int);
}
}
}
},
else => @compileError("expected int type, found '" ++ @typeName(Int) ++ "'"),
}
}
/// Convert self to `Float`.
pub fn toFloat(self: Const, comptime Float: type, round: Round) struct { Float, Exactness } {
if (Float == comptime_float) return self.toFloat(f128, round);
const normalized_abs: Const = .{
.limbs = self.limbs[0..llnormalize(self.limbs)],
.positive = true,
};
if (normalized_abs.eqlZero()) return .{ if (self.positive) 0.0 else -0.0, .exact };
const Repr = std.math.FloatRepr(Float);
var mantissa_limbs: [calcNonZeroTwosCompLimbCount(1 + @bitSizeOf(Repr.Mantissa))]Limb = undefined;
var mantissa: Mutable = .{
.limbs = &mantissa_limbs,
.positive = undefined,
.len = undefined,
};
var exponent = normalized_abs.bitCountAbs() - 1;
const exactness: Exactness = exactness: {
if (exponent <= @bitSizeOf(Repr.Normalized.Fraction)) {
mantissa.shiftLeft(normalized_abs, @intCast(@bitSizeOf(Repr.Normalized.Fraction) - exponent));
break :exactness .exact;
}
const shift: usize = @intCast(exponent - @bitSizeOf(Repr.Normalized.Fraction));
mantissa.shiftRight(normalized_abs, shift);
const final_limb_index = (shift - 1) / limb_bits;
const round_bits = normalized_abs.limbs[final_limb_index] << @truncate(-%shift) |
@intFromBool(!std.mem.allEqual(Limb, normalized_abs.limbs[0..final_limb_index], 0));
if (round_bits == 0) break :exactness .exact;
round: switch (round) {
.nearest_even => {
const half: Limb = 1 << (limb_bits - 1);
if (round_bits >= half) mantissa.addScalar(mantissa.toConst(), 1);
if (round_bits == half) mantissa.limbs[0] &= ~@as(Limb, 1);
},
.away => mantissa.addScalar(mantissa.toConst(), 1),
.trunc => {},
.floor => if (!self.positive) continue :round .away,
.ceil => if (self.positive) continue :round .away,
}
break :exactness .inexact;
};
const normalized_res: Repr.Normalized = .{
.fraction = @truncate(mantissa.toInt(Repr.Mantissa) catch |err| switch (err) {
error.NegativeIntoUnsigned => unreachable,
error.TargetTooSmall => fraction: {
assert(mantissa.toConst().orderAgainstScalar(1 << @bitSizeOf(Repr.Mantissa)).compare(.eq));
exponent += 1;
break :fraction 1 << (@bitSizeOf(Repr.Mantissa) - 1);
},
}),
.exponent = std.math.lossyCast(Repr.Normalized.Exponent, exponent),
};
return .{ normalized_res.reconstruct(if (self.positive) .positive else .negative), exactness };
}
pub fn format(self: Const, w: *std.Io.Writer) std.Io.Writer.Error!void {
return self.formatNumber(w, .{});
}
/// If the absolute value of integer is greater than or equal to `pow(2, 64 * @sizeOf(usize) * 8)`,
/// this function will fail to print the string, printing "(BigInt)" instead of a number.
/// This is because the rendering algorithm requires reversing a string, which requires O(N) memory.
/// See `toString` and `toStringAlloc` for a way to print big integers without failure.
pub fn formatNumber(self: Const, w: *std.Io.Writer, number: std.fmt.Number) std.Io.Writer.Error!void {
const available_len = 64;
if (self.limbs.len > available_len)
return w.writeAll("(BigInt)");
var limbs: [calcToStringLimbsBufferLen(available_len, 10)]Limb = undefined;
const biggest: Const = .{
.limbs = &@as([available_len]Limb, @splat(comptime math.maxInt(Limb))),
.positive = false,
};
var buf: [biggest.sizeInBaseUpperBound(2)]u8 = undefined;
const base: u8 = number.mode.base() orelse @panic("TODO print big int in scientific form");
const len = self.toString(&buf, base, number.case, &limbs);
return w.writeAll(buf[0..len]);
}
/// Converts self to a string in the requested base.
/// Caller owns returned memory.
/// Asserts that `base` is in the range [2, 36].
/// See also `toString`, a lower level function than this.
pub fn toStringAlloc(self: Const, allocator: Allocator, base: u8, case: std.fmt.Case) Allocator.Error![]u8 {
assert(base >= 2);
assert(base <= 36);
if (self.eqlZero()) {
return allocator.dupe(u8, "0");
}
const string = try allocator.alloc(u8, self.sizeInBaseUpperBound(base));
errdefer allocator.free(string);
const limbs = try allocator.alloc(Limb, calcToStringLimbsBufferLen(self.limbs.len, base));
defer allocator.free(limbs);
return allocator.realloc(string, self.toString(string, base, case, limbs));
}
/// Converts self to a string in the requested base.
/// Asserts that `base` is in the range [2, 36].
/// `string` is a caller-provided slice of at least `sizeInBaseUpperBound` bytes,
/// where the result is written to.
/// Returns the length of the string.
/// `limbs_buffer` is caller-provided memory for `toString` to use as a working area. It must have
/// length of at least `calcToStringLimbsBufferLen`.
/// In the case of power-of-two base, `limbs_buffer` is ignored.
/// See also `toStringAlloc`, a higher level function than this.
pub fn toString(self: Const, string: []u8, base: u8, case: std.fmt.Case, limbs_buffer: []Limb) usize {
assert(base >= 2);
assert(base <= 36);
if (self.eqlZero()) {
string[0] = '0';
return 1;
}
var digits_len: usize = 0;
// Power of two: can do a single pass and use masks to extract digits.
if (math.isPowerOfTwo(base)) {
const base_shift = math.log2_int(Limb, base);
outer: for (self.limbs[0..self.limbs.len]) |limb| {
var shift: usize = 0;
while (shift < limb_bits) : (shift += base_shift) {
const r = @as(u8, @intCast((limb >> @as(Log2Limb, @intCast(shift))) & @as(Limb, base - 1)));
const ch = std.fmt.digitToChar(r, case);
string[digits_len] = ch;
digits_len += 1;
// If we hit the end, it must be all zeroes from here.
if (digits_len == string.len) break :outer;
}
}
// Always will have a non-zero digit somewhere.
while (string[digits_len - 1] == '0') {
digits_len -= 1;
}
} else {
// Non power-of-two: batch divisions per word size.
// We use a HalfLimb here so the division uses the faster lldiv0p5 over lldiv1 codepath.
const digits_per_limb = math.log(HalfLimb, base, maxInt(HalfLimb));
var limb_base: Limb = 1;
var j: usize = 0;
while (j < digits_per_limb) : (j += 1) {
limb_base *= base;
}
const b: Const = .{ .limbs = &[_]Limb{limb_base}, .positive = true };
var q: Mutable = .{
.limbs = limbs_buffer[0 .. self.limbs.len + 2],
.positive = true, // Make absolute by ignoring self.positive.
.len = self.limbs.len,
};
@memcpy(q.limbs[0..self.limbs.len], self.limbs);
var r: Mutable = .{
.limbs = limbs_buffer[q.limbs.len..][0..self.limbs.len],
.positive = true,
.len = 1,
};
r.limbs[0] = 0;
const rest_of_the_limbs_buf = limbs_buffer[q.limbs.len + r.limbs.len ..];
while (q.len >= 2) {
// Passing an allocator here would not be helpful since this division is destroying
// information, not creating it. [TODO citation needed]
q.divTrunc(&r, q.toConst(), b, rest_of_the_limbs_buf);
var r_word = r.limbs[0];
var i: usize = 0;
while (i < digits_per_limb) : (i += 1) {
const ch = std.fmt.digitToChar(@as(u8, @intCast(r_word % base)), case);
r_word /= base;
string[digits_len] = ch;
digits_len += 1;
}
}
{
assert(q.len == 1);
var r_word = q.limbs[0];
while (r_word != 0) {
const ch = std.fmt.digitToChar(@as(u8, @intCast(r_word % base)), case);
r_word /= base;
string[digits_len] = ch;
digits_len += 1;
}
}
}
if (!self.positive) {
string[digits_len] = '-';
digits_len += 1;
}
const s = string[0..digits_len];
mem.reverse(u8, s);
return s.len;
}
/// Write the value of `x` into `buffer`
/// Asserts that `buffer` is large enough to store the value.
///
/// `buffer` is filled so that its contents match what would be observed via
/// @ptrCast(*[buffer.len]const u8, &x). Byte ordering is determined by `endian`,
/// and any required padding bits are added on the MSB end.
pub fn writeTwosComplement(x: Const, buffer: []u8, endian: Endian) void {
return writePackedTwosComplement(x, buffer, 0, 8 * buffer.len, endian);
}
/// Write the value of `x` to a packed memory `buffer`.
/// Asserts that `buffer` is large enough to contain a value of bit-size `bit_count`
/// at offset `bit_offset`.
///
/// This is equivalent to storing the value of an integer with `bit_count` bits as
/// if it were a field in packed memory at the provided bit offset.
pub fn writePackedTwosComplement(x: Const, buffer: []u8, bit_offset: usize, bit_count: usize, endian: Endian) void {
assert(x.fitsInTwosComp(if (x.positive) .unsigned else .signed, bit_count));
// Copy all complete limbs
var carry: u1 = 1;
var limb_index: usize = 0;
var bit_index: usize = 0;
while (limb_index < bit_count / @bitSizeOf(Limb)) : (limb_index += 1) {
var limb: Limb = if (limb_index < x.limbs.len) x.limbs[limb_index] else 0;
// 2's complement (bitwise not, then add carry bit)
if (!x.positive) {
const ov = @addWithOverflow(~limb, carry);
limb = ov[0];
carry = ov[1];
}
// Write one Limb of bits
mem.writePackedInt(Limb, buffer, bit_index + bit_offset, limb, endian);
bit_index += @bitSizeOf(Limb);
}
// Copy the remaining bits
if (bit_count != bit_index) {
var limb: Limb = if (limb_index < x.limbs.len) x.limbs[limb_index] else 0;
// 2's complement (bitwise not, then add carry bit)
if (!x.positive) limb = ~limb +% carry;
// Write all remaining bits
mem.writeVarPackedInt(buffer, bit_index + bit_offset, bit_count - bit_index, limb, endian);
}
}
/// Returns `math.Order.lt`, `math.Order.eq`, `math.Order.gt` if
/// `|a| < |b|`, `|a| == |b|`, or `|a| > |b|` respectively.
pub fn orderAbs(a: Const, b: Const) math.Order {
if (a.limbs.len < b.limbs.len) {
return .lt;
}
if (a.limbs.len > b.limbs.len) {
return .gt;
}
var i: usize = a.limbs.len - 1;
while (i != 0) : (i -= 1) {
if (a.limbs[i] != b.limbs[i]) {
break;
}
}
if (a.limbs[i] < b.limbs[i]) {
return .lt;
} else if (a.limbs[i] > b.limbs[i]) {
return .gt;
} else {
return .eq;
}
}
/// Returns `math.Order.lt`, `math.Order.eq`, `math.Order.gt` if `a < b`, `a == b` or `a > b` respectively.
pub fn order(a: Const, b: Const) math.Order {
if (a.positive != b.positive) {
if (eqlZero(a) and eqlZero(b)) {
return .eq;
} else {
return if (a.positive) .gt else .lt;
}
} else {
const r = orderAbs(a, b);
return if (a.positive) r else switch (r) {
.lt => math.Order.gt,
.eq => math.Order.eq,
.gt => math.Order.lt,
};
}
}
/// Same as `order` but the right-hand operand is a primitive integer.
pub fn orderAgainstScalar(lhs: Const, scalar: anytype) math.Order {
// Normally we could just determine the number of limbs needed with calcLimbLen,
// but that is not comptime-known when scalar is not a comptime_int. Instead, we
// use calcTwosCompLimbCount for a non-comptime_int scalar, which can be pessimistic
// in the case that scalar happens to be small in magnitude within its type, but it
// is well worth being able to use the stack and not needing an allocator passed in.
// Note that Mutable.init still sets len to calcLimbLen(scalar) in any case.
const limbs_len = comptime switch (@typeInfo(@TypeOf(scalar))) {
.comptime_int => calcLimbLen(scalar),
.int => |info| calcTwosCompLimbCount(info.bits),
else => @compileError("expected scalar to be an int"),
};
var limbs: [limbs_len]Limb = undefined;
const rhs = Mutable.init(&limbs, scalar);
return order(lhs, rhs.toConst());
}
/// Returns true if `a == 0`.
pub fn eqlZero(a: Const) bool {
var d: Limb = 0;
for (a.limbs) |limb| d |= limb;
return d == 0;
}
/// Returns true if `|a| == |b|`.
pub fn eqlAbs(a: Const, b: Const) bool {
return orderAbs(a, b) == .eq;
}
/// Returns true if `a == b`.
pub fn eql(a: Const, b: Const) bool {
return order(a, b) == .eq;
}
/// Returns the number of leading zeros in twos-complement form.
pub fn clz(a: Const, bits: Limb) Limb {
// Limbs are stored in little-endian order but we need to iterate big-endian.
if (!a.positive and !a.eqlZero()) return 0;
var total_limb_lz: Limb = 0;
var i: usize = a.limbs.len;
const bits_per_limb = @bitSizeOf(Limb);
while (i != 0) {
i -= 1;
const this_limb_lz = @clz(a.limbs[i]);
total_limb_lz += this_limb_lz;
if (this_limb_lz != bits_per_limb) break;
}
const total_limb_bits = a.limbs.len * bits_per_limb;
return total_limb_lz + bits - total_limb_bits;
}
/// Returns the number of trailing zeros in twos-complement form.
pub fn ctz(a: Const, bits: Limb) Limb {
// Limbs are stored in little-endian order. Converting a negative number to twos-complement
// flips all bits above the lowest set bit, which does not affect the trailing zero count.
if (a.eqlZero()) return bits;
var result: Limb = 0;
for (a.limbs) |limb| {
const limb_tz = @ctz(limb);
result += limb_tz;
if (limb_tz != @bitSizeOf(Limb)) break;
}
return @min(result, bits);
}
/// Calculate the base 2 logarithm, rounded down.
pub fn log2(a: Const) Limb {
assert(a.positive);
assert(!a.eqlZero());
return a.bitCountAbs() - 1;
}
/// Calculate the base 10 logarithm, rounded down.
///
/// The allocator is used to allocate a temporary buffer.
pub fn log10Alloc(a: Const, allocator: Allocator) Allocator.Error!Limb {
const limbs_buffer = try allocator.alloc(Limb, calcLog10LimbsBufferLen(a.limbs.len));
defer allocator.free(limbs_buffer);
return a.log10(limbs_buffer);
}
/// Calculate the base 10 logarithm, rounded down.
///
/// `limbs_buffer` is used for temporary storage. The amount required is given by `calcLog10LimbsBufferLen`.
pub fn log10(a: Const, limbs_buffer: []Limb) Limb {
assert(a.positive);
assert(!a.eqlZero());
const limb_base_as_bigint: Const = .{ .limbs = &.{constants.big_bases[10]}, .positive = true };
var q: Mutable = .{
.limbs = limbs_buffer[0 .. a.limbs.len + 2],
.positive = true,
.len = a.limbs.len,
};
@memcpy(q.limbs[0..a.limbs.len], a.limbs);
var remainder: Mutable = .{
.limbs = limbs_buffer[q.limbs.len..][0..a.limbs.len],
.positive = true,
.len = 1,
};
const division_buf = limbs_buffer[q.limbs.len + remainder.limbs.len ..];
var num_digits: Limb = 0;
while (q.len >= 2) {
q.divTrunc(&remainder, q.toConst(), limb_base_as_bigint, division_buf);
num_digits += constants.digits_per_limb[10];
}
var remaining_limb = q.limbs[0];
while (remaining_limb != 0) {
remaining_limb /= 10;
num_digits += 1;
}
return num_digits - 1;
}
}