mirrors/Odin
1
0
Fork 0
mirror of https://github.com/odin-lang/Odin.git synced 2026-02-26 12:54:57 +00:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #1113 from Kelimion/bigint

Browse Source 
big: Add `expt_mod`, new comparison helpers, etc.
This commit is contained in:
Jeroen van Rijn
2021-09-01 19:25:52 +02:00
committed by GitHub
parent b2cf0755f2 ae354731ed
commit e2f035d6ee
14 changed files with 1729 additions and 365 deletions
Download Patch File Download Diff File Expand all files Collapse all files

7
core/math/big/api.odin
View File

@@ -1,14 +1,15 @@
package math_big
/*
Copyright 2021 Jeroen van Rijn <nom@duclavier.com>.
Made available under Odin's BSD-2 license.
Made available under Odin's BSD-3 license.
An arbitrary precision mathematics 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.
This file collects public proc maps and their aliases.
*/
package math_big
/*
=== === === === === === === === === === === === === === === === === === === === === === === ===
Basic arithmetic.

1
core/math/big/build.bat
View File

@@ -2,6 +2,7 @@
:odin run . -vet
set TEST_ARGS=-fast-tests
:set TEST_ARGS=
:odin build . -build-mode:shared -show-timings -o:minimal -no-bounds-check -define:MATH_BIG_EXE=false && python test.py %TEST_ARGS%
odin build . -build-mode:shared -show-timings -o:size -no-bounds-check -define:MATH_BIG_EXE=false && python test.py %TEST_ARGS%
:odin build . -build-mode:shared -show-timings -o:size -define:MATH_BIG_EXE=false && python test.py %TEST_ARGS%

28
core/math/big/common.odin
View File

@@ -1,5 +1,3 @@
package math_big
/*
Copyright 2021 Jeroen van Rijn <nom@duclavier.com>.
Made available under Odin's BSD-3 license.
@@ -8,6 +6,7 @@ package math_big
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.
*/
package math_big
import "core:intrinsics"
@@ -57,10 +56,10 @@ when #config(MATH_BIG_EXE, true) {
debugged where necessary.
*/
_DEFAULT_MUL_KARATSUBA_CUTOFF :: #config(MUL_KARATSUBA_CUTOFF, 80);
_DEFAULT_SQR_KARATSUBA_CUTOFF :: #config(SQR_KARATSUBA_CUTOFF, 120);
_DEFAULT_MUL_TOOM_CUTOFF :: #config(MUL_TOOM_CUTOFF, 350);
_DEFAULT_SQR_TOOM_CUTOFF :: #config(SQR_TOOM_CUTOFF, 400);
_DEFAULT_MUL_KARATSUBA_CUTOFF :: #config(MATH_BIG_MUL_KARATSUBA_CUTOFF, 80);
_DEFAULT_SQR_KARATSUBA_CUTOFF :: #config(MATH_BIG_SQR_KARATSUBA_CUTOFF, 120);
_DEFAULT_MUL_TOOM_CUTOFF :: #config(MATH_BIG_MUL_TOOM_CUTOFF, 350);
_DEFAULT_SQR_TOOM_CUTOFF :: #config(MATH_BIG_SQR_TOOM_CUTOFF, 400);
MAX_ITERATIONS_ROOT_N := 500;
@@ -85,15 +84,22 @@ FACTORIAL_BINARY_SPLIT_MAX_RECURSIONS := 100;
2) Optimizations thanks to precomputed masks wouldn't work.
*/
MATH_BIG_FORCE_64_BIT :: #config(MATH_BIG_FORCE_64_BIT, false);
MATH_BIG_FORCE_32_BIT :: #config(MATH_BIG_FORCE_32_BIT, false);
MATH_BIG_FORCE_64_BIT :: #config(MATH_BIG_FORCE_64_BIT, false);
MATH_BIG_FORCE_32_BIT :: #config(MATH_BIG_FORCE_32_BIT, false);
when (MATH_BIG_FORCE_32_BIT && MATH_BIG_FORCE_64_BIT) { #panic("Cannot force 32-bit and 64-bit big backend simultaneously."); };
_LOW_MEMORY :: #config(BIGINT_SMALL_MEMORY, false);
/*
Trade a smaller memory footprint for more processing overhead?
*/
_LOW_MEMORY :: #config(MATH_BIG_SMALL_MEMORY, false);
when _LOW_MEMORY {
_DEFAULT_DIGIT_COUNT :: 8;
_DEFAULT_DIGIT_COUNT :: 8;
_TAB_SIZE :: 32;
_MAX_WIN_SIZE :: 5;
} else {
_DEFAULT_DIGIT_COUNT :: 32;
_DEFAULT_DIGIT_COUNT :: 32;
_TAB_SIZE :: 256;
_MAX_WIN_SIZE :: 0;
}
/*

143
core/math/big/example.odin
View File

@@ -1,6 +1,4 @@
//+ignore
package math_big
/*
Copyright 2021 Jeroen van Rijn <nom@duclavier.com>.
Made available under Odin's BSD-3 license.
@@ -9,6 +7,8 @@ package math_big
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.
*/
package math_big
import "core:fmt"
import "core:mem"
@@ -18,11 +18,14 @@ print_configation :: proc() {
`
Configuration:
_DIGIT_BITS %v
_SMALL_MEMORY %v
_MIN_DIGIT_COUNT %v
_MAX_DIGIT_COUNT %v
_DEFAULT_DIGIT_COUNT %v
_MAX_COMBA %v
_WARRAY %v
_TAB_SIZE %v
_MAX_WIN_SIZE %v
Runtime tunable:
MUL_KARATSUBA_CUTOFF %v
SQR_KARATSUBA_CUTOFF %v
@@ -34,11 +37,14 @@ Runtime tunable:
FACTORIAL_BINARY_SPLIT_MAX_RECURSIONS %v
`, _DIGIT_BITS,
_LOW_MEMORY,
_MIN_DIGIT_COUNT,
_MAX_DIGIT_COUNT,
_DEFAULT_DIGIT_COUNT,
_MAX_COMBA,
_WARRAY,
_TAB_SIZE,
_MAX_WIN_SIZE,
MUL_KARATSUBA_CUTOFF,
SQR_KARATSUBA_CUTOFF,
MUL_TOOM_CUTOFF,
@@ -73,138 +79,11 @@ print :: proc(name: string, a: ^Int, base := i8(10), print_name := true, newline
}
}
int_to_byte :: proc(v: ^Int) {
err: Error;
size: int;
print("v: ", v);
fmt.println();
t := &Int{};
defer destroy(t);
if size, err = int_to_bytes_size(v); err != nil {
fmt.printf("int_to_bytes_size returned: %v\n", err);
return;
}
b1 := make([]u8, size, context.temp_allocator);
err = int_to_bytes_big(v, b1);
int_from_bytes_big(t, b1);
fmt.printf("big: %v | err: %v\n", b1, err);
int_from_bytes_big(t, b1);
if internal_cmp_mag(t, v) != 0 {
print("\tError parsing t: ", t);
}
if size, err = int_to_bytes_size(v); err != nil {
fmt.printf("int_to_bytes_size returned: %v\n", err);
return;
}
b2 := make([]u8, size, context.temp_allocator);
err = int_to_bytes_big_python(v, b2);
fmt.printf("big python: %v | err: %v\n", b2, err);
if err == nil {
int_from_bytes_big_python(t, b2);
if internal_cmp_mag(t, v) != 0 {
print("\tError parsing t: ", t);
}
}
if size, err = int_to_bytes_size(v, true); err != nil {
fmt.printf("int_to_bytes_size returned: %v\n", err);
return;
}
b3 := make([]u8, size, context.temp_allocator);
err = int_to_bytes_big(v, b3, true);
fmt.printf("big signed: %v | err: %v\n", b3, err);
int_from_bytes_big(t, b3, true);
if internal_cmp(t, v) != 0 {
print("\tError parsing t: ", t);
}
if size, err = int_to_bytes_size(v, true); err != nil {
fmt.printf("int_to_bytes_size returned: %v\n", err);
return;
}
b4 := make([]u8, size, context.temp_allocator);
err = int_to_bytes_big_python(v, b4, true);
fmt.printf("big signed python: %v | err: %v\n", b4, err);
int_from_bytes_big_python(t, b4, true);
if internal_cmp(t, v) != 0 {
print("\tError parsing t: ", t);
}
}
int_to_byte_little :: proc(v: ^Int) {
err: Error;
size: int;
print("v: ", v);
fmt.println();
t := &Int{};
defer destroy(t);
if size, err = int_to_bytes_size(v); err != nil {
fmt.printf("int_to_bytes_size returned: %v\n", err);
return;
}
b1 := make([]u8, size, context.temp_allocator);
err = int_to_bytes_little(v, b1);
fmt.printf("little: %v | err: %v\n", b1, err);
int_from_bytes_little(t, b1);
if internal_cmp_mag(t, v) != 0 {
print("\tError parsing t: ", t);
}
if size, err = int_to_bytes_size(v); err != nil {
fmt.printf("int_to_bytes_size returned: %v\n", err);
return;
}
b2 := make([]u8, size, context.temp_allocator);
err = int_to_bytes_little_python(v, b2);
fmt.printf("little python: %v | err: %v\n", b2, err);
if err == nil {
int_from_bytes_little_python(t, b2);
if internal_cmp_mag(t, v) != 0 {
print("\tError parsing t: ", t);
}
}
if size, err = int_to_bytes_size(v, true); err != nil {
fmt.printf("int_to_bytes_size returned: %v\n", err);
return;
}
b3 := make([]u8, size, context.temp_allocator);
err = int_to_bytes_little(v, b3, true);
fmt.printf("little signed: %v | err: %v\n", b3, err);
int_from_bytes_little(t, b3, true);
if internal_cmp(t, v) != 0 {
print("\tError parsing t: ", t);
}
if size, err = int_to_bytes_size(v, true); err != nil {
fmt.printf("int_to_bytes_size returned: %v\n", err);
return;
}
b4 := make([]u8, size, context.temp_allocator);
err = int_to_bytes_little_python(v, b4, true);
fmt.printf("little signed python: %v | err: %v\n", b4, err);
int_from_bytes_little_python(t, b4, true);
if internal_cmp(t, v) != 0 {
print("\tError parsing t: ", t);
}
}
// printf :: fmt.printf;
demo :: proc() {
a, b, c, d, e, f := &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{};
defer destroy(a, b, c, d, e, f);
a, b, c, d, e, f, res := &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{};
defer destroy(a, b, c, d, e, f, res);
}
main :: proc() {

3
core/math/big/helpers.odin
View File

@@ -1,5 +1,3 @@
package math_big
/*
Copyright 2021 Jeroen van Rijn <nom@duclavier.com>.
Made available under Odin's BSD-3 license.
@@ -8,6 +6,7 @@ package math_big
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.
*/
package math_big
import "core:intrinsics"
import rnd "core:math/rand"

223
core/math/big/internal.odin
View File

@@ -1,6 +1,4 @@
//+ignore
package math_big
/*
Copyright 2021 Jeroen van Rijn <nom@duclavier.com>.
Made available under Odin's BSD-3 license.
@@ -31,6 +29,7 @@ package math_big
TODO: Handle +/- Infinity and NaN.
*/
package math_big
import "core:mem"
import "core:intrinsics"
@@ -137,7 +136,7 @@ internal_int_add_signed :: proc(dest, a, b: ^Int, allocator := context.allocator
Subtract the one with the greater magnitude from the other.
The result gets the sign of the one with the greater magnitude.
*/
if #force_inline internal_cmp_mag(a, b) == -1 {
if #force_inline internal_lt_abs(a, b) {
x, y = y, x;
}
@@ -359,7 +358,7 @@ internal_int_sub_signed :: proc(dest, number, decrease: ^Int, allocator := conte
Subtract a positive from a positive, OR negative from a negative.
First, take the difference between their magnitudes, then...
*/
if #force_inline internal_cmp_mag(number, decrease) == -1 {
if #force_inline internal_lt_abs(number, decrease) {
/*
The second has a larger magnitude.
The result has the *opposite* sign from the first number.
@@ -719,7 +718,7 @@ internal_int_divmod :: proc(quotient, remainder, numerator, denominator: ^Int, a
/*
If numerator < denominator then quotient = 0, remainder = numerator.
*/
if #force_inline internal_cmp_mag(numerator, denominator) == -1 {
if #force_inline internal_lt_abs(numerator, denominator) {
if remainder != nil {
internal_copy(remainder, numerator) or_return;
}
@@ -732,7 +731,6 @@ internal_int_divmod :: proc(quotient, remainder, numerator, denominator: ^Int, a
if (denominator.used > 2 * MUL_KARATSUBA_CUTOFF) && (denominator.used <= (numerator.used / 3) * 2) {
assert(denominator.used >= 160 && numerator.used >= 240, "MUL_KARATSUBA_CUTOFF global not properly set.");
err = _private_int_div_recursive(quotient, remainder, numerator, denominator);
// err = #force_inline _private_int_div_school(quotient, remainder, numerator, denominator);
} else {
when true {
err = #force_inline _private_int_div_school(quotient, remainder, numerator, denominator);
@@ -992,13 +990,21 @@ internal_int_mod_bits :: proc(remainder, numerator: ^Int, bits: int, allocator :
public ones that have already satisfied these constraints.
*/
/*
This procedure returns the allocated capacity of an Int.
Assumes `a` not to be `nil`.
*/
internal_int_allocated_cap :: #force_inline proc(a: ^Int) -> (cap: int) {
raw := transmute(mem.Raw_Dynamic_Array)a.digit;
return raw.cap;
}
/*
This procedure will return `true` if the `Int` is initialized, `false` if not.
Assumes `a` not to be `nil`.
*/
internal_int_is_initialized :: #force_inline proc(a: ^Int) -> (initialized: bool) {
raw := transmute(mem.Raw_Dynamic_Array)a.digit;
return raw.cap >= _MIN_DIGIT_COUNT;
return internal_int_allocated_cap(a) >= _MIN_DIGIT_COUNT;
}
internal_is_initialized :: proc { internal_int_is_initialized, };
@@ -1091,6 +1097,7 @@ internal_is_power_of_two :: proc { internal_int_is_power_of_two, };
Expects `a` and `b` both to be valid `Int`s, i.e. initialized and not `nil`.
*/
internal_int_compare :: #force_inline proc(a, b: ^Int) -> (comparison: int) {
assert_if_nil(a, b);
a_is_negative := #force_inline internal_is_negative(a);
/*
@@ -1114,6 +1121,7 @@ internal_cmp :: internal_compare;
Expects: `a` and `b` both to be valid `Int`s, i.e. initialized and not `nil`.
*/
internal_int_compare_digit :: #force_inline proc(a: ^Int, b: DIGIT) -> (comparison: int) {
assert_if_nil(a);
a_is_negative := #force_inline internal_is_negative(a);
switch {
@@ -1145,6 +1153,7 @@ internal_cmp_digit :: internal_compare_digit;
Compare the magnitude of two `Int`s, unsigned.
*/
internal_int_compare_magnitude :: #force_inline proc(a, b: ^Int) -> (comparison: int) {
assert_if_nil(a, b);
/*
Compare based on used digits.
*/
@@ -1172,6 +1181,177 @@ internal_int_compare_magnitude :: #force_inline proc(a, b: ^Int) -> (comparison:
internal_compare_magnitude :: proc { internal_int_compare_magnitude, };
internal_cmp_mag :: internal_compare_magnitude;
/*
bool := a < b
*/
internal_int_less_than :: #force_inline proc(a, b: ^Int) -> (less_than: bool) {
return internal_cmp(a, b) == -1;
}
/*
bool := a < b
*/
internal_int_less_than_digit :: #force_inline proc(a: ^Int, b: DIGIT) -> (less_than: bool) {
return internal_cmp_digit(a, b) == -1;
}
/*
bool := |a| < |b|
Compares the magnitudes only, ignores the sign.
*/
internal_int_less_than_abs :: #force_inline proc(a, b: ^Int) -> (less_than: bool) {
return internal_cmp_mag(a, b) == -1;
}
internal_less_than :: proc {
internal_int_less_than,
internal_int_less_than_digit,
};
internal_lt :: internal_less_than;
internal_less_than_abs :: proc {
internal_int_less_than_abs,
};
internal_lt_abs :: internal_less_than_abs;
/*
bool := a <= b
*/
internal_int_less_than_or_equal :: #force_inline proc(a, b: ^Int) -> (less_than_or_equal: bool) {
return internal_cmp(a, b) <= 0;
}
/*
bool := a <= b
*/
internal_int_less_than_or_equal_digit :: #force_inline proc(a: ^Int, b: DIGIT) -> (less_than_or_equal: bool) {
return internal_cmp_digit(a, b) <= 0;
}
/*
bool := |a| <= |b|
Compares the magnitudes only, ignores the sign.
*/
internal_int_less_than_or_equal_abs :: #force_inline proc(a, b: ^Int) -> (less_than_or_equal: bool) {
return internal_cmp_mag(a, b) <= 0;
}
internal_less_than_or_equal :: proc {
internal_int_less_than_or_equal,
internal_int_less_than_or_equal_digit,
};
internal_lte :: internal_less_than_or_equal;
internal_less_than_or_equal_abs :: proc {
internal_int_less_than_or_equal_abs,
};
internal_lte_abs :: internal_less_than_or_equal_abs;
/*
bool := a == b
*/
internal_int_equals :: #force_inline proc(a, b: ^Int) -> (equals: bool) {
return internal_cmp(a, b) == 0;
}
/*
bool := a == b
*/
internal_int_equals_digit :: #force_inline proc(a: ^Int, b: DIGIT) -> (equals: bool) {
return internal_cmp_digit(a, b) == 0;
}
/*
bool := |a| == |b|
Compares the magnitudes only, ignores the sign.
*/
internal_int_equals_abs :: #force_inline proc(a, b: ^Int) -> (equals: bool) {
return internal_cmp_mag(a, b) == 0;
}
internal_equals :: proc {
internal_int_equals,
internal_int_equals_digit,
};
internal_eq :: internal_equals;
internal_equals_abs :: proc {
internal_int_equals_abs,
};
internal_eq_abs :: internal_equals_abs;
/*
bool := a >= b
*/
internal_int_greater_than_or_equal :: #force_inline proc(a, b: ^Int) -> (greater_than_or_equal: bool) {
return internal_cmp(a, b) >= 0;
}
/*
bool := a >= b
*/
internal_int_greater_than_or_equal_digit :: #force_inline proc(a: ^Int, b: DIGIT) -> (greater_than_or_equal: bool) {
return internal_cmp_digit(a, b) >= 0;
}
/*
bool := |a| >= |b|
Compares the magnitudes only, ignores the sign.
*/
internal_int_greater_than_or_equal_abs :: #force_inline proc(a, b: ^Int) -> (greater_than_or_equal: bool) {
return internal_cmp_mag(a, b) >= 0;
}
internal_greater_than_or_equal :: proc {
internal_int_greater_than_or_equal,
internal_int_greater_than_or_equal_digit,
};
internal_gte :: internal_greater_than_or_equal;
internal_greater_than_or_equal_abs :: proc {
internal_int_greater_than_or_equal_abs,
};
internal_gte_abs :: internal_greater_than_or_equal_abs;
/*
bool := a > b
*/
internal_int_greater_than :: #force_inline proc(a, b: ^Int) -> (greater_than: bool) {
return internal_cmp(a, b) == 1;
}
/*
bool := a > b
*/
internal_int_greater_than_digit :: #force_inline proc(a: ^Int, b: DIGIT) -> (greater_than: bool) {
return internal_cmp_digit(a, b) == 1;
}
/*
bool := |a| > |b|
Compares the magnitudes only, ignores the sign.
*/
internal_int_greater_than_abs :: #force_inline proc(a, b: ^Int) -> (greater_than: bool) {
return internal_cmp_mag(a, b) == 1;
}
internal_greater_than :: proc {
internal_int_greater_than,
internal_int_greater_than_digit,
};
internal_gt :: internal_greater_than;
internal_greater_than_abs :: proc {
internal_int_greater_than_abs,
};
internal_gt_abs :: internal_greater_than_abs;
/*
Check if remainders are possible squares - fast exclude non-squares.
@@ -1229,7 +1409,7 @@ internal_int_is_square :: proc(a: ^Int, allocator := context.allocator) -> (squa
sqrt(t, a) or_return;
sqr(t, t) or_return;
square = internal_cmp_mag(t, a) == 0;
square = internal_eq_abs(t, a);
return;
}
@@ -1461,7 +1641,7 @@ internal_int_sqrt :: proc(dest, src: ^Int, allocator := context.allocator) -> (e
internal_add(t2, t1, x) or_return;
internal_shr(y, t2, 1) or_return;
if c := internal_cmp(y, x); c == 0 || c == 1 {
if internal_gte(y, x) {
internal_swap(dest, x);
return nil;
}
@@ -1576,8 +1756,8 @@ internal_int_root_n :: proc(dest, src: ^Int, n: int, allocator := context.alloca
Number of rounds is at most log_2(root). If it is more it
got stuck, so break out of the loop and do the rest manually.
*/
if ilog2 -= 1; ilog2 == 0 { break; }
if internal_cmp(t1, t2) == 0 { break; }
if ilog2 -= 1; ilog2 == 0 { break; }
if internal_eq(t1, t2) { break; }
iterations += 1;
if iterations == MAX_ITERATIONS_ROOT_N {
@@ -1615,7 +1795,7 @@ internal_int_root_n :: proc(dest, src: ^Int, n: int, allocator := context.alloca
for {
internal_pow(t2, t1, n) or_return;
if internal_cmp(t2, a) != 1 { break; }
if internal_lt(t2, a) { break; }
internal_sub(t1, t1, DIGIT(1)) or_return;
@@ -1651,8 +1831,7 @@ internal_int_destroy :: proc(integers: ..^Int) {
integers := integers;
for a in &integers {
raw := transmute(mem.Raw_Dynamic_Array)a.digit;
if raw.cap > 0 {
if internal_int_allocated_cap(a) > 0 {
mem.zero_slice(a.digit[:]);
free(&a.digit[0]);
}
@@ -1821,12 +2000,12 @@ internal_int_inverse_modulo :: proc(dest, a, b: ^Int, allocator := context.alloc
/*
For all n in N and n > 0, n = 0 mod 1.
*/
if internal_is_positive(a) && internal_cmp(b, 1) == 0 { return internal_zero(dest); }
if internal_is_positive(a) && internal_eq(b, 1) { return internal_zero(dest); }
/*
`b` cannot be negative and has to be > 1
*/
if internal_is_negative(b) && internal_cmp(b, 1) != 1 { return .Invalid_Argument; }
if internal_is_negative(b) || internal_gt(b, 1) { return .Invalid_Argument; }
/*
If the modulus is odd we can use a faster routine instead.
@@ -1914,23 +2093,23 @@ internal_int_shrink :: proc(a: ^Int) -> (err: Error) {
internal_shrink :: proc { internal_int_shrink, };
internal_int_grow :: proc(a: ^Int, digits: int, allow_shrink := false, allocator := context.allocator) -> (err: Error) {
raw := transmute(mem.Raw_Dynamic_Array)a.digit;
/*
We need at least _MIN_DIGIT_COUNT or a.used digits, whichever is bigger.
The caller is asking for `digits`. Let's be accomodating.
*/
cap := internal_int_allocated_cap(a);
needed := max(_MIN_DIGIT_COUNT, a.used, digits);
if !allow_shrink {
needed = max(needed, raw.cap);
needed = max(needed, cap);
}
/*
If not yet iniialized, initialize the `digit` backing with the allocator we were passed.
*/
if raw.cap == 0 {
if cap == 0 {
a.digit = make([dynamic]DIGIT, needed, allocator);
} else if raw.cap != needed {
} else if cap != needed {
/*
`[dynamic]DIGIT` already knows what allocator was used for it, so resize will do the right thing.
*/

3
core/math/big/logical.odin
View File

@@ -1,5 +1,3 @@
package math_big
/*
Copyright 2021 Jeroen van Rijn <nom@duclavier.com>.
Made available under Odin's BSD-3 license.
@@ -10,6 +8,7 @@ package math_big
This file contains logical operations like `and`, `or` and `xor`.
*/
package math_big
/*
The `and`, `or` and `xor` binops differ in two lines only.

311
core/math/big/prime.odin
View File

@@ -1,5 +1,3 @@
package math_big
/*
Copyright 2021 Jeroen van Rijn <nom@duclavier.com>.
Made available under Odin's BSD-3 license.
@@ -10,12 +8,13 @@ package math_big
This file contains prime finding operations.
*/
package math_big
/*
Determines if an Integer is divisible by one of the _PRIME_TABLE primes.
Returns true if it is, false if not.
*/
int_prime_is_divisible :: proc(a: ^Int, allocator := context.allocator) -> (res: bool, err: Error) {
internal_int_prime_is_divisible :: proc(a: ^Int, allocator := context.allocator) -> (res: bool, err: Error) {
assert_if_nil(a);
context.allocator = allocator;
@@ -34,177 +33,175 @@ int_prime_is_divisible :: proc(a: ^Int, allocator := context.allocator) -> (res:
}
/*
Computes xR**-1 == x (mod N) via Montgomery Reduction.
This is a shell function that calls either the normal or Montgomery exptmod functions.
Originally the call to the Montgomery code was embedded in the normal function but that
wasted alot of stack space for nothing (since 99% of the time the Montgomery code would be called).
Computes res == G**X mod P.
Assumes `res`, `G`, `X` and `P` to not be `nil` and for `G`, `X` and `P` to have been initialized.
*/
internal_int_montgomery_reduce :: proc(x, n: ^Int, rho: DIGIT, allocator := context.allocator) -> (err: Error) {
context.allocator = allocator;
/*
Can the fast reduction [comba] method be used?
Note that unlike in mul, you're safely allowed *less* than the available columns [255 per default],
since carries are fixed up in the inner loop.
*/
digs := (n.used * 2) + 1;
if digs < _WARRAY && x.used <= _WARRAY && n.used < _MAX_COMBA {
return _private_montgomery_reduce_comba(x, n, rho);
}
/*
Grow the input as required
*/
internal_grow(x, digs) or_return;
x.used = digs;
for ix := 0; ix < n.used; ix += 1 {
/*
`mu = ai * rho mod b`
The value of rho must be precalculated via `int_montgomery_setup()`,
such that it equals -1/n0 mod b this allows the following inner loop
to reduce the input one digit at a time.
*/
mu := DIGIT((_WORD(x.digit[ix]) * _WORD(rho)) & _WORD(_MASK));
/*
a = a + mu * m * b**i
Multiply and add in place.
*/
u := DIGIT(0);
iy := int(0);
for ; iy < n.used; iy += 1 {
/*
Compute product and sum.
*/
r := (_WORD(mu) * _WORD(n.digit[iy]) + _WORD(u) + _WORD(x.digit[ix + iy]));
/*
Get carry.
*/
u = DIGIT(r >> _DIGIT_BITS);
/*
Fix digit.
*/
x.digit[ix + iy] = DIGIT(r & _WORD(_MASK));
}
/*
At this point the ix'th digit of x should be zero.
Propagate carries upwards as required.
*/
for u != 0 {
x.digit[ix + iy] += u;
u = x.digit[ix + iy] >> _DIGIT_BITS;
x.digit[ix + iy] &= _MASK;
iy += 1;
}
}
/*
At this point the n.used'th least significant digits of x are all zero,
which means we can shift x to the right by n.used digits and the
residue is unchanged.
x = x/b**n.used.
*/
internal_clamp(x);
internal_shr_digit(x, n.used);
/*
if x >= n then x = x - n
*/
if internal_cmp_mag(x, n) != -1 {
return internal_sub(x, x, n);
}
return nil;
}
int_montgomery_reduce :: proc(x, n: ^Int, rho: DIGIT, allocator := context.allocator) -> (err: Error) {
assert_if_nil(x, n);
internal_int_exponent_mod :: proc(res, G, X, P: ^Int, allocator := context.allocator) -> (err: Error) {
context.allocator = allocator;
internal_clear_if_uninitialized(x, n) or_return;
dr: int;
return #force_inline internal_int_montgomery_reduce(x, n, rho);
/*
Modulus P must be positive.
*/
if internal_is_negative(P) { return .Invalid_Argument; }
/*
If exponent X is negative we have to recurse.
*/
if internal_is_negative(X) {
tmpG, tmpX := &Int{}, &Int{};
defer internal_destroy(tmpG, tmpX);
internal_init_multi(tmpG, tmpX) or_return;
/*
First compute 1/G mod P.
*/
internal_invmod(tmpG, G, P) or_return;
/*
now get |X|.
*/
internal_abs(tmpX, X) or_return;
/*
And now compute (1/G)**|X| instead of G**X [X < 0].
*/
return internal_int_exponent_mod(res, tmpG, tmpX, P);
}
/*
Modified diminished radix reduction.
*/
can_reduce_2k_l := _private_int_reduce_is_2k_l(P) or_return;
if can_reduce_2k_l {
return _private_int_exponent_mod(res, G, X, P, 1);
}
/*
Is it a DR modulus? default to no.
*/
dr = 1 if _private_dr_is_modulus(P) else 0;
/*
If not, is it a unrestricted DR modulus?
*/
if dr == 0 {
reduce_is_2k := _private_int_reduce_is_2k(P) or_return;
dr = 2 if reduce_is_2k else 0;
}
/*
If the modulus is odd or dr != 0 use the montgomery method.
*/
if internal_int_is_odd(P) || dr != 0 {
return _private_int_exponent_mod(res, G, X, P, dr);
}
/*
Otherwise use the generic Barrett reduction technique.
*/
return _private_int_exponent_mod(res, G, X, P, 0);
}
/*
Shifts with subtractions when the result is greater than b.
Kronecker symbol (a|p)
Straightforward implementation of algorithm 1.4.10 in
Henri Cohen: "A Course in Computational Algebraic Number Theory"
The method is slightly modified to shift B unconditionally upto just under
the leading bit of b. This saves alot of multiple precision shifting.
*/
internal_int_montgomery_calc_normalization :: proc(a, b: ^Int, allocator := context.allocator) -> (err: Error) {
context.allocator = allocator;
/*
How many bits of last digit does b use.
*/
bits := internal_count_bits(b) % _DIGIT_BITS;
if b.used > 1 {
power := ((b.used - 1) * _DIGIT_BITS) + bits - 1;
internal_int_power_of_two(a, power) or_return;
} else {
internal_one(a);
bits = 1;
@book{cohen2013course,
title={A course in computational algebraic number theory},
author={Cohen, Henri},
volume={138},
year={2013},
publisher={Springer Science \& Business Media}
}
/*
Now compute C = A * B mod b.
*/
for x := bits - 1; x < _DIGIT_BITS; x += 1 {
internal_int_shl1(a, a) or_return;
if internal_cmp_mag(a, b) != -1 {
internal_sub(a, a, b) or_return;
Assumes `a` and `p` to not be `nil` and to have been initialized.
*/
internal_int_kronecker :: proc(a, p: ^Int, allocator := context.allocator) -> (kronecker: int, err: Error) {
context.allocator = allocator;
a1, p1, r := &Int{}, &Int{}, &Int{};
defer internal_destroy(a1, p1, r);
table := []int{0, 1, 0, -1, 0, -1, 0, 1};
if internal_int_is_zero(p) {
if a.used == 1 && a.digit[0] == 1 {
return 1, nil;
} else {
return 0, nil;
}
}
return nil;
}
int_montgomery_calc_normalization :: proc(a, b: ^Int, allocator := context.allocator) -> (err: Error) {
assert_if_nil(a, b);
context.allocator = allocator;
internal_clear_if_uninitialized(a, b) or_return;
return #force_inline internal_int_montgomery_calc_normalization(a, b);
}
/*
Sets up the Montgomery reduction stuff.
*/
internal_int_montgomery_setup :: proc(n: ^Int) -> (rho: DIGIT, err: Error) {
/*
Fast inversion mod 2**k
Based on the fact that:
XA = 1 (mod 2**n) => (X(2-XA)) A = 1 (mod 2**2n)
=> 2*X*A - X*X*A*A = 1
=> 2*(1) - (1) = 1
*/
b := n.digit[0];
if b & 1 == 0 { return 0, .Invalid_Argument; }
x := (((b + 2) & 4) << 1) + b; /* here x*a==1 mod 2**4 */
x *= 2 - (b * x); /* here x*a==1 mod 2**8 */
x *= 2 - (b * x); /* here x*a==1 mod 2**16 */
when _WORD_TYPE_BITS == 64 {
x *= 2 - (b * x); /* here x*a==1 mod 2**32 */
x *= 2 - (b * x); /* here x*a==1 mod 2**64 */
if internal_is_even(a) && internal_is_even(p) {
return 0, nil;
}
/*
rho = -1/m mod b
*/
rho = DIGIT(((_WORD(1) << _WORD(_DIGIT_BITS)) - _WORD(x)) & _WORD(_MASK));
return rho, nil;
}
internal_copy(a1, a) or_return;
internal_copy(p1, p) or_return;
int_montgomery_setup :: proc(n: ^Int, allocator := context.allocator) -> (rho: DIGIT, err: Error) {
assert_if_nil(n);
internal_clear_if_uninitialized(n, allocator) or_return;
v := internal_count_lsb(p1) or_return;
internal_shr(p1, p1, v) or_return;
return #force_inline internal_int_montgomery_setup(n);
k := 1 if v & 1 == 0 else table[a.digit[0] & 7];
if internal_is_negative(p1) {
p1.sign = .Zero_or_Positive;
if internal_is_negative(a1) {
k = -k;
}
}
internal_zero(r) or_return;
for {
if internal_is_zero(a1) {
if internal_eq(p1, 1) {
return k, nil;
} else {
return 0, nil;
}
}
v = internal_count_lsb(a1) or_return;
internal_shr(a1, a1, v) or_return;
if v & 1 == 1 {
k = k * table[p1.digit[0] & 7];
}
if internal_is_negative(a1) {
/*
Compute k = (-1)^((a1)*(p1-1)/4) * k.
a1.digit[0] + 1 cannot overflow because the MSB
of the DIGIT type is not set by definition.
*/
if a1.digit[0] + 1 & p1.digit[0] & 2 != 0 {
k = -k;
}
} else {
/*
Compute k = (-1)^((a1-1)*(p1-1)/4) * k.
*/
if a1.digit[0] & p1.digit[0] & 2 != 0 {
k = -k;
}
}
internal_copy(r, a1) or_return;
r.sign = .Zero_or_Positive;
internal_mod(a1, p1, r) or_return;
internal_copy(p1, r) or_return;
}
return;
}
/*

1055
core/math/big/private.odin
View File

File diff suppressed because it is too large Load Diff

301
core/math/big/public.odin
View File

@@ -1,5 +1,3 @@
package math_big
/*
Copyright 2021 Jeroen van Rijn <nom@duclavier.com>.
Made available under Odin's BSD-3 license.
@@ -10,6 +8,9 @@ package math_big
This file contains basic arithmetic operations like `add`, `sub`, `mul`, `div`, ...
*/
package math_big
import "core:intrinsics"
/*
===========================
@@ -384,6 +385,10 @@ digit_log :: proc(a: DIGIT, base: DIGIT) -> (log: int, err: Error) {
}
log :: proc { int_log, digit_log, };
ilog2 :: proc(value: $T) -> (log2: T) {
return (size_of(T) * 8) - intrinsics.count_leading_zeros(value);
}
/*
Calculate `dest = base^power` using a square-multiply algorithm.
*/
@@ -556,6 +561,298 @@ int_compare_magnitude :: proc(a, b: ^Int, allocator := context.allocator) -> (re
return #force_inline internal_cmp_mag(a, b), nil;
}
int_cmp_mag :: int_compare_magnitude;
/*
bool := a < b
*/
int_less_than :: #force_inline proc(a, b: ^Int, allocator := context.allocator) -> (less_than: bool, err: Error) {
assert_if_nil(a, b);
context.allocator = allocator;
internal_clear_if_uninitialized(a, b) or_return;
c: int;
c, err = cmp(a, b);
return c == -1, err;
}
/*
bool := a < b
*/
int_less_than_digit :: #force_inline proc(a: ^Int, b: DIGIT, allocator := context.allocator) -> (less_than: bool, err: Error) {
assert_if_nil(a);
context.allocator = allocator;
internal_clear_if_uninitialized(a) or_return;
c: int;
c, err = cmp(a, b);
return c == -1, err;
}
/*
bool := |a| < |b|
Compares the magnitudes only, ignores the sign.
*/
int_less_than_abs :: #force_inline proc(a, b: ^Int, allocator := context.allocator) -> (less_than: bool, err: Error) {
assert_if_nil(a, b);
context.allocator = allocator;
internal_clear_if_uninitialized(a, b) or_return;
c: int;
c, err = cmp_mag(a, b);
return c == -1, err;
}
less_than :: proc {
int_less_than,
int_less_than_digit,
};
lt :: less_than;
less_than_abs :: proc {
int_less_than_abs,
};
lt_abs :: less_than_abs;
/*
bool := a <= b
*/
int_less_than_or_equal :: #force_inline proc(a, b: ^Int, allocator := context.allocator) -> (less_than_or_equal: bool, err: Error) {
assert_if_nil(a, b);
context.allocator = allocator;
internal_clear_if_uninitialized(a, b) or_return;
c: int;
c, err = cmp(a, b);
return c <= 0, err;
}
/*
bool := a <= b
*/
int_less_than_or_equal_digit :: #force_inline proc(a: ^Int, b: DIGIT, allocator := context.allocator) -> (less_than_or_equal: bool, err: Error) {
assert_if_nil(a);
context.allocator = allocator;
internal_clear_if_uninitialized(a) or_return;
c: int;
c, err = cmp(a, b);
return c <= 0, err;
}
/*
bool := |a| <= |b|
Compares the magnitudes only, ignores the sign.
*/
int_less_than_or_equal_abs :: #force_inline proc(a, b: ^Int, allocator := context.allocator) -> (less_than_or_equal: bool, err: Error) {
assert_if_nil(a, b);
context.allocator = allocator;
internal_clear_if_uninitialized(a, b) or_return;
c: int;
c, err = cmp_mag(a, b);
return c <= 0, err;
}
less_than_or_equal :: proc {
int_less_than_or_equal,
int_less_than_or_equal_digit,
};
lteq :: less_than_or_equal;
less_than_or_equal_abs :: proc {
int_less_than_or_equal_abs,
};
lteq_abs :: less_than_or_equal_abs;
/*
bool := a == b
*/
int_equals :: #force_inline proc(a, b: ^Int, allocator := context.allocator) -> (equals: bool, err: Error) {
assert_if_nil(a, b);
context.allocator = allocator;
internal_clear_if_uninitialized(a, b) or_return;
c: int;
c, err = cmp(a, b);
return c == 0, err;
}
/*
bool := a == b
*/
int_equals_digit :: #force_inline proc(a: ^Int, b: DIGIT, allocator := context.allocator) -> (equals: bool, err: Error) {
assert_if_nil(a);
context.allocator = allocator;
internal_clear_if_uninitialized(a) or_return;
c: int;
c, err = cmp(a, b);
return c == 0, err;
}
/*
bool := |a| == |b|
Compares the magnitudes only, ignores the sign.
*/
int_equals_abs :: #force_inline proc(a, b: ^Int, allocator := context.allocator) -> (equals: bool, err: Error) {
assert_if_nil(a, b);
context.allocator = allocator;
internal_clear_if_uninitialized(a, b) or_return;
c: int;
c, err = cmp_mag(a, b);
return c == 0, err;
}
equals :: proc {
int_equals,
int_equals_digit,
};
eq :: equals;
equals_abs :: proc {
int_equals_abs,
};
eq_abs :: equals_abs;
/*
bool := a >= b
*/
int_greater_than_or_equal :: #force_inline proc(a, b: ^Int, allocator := context.allocator) -> (greater_than_or_equal: bool, err: Error) {
assert_if_nil(a, b);
context.allocator = allocator;
internal_clear_if_uninitialized(a, b) or_return;
c: int;
c, err = cmp(a, b);
return c >= 0, err;
}
/*
bool := a >= b
*/
int_greater_than_or_equal_digit :: #force_inline proc(a: ^Int, b: DIGIT, allocator := context.allocator) -> (greater_than_or_equal: bool, err: Error) {
assert_if_nil(a);
context.allocator = allocator;
internal_clear_if_uninitialized(a) or_return;
c: int;
c, err = cmp(a, b);
return c >= 0, err;
}
/*
bool := |a| >= |b|
Compares the magnitudes only, ignores the sign.
*/
int_greater_than_or_equal_abs :: #force_inline proc(a, b: ^Int, allocator := context.allocator) -> (greater_than_or_equal: bool, err: Error) {
assert_if_nil(a, b);
context.allocator = allocator;
internal_clear_if_uninitialized(a, b) or_return;
c: int;
c, err = cmp_mag(a, b);
return c >= 0, err;
}
greater_than_or_equal :: proc {
int_greater_than_or_equal,
int_greater_than_or_equal_digit,
};
gteq :: greater_than_or_equal;
greater_than_or_equal_abs :: proc {
int_greater_than_or_equal_abs,
};
gteq_abs :: greater_than_or_equal_abs;
/*
bool := a > b
*/
int_greater_than :: #force_inline proc(a, b: ^Int, allocator := context.allocator) -> (greater_than: bool, err: Error) {
assert_if_nil(a, b);
context.allocator = allocator;
internal_clear_if_uninitialized(a, b) or_return;
c: int;
c, err = cmp(a, b);
return c > 0, err;
}
/*
bool := a > b
*/
int_greater_than_digit :: #force_inline proc(a: ^Int, b: DIGIT, allocator := context.allocator) -> (greater_than: bool, err: Error) {
assert_if_nil(a);
context.allocator = allocator;
internal_clear_if_uninitialized(a) or_return;
c: int;
c, err = cmp(a, b);
return c > 0, err;
}
/*
bool := |a| > |b|
Compares the magnitudes only, ignores the sign.
*/
int_greater_than_abs :: #force_inline proc(a, b: ^Int, allocator := context.allocator) -> (greater_than: bool, err: Error) {
assert_if_nil(a, b);
context.allocator = allocator;
internal_clear_if_uninitialized(a, b) or_return;
c: int;
c, err = cmp_mag(a, b);
return c > 0, err;
}
greater_than :: proc {
int_greater_than,
int_greater_than_digit,
};
gt :: greater_than;
greater_than_abs :: proc {
int_greater_than_abs,
};
gt_abs :: greater_than_abs;
/*
Check if remainders are possible squares - fast exclude non-squares.

3
core/math/big/radix.odin
View File

@@ -1,5 +1,3 @@
package math_big
/*
Copyright 2021 Jeroen van Rijn <nom@duclavier.com>.
Made available under Odin's BSD-3 license.
@@ -14,6 +12,7 @@ package math_big
- Use Barrett reduction for non-powers-of-two.
- Also look at extracting and splatting several digits at once.
*/
package math_big
import "core:intrinsics"
import "core:mem"

3
core/math/big/test.odin
View File

@@ -1,6 +1,4 @@
//+ignore
package math_big
/*
Copyright 2021 Jeroen van Rijn <nom@duclavier.com>.
Made available under Odin's BSD-3 license.
@@ -11,6 +9,7 @@ package math_big
This file exports procedures for use with the test.py test suite.
*/
package math_big
/*
TODO: Write tests for `internal_*` and test reusing parameters with the public implementations.

10
core/math/big/test.py
View File

@@ -1,3 +1,12 @@
#
# Copyright 2021 Jeroen van Rijn <nom@duclavier.com>.
# Made available under Odin's BSD-3 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.
#
from ctypes import *
from random import *
import math
@@ -404,7 +413,6 @@ def test_shr_signed(a = 0, bits = 0, expected_error = Error.Okay):
return test("test_shr_signed", res, [a, bits], expected_error, expected_result)
def test_factorial(number = 0, expected_error = Error.Okay):
print("Factorial:", number)
args = [number]
try:
res = int_factorial(*args)

3
core/math/big/tune.odin
View File

@@ -1,6 +1,4 @@
//+ignore
package math_big
/*
Copyright 2021 Jeroen van Rijn <nom@duclavier.com>.
Made available under Odin's BSD-3 license.
@@ -9,6 +7,7 @@ package math_big
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.
*/
package math_big
import "core:fmt"
import "core:time"