big: Add binary split factorial.

core/math/big/basic.odin

@@ -751,12 +751,15 @@ sqrmod :: proc { int_sqrmod, };
 
 
 int_factorial :: proc(res: ^Int, n: DIGIT) -> (err: Error) {
-	if n < 0 { return .Invalid_Argument; }
+	if n < 0 || n > _FACTORIAL_MAX_N || res == nil { return .Invalid_Argument; }
 
 	i := DIGIT(len(_factorial_table));
 	if n < i {
 		return set(res, _factorial_table[n]);
 	}
+	if n >= _FACTORIAL_BINARY_SPLIT_CUTOFF {
+		return int_factorial_binary_split(res, n);
+	}
 
 	a := &Int{};
 	defer destroy(a);
@@ -771,6 +774,58 @@ int_factorial :: proc(res: ^Int, n: DIGIT) -> (err: Error) {
 
 	return .None;
 }
+
+_int_recursive_product :: proc(res: ^Int, start, stop: DIGIT, level := int(0)) -> (err: Error) {
+	t1, t2 := &Int{}, &Int{};
+	defer destroy(t1, t2);
+
+	if level > _FACTORIAL_BINARY_SPLIT_MAX_RECURSIONS { return .Max_Iterations_Reached; }
+
+	num_factors := (stop - start) >> 1;
+	if num_factors == 2 {
+		if err = set(t1, start); err != .None { return err; }
+		if err = add(t2, t1, 2); err != .None { return err; }
+		return mul(res, t1, t2);
+	}
+
+	if num_factors > 1 {
+		mid := (start + num_factors) | 1;
+		if err = _int_recursive_product(t1, start,  mid, level + 1); err != .None { return err; }
+		if err = _int_recursive_product(t2,   mid, stop, level + 1); err != .None { return err; }
+		return mul(res, t1, t2);
+	}
+
+	if num_factors == 1 { return set(res, start); }
+
+	return one(res);
+}
+
+/*
+	Binary split factorial algo due to: http://www.luschny.de/math/factorial/binarysplitfact.html
+*/
+int_factorial_binary_split :: proc(res: ^Int, n: DIGIT) -> (err: Error) {
+	if n < 0 || n > _FACTORIAL_MAX_N || res == nil { return .Invalid_Argument; }
+
+	inner, outer, start, stop, temp := &Int{}, &Int{}, &Int{}, &Int{}, &Int{};
+	defer destroy(inner, outer, start, stop, temp);
+
+	if err = one(inner); err != .None { return err; }
+	if err = one(outer); err != .None { return err; }
+
+	bits_used := int(_DIGIT_TYPE_BITS - intrinsics.count_leading_zeros(n));
+
+	for i := bits_used; i >= 0; i -= 1 {
+		start := (n >> (uint(i) + 1)) + 1 | 1;
+		stop  := (n >> uint(i)) + 1 | 1;
+		if err = _int_recursive_product(temp, start, stop); err != .None { return err; }
+		if err = mul(inner, inner, temp);                   err != .None { return err; }
+		if err = mul(outer, outer, inner);                  err != .None { return err; }
+	}
+	shift := n - intrinsics.count_ones(n);
+
+	return shl(res, outer, int(shift));
+}
+
 factorial :: proc { int_factorial, };
 
 /*

core/math/big/build.bat

@@ -1,5 +1,6 @@
 @echo off
 odin run . -vet
+: -o:size -no-bounds-check
 :odin build . -build-mode:shared -show-timings -o:minimal -use-separate-modules
 :odin build . -build-mode:shared -show-timings -o:size -use-separate-modules -no-bounds-check
 :odin build . -build-mode:shared -show-timings -o:size -use-separate-modules

core/math/big/common.odin

@@ -26,7 +26,6 @@ when _LOW_MEMORY {
 	_DEFAULT_DIGIT_COUNT :: 32;
 }
 
-
 _MUL_KARATSUBA_CUTOFF :: #config(MUL_KARATSUBA_CUTOFF, _DEFAULT_MUL_KARATSUBA_CUTOFF);
 _SQR_KARATSUBA_CUTOFF :: #config(SQR_KARATSUBA_CUTOFF, _DEFAULT_SQR_KARATSUBA_CUTOFF);
 _MUL_TOOM_CUTOFF      :: #config(MUL_TOOM_CUTOFF,      _DEFAULT_MUL_TOOM_CUTOFF);
@@ -43,6 +42,17 @@ _DEFAULT_SQR_TOOM_CUTOFF      :: 400;
 
 _MAX_ITERATIONS_ROOT_N        :: 500;
 
+/*
+	Largest `N` for which we'll compute `N!`
+*/
+_FACTORIAL_MAX_N              :: 100_000;
+
+/*
+	Cutoff to switch to int_factorial_binary_split, and its max recursion level.
+*/
+_FACTORIAL_BINARY_SPLIT_CUTOFF         :: 6100;
+_FACTORIAL_BINARY_SPLIT_MAX_RECURSIONS :: 100;
+
 Sign :: enum u8 {
 	Zero_or_Positive = 0,
 	Negative         = 1,

core/math/big/example.odin

@@ -42,7 +42,7 @@ _SQR_TOOM_CUTOFF,
 }
 
 print_timings :: proc() {
-	fmt.printf("\nTimings:\n");
+	fmt.printf("Timings:\n");
 	for v, i in Timings {
 		if v.c > 0 {
 			avg   := time.Duration(f64(v.t) / f64(v.c));
@@ -76,6 +76,7 @@ Category :: enum {
 	itoa,
 	atoi,
 	factorial,
+	factorial_bin,
 	choose,
 	lsb,
 	ctz,
@@ -116,10 +117,11 @@ demo :: proc() {
 	defer destroy(a, b, c, d, e, f);
 
 	s := time.tick_now();
-	err = choose(a, 65535, 255);
+	err = choose(a, 1024, 255);
 	Timings[.choose].t += time.tick_since(s); Timings[.choose].c += 1;
-	print("choose", a);
+
-	fmt.println(err);
+	print("1024 choose 255", a, 10, true, true, true);
+	fmt.printf("Error: %v\n", err);
 }
 
 main :: proc() {

core/math/big/helpers.odin

@@ -337,9 +337,7 @@ zero  :: clear;
 	Set the `Int` to 1 and optionally shrink it to the minimum backing size.
 */
 int_one :: proc(a: ^Int, minimize := false, allocator := context.allocator) -> (err: Error) {
-	if err = clear(a, minimize, allocator); err != .None {
+	if err = clear(a, minimize, allocator); err != .None { return err; }
-		return err;
-	}
 
 	a.used     = 1;
 	a.digit[0] = 1;
@@ -429,9 +427,7 @@ get_i32 :: proc { int_get_i32, };
 	and maybe return max(T), .Integer_Overflow if not?
 */
 int_get :: proc(a: ^Int, $T: typeid) -> (res: T, err: Error) where intrinsics.type_is_integer(T) {
-	if err = clear_if_uninitialized(a); err != .None {
+	if err = clear_if_uninitialized(a); err != .None { return 0, err; }
-		return 0, err;
-	}
 
 	size_in_bits := int(size_of(T) * 8);
 	i := int((size_in_bits + _DIGIT_BITS - 1) / _DIGIT_BITS);
@@ -657,3 +653,9 @@ clamp :: proc(a: ^Int) -> (err: Error) {
 	}
 	return .None;
 }
+
+
+_STATIC_ZERO := &Int{
+	used = 0,
+	sign = .Zero_or_Positive,
+};

core/math/big/logical.odin

@@ -16,25 +16,21 @@ import "core:mem"
 /*
 	The `and`, `or` and `xor` binops differ in two lines only.
 	We could handle those with a switch, but that adds overhead.
+
+	TODO: Implement versions that take a DIGIT immediate.
 */
 
 /*
 	2's complement `and`, returns `dest = a & b;`
 */
-and :: proc(dest, a, b: ^Int) -> (err: Error) {
+int_and :: proc(dest, a, b: ^Int) -> (err: Error) {
-	if err = clear_if_uninitialized(a); err != .None {
+	if err = clear_if_uninitialized(a, b); err != .None { return err; }
-		return err;
+
-	}
-	if err = clear_if_uninitialized(b); err != .None {
-		return err;
-	}
 	used := max(a.used, b.used) + 1;
 	/*
 		Grow the destination to accomodate the result.
 	*/
-	if err = grow(dest, used); err != .None {
+	if err = grow(dest, used); err != .None { return err; }
-		return err;
-	}
 
 	neg_a, _ := is_neg(a);
 	neg_b, _ := is_neg(b);
@@ -83,24 +79,18 @@ and :: proc(dest, a, b: ^Int) -> (err: Error) {
 	dest.sign = .Negative if neg else .Zero_or_Positive;
 	return clamp(dest);
 }
+and :: proc { int_add, };
 
 /*
 	2's complement `or`, returns `dest = a | b;`
 */
-or :: proc(dest, a, b: ^Int) -> (err: Error) {
+int_or :: proc(dest, a, b: ^Int) -> (err: Error) {
-	if err = clear_if_uninitialized(a); err != .None {
+	if err = clear_if_uninitialized(a, b); err != .None { return err; }
-		return err;
-	}
-	if err = clear_if_uninitialized(b); err != .None {
-		return err;
-	}
 	used := max(a.used, b.used) + 1;
 	/*
 		Grow the destination to accomodate the result.
 	*/
-	if err = grow(dest, used); err != .None {
+	if err = grow(dest, used); err != .None { return err; }
-		return err;
-	}
 
 	neg_a, _ := is_neg(a);
 	neg_b, _ := is_neg(b);
@@ -149,24 +139,19 @@ or :: proc(dest, a, b: ^Int) -> (err: Error) {
 	dest.sign = .Negative if neg else .Zero_or_Positive;
 	return clamp(dest);
 }
+or :: proc { int_or, };
 
 /*
 	2's complement `xor`, returns `dest = a ~ b;`
 */
-xor :: proc(dest, a, b: ^Int) -> (err: Error) {
+int_xor :: proc(dest, a, b: ^Int) -> (err: Error) {
-	if err = clear_if_uninitialized(a); err != .None {
+	if err = clear_if_uninitialized(a, b); err != .None { return err; }
-		return err;
+
-	}
-	if err = clear_if_uninitialized(b); err != .None {
-		return err;
-	}
 	used := max(a.used, b.used) + 1;
 	/*
 		Grow the destination to accomodate the result.
 	*/
-	if err = grow(dest, used); err != .None {
+	if err = grow(dest, used); err != .None { return err; }
-		return err;
-	}
 
 	neg_a, _ := is_neg(a);
 	neg_b, _ := is_neg(b);
@@ -215,20 +200,16 @@ xor :: proc(dest, a, b: ^Int) -> (err: Error) {
 	dest.sign = .Negative if neg else .Zero_or_Positive;
 	return clamp(dest);
 }
+xor :: proc { int_xor, };
 
 /*
 	dest = ~src
 */
 int_complement :: proc(dest, src: ^Int) -> (err: Error) {
 	/*
-		Check that src and dest are usable.
+		Check that src is usable. Dest will get checked by `sub`.
 	*/
-	if err = clear_if_uninitialized(src); err != .None {
+	if err = clear_if_uninitialized(src); err != .None { return err; }
-		return err;
-	}
-	if err = clear_if_uninitialized(dest); err != .None {
-		return err;
-	}
 
 	/*
 		Temporarily fix sign.
@@ -254,8 +235,7 @@ complement :: proc { int_complement, };
 */
 int_shrmod :: proc(quotient, remainder, numerator: ^Int, bits: int) -> (err: Error) {
 	bits := bits;
-	if err = clear_if_uninitialized(quotient);  err != .None { return err; }
+	if err = clear_if_uninitialized(quotient, numerator);  err != .None { return err; }
-	if err = clear_if_uninitialized(numerator); err != .None { return err; }
 
 	if bits < 0 { return .Invalid_Argument; }
 
@@ -371,8 +351,7 @@ shr_signed :: proc { int_shr_signed, };
 */
 int_shl :: proc(dest, src: ^Int, bits: int) -> (err: Error) {
 	bits := bits;
-	if err = clear_if_uninitialized(src);  err != .None { return err; }
+	if err = clear_if_uninitialized(src, dest);  err != .None { return err; }
-	if err = clear_if_uninitialized(dest); err != .None { return err; }
 
 	if bits < 0 {
 		return .Invalid_Argument;