big: `Add `_private_int_mul_toom`.

core/math/big/build.bat

@@ -1,8 +1,8 @@
 @echo off
-odin run . -vet
+:odin run . -vet
 : -o:size
 :odin build . -build-mode:shared -show-timings -o:minimal -no-bounds-check -define:MATH_BIG_EXE=false && python test.py -fast-tests
 :odin build . -build-mode:shared -show-timings -o:size -no-bounds-check -define:MATH_BIG_EXE=false && python test.py -fast-tests
 :odin build . -build-mode:shared -show-timings -o:size -define:MATH_BIG_EXE=false && python test.py -fast-tests
-:odin build . -build-mode:shared -show-timings -o:speed -no-bounds-check -define:MATH_BIG_EXE=false && python test.py -fast-tests
+odin build . -build-mode:shared -show-timings -o:speed -no-bounds-check -define:MATH_BIG_EXE=false && python test.py -fast-tests
 :odin build . -build-mode:shared -show-timings -o:speed -define:MATH_BIG_EXE=false && python test.py -fast-tests

core/math/big/example.odin

@@ -206,11 +206,14 @@ demo :: proc() {
 	a, b, c, d, e, f := &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{};
 	defer destroy(a, b, c, d, e, f);
 
-	power_of_two(a, 312);
+	atoi(a, "615037959146039477924633848896619112832171971562900618409305032006863881436080", 10);
 	print("a: ", a, 10, true, true, true);
-	power_of_two(b, 314);
+	atoi(b, "378271691190525325893712245607881659587045836991909505715443874842659307597325888631898626653926188084180707310543535657996185416604973577488563643125766400", 10);
 	print("b: ", b, 10, true, true, true);
-	_private_mul_karatsuba(c, a, b);
+
+	// _private_mul_karatsuba(c, a, b);
+	_private_int_mul_toom(c, a, b);
+	// 232651448952541576870611266174879305550351371288854695862580414333123414997160350830885091499735909790287667499899722495800734048928379224433901855785208987458832826418636718381316545267329375006999278984386253755079362097682611712000
 	print("c: ", c, 10, true, true, true);
 }
 

core/math/big/internal.odin

@@ -674,10 +674,10 @@ internal_int_mul :: proc(dest, src, multiplier: ^Int, allocator := context.alloc
 			*/
 							max_used     >= 2 * min_used {
 			// err = s_mp_mul_balance(a,b,c);
-		} else if false && min_used >= MUL_TOOM_CUTOFF {
-			// err = s_mp_mul_toom(a, b, c);
+		} else if min_used >= MUL_TOOM_CUTOFF {
+			err = #force_inline _private_int_mul_toom(dest, src, multiplier);
 		} else if min_used >= MUL_KARATSUBA_CUTOFF {
-			err = #force_inline _private_mul_karatsuba(dest, src, multiplier);
+			err = #force_inline _private_int_mul_karatsuba(dest, src, multiplier);
 		} else if digits < _WARRAY && min_used <= _MAX_COMBA {
 			/*
 				Can we use the fast multiplier?

core/math/big/private.odin

@@ -89,6 +89,143 @@ _private_int_mul :: proc(dest, a, b: ^Int, digits: int, allocator := context.all
 	return internal_clamp(dest);
 }
 
+
+/*
+	Multiplication using the Toom-Cook 3-way algorithm.
+
+	Much more complicated than Karatsuba but has a lower asymptotic running time of O(N**1.464).
+	This algorithm is only particularly useful on VERY large inputs.
+	(We're talking 1000s of digits here...).
+
+	This file contains code from J. Arndt's book  "Matters Computational"
+	and the accompanying FXT-library with permission of the author.
+
+	Setup from:
+		Chung, Jaewook, and M. Anwar Hasan. "Asymmetric squaring formulae."
+		18th IEEE Symposium on Computer Arithmetic (ARITH'07). IEEE, 2007.
+
+	The interpolation from above needed one temporary variable more than the interpolation here:
+
+		Bodrato, Marco, and Alberto Zanoni. "What about Toom-Cook matrices optimality."
+		Centro Vito Volterra Universita di Roma Tor Vergata (2006)
+*/
+_private_int_mul_toom :: proc(dest, a, b: ^Int, allocator := context.allocator) -> (err: Error) {
+	context.allocator = allocator;
+
+	S1, S2, T1, a0, a1, a2, b0, b1, b2 := &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{};
+	defer destroy(S1, S2, T1, a0, a1, a2, b0, b1, b2);
+
+	/*
+		Init temps.
+	*/
+	internal_init_multi(S1, S2, T1)             or_return;
+
+	/*
+		B
+	*/
+	B := min(a.used, b.used) / 3;
+
+	/*
+		a = a2 * x^2 + a1 * x + a0;
+	*/
+	internal_grow(a0, B)                        or_return;
+	internal_grow(a1, B)                        or_return;
+	internal_grow(a2, a.used - 2 * B)           or_return;
+
+	a0.used, a1.used = B, B;
+	a2.used = a.used - 2 * B;
+
+	internal_copy_digits(a0, a, a0.used)        or_return;
+	internal_copy_digits(a1, a, a1.used, B)     or_return;
+	internal_copy_digits(a2, a, a2.used, 2 * B) or_return;
+
+	internal_clamp(a0);
+	internal_clamp(a1);
+	internal_clamp(a2);
+
+	/*
+		b = b2 * x^2 + b1 * x + b0;
+	*/
+	internal_grow(b0, B)                        or_return;
+	internal_grow(b1, B)                        or_return;
+	internal_grow(b2, b.used - 2 * B)           or_return;
+
+	b0.used, b1.used = B, B;
+	b2.used = b.used - 2 * B;
+
+	internal_copy_digits(b0, b, b0.used)        or_return;
+	internal_copy_digits(b1, b, b1.used, B)     or_return;
+	internal_copy_digits(b2, b, b2.used, 2 * B) or_return;
+
+	internal_clamp(b0);
+	internal_clamp(b1);
+	internal_clamp(b2);
+
+	/*
+		\\ S1 = (a2+a1+a0) * (b2+b1+b0);
+	*/
+	internal_add(T1, a2, a1)                    or_return; /*   T1 = a2 + a1; */
+	internal_add(S2, T1, a0)                    or_return; /*   S2 = T1 + a0; */
+	internal_add(dest, b2, b1)                  or_return; /* dest = b2 + b1; */
+	internal_add(S1, dest, b0)                  or_return; /*   S1 =  c + b0; */
+	internal_mul(S1, S1, S2)                    or_return; /*   S1 = S1 * S2; */
+
+	/*
+		\\S2 = (4*a2+2*a1+a0) * (4*b2+2*b1+b0);
+	*/
+	internal_add(T1, T1, a2)                    or_return; /*   T1 = T1 + a2; */
+	internal_int_shl1(T1, T1)                   or_return; /*   T1 = T1 << 1; */
+	internal_add(T1, T1, a0)                    or_return; /*   T1 = T1 + a0; */
+	internal_add(dest, dest, b2)                or_return; /*    c =  c + b2; */
+	internal_int_shl1(dest, dest)               or_return; /*    c =  c << 1; */
+	internal_add(dest, dest, b0)                or_return; /*    c =  c + b0; */
+	internal_mul(S2, T1, dest)                  or_return; /*   S2 = T1 *  c; */
+
+	/*
+		\\S3 = (a2-a1+a0) * (b2-b1+b0);
+	*/
+	internal_sub(a1, a2, a1)                    or_return; /*   a1 = a2 - a1; */
+	internal_add(a1, a1, a0)                    or_return; /*   a1 = a1 + a0; */
+	internal_sub(b1, b2, b1)                    or_return; /*   b1 = b2 - b1; */
+	internal_add(b1, b1, b0)                    or_return; /*   b1 = b1 + b0; */
+	internal_mul(a1, a1, b1)                    or_return; /*   a1 = a1 * b1; */
+	internal_mul(b1, a2, b2)                    or_return; /*   b1 = a2 * b2; */
+
+	/*
+		\\S2 = (S2 - S3) / 3;
+	*/
+	internal_sub(S2, S2, a1)                    or_return; /*   S2 = S2 - a1; */
+	_private_int_div_3(S2, S2)                  or_return; /*   S2 = S2 / 3; \\ this is an exact division  */
+	internal_sub(a1, S1, a1)                    or_return; /*   a1 = S1 - a1; */
+	internal_int_shr1(a1, a1)                   or_return; /*   a1 = a1 >> 1; */
+	internal_mul(a0, a0, b0)                    or_return; /*   a0 = a0 * b0; */
+	internal_sub(S1, S1, a0)                    or_return; /*   S1 = S1 - a0; */
+	internal_sub(S2, S2, S1)                    or_return; /*   S2 = S2 - S1; */
+	internal_int_shr1(S2, S2)                   or_return; /*   S2 = S2 >> 1; */
+	internal_sub(S1, S1, a1)                    or_return; /*   S1 = S1 - a1; */
+	internal_sub(S1, S1, b1)                    or_return; /*   S1 = S1 - b1; */
+	internal_int_shl1(T1, b1)                   or_return; /*   T1 = b1 << 1; */
+	internal_sub(S2, S2, T1)                    or_return; /*   S2 = S2 - T1; */
+	internal_sub(a1, a1, S2)                    or_return; /*   a1 = a1 - S2; */
+
+	/*
+		P = b1*x^4+ S2*x^3+ S1*x^2+ a1*x + a0;
+	*/
+	internal_shl_digit(b1, 4 * B)               or_return;
+	internal_shl_digit(S2, 3 * B)               or_return;
+	internal_add(b1, b1, S2)                    or_return;
+	internal_shl_digit(S1, 2 * B)               or_return;
+	internal_add(b1, b1, S1)                    or_return;
+	internal_shl_digit(a1, 1 * B)               or_return;
+	internal_add(b1, b1, a1)                    or_return;
+	internal_add(dest, b1, a0)                  or_return;
+
+	/*
+		a * b - P
+	*/
+	return nil;
+}
+
 /*
 	product = |a| * |b| using Karatsuba Multiplication using three half size multiplications.
 
@@ -116,7 +253,7 @@ _private_int_mul :: proc(dest, a, b: ^Int, digits: int, allocator := context.all
 	baseline/comba methods use. Generally though, the overhead of this method doesn't pay off
 	until a certain size is reached, of around 80 used DIGITs.
 */
-_private_mul_karatsuba :: proc(dest, a, b: ^Int, allocator := context.allocator) -> (err: Error) {
+_private_int_mul_karatsuba :: proc(dest, a, b: ^Int, allocator := context.allocator) -> (err: Error) {
 	context.allocator = allocator;
 
 	x0, x1, y0, y1, t1, x0y0, x1y1 := &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{};

core/math/big/test.py

@@ -446,6 +446,7 @@ TESTS = {
 	test_mul: [
 		[ 1234,   5432],
 		[ 0xd3b4e926aaba3040e1c12b5ea553b5, 0x1a821e41257ed9281bee5bc7789ea7 ],
+		[ 1 << 21_105, 1 << 21_501 ],
 	],
 	test_sqr: [
 		[ 5432],