Start of core:math/bigint

We have:
- `init` to create a new `Int`
- `init(from_integer)` to create a new `Int` and set it to `from_integer`.
- `set(Int, from_integer)` to set an `Int` to `from_integer`
- `add(dest, a, b)` to add `a` and `b` into `dest`.
- `sub(dest, a, b)` to subtract `b` from `a` and put the result in `dest`.

And a few helper functions, like:
- `is_zero`, `is_negative`, ...
- `grow`, `shrink`, `clear`, `zero`
This commit is contained in:
Jeroen van Rijn
2021-07-15 23:34:02 +02:00
parent 7afc367275
commit 18dda6ff9d
6 changed files with 693 additions and 0 deletions

196
core/math/bigint/basic.odin Normal file
View File

@@ -0,0 +1,196 @@
package bigint
/*
Copyright 2021 Jeroen van Rijn <nom@duclavier.com>.
Made available under Odin's BSD-2 license.
A BigInt implementation in Odin.
For the theoretical underpinnings, see Knuth's The Art of Computer Programming, Volume 2, section 4.3.
The code started out as an idiomatic source port of libTomMath, which is in the public domain, with thanks.
*/
import "core:mem"
import "core:intrinsics"
/*
High-level addition. Handles sign.
*/
add :: proc(dest, a, b: ^Int) -> (err: Error) {
dest := dest; x := a; y := b;
_panic_if_uninitialized(a); _panic_if_uninitialized(b); _panic_if_uninitialized(dest);
/*
Handle both negative or both positive.
*/
if x.sign == y.sign {
dest.sign = x.sign;
return _add(dest, x, y);
}
/*
One positive, the other negative.
Subtract the one with the greater magnitude from the other.
The result gets the sign of the one with the greater magnitude.
*/
if cmp_mag(x, y) == .Less_Than {
x, y = y, x;
}
dest.sign = x.sign;
return _sub(dest, x, y);
}
/*
Low-level addition, unsigned.
Handbook of Applied Cryptography, algorithm 14.7.
*/
_add :: proc(dest, a, b: ^Int) -> (err: Error) {
dest := dest; x := a; y := b;
_panic_if_uninitialized(a); _panic_if_uninitialized(b); _panic_if_uninitialized(dest);
old_used, min_used, max_used, i: int;
if x.used < y.used {
x, y = y, x;
}
min_used = x.used;
max_used = y.used;
old_used = dest.used;
err = grow(dest, max(max_used + 1, _DEFAULT_DIGIT_COUNT));
if err != .OK {
return err;
}
dest.used = max_used + 1;
/* Zero the carry */
carry := DIGIT(0);
for i = 0; i < min_used; i += 1 {
/*
Compute the sum one _DIGIT at a time.
dest[i] = a[i] + b[i] + carry;
*/
dest.digit[i] = x.digit[i] + y.digit[i] + carry;
/*
Compute carry
*/
carry = dest.digit[i] >> _DIGIT_BITS;
/*
Mask away carry from result digit.
*/
dest.digit[i] &= _MASK;
}
if min_used != max_used {
/*
Now copy higher words, if any, in A+B.
If A or B has more digits, add those in.
*/
for ; i < max_used; i += 1 {
dest.digit[i] = x.digit[i] + carry;
/*
Compute carry
*/
carry = dest.digit[i] >> _DIGIT_BITS;
/*
Mask away carry from result digit.
*/
dest.digit[i] &= _MASK;
}
}
/*
Add remaining carry.
*/
dest.digit[i] = carry;
zero_count := old_used - dest.used;
/*
Zero remainder.
*/
if zero_count > 0 {
mem.zero_slice(dest.digit[dest.used:][:zero_count]);
}
/*
Adjust dest.used based on leading zeroes.
*/
clamp(dest);
return .OK;
}
/*
Low-level subtraction. Assumes |a| > |b|.
Handbook of Applied Cryptography, algorithm 14.9.
*/
_sub :: proc(dest, a, b: ^Int) -> (err: Error) {
dest := dest; x := a; y := b;
_panic_if_uninitialized(a); _panic_if_uninitialized(b); _panic_if_uninitialized(dest);
for n in 0..=12 {
dest.digit[n] = DIGIT(n);
dest.used = n+1;
}
old_used := dest.used;
min_used := y.used;
max_used := a.used;
i: int;
err = grow(dest, max(max_used, _DEFAULT_DIGIT_COUNT));
if err != .OK {
return err;
}
dest.used = max_used;
borrow := DIGIT(0);
for i = 0; i < min_used; i += 1 {
dest.digit[i] = (x.digit[i] - y.digit[i] - borrow);
/*
borrow = carry bit of dest[i]
Note this saves performing an AND operation since if a carry does occur,
it will propagate all the way to the MSB.
As a result a single shift is enough to get the carry.
*/
borrow = dest.digit[i] >> (_DIGIT_BITS - 1);
/*
Clear borrow from dest[i].
*/
dest.digit[i] &= _MASK;
}
/*
Now copy higher words if any, e.g. if A has more digits than B
*/
for ; i < max_used; i += 1 {
dest.digit[i] = x.digit[i] - borrow;
/*
borrow = carry bit of dest[i]
Note this saves performing an AND operation since if a carry does occur,
it will propagate all the way to the MSB.
As a result a single shift is enough to get the carry.
*/
borrow = dest.digit[i] >> (_DIGIT_BITS - 1);
/*
Clear borrow from dest[i].
*/
dest.digit[i] &= _MASK;
}
zero_count := old_used - dest.used;
/*
Zero remainder.
*/
if zero_count > 0 {
mem.zero_slice(dest.digit[dest.used:][:zero_count]);
}
/*
Adjust dest.used based on leading zeroes.
*/
clamp(dest);
return .OK;
}

