From 18dda6ff9d8f67f259bafbcdfe5fd5286f5a45fa Mon Sep 17 00:00:00 2001 From: Jeroen van Rijn Date: Thu, 15 Jul 2021 23:34:02 +0200 Subject: [PATCH] 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` --- core/math/bigint/basic.odin | 196 ++++++++++++++++++++++++++++ core/math/bigint/bigint.odin | 115 ++++++++++++++++ core/math/bigint/build.bat | 2 + core/math/bigint/compares.odin | 98 ++++++++++++++ core/math/bigint/example.odin | 50 +++++++ core/math/bigint/helpers.odin | 232 +++++++++++++++++++++++++++++++++ 6 files changed, 693 insertions(+) create mode 100644 core/math/bigint/basic.odin create mode 100644 core/math/bigint/bigint.odin create mode 100644 core/math/bigint/build.bat create mode 100644 core/math/bigint/compares.odin create mode 100644 core/math/bigint/example.odin create mode 100644 core/math/bigint/helpers.odin diff --git a/core/math/bigint/basic.odin b/core/math/bigint/basic.odin new file mode 100644 index 000000000..56ac692ff --- /dev/null +++ b/core/math/bigint/basic.odin @@ -0,0 +1,196 @@ +package bigint + +/* + Copyright 2021 Jeroen van Rijn . + 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; +} \ No newline at end of file diff --git a/core/math/bigint/bigint.odin b/core/math/bigint/bigint.odin new file mode 100644 index 000000000..1445bc69d --- /dev/null +++ b/core/math/bigint/bigint.odin @@ -0,0 +1,115 @@ +package bigint + +/* + Copyright 2021 Jeroen van Rijn . + 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, +}; \ No newline at end of file diff --git a/core/math/bigint/build.bat b/core/math/bigint/build.bat new file mode 100644 index 000000000..df9cd1e85 --- /dev/null +++ b/core/math/bigint/build.bat @@ -0,0 +1,2 @@ +@echo off +odin run . -vet \ No newline at end of file diff --git a/core/math/bigint/compares.odin b/core/math/bigint/compares.odin new file mode 100644 index 000000000..1838d901e --- /dev/null +++ b/core/math/bigint/compares.odin @@ -0,0 +1,98 @@ +package bigint + +/* + Copyright 2021 Jeroen van Rijn . + 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; +} \ No newline at end of file diff --git a/core/math/bigint/example.odin b/core/math/bigint/example.odin new file mode 100644 index 000000000..4b89b193d --- /dev/null +++ b/core/math/bigint/example.odin @@ -0,0 +1,50 @@ +//+ignore +package bigint + +/* + Copyright 2021 Jeroen van Rijn . + 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); + } + } +} \ No newline at end of file diff --git a/core/math/bigint/helpers.odin b/core/math/bigint/helpers.odin new file mode 100644 index 000000000..f96942954 --- /dev/null +++ b/core/math/bigint/helpers.odin @@ -0,0 +1,232 @@ +package bigint + +/* + Copyright 2021 Jeroen van Rijn . + 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; + } +} \ No newline at end of file