View File

@@ -0,0 +1,115 @@
package bigint
/*
Copyright 2021 Jeroen van Rijn <nom@duclavier.com>.
Made available under Odin's BSD-2 license.
A BigInt implementation in Odin.
For the theoretical underpinnings, see Knuth's The Art of Computer Programming, Volume 2, section 4.3.
The code started out as an idiomatic source port of libTomMath, which is in the public domain, with thanks.
*/
/*
Tunables
*/
_LOW_MEMORY :: #config(BIGINT_SMALL_MEMORY, false);
when _LOW_MEMORY {
_DEFAULT_DIGIT_COUNT :: 8;
} else {
_DEFAULT_DIGIT_COUNT :: 32;
}
// /* tunable cutoffs */
// #ifndef MP_FIXED_CUTOFFS
// extern int
// MP_MUL_KARATSUBA_CUTOFF,
// MP_SQR_KARATSUBA_CUTOFF,
// MP_MUL_TOOM_CUTOFF,
// MP_SQR_TOOM_CUTOFF;
// #endif
Sign :: enum u8 {
Zero_or_Positive = 0,
Negative = 1,
};
Int :: struct {
used: int,
allocated: int,
sign: Sign,
digit: [dynamic]DIGIT,
};
Comparison_Flag :: enum i8 {
Less_Than = -1,
Equal = 0,
Greater_Than = 1,
/* One of the numbers was uninitialized */
Uninitialized = -127,
};
Error :: enum i8 {
OK = 0,
Unknown_Error = -1,
Out_of_Memory = -2,
Invalid_Input = -3,
Max_Iterations_Reached = -4,
Buffer_Overflow = -5,
Integer_Overflow = -6,
Unimplemented = -127,
};
Primality_Flag :: enum u8 {
Blum_Blum_Shub = 0, /* BBS style prime */
Safe = 1, /* Safe prime (p-1)/2 == prime */
Second_MSB_On = 3, /* force 2nd MSB to 1 */
};
Primality_Flags :: bit_set[Primality_Flag; u8];
/*
How do we store the Ints?
Minimum number of available digits in `Int`, `_DEFAULT_DIGIT_COUNT` >= `_MIN_DIGIT_COUNT`
- Must be at least 3 for `_div_school`.
- Must be large enough such that `init_integer` can store `u128` in the `Int` without growing.
*/
_MIN_DIGIT_COUNT :: max(3, ((size_of(u128) + _DIGIT_BITS) - 1) / _DIGIT_BITS);
#assert(_DEFAULT_DIGIT_COUNT >= _MIN_DIGIT_COUNT);
/*
Maximum number of digits.
- Must be small enough such that `_bit_count` does not overflow.
- Must be small enough such that `_radix_size` for base 2 does not overflow.
`_radix_size` needs two additional bytes for zero termination and sign.
*/
_MAX_DIGIT_COUNT :: (max(int) - 2) / _DIGIT_BITS;
when size_of(rawptr) == 8 {
/*
We can use u128 as an intermediary.
*/
DIGIT :: distinct(u64);
_DIGIT_BITS :: 60;
_WORD :: u128;
} else {
DIGIT :: distinct(u32);
_DIGIT_BITS :: 28;
_WORD :: u64;
}
#assert(size_of(_WORD) == 2 * size_of(DIGIT));
_MASK :: (DIGIT(1) << DIGIT(_DIGIT_BITS)) - DIGIT(1);
_DIGIT_MAX :: _MASK;
Order :: enum i8 {
LSB_First = -1,
MSB_First = 1,
};
Endianness :: enum i8 {
Little = -1,
Platform = 0,
Big = 1,
};

View File

@@ -0,0 +1,2 @@
@echo off
odin run . -vet

View File

@@ -0,0 +1,98 @@
package bigint
/*
Copyright 2021 Jeroen van Rijn <nom@duclavier.com>.
Made available under Odin's BSD-2 license.
A BigInt implementation in Odin.
For the theoretical underpinnings, see Knuth's The Art of Computer Programming, Volume 2, section 4.3.
The code started out as an idiomatic source port of libTomMath, which is in the public domain, with thanks.
*/
import "core:intrinsics"
is_initialized :: proc(a: ^Int) -> bool {
return a != rawptr(uintptr(0));
}
is_zero :: proc(a: ^Int) -> bool {
return is_initialized(a) && a.used == 0;
}
is_positive :: proc(a: ^Int) -> bool {
return is_initialized(a) && a.sign == .Zero_or_Positive;
}
is_pos :: is_positive;;
is_negative :: proc(a: ^Int) -> bool {
return is_initialized(a) && a.sign == .Negative;
}
is_neg :: is_negative;
/*
Compare two `Int`s, signed.
*/
compare :: proc(a, b: ^Int) -> Comparison_Flag {
if !is_initialized(a) { return .Uninitialized; }
if !is_initialized(b) { return .Uninitialized; }
/* Compare based on sign */
if a.sign != b.sign {
return .Less_Than if is_negative(a) else .Greater_Than;
}
x, y := a, b;
/* If negative, compare in the opposite direction */
if is_neg(a) {
x, y = b, a;
}
return cmp_mag(x, y);
}
cmp :: compare;
/*
Compare the magnitude of two `Int`s, unsigned.
*/
compare_magnitude :: proc(a, b: ^Int) -> Comparison_Flag {
if !is_initialized(a) { return .Uninitialized; }
if !is_initialized(b) { return .Uninitialized; }
/* Compare based on used digits */
if a.used != b.used {
return .Greater_Than if a.used > b.used else .Less_Than;
}
/* Same number of used digits, compare based on their value */
for n := a.used - 1; n >= 0; n -= 1 {
if a.digit[n] != b.digit[n] {
return .Greater_Than if a.digit[n] > b.digit[n] else .Less_Than;
}
}
return .Equal;
}
cmp_mag :: compare_magnitude;
/*
Compare an `Int` to an unsigned number upto the size of the backing type.
*/
compare_digit :: proc(a: ^Int, u: DIGIT) -> Comparison_Flag {
if !is_initialized(a) { return .Uninitialized; }
/* Compare based on sign */
if is_neg(a) {
return .Less_Than;
}
/* Compare based on magnitude */
if a.used > 1 {
return .Greater_Than;
}
/* Compare the only digit in `a` to `u`. */
if a.digit[0] != u {
return .Greater_Than if a.digit[0] > u else .Less_Than;
}
return .Equal;
}

View File

@@ -0,0 +1,50 @@
//+ignore
package bigint
/*
Copyright 2021 Jeroen van Rijn <nom@duclavier.com>.
Made available under Odin's BSD-2 license.
A BigInt implementation in Odin.
For the theoretical underpinnings, see Knuth's The Art of Computer Programming, Volume 2, section 4.3.
The code started out as an idiomatic source port of libTomMath, which is in the public domain, with thanks.
*/
import "core:fmt"
import "core:mem"
demo :: proc() {
a, b, c: ^Int;
err: Error;
a, err = init(-21);
defer destroy(a);
fmt.printf("a: %v, err: %v\n\n", a, err);
b, err = init(42);
defer destroy(b);
fmt.printf("b: %v, err: %v\n\n", b, err);
c, err = init();
defer destroy(c);
err = add(c, a, b);
fmt.printf("c: %v, err: %v\n\n", c, err);
}
main :: proc() {
ta := mem.Tracking_Allocator{};
mem.tracking_allocator_init(&ta, context.allocator);
context.allocator = mem.tracking_allocator(&ta);
fmt.printf("_DIGIT_BITS: %v\n_MIN_DIGIT_COUNT: %v\n_MAX_DIGIT_COUNT: %v\n_DEFAULT_DIGIT_COUNT: %v\n\n", _DIGIT_BITS, _MIN_DIGIT_COUNT, _MAX_DIGIT_COUNT, _DEFAULT_DIGIT_COUNT);
demo();
if len(ta.allocation_map) > 0 {
for _, v in ta.allocation_map {
fmt.printf("Leaked %v bytes @ %v\n", v.size, v.location);
}
}
}

View File

@@ -0,0 +1,232 @@
package bigint
/*
Copyright 2021 Jeroen van Rijn <nom@duclavier.com>.
Made available under Odin's BSD-2 license.
A BigInt implementation in Odin.
For the theoretical underpinnings, see Knuth's The Art of Computer Programming, Volume 2, section 4.3.
The code started out as an idiomatic source port of libTomMath, which is in the public domain, with thanks.
*/
import "core:mem"
import "core:intrinsics"
/*
Deallocates the backing memory of an Int.
*/
destroy :: proc(a: ^Int, allocator_zeroes := false, free_int := true, loc := #caller_location) {
if !is_initialized(a) {
// Nothing to do.
return;
}
if !allocator_zeroes {
mem.zero_slice(a.digit[:]);
}
free(&a.digit[0]);
a.used = 0;
a.allocated = 0;
if free_int {
free(a);
}
}
/*
Creates and returns a new `Int`.
*/
init_new :: proc(allocator_zeroes := true, allocator := context.allocator, size := _DEFAULT_DIGIT_COUNT) -> (a: ^Int, err: Error) {
/*
Allocating a new variable.
*/
a = new(Int, allocator);
a.digit = mem.make_dynamic_array_len_cap([dynamic]DIGIT, size, size, allocator);
a.allocated = 0;
a.used = 0;
a.sign = .Zero_or_Positive;
if len(a.digit) != size {
return a, .Out_of_Memory;
}
a.allocated = size;
if !allocator_zeroes {
_zero_unused(a);
}
return a, .OK;
}
/*
Initialize from a signed or unsigned integer.
Inits a new `Int` and then calls the appropriate `set` routine.
*/
init_new_integer :: proc(u: $T, minimize := false, allocator_zeroes := true, allocator := context.allocator) -> (a: ^Int, err: Error) where intrinsics.type_is_integer(T) {
n := _DEFAULT_DIGIT_COUNT;
if minimize {
n = _MIN_DIGIT_COUNT;
}
a, err = init_new(allocator_zeroes, allocator, n);
if err == .OK {
set(a, u, minimize);
}
return;
}
init :: proc{init_new, init_new_integer};
/*
Helpers to set an `Int` to a specific value.
*/
set_integer :: proc(a: ^Int, n: $T, minimize := false) where intrinsics.type_is_integer(T) {
n := n;
_panic_if_uninitialized(a);
a.used = 0;
a.sign = .Zero_or_Positive if n >= 0 else .Negative;
n = abs(n);
for n != 0 {
a.digit[a.used] = DIGIT(n) & _MASK;
a.used += 1;
n >>= _DIGIT_BITS;
}
if minimize {
shrink(a);
}
_zero_unused(a);
}
set :: proc{set_integer};
/*
Resize backing store.
*/
shrink :: proc(a: ^Int) -> (err: Error) {
needed := max(_MIN_DIGIT_COUNT, a.used);
if a.used != needed {
return grow(a, needed);
}
return .OK;
}
grow :: proc(a: ^Int, n: int) -> (err: Error) {
_panic_if_uninitialized(a);
resize(&a.digit, n);
if len(a.digit) != n {
return .Out_of_Memory;
}
a.used = min(n, a.used);
a.allocated = n;
return .OK;
}
/*
Clear `Int` and resize it to the default size.
*/
clear :: proc(a: ^Int) -> (err: Error) {
_panic_if_uninitialized(a);
mem.zero_slice(a.digit[:]);
a.sign = .Zero_or_Positive;
a.used = 0;
grow(a, _DEFAULT_DIGIT_COUNT);
return .OK;
}
/*
Set the `Int` to 0 and optionally shrink it to the minimum backing size.
*/
zero :: proc(a: ^Int, minimize := false) -> (err: Error) {
_panic_if_uninitialized(a);
mem.zero_slice(a.digit[:]);
a.sign = .Zero_or_Positive;
a.used = 0;
if minimize {
return shrink(a);
}
return .OK;
}
/*
Set the `Int` to 1 and optionally shrink it to the minimum backing size.
*/
one :: proc(a: ^Int, minimize := false) -> (err: Error) {
_panic_if_uninitialized(a);
mem.zero_slice(a.digit[:]);
a.sign = .Zero_or_Positive;
a.used = 1;
a.digit[0] = 1;
if minimize {
return shrink(a);
}
return .OK;
}
/*
Set the `Int` to -1 and optionally shrink it to the minimum backing size.
*/
minus_one :: proc(a: ^Int, minimize := false) -> (err: Error) {
_panic_if_uninitialized(a);
mem.zero_slice(a.digit[:]);
a.sign = .Negative;
a.used = 1;
a.digit[0] = 1;
if minimize {
return shrink(a);
}
return .OK;
}
/*
Internal helpers.
*/
_panic_if_uninitialized :: proc(a: ^Int, loc := #caller_location) {
if !is_initialized(a) {
panic("Int was not properly initialized.", loc);
}
}
_zero_unused :: proc(a: ^Int) {
_panic_if_uninitialized(a);
if a.used < a.allocated {
mem.zero_slice(a.digit[a.used:]);
}
}
clamp :: proc(a: ^Int) {
_panic_if_uninitialized(a);
/*
Trim unused digits
This is used to ensure that leading zero digits are
trimmed and the leading "used" digit will be non-zero.
Typically very fast. Also fixes the sign if there
are no more leading digits.
*/
for a.used > 0 && a.digit[a.used - 1] == 0 {
a.used -= 1;
}
if is_zero(a) {
a.sign = .Zero_or_Positive;
}
}