Fix typo in big.internal_invmod

Fixes #3253

2 files changed, 3359 insertions, 3416 deletions

diff --git a/core/math/big/internal.odin b/core/math/big/internal.odin
index 35c95f465..03623e7f2 100644
--- a/core/math/big/internal.odin
+++ b/core/math/big/internal.odin
@@ -1181,28 +1181,18 @@ internal_cmp_digit :: internal_compare_digit
 */
 internal_int_compare_magnitude :: #force_inline proc(a, b: ^Int) -> (comparison: int) {
 	assert_if_nil(a, b)
-	/*
-		Compare based on used digits.
-	*/
+
+	// Compare based on used digits.
 	if a.used != b.used {
-		if a.used > b.used {
-			return +1
-		}
-		return -1
+		return +1 if a.used > b.used else -1
 	}
 
-	/*
-		Same number of used digits, compare based on their value.
-	*/
+	// Same number of used digits, compare based on their value.
 	#no_bounds_check for n := a.used - 1; n >= 0; n -= 1 {
 		if a.digit[n] != b.digit[n] {
-			if a.digit[n] > b.digit[n] {
-				return +1
-			}
-			return -1
+			return +1 if a.digit[n] > b.digit[n] else -1
 		}
 	}
-
 	return 0
 }
 internal_compare_magnitude :: proc { internal_int_compare_magnitude, }
diff --git a/core/math/big/private.odin b/core/math/big/private.odin
index d045b4239..2ee6cfafa 100644
--- a/core/math/big/private.odin
+++ b/core/math/big/private.odin
@@ -1,3402 +1,3355 @@
-/*

-	Copyright 2021 Jeroen van Rijn <nom@duclavier.com>.

-	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.

-

-	=============================    Private procedures    =============================

-

-	Private procedures used by the above low-level routines follow.

-

-	Don't call these yourself unless you really know what you're doing.

-	They include implementations that are optimimal for certain ranges of input only.

-

-	These aren't exported for the same reasons.

-*/

-

-

-package math_big

-

-import "base:intrinsics"

-import "core:mem"

-

-/*

-	Multiplies |a| * |b| and only computes upto digs digits of result.

-	HAC pp. 595, Algorithm 14.12  Modified so you can control how

-	many digits of output are created.

-*/

-_private_int_mul :: proc(dest, a, b: ^Int, digits: int, allocator := context.allocator) -> (err: Error) {

-	context.allocator = allocator

-

-	/*

-		Can we use the fast multiplier?

-	*/

-	if digits < _WARRAY && min(a.used, b.used) < _MAX_COMBA {

-		return #force_inline _private_int_mul_comba(dest, a, b, digits)

-	}

-

-	/*

-		Set up temporary output `Int`, which we'll swap for `dest` when done.

-	*/

-

-	t := &Int{}

-

-	internal_grow(t, max(digits, _DEFAULT_DIGIT_COUNT)) or_return

-	t.used = digits

-

-	/*

-		Compute the digits of the product directly.

-	*/

-	pa := a.used

-	for ix := 0; ix < pa; ix += 1 {

-		/*

-			Limit ourselves to `digits` DIGITs of output.

-		*/

-		pb    := min(b.used, digits - ix)

-		carry := _WORD(0)

-		iy    := 0

-

-		/*

-			Compute the column of the output and propagate the carry.

-		*/

-		#no_bounds_check for iy = 0; iy < pb; iy += 1 {

-			/*

-				Compute the column as a _WORD.

-			*/

-			column := _WORD(t.digit[ix + iy]) + _WORD(a.digit[ix]) * _WORD(b.digit[iy]) + carry

-

-			/*

-				The new column is the lower part of the result.

-			*/

-			t.digit[ix + iy] = DIGIT(column & _WORD(_MASK))

-

-			/*

-				Get the carry word from the result.

-			*/

-			carry = column >> _DIGIT_BITS

-		}

-		/*

-			Set carry if it is placed below digits

-		*/

-		if ix + iy < digits {

-			t.digit[ix + pb] = DIGIT(carry)

-		}

-	}

-

-	internal_swap(dest, t)

-	internal_destroy(t)

-	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 internal_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;

-	*/

-	_private_int_shl_leg(b1, 4 * B)             or_return

-	_private_int_shl_leg(S2, 3 * B)             or_return

-	internal_add(b1, b1, S2)                    or_return

-	_private_int_shl_leg(S1, 2 * B)             or_return

-	internal_add(b1, b1, S1)                    or_return

-	_private_int_shl_leg(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.

-

-	Let `B` represent the radix [e.g. 2**_DIGIT_BITS] and let `n` represent

-	half of the number of digits in the min(a,b)

-

-	`a` = `a1` * `B`**`n` + `a0`

-	`b` = `b`1 * `B`**`n` + `b0`

-

-	Then, a * b => 1b1 * B**2n + ((a1 + a0)(b1 + b0) - (a0b0 + a1b1)) * B + a0b0

-

-	Note that a1b1 and a0b0 are used twice and only need to be computed once.

-	So in total three half size (half # of digit) multiplications are performed,

-		a0b0, a1b1 and (a1+b1)(a0+b0)

-

-	Note that a multiplication of half the digits requires 1/4th the number of

-	single precision multiplications, so in total after one call 25% of the

-	single precision multiplications are saved.

-

-	Note also that the call to `internal_mul` can end up back in this function

-	if the a0, a1, b0, or b1 are above the threshold.

-

-	This is known as divide-and-conquer and leads to the famous O(N**lg(3)) or O(N**1.584)

-	work which is asymptopically lower than the standard O(N**2) that the

-	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_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{}

-	defer internal_destroy(x0, x1, y0, y1, t1, x0y0, x1y1)

-

-	/*

-		min # of digits, divided by two.

-	*/

-	B := min(a.used, b.used) >> 1

-

-	/*

-		Init all the temps.

-	*/

-	internal_grow(x0, B)          or_return

-	internal_grow(x1, a.used - B) or_return

-	internal_grow(y0, B)          or_return

-	internal_grow(y1, b.used - B) or_return

-	internal_grow(t1, B * 2)      or_return

-	internal_grow(x0y0, B * 2)    or_return

-	internal_grow(x1y1, B * 2)    or_return

-

-	/*

-		Now shift the digits.

-	*/

-	x0.used, y0.used = B, B

-	x1.used = a.used - B

-	y1.used = b.used - B

-

-	/*

-		We copy the digits directly instead of using higher level functions

-		since we also need to shift the digits.

-	*/

-	internal_copy_digits(x0, a, x0.used)

-	internal_copy_digits(y0, b, y0.used)

-	internal_copy_digits(x1, a, x1.used, B)

-	internal_copy_digits(y1, b, y1.used, B)

-

-	/*

-		Only need to clamp the lower words since by definition the

-		upper words x1/y1 must have a known number of digits.

-	*/

-	clamp(x0)

-	clamp(y0)

-

-	/*

-		Now calc the products x0y0 and x1y1,

-		after this x0 is no longer required, free temp [x0==t2]!

-	*/

-	internal_mul(x0y0, x0, y0)      or_return /* x0y0 = x0*y0 */

-	internal_mul(x1y1, x1, y1)      or_return /* x1y1 = x1*y1 */

-	internal_add(t1,   x1, x0)      or_return /* now calc x1+x0 and */

-	internal_add(x0,   y1, y0)      or_return /* t2 = y1 + y0 */

-	internal_mul(t1,   t1, x0)      or_return /* t1 = (x1 + x0) * (y1 + y0) */

-

-	/*

-		Add x0y0.

-	*/

-	internal_add(x0, x0y0, x1y1)    or_return /* t2 = x0y0 + x1y1 */

-	internal_sub(t1,   t1,   x0)    or_return /* t1 = (x1+x0)*(y1+y0) - (x1y1 + x0y0) */

-

-	/*

-		shift by B.

-	*/

-	_private_int_shl_leg(t1, B)       or_return /* t1 = (x0y0 + x1y1 - (x1-x0)*(y1-y0))<<B */

-	_private_int_shl_leg(x1y1, B * 2) or_return /* x1y1 = x1y1 << 2*B */

-

-	internal_add(t1, x0y0, t1)      or_return /* t1 = x0y0 + t1 */

-	internal_add(dest, t1, x1y1)    or_return /* t1 = x0y0 + t1 + x1y1 */

-

-	return nil

-}

-

-

-

-/*

-	Fast (comba) multiplier

-

-	This is the fast column-array [comba] multiplier.  It is

-	designed to compute the columns of the product first

-	then handle the carries afterwards.  This has the effect

-	of making the nested loops that compute the columns very

-	simple and schedulable on super-scalar processors.

-

-	This has been modified to produce a variable number of

-	digits of output so if say only a half-product is required

-	you don't have to compute the upper half (a feature

-	required for fast Barrett reduction).

-

-	Based on Algorithm 14.12 on pp.595 of HAC.

-*/

-_private_int_mul_comba :: proc(dest, a, b: ^Int, digits: int, allocator := context.allocator) -> (err: Error) {

-	context.allocator = allocator

-

-	/*

-		Set up array.

-	*/

-	W: [_WARRAY]DIGIT = ---

-

-	/*

-		Grow the destination as required.

-	*/

-	internal_grow(dest, digits) or_return

-

-	/*

-		Number of output digits to produce.

-	*/

-	pa := min(digits, a.used + b.used)

-

-	/*

-		Clear the carry

-	*/

-	_W := _WORD(0)

-

-	ix: int

-	for ix = 0; ix < pa; ix += 1 {

-		tx, ty, iy, iz: int

-

-		/*

-			Get offsets into the two bignums.

-		*/

-		ty = min(b.used - 1, ix)

-		tx = ix - ty

-

-		/*

-			This is the number of times the loop will iterate, essentially.

-			while (tx++ < a->used && ty-- >= 0) { ... }

-		*/

-		 

-		iy = min(a.used - tx, ty + 1)

-

-		/*

-			Execute loop.

-		*/

-		#no_bounds_check for iz = 0; iz < iy; iz += 1 {

-			_W += _WORD(a.digit[tx + iz]) * _WORD(b.digit[ty - iz])

-		}

-

-		/*

-			Store term.

-		*/

-		W[ix] = DIGIT(_W) & _MASK

-

-		/*

-			Make next carry.

-		*/

-		_W = _W >> _WORD(_DIGIT_BITS)

-	}

-

-	/*

-		Setup dest.

-	*/

-	old_used := dest.used

-	dest.used = pa

-

-	/*

-		Now extract the previous digit [below the carry].

-	*/

-	copy_slice(dest.digit[0:], W[:pa])	

-

-	/*

-		Clear unused digits [that existed in the old copy of dest].

-	*/

-	internal_zero_unused(dest, old_used)

-

-	/*

-		Adjust dest.used based on leading zeroes.

-	*/

-

-	return internal_clamp(dest)

-}

-

-/*

-	Multiplies |a| * |b| and does not compute the lower digs digits

-	[meant to get the higher part of the product]

-*/

-_private_int_mul_high :: proc(dest, a, b: ^Int, digits: int, allocator := context.allocator) -> (err: Error) {

-	context.allocator = allocator

-

-	/*

-		Can we use the fast multiplier?

-	*/

-	if a.used + b.used + 1 < _WARRAY && min(a.used, b.used) < _MAX_COMBA {

-		return _private_int_mul_high_comba(dest, a, b, digits)

-	}

-

-	internal_grow(dest, a.used + b.used + 1) or_return

-	dest.used = a.used + b.used + 1

-

-	pa := a.used

-	pb := b.used

-	for ix := 0; ix < pa; ix += 1 {

-		carry := DIGIT(0)

-

-		for iy := digits - ix; iy < pb; iy += 1 {

-			/*

-				Calculate the double precision result.

-			*/

-			r := _WORD(dest.digit[ix + iy]) + _WORD(a.digit[ix]) * _WORD(b.digit[iy]) + _WORD(carry)

-

-			/*

-				Get the lower part.

-			*/

-			dest.digit[ix + iy] = DIGIT(r & _WORD(_MASK))

-

-			/*

-				Carry the carry.

-			*/

-			carry = DIGIT(r >> _WORD(_DIGIT_BITS))

-		}

-		dest.digit[ix + pb] = carry

-	}

-	return internal_clamp(dest)

-}

-

-/*

-	This is a modified version of `_private_int_mul_comba` that only produces output digits *above* `digits`.

-	See the comments for `_private_int_mul_comba` to see how it works.

-

-	This is used in the Barrett reduction since for one of the multiplications

-	only the higher digits were needed.  This essentially halves the work.

-

-	Based on Algorithm 14.12 on pp.595 of HAC.

-*/

-_private_int_mul_high_comba :: proc(dest, a, b: ^Int, digits: int, allocator := context.allocator) -> (err: Error) {

-	context.allocator = allocator

-

-	W: [_WARRAY]DIGIT = ---

-	_W: _WORD = 0

-

-	/*

-		Number of output digits to produce. Grow the destination as required.

-	*/

-	pa := a.used + b.used

-	internal_grow(dest, pa) or_return

-

-	ix: int

-	for ix = digits; ix < pa; ix += 1 {

-		/*

-			Get offsets into the two bignums.

-		*/

-		ty := min(b.used - 1, ix)

-		tx := ix - ty

-

-		/*

-			This is the number of times the loop will iterrate, essentially it's

-			while (tx++ < a->used && ty-- >= 0) { ... }

-		*/

-		iy := min(a.used - tx, ty + 1)

-

-		/*

-			Execute loop.

-		*/

-		for iz := 0; iz < iy; iz += 1 {

-			_W += _WORD(a.digit[tx + iz]) * _WORD(b.digit[ty - iz])

-		}

-

-		/*

-			Store term.

-		*/

-		W[ix] = DIGIT(_W) & DIGIT(_MASK)

-

-		/*

-			Make next carry.

-		*/

-		_W = _W >> _WORD(_DIGIT_BITS)

-	}

-

-	/*

-		Setup dest

-	*/

-	old_used := dest.used

-	dest.used = pa

-

-	for ix = digits; ix < pa; ix += 1 {

-		/*

-			Now extract the previous digit [below the carry].

-		*/

-		dest.digit[ix] = W[ix]

-	}

-

-	/*

-		Zero remainder.

-	*/

-	internal_zero_unused(dest, old_used)

-

-	/*

-		Adjust dest.used based on leading zeroes.

-	*/

-	return internal_clamp(dest)

-}

-

-/*

-	Single-digit multiplication with the smaller number as the single-digit.

-*/

-_private_int_mul_balance :: proc(dest, a, b: ^Int, allocator := context.allocator) -> (err: Error) {

-	context.allocator = allocator

-	a, b := a, b

-

-	a0, tmp, r := &Int{}, &Int{}, &Int{}

-	defer internal_destroy(a0, tmp, r)

-

-	b_size   := min(a.used, b.used)

-	n_blocks := max(a.used, b.used) / b_size

-

-	internal_grow(a0, b_size + 2) or_return

-	internal_init_multi(tmp, r)   or_return

-

-	/*

-		Make sure that `a` is the larger one.

-	*/

-	if a.used < b.used {

-		a, b = b, a

-	}

-	assert(a.used >= b.used)

-

-	i, j := 0, 0

-	for ; i < n_blocks; i += 1 {

-		/*

-			Cut a slice off of `a`.

-		*/

-

-		a0.used = b_size

-		internal_copy_digits(a0, a, a0.used, j)

-		j += a0.used

-		internal_clamp(a0)

-

-		/*

-			Multiply with `b`.

-		*/

-		internal_mul(tmp, a0, b)                                     or_return

-

-		/*

-			Shift `tmp` to the correct position.

-		*/

-		_private_int_shl_leg(tmp, b_size * i)                          or_return

-

-		/*

-			Add to output. No carry needed.

-		*/

-		internal_add(r, r, tmp)                                      or_return

-	}

-

-	/*

-		The left-overs; there are always left-overs.

-	*/

-	if j < a.used {

-		a0.used = a.used - j

-		internal_copy_digits(a0, a, a0.used, j)

-		j += a0.used

-		internal_clamp(a0)

-

-		internal_mul(tmp, a0, b)                                     or_return

-		_private_int_shl_leg(tmp, b_size * i)                          or_return

-		internal_add(r, r, tmp)                                      or_return

-	}

-

-	internal_swap(dest, r)

-	return

-}

-

-/*

-	Low level squaring, b = a*a, HAC pp.596-597, Algorithm 14.16

-	Assumes `dest` and `src` to not be `nil`, and `src` to have been initialized.

-*/

-_private_int_sqr :: proc(dest, src: ^Int, allocator := context.allocator) -> (err: Error) {

-	context.allocator = allocator

-	pa := src.used

-

-	t := &Int{}; ix, iy: int

-	/*

-		Grow `t` to maximum needed size, or `_DEFAULT_DIGIT_COUNT`, whichever is bigger.

-	*/

-	internal_grow(t, max((2 * pa) + 1, _DEFAULT_DIGIT_COUNT)) or_return

-	t.used = (2 * pa) + 1

-

-	#no_bounds_check for ix = 0; ix < pa; ix += 1 {

-		carry := DIGIT(0)

-		/*

-			First calculate the digit at 2*ix; calculate double precision result.

-		*/

-		r := _WORD(t.digit[ix+ix]) + (_WORD(src.digit[ix]) * _WORD(src.digit[ix]))

-

-		/*

-			Store lower part in result.

-		*/

-		t.digit[ix+ix] = DIGIT(r & _WORD(_MASK))

-		/*

-			Get the carry.

-		*/

-		carry = DIGIT(r >> _DIGIT_BITS)

-

-		#no_bounds_check for iy = ix + 1; iy < pa; iy += 1 {

-			/*

-				First calculate the product.

-			*/

-			r = _WORD(src.digit[ix]) * _WORD(src.digit[iy])

-

-			/* Now calculate the double precision result. Nóte we use

-			 * addition instead of *2 since it's easier to optimize

-			 */

-			r = _WORD(t.digit[ix+iy]) + r + r + _WORD(carry)

-

-			/*

-				Store lower part.

-			*/

-			t.digit[ix+iy] = DIGIT(r & _WORD(_MASK))

-

-			/*

-				Get carry.

-			*/

-			carry = DIGIT(r >> _DIGIT_BITS)

-		}

-		/*

-			Propagate upwards.

-		*/

-		#no_bounds_check for carry != 0 {

-			r     = _WORD(t.digit[ix+iy]) + _WORD(carry)

-			t.digit[ix+iy] = DIGIT(r & _WORD(_MASK))

-			carry = DIGIT(r >> _WORD(_DIGIT_BITS))

-			iy += 1

-		}

-	}

-

-	err = internal_clamp(t)

-	internal_swap(dest, t)

-	internal_destroy(t)

-	return err

-}

-

-/*

-	The jist of squaring...

-	You do like mult except the offset of the tmpx [one that starts closer to zero] can't equal the offset of tmpy.

-	So basically you set up iy like before then you min it with (ty-tx) so that it never happens.

-	You double all those you add in the inner loop. After that loop you do the squares and add them in.

-

-	Assumes `dest` and `src` not to be `nil` and `src` to have been initialized.	

-*/

-_private_int_sqr_comba :: proc(dest, src: ^Int, allocator := context.allocator) -> (err: Error) {

-	context.allocator = allocator

-

-	W: [_WARRAY]DIGIT = ---

-

-	/*

-		Grow the destination as required.

-	*/

-	pa := uint(src.used) + uint(src.used)

-	internal_grow(dest, int(pa)) or_return

-

-	/*

-		Number of output digits to produce.

-	*/

-	W1 := _WORD(0)

-	_W  : _WORD = ---

-	ix := uint(0)

-

-	#no_bounds_check for ; ix < pa; ix += 1 {

-		/*

-			Clear counter.

-		*/

-		_W = {}

-

-		/*

-			Get offsets into the two bignums.

-		*/

-		ty := min(uint(src.used) - 1, ix)

-		tx := ix - ty

-

-		/*

-			This is the number of times the loop will iterate,

-			essentially while (tx++ < a->used && ty-- >= 0) { ... }

-		*/

-		iy := min(uint(src.used) - tx, ty + 1)

-

-		/*

-			Now for squaring, tx can never equal ty.

-			We halve the distance since they approach at a rate of 2x,

-			and we have to round because odd cases need to be executed.

-		*/

-		iy = min(iy, ((ty - tx) + 1) >> 1 )

-

-		/*

-			Execute loop.

-		*/

-		#no_bounds_check for iz := uint(0); iz < iy; iz += 1 {

-			_W += _WORD(src.digit[tx + iz]) * _WORD(src.digit[ty - iz])

-		}

-

-		/*

-			Double the inner product and add carry.

-		*/

-		_W = _W + _W + W1

-

-		/*

-			Even columns have the square term in them.

-		*/

-		if ix & 1 == 0 {

-			_W += _WORD(src.digit[ix >> 1]) * _WORD(src.digit[ix >> 1])

-		}

-

-		/*

-			Store it.

-		*/

-		W[ix] = DIGIT(_W & _WORD(_MASK))

-

-		/*

-			Make next carry.

-		*/

-		W1 = _W >> _DIGIT_BITS

-	}

-

-	/*

-		Setup dest.

-	*/

-	old_used := dest.used

-	dest.used = src.used + src.used

-

-	#no_bounds_check for ix = 0; ix < pa; ix += 1 {

-		dest.digit[ix] = W[ix] & _MASK

-	}

-

-	/*

-		Clear unused digits [that existed in the old copy of dest].

-	*/

-	internal_zero_unused(dest, old_used)

-

-	return internal_clamp(dest)

-}

-

-/*

-	Karatsuba squaring, computes `dest` = `src` * `src` using three half-size squarings.

- 

- 	See comments of `_private_int_mul_karatsuba` for details.

- 	It is essentially the same algorithm but merely tuned to perform recursive squarings.

-*/

-_private_int_sqr_karatsuba :: proc(dest, src: ^Int, allocator := context.allocator) -> (err: Error) {

-	context.allocator = allocator

-

-	x0, x1, t1, t2, x0x0, x1x1 := &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}

-	defer internal_destroy(x0, x1, t1, t2, x0x0, x1x1)

-

-	/*

-		Min # of digits, divided by two.

-	*/

-	B := src.used >> 1

-

-	/*

-		Init temps.

-	*/

-	internal_grow(x0,   B) or_return

-	internal_grow(x1,   src.used - B) or_return

-	internal_grow(t1,   src.used * 2) or_return

-	internal_grow(t2,   src.used * 2) or_return

-	internal_grow(x0x0, B * 2       ) or_return

-	internal_grow(x1x1, (src.used - B) * 2) or_return

-

-	/*

-		Now shift the digits.

-	*/

-	x0.used = B

-	x1.used = src.used - B

-

-	#force_inline internal_copy_digits(x0, src, x0.used)

-	#force_inline mem.copy_non_overlapping(&x1.digit[0], &src.digit[B], size_of(DIGIT) * x1.used)

-	#force_inline internal_clamp(x0)

-

-	/*

-		Now calc the products x0*x0 and x1*x1.

-	*/

-	internal_sqr(x0x0, x0) or_return

-	internal_sqr(x1x1, x1) or_return

-

-	/*

-		Now calc (x1+x0)^2

-	*/

-	internal_add(t1, x0, x1) or_return

-	internal_sqr(t1, t1) or_return

-

-	/*

-		Add x0y0

-	*/

-	internal_add(t2, x0x0, x1x1) or_return

-	internal_sub(t1, t1, t2) or_return

-

-	/*

-		Shift by B.

-	*/

-	_private_int_shl_leg(t1, B) or_return

-	_private_int_shl_leg(x1x1, B * 2) or_return

-	internal_add(t1, t1, x0x0) or_return

-	internal_add(dest, t1, x1x1) or_return

-

-	return #force_inline internal_clamp(dest)

-}

-

-/*

-	Squaring using Toom-Cook 3-way algorithm.

-

-	Setup and interpolation from algorithm SQR_3 in Chung, Jaewook, and M. Anwar Hasan. "Asymmetric squaring formulae."

-	  18th IEEE Symposium on Computer Arithmetic (ARITH'07). IEEE, 2007.

-*/

-_private_int_sqr_toom :: proc(dest, src: ^Int, allocator := context.allocator) -> (err: Error) {

-	context.allocator = allocator

-

-	S0, a0, a1, a2 := &Int{}, &Int{}, &Int{}, &Int{}

-	defer internal_destroy(S0, a0, a1, a2)

-

-	/*

-		Init temps.

-	*/

-	internal_zero(S0) or_return

-

-	/*

-		B

-	*/

-	B := src.used / 3

-

-	/*

-		a = a2 * x^2 + a1 * x + a0;

-	*/

-	internal_grow(a0, B) or_return

-	internal_grow(a1, B) or_return

-	internal_grow(a2, src.used - (2 * B)) or_return

-

-	a0.used = B

-	a1.used = B

-	a2.used = src.used - 2 * B

-

-	#force_inline mem.copy_non_overlapping(&a0.digit[0], &src.digit[    0], size_of(DIGIT) * a0.used)

-	#force_inline mem.copy_non_overlapping(&a1.digit[0], &src.digit[    B], size_of(DIGIT) * a1.used)

-	#force_inline mem.copy_non_overlapping(&a2.digit[0], &src.digit[2 * B], size_of(DIGIT) * a2.used)

-

-	internal_clamp(a0)

-	internal_clamp(a1)

-	internal_clamp(a2)

-

-	/** S0 = a0^2;  */

-	internal_sqr(S0, a0) or_return

-

-	/** \\S1 = (a2 + a1 + a0)^2 */

-	/** \\S2 = (a2 - a1 + a0)^2  */

-	/** \\S1 = a0 + a2; */

-	/** a0 = a0 + a2; */

-	internal_add(a0, a0, a2) or_return

-	/** \\S2 = S1 - a1; */

-	/** b = a0 - a1; */

-	internal_sub(dest, a0, a1) or_return

-	/** \\S1 = S1 + a1; */

-	/** a0 = a0 + a1; */

-	internal_add(a0, a0, a1) or_return

-	/** \\S1 = S1^2;  */

-	/** a0 = a0^2; */

-	internal_sqr(a0, a0) or_return

-	/** \\S2 = S2^2;  */

-	/** b = b^2; */

-	internal_sqr(dest, dest) or_return

-	/** \\ S3 = 2 * a1 * a2  */

-	/** \\S3 = a1 * a2;  */

-	/** a1 = a1 * a2; */

-	internal_mul(a1, a1, a2) or_return

-	/** \\S3 = S3 << 1;  */

-	/** a1 = a1 << 1; */

-	internal_shl(a1, a1, 1) or_return

-	/** \\S4 = a2^2;  */

-	/** a2 = a2^2; */

-	internal_sqr(a2, a2) or_return

-	/** \\ tmp = (S1 + S2)/2  */

-	/** \\tmp = S1 + S2; */

-	/** b = a0 + b; */

-	internal_add(dest, a0, dest) or_return

-	/** \\tmp = tmp >> 1; */

-	/** b = b >> 1; */

-	internal_shr(dest, dest, 1) or_return

-	/** \\ S1 = S1 - tmp - S3  */

-	/** \\S1 = S1 - tmp; */

-	/** a0 = a0 - b; */

-	internal_sub(a0, a0, dest) or_return

-	/** \\S1 = S1 - S3;  */

-	/** a0 = a0 - a1; */

-	internal_sub(a0, a0, a1) or_return

-	/** \\S2 = tmp - S4 -S0  */

-	/** \\S2 = tmp - S4;  */

-	/** b = b - a2; */

-	internal_sub(dest, dest, a2) or_return

-	/** \\S2 = S2 - S0;  */

-	/** b = b - S0; */

-	internal_sub(dest, dest, S0) or_return

-	/** \\P = S4*x^4 + S3*x^3 + S2*x^2 + S1*x + S0; */

-	/** P = a2*x^4 + a1*x^3 + b*x^2 + a0*x + S0; */

-	_private_int_shl_leg(  a2, 4 * B) or_return

-	_private_int_shl_leg(  a1, 3 * B) or_return

-	_private_int_shl_leg(dest, 2 * B) or_return

-	_private_int_shl_leg(  a0, 1 * B) or_return

-

-	internal_add(a2, a2, a1) or_return

-	internal_add(dest, dest, a2) or_return

-	internal_add(dest, dest, a0) or_return

-	internal_add(dest, dest, S0) or_return

-	/** a^2 - P  */

-

-	return #force_inline internal_clamp(dest)

-}

-

-/*

-	Divide by three (based on routine from MPI and the GMP manual).

-*/

-_private_int_div_3 :: proc(quotient, numerator: ^Int, allocator := context.allocator) -> (remainder: DIGIT, err: Error) {

-	context.allocator = allocator

-

-	/*

-		b = 2^_DIGIT_BITS / 3

-	*/

- 	b := _WORD(1) << _WORD(_DIGIT_BITS) / _WORD(3)

-

-	q := &Int{}

-	internal_grow(q, numerator.used) or_return

-	q.used = numerator.used

-	q.sign = numerator.sign

-

-	w, t: _WORD

-	#no_bounds_check for ix := numerator.used; ix >= 0; ix -= 1 {

-		w = (w << _WORD(_DIGIT_BITS)) | _WORD(numerator.digit[ix])

-		if w >= 3 {

-			/*

-				Multiply w by [1/3].

-			*/

-			t = (w * b) >> _WORD(_DIGIT_BITS)

-

-			/*

-				Now subtract 3 * [w/3] from w, to get the remainder.

-			*/

-			w -= t+t+t

-

-			/*

-				Fixup the remainder as required since the optimization is not exact.

-			*/

-			for w >= 3 {

-				t += 1

-				w -= 3

-			}

-		} else {

-			t = 0

-		}

-		q.digit[ix] = DIGIT(t)

-	}

-	remainder = DIGIT(w)

-

-	/*

-		[optional] store the quotient.

-	*/

-	if quotient != nil {

-		err = clamp(q)

- 		internal_swap(q, quotient)

- 	}

-	internal_destroy(q)

-	return remainder, nil

-}

-

-/*

-	Signed Integer Division

-

-	c*b + d == a [i.e. a/b, c=quotient, d=remainder], HAC pp.598 Algorithm 14.20

-

-	Note that the description in HAC is horribly incomplete.

-	For example, it doesn't consider the case where digits are removed from 'x' in

-	the inner loop.

-

-	It also doesn't consider the case that y has fewer than three digits, etc.

-	The overall algorithm is as described as 14.20 from HAC but fixed to treat these cases.

-*/

-_private_int_div_school :: proc(quotient, remainder, numerator, denominator: ^Int, allocator := context.allocator) -> (err: Error) {

-	context.allocator = allocator

-

-	error_if_immutable(quotient, remainder) or_return

-

-	q, x, y, t1, t2 := &Int{}, &Int{}, &Int{}, &Int{}, &Int{}

-	defer internal_destroy(q, x, y, t1, t2)

-

-	internal_grow(q, numerator.used + 2) or_return

-	q.used = numerator.used + 2

-

-	internal_init_multi(t1, t2) or_return

-	internal_copy(x, numerator) or_return

-	internal_copy(y, denominator) or_return

-

-	/*

-		Fix the sign.

-	*/

-	neg   := numerator.sign != denominator.sign

-	x.sign = .Zero_or_Positive

-	y.sign = .Zero_or_Positive

-

-	/*

-		Normalize both x and y, ensure that y >= b/2, [b == 2**MP_DIGIT_BIT]

-	*/

-	norm := internal_count_bits(y) % _DIGIT_BITS

-

-	if norm < _DIGIT_BITS - 1 {

-		norm = (_DIGIT_BITS - 1) - norm

-		internal_shl(x, x, norm) or_return

-		internal_shl(y, y, norm) or_return

-	} else {

-		norm = 0

-	}

-

-	/*

-		Note: HAC does 0 based, so if used==5 then it's 0,1,2,3,4, i.e. use 4

-	*/

-	n := x.used - 1

-	t := y.used - 1

-

-	/*

-		while (x >= y*b**n-t) do { q[n-t] += 1; x -= y*b**{n-t} }

-		y = y*b**{n-t}

-	*/

-

-	_private_int_shl_leg(y, n - t) or_return

-

-	gte := internal_gte(x, y)

-	for gte {

-		q.digit[n - t] += 1

-		internal_sub(x, x, y) or_return

-		gte = internal_gte(x, y)

-	}

-

-	/*

-		Reset y by shifting it back down.

-	*/

-	_private_int_shr_leg(y, n - t)

-

-	/*

-		Step 3. for i from n down to (t + 1).

-	*/

-	#no_bounds_check for i := n; i >= (t + 1); i -= 1 {

-		if i > x.used { continue }

-

-		/*

-			step 3.1 if xi == yt then set q{i-t-1} to b-1, otherwise set q{i-t-1} to (xi*b + x{i-1})/yt

-		*/

-		if x.digit[i] == y.digit[t] {

-			q.digit[(i - t) - 1] = 1 << (_DIGIT_BITS - 1)

-		} else {

-

-			tmp := _WORD(x.digit[i]) << _DIGIT_BITS

-			tmp |= _WORD(x.digit[i - 1])

-			tmp /= _WORD(y.digit[t])

-			if tmp > _WORD(_MASK) {

-				tmp = _WORD(_MASK)

-			}

-			q.digit[(i - t) - 1] = DIGIT(tmp & _WORD(_MASK))

-		}

-

-		/* while (q{i-t-1} * (yt * b + y{t-1})) >

-					xi * b**2 + xi-1 * b + xi-2

-

-			do q{i-t-1} -= 1;

-		*/

-

-		iter := 0

-

-		q.digit[(i - t) - 1] = (q.digit[(i - t) - 1] + 1) & _MASK

-		#no_bounds_check for {

-			q.digit[(i - t) - 1] = (q.digit[(i - t) - 1] - 1) & _MASK

-

-			/*

-				Find left hand.

-			*/

-			internal_zero(t1)

-			t1.digit[0] = ((t - 1) < 0) ? 0 : y.digit[t - 1]

-			t1.digit[1] = y.digit[t]

-			t1.used = 2

-			internal_mul(t1, t1, q.digit[(i - t) - 1]) or_return

-

-			/*

-				Find right hand.

-			*/

-			t2.digit[0] = ((i - 2) < 0) ? 0 : x.digit[i - 2]

-			t2.digit[1] = x.digit[i - 1] /* i >= 1 always holds */

-			t2.digit[2] = x.digit[i]

-			t2.used = 3

-

-			if internal_lte(t1, t2) {

-				break

-			}

-			iter += 1; if iter > 100 {

-				return .Max_Iterations_Reached

-			}

-		}

-

-		/*

-			Step 3.3 x = x - q{i-t-1} * y * b**{i-t-1}

-		*/

-		int_mul_digit(t1, y, q.digit[(i - t) - 1]) or_return

-		_private_int_shl_leg(t1, (i - t) - 1) or_return

-		internal_sub(x, x, t1) or_return

-

-		/*

-			if x < 0 then { x = x + y*b**{i-t-1}; q{i-t-1} -= 1; }

-		*/

-		if x.sign == .Negative {

-			internal_copy(t1, y) or_return

-			_private_int_shl_leg(t1, (i - t) - 1) or_return

-			internal_add(x, x, t1) or_return

-

-			q.digit[(i - t) - 1] = (q.digit[(i - t) - 1] - 1) & _MASK

-		}

-	}

-

-	/*

-		Now q is the quotient and x is the remainder, [which we have to normalize]

-		Get sign before writing to c.

-	*/

-	z, _ := is_zero(x)

-	x.sign = .Zero_or_Positive if z else numerator.sign

-

-	if quotient != nil {

-		internal_clamp(q)

-		internal_swap(q, quotient)

-		quotient.sign = .Negative if neg else .Zero_or_Positive

-	}

-

-	if remainder != nil {

-		internal_shr(x, x, norm) or_return

-		internal_swap(x, remainder)

-	}

-

-	return nil

-}

-

-/*

-	Direct implementation of algorithms 1.8 "RecursiveDivRem" and 1.9 "UnbalancedDivision" from:

-

-		Brent, Richard P., and Paul Zimmermann. "Modern computer arithmetic"

-		Vol. 18. Cambridge University Press, 2010

-		Available online at https://arxiv.org/pdf/1004.4710

-

-	pages 19ff. in the above online document.

-*/

-_private_div_recursion :: proc(quotient, remainder, a, b: ^Int, allocator := context.allocator) -> (err: Error) {

-	context.allocator = allocator

-

-	A1, A2, B1, B0, Q1, Q0, R1, R0, t := &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}

-	defer internal_destroy(A1, A2, B1, B0, Q1, Q0, R1, R0, t)

-

-	m := a.used - b.used

-	k := m / 2

-

-	if m < MUL_KARATSUBA_CUTOFF {

-		return _private_int_div_school(quotient, remainder, a, b)

-	}

-

-	internal_init_multi(A1, A2, B1, B0, Q1, Q0, R1, R0, t) or_return

-

-	/*

-		`B1` = `b` / `beta`^`k`, `B0` = `b` % `beta`^`k`

-	*/

-	internal_shrmod(B1, B0, b, k * _DIGIT_BITS) or_return

-

-	/*

-		(Q1, R1) =  RecursiveDivRem(A / beta^(2k), B1)

-	*/

-	internal_shrmod(A1, t, a, 2 * k * _DIGIT_BITS) or_return

-	_private_div_recursion(Q1, R1, A1, B1) or_return

-

-	/*

-		A1 = (R1 * beta^(2k)) + (A % beta^(2k)) - (Q1 * B0 * beta^k)

-	*/

-	_private_int_shl_leg(R1, 2 * k) or_return

-	internal_add(A1, R1, t) or_return

-	internal_mul(t, Q1, B0) or_return

-

-	/*

-		While A1 < 0 do Q1 = Q1 - 1, A1 = A1 + (beta^k * B)

-	*/

-	if internal_lt(A1, 0) {

-		internal_shl(t, b, k * _DIGIT_BITS) or_return

-

-		for {

-			internal_decr(Q1) or_return

-			internal_add(A1, A1, t) or_return

-			if internal_gte(A1, 0) { break }

-		}

-	}

-

-	/*

-		(Q0, R0) =  RecursiveDivRem(A1 / beta^(k), B1)

-	*/

-	internal_shrmod(A1, t, A1, k * _DIGIT_BITS) or_return

-	_private_div_recursion(Q0, R0, A1, B1) or_return

-

-	/*

-		A2 = (R0*beta^k) +  (A1 % beta^k) - (Q0*B0)

-	*/

-	_private_int_shl_leg(R0, k) or_return

-	internal_add(A2, R0, t) or_return

-	internal_mul(t, Q0, B0) or_return

-	internal_sub(A2, A2, t) or_return

-

-	/*

-		While A2 < 0 do Q0 = Q0 - 1, A2 = A2 + B.

-	*/

-	for internal_is_negative(A2) { // internal_lt(A2, 0) {

-		internal_decr(Q0) or_return

-		internal_add(A2, A2, b) or_return

-	}

-

-	/*

-		Return q = (Q1*beta^k) + Q0, r = A2.

-	*/

-	_private_int_shl_leg(Q1, k) or_return

-	internal_add(quotient, Q1, Q0) or_return

-

-	return internal_copy(remainder, A2)

-}

-

-_private_int_div_recursive :: proc(quotient, remainder, a, b: ^Int, allocator := context.allocator) -> (err: Error) {

-	context.allocator = allocator

-

-	A, B, Q, Q1, R, A_div, A_mod := &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}

-	defer internal_destroy(A, B, Q, Q1, R, A_div, A_mod)

-

-	internal_init_multi(A, B, Q, Q1, R, A_div, A_mod) or_return

-

-	/*

-		Most significant bit of a limb.

-		Assumes  _DIGIT_MAX < (sizeof(DIGIT) * sizeof(u8)).

-	*/

-	msb := (_DIGIT_MAX + DIGIT(1)) >> 1

-	sigma := 0

-	msb_b := b.digit[b.used - 1]

-	for msb_b < msb {

-		sigma += 1

-		msb_b <<= 1

-	}

-

-	/*

-		Use that sigma to normalize B.

-	*/

-	internal_shl(B, b, sigma) or_return

-	internal_shl(A, a, sigma) or_return

-

-	/*

-		Fix the sign.

-	*/

-	neg := a.sign != b.sign

-	A.sign = .Zero_or_Positive; B.sign = .Zero_or_Positive

-

-	/*

-		If the magnitude of "A" is not more more than twice that of "B" we can work

-		on them directly, otherwise we need to work at "A" in chunks.

-	*/

-	n := B.used

-	m := A.used - B.used

-

-	/*

-		Q = 0. We already ensured that when we called `internal_init_multi`.

-	*/

-	for m > n {

-		/*

-			(q, r) = RecursiveDivRem(A / (beta^(m-n)), B)

-		*/

-		j := (m - n) * _DIGIT_BITS

-		internal_shrmod(A_div, A_mod, A, j) or_return

-		_private_div_recursion(Q1, R, A_div, B) or_return

-

-		/*

-			Q = (Q*beta!(n)) + q

-		*/

-		internal_shl(Q, Q, n * _DIGIT_BITS) or_return

-		internal_add(Q, Q, Q1) or_return

-

-		/*

-			A = (r * beta^(m-n)) + (A % beta^(m-n))

-		*/

-		internal_shl(R, R, (m - n) * _DIGIT_BITS) or_return

-		internal_add(A, R, A_mod) or_return

-

-		/*

-			m = m - n

-		*/

-		m -= n

-	}

-

-	/*

-		(q, r) = RecursiveDivRem(A, B)

-	*/

-	_private_div_recursion(Q1, R, A, B) or_return

-

-	/*

-		Q = (Q * beta^m) + q, R = r

-	*/

-	internal_shl(Q, Q, m * _DIGIT_BITS) or_return

-	internal_add(Q, Q, Q1) or_return

-

-	/*

-		Get sign before writing to dest.

-	*/

-	R.sign = .Zero_or_Positive if internal_is_zero(Q) else a.sign

-

-	if quotient != nil {

-		swap(quotient, Q)

-		quotient.sign = .Negative if neg else .Zero_or_Positive

-	}

-	if remainder != nil {

-		/*

-			De-normalize the remainder.

-		*/

-		internal_shrmod(R, nil, R, sigma) or_return

-		swap(remainder, R)

-	}

-	return nil

-}

-

-/*

-	Slower bit-bang division... also smaller.

-*/

-@(deprecated="Use `_int_div_school`, it's 3.5x faster.")

-_private_int_div_small :: proc(quotient, remainder, numerator, denominator: ^Int) -> (err: Error) {

-

-	ta, tb, tq, q := &Int{}, &Int{}, &Int{}, &Int{}

-

-	defer internal_destroy(ta, tb, tq, q)

-

-	for {

-		internal_one(tq) or_return

-

-		num_bits, _ := count_bits(numerator)

-		den_bits, _ := count_bits(denominator)

-		n := num_bits - den_bits

-

-		abs(ta, numerator)   or_return

-		abs(tb, denominator) or_return

-		shl(tb, tb, n)       or_return

-		shl(tq, tq, n)       or_return

-

-		for n >= 0 {

-			if internal_gte(ta, tb) {

-				// ta -= tb

-				sub(ta, ta, tb) or_return

-				//  q += tq

-				add( q, q,  tq) or_return

-			}

-			shr1(tb, tb) or_return

-			shr1(tq, tq) or_return

-

-			n -= 1

-		}

-

-		/*

-			Now q == quotient and ta == remainder.

-		*/

-		neg := numerator.sign != denominator.sign

-		if quotient != nil {

-			swap(quotient, q)

-			z, _ := is_zero(quotient)

-			quotient.sign = .Negative if neg && !z else .Zero_or_Positive

-		}

-		if remainder != nil {

-			swap(remainder, ta)

-			z, _ := is_zero(numerator)

-			remainder.sign = .Zero_or_Positive if z else numerator.sign

-		}

-

-		break

-	}

-	return err

-}

-

-

-

-/*

-	Binary split factorial algo due to: http://www.luschny.de/math/factorial/binarysplitfact.html

-*/

-_private_int_factorial_binary_split :: proc(res: ^Int, n: int, allocator := context.allocator) -> (err: Error) {

-	context.allocator = allocator

-

-	inner, outer, start, stop, temp := &Int{}, &Int{}, &Int{}, &Int{}, &Int{}

-	defer internal_destroy(inner, outer, start, stop, temp)

-

-	internal_one(inner, false)                                       or_return

-	internal_one(outer, false)                                       or_return

-

-	bits_used := ilog2(n)

-

-	for i := bits_used; i >= 0; i -= 1 {

-		start := (n >> (uint(i) + 1)) + 1 | 1

-		stop  := (n >> uint(i)) + 1 | 1

-		_private_int_recursive_product(temp, start, stop, 0)         or_return

-		internal_mul(inner, inner, temp)                             or_return

-		internal_mul(outer, outer, inner)                            or_return

-	}

-	shift := n - intrinsics.count_ones(n)

-

-	return internal_shl(res, outer, int(shift))

-}

-

-/*

-	Recursive product used by binary split factorial algorithm.

-*/

-_private_int_recursive_product :: proc(res: ^Int, start, stop: int, level := int(0), allocator := context.allocator) -> (err: Error) {

-	context.allocator = allocator

-

-	t1, t2 := &Int{}, &Int{}

-	defer internal_destroy(t1, t2)

-

-	if level > FACTORIAL_BINARY_SPLIT_MAX_RECURSIONS {

-		return .Max_Iterations_Reached

-	}

-

-	num_factors := (stop - start) >> 1

-	if num_factors == 2 {

-		internal_set(t1, start, false)                               or_return

-		when true {

-			internal_grow(t2, t1.used + 1, false)                    or_return

-			internal_add(t2, t1, 2)                                  or_return

-		} else {

-			internal_add(t2, t1, 2)                                  or_return

-		}

-		return internal_mul(res, t1, t2)

-	}

-

-	if num_factors > 1 {

-		mid := (start + num_factors) | 1

-		_private_int_recursive_product(t1, start,  mid, level + 1)   or_return

-		_private_int_recursive_product(t2,   mid, stop, level + 1)   or_return

-		return internal_mul(res, t1, t2)

-	}

-

-	if num_factors == 1 {

-		return #force_inline internal_set(res, start, true)

-	}

-

-	return #force_inline internal_one(res, true)

-}

-

-/*

-	Internal function computing both GCD using the binary method,

-		and, if target isn't `nil`, also LCM.

-

-	Expects the `a` and `b` to have been initialized

-		and one or both of `res_gcd` or `res_lcm` not to be `nil`.

-

-	If both `a` and `b` are zero, return zero.

-	If either `a` or `b`, return the other one.

-

-	The `gcd` and `lcm` wrappers have already done this test,

-	but `gcd_lcm` wouldn't have, so we still need to perform it.

-

-	If neither result is wanted, we have nothing to do.

-*/

-_private_int_gcd_lcm :: proc(res_gcd, res_lcm, a, b: ^Int, allocator := context.allocator) -> (err: Error) {

-	context.allocator = allocator

-

-	if res_gcd == nil && res_lcm == nil {

-		return nil

-	}

-

-	/*

-		We need a temporary because `res_gcd` is allowed to be `nil`.

-	*/

-	if a.used == 0 && b.used == 0 {

-		/*

-			GCD(0, 0) and LCM(0, 0) are both 0.

-		*/

-		if res_gcd != nil {

-			internal_zero(res_gcd) or_return

-		}

-		if res_lcm != nil {

-			internal_zero(res_lcm) or_return

-		}

-		return nil

-	} else if a.used == 0 {

-		/*

-			We can early out with GCD = B and LCM = 0

-		*/

-		if res_gcd != nil {

-			internal_abs(res_gcd, b) or_return

-		}

-		if res_lcm != nil {

-			internal_zero(res_lcm) or_return

-		}

-		return nil

-	} else if b.used == 0 {

-		/*

-			We can early out with GCD = A and LCM = 0

-		*/

-		if res_gcd != nil {

-			internal_abs(res_gcd, a) or_return

-		}

-		if res_lcm != nil {

-			internal_zero(res_lcm) or_return

-		}

-		return nil

-	}

-

-	temp_gcd_res := &Int{}

-	defer internal_destroy(temp_gcd_res)

-

-	/*

-		If neither `a` or `b` was zero, we need to compute `gcd`.

- 		Get copies of `a` and `b` we can modify.

- 	*/

-	u, v := &Int{}, &Int{}

-	defer internal_destroy(u, v)

-	internal_copy(u, a) or_return

-	internal_copy(v, b) or_return

-

- 	/*

- 		Must be positive for the remainder of the algorithm.

- 	*/

-	u.sign = .Zero_or_Positive; v.sign = .Zero_or_Positive

-

- 	/*

- 		B1.  Find the common power of two for `u` and `v`.

- 	*/

- 	u_lsb, _ := internal_count_lsb(u)

- 	v_lsb, _ := internal_count_lsb(v)

- 	k        := min(u_lsb, v_lsb)

-

-	if k > 0 {

-		/*

-			Divide the power of two out.

-		*/

-		internal_shr(u, u, k) or_return

-		internal_shr(v, v, k) or_return

-	}

-

-	/*

-		Divide any remaining factors of two out.

-	*/

-	if u_lsb != k {

-		internal_shr(u, u, u_lsb - k) or_return

-	}

-	if v_lsb != k {

-		internal_shr(v, v, v_lsb - k) or_return

-	}

-

-	for v.used != 0 {

-		/*

-			Make sure `v` is the largest.

-		*/

-		if internal_gt(u, v) {

-			/*

-				Swap `u` and `v` to make sure `v` is >= `u`.

-			*/

-			internal_swap(u, v)

-		}

-

-		/*

-			Subtract smallest from largest.

-		*/

-		internal_sub(v, v, u) or_return

-

-		/*

-			Divide out all factors of two.

-		*/

-		b, _ := internal_count_lsb(v)

-		internal_shr(v, v, b) or_return

-	}

-

- 	/*

- 		Multiply by 2**k which we divided out at the beginning.

- 	*/

- 	internal_shl(temp_gcd_res, u, k) or_return

- 	temp_gcd_res.sign = .Zero_or_Positive

-

-	/*

-		We've computed `gcd`, either the long way, or because one of the inputs was zero.

-		If we don't want `lcm`, we're done.

-	*/

-	if res_lcm == nil {

-		internal_swap(temp_gcd_res, res_gcd)

-		return nil

-	}

-

-	/*

-		Computes least common multiple as `|a*b|/gcd(a,b)`

-		Divide the smallest by the GCD.

-	*/

-	if internal_lt_abs(a, b) {

-		/*

-			Store quotient in `t2` such that `t2 * b` is the LCM.

-		*/

-		internal_div(res_lcm, a, temp_gcd_res) or_return

-		err = internal_mul(res_lcm, res_lcm, b)

-	} else {

-		/*

-			Store quotient in `t2` such that `t2 * a` is the LCM.

-		*/

-		internal_div(res_lcm, b, temp_gcd_res) or_return

-		err = internal_mul(res_lcm, res_lcm, a)

-	}

-

-	if res_gcd != nil {

-		internal_swap(temp_gcd_res, res_gcd)

-	}

-

-	/*

-		Fix the sign to positive and return.

-	*/

-	res_lcm.sign = .Zero_or_Positive

-	return err

-}

-

-/*

-	Internal implementation of log.

-	Assumes `a` not to be `nil` and to have been initialized.

-*/

-_private_int_log :: proc(a: ^Int, base: DIGIT, allocator := context.allocator) -> (res: int, err: Error) {

-	bracket_low, bracket_high, bracket_mid, t, bi_base := &Int{}, &Int{}, &Int{}, &Int{}, &Int{}

-	defer internal_destroy(bracket_low, bracket_high, bracket_mid, t, bi_base)

-

-	ic := #force_inline internal_cmp(a, base)

-	if ic == -1 || ic == 0 {

-		return 1 if ic == 0 else 0, nil

-	}

-	defer if err != nil {

-		res = -1

-	}

-

-	internal_set(bi_base, base, true, allocator)       or_return

-	internal_clear(bracket_mid, false, allocator)      or_return

-	internal_clear(t, false, allocator)                or_return

-	internal_one(bracket_low, false, allocator)        or_return

-	internal_set(bracket_high, base, false, allocator) or_return

-

-	low := 0; high := 1

-

-	/*

-		A kind of Giant-step/baby-step algorithm.

-		Idea shamelessly stolen from https://programmingpraxis.com/2010/05/07/integer-logarithms/2/

-		The effect is asymptotic, hence needs benchmarks to test if the Giant-step should be skipped

-		for small n.

-	*/

-

-	for {

-		/*

-			Iterate until `a` is bracketed between low + high.

-		*/

-		if #force_inline internal_gte(bracket_high, a) { break }

-

-		low = high

-		#force_inline internal_copy(bracket_low, bracket_high) or_return

-		high <<= 1

-		#force_inline internal_sqr(bracket_high, bracket_high) or_return

-	}

-

-	for (high - low) > 1 {

-		mid := (high + low) >> 1

-

-		#force_inline internal_pow(t, bi_base, mid - low) or_return

-

-		#force_inline internal_mul(bracket_mid, bracket_low, t) or_return

-

-		mc := #force_inline internal_cmp(a, bracket_mid)

-		switch mc {

-		case -1:

-			high = mid

-			internal_swap(bracket_mid, bracket_high)

-		case  0:

-			return mid, nil

-		case  1:

-			low = mid

-			internal_swap(bracket_mid, bracket_low)

-		}

-	}

-

-	fc := #force_inline internal_cmp(bracket_high, a)

-	res = high if fc == 0 else low

-

-	return

-}

-

-/*

-	Computes xR**-1 == x (mod N) via Montgomery Reduction.

-	This is an optimized implementation of `internal_montgomery_reduce`

-	which uses the comba method to quickly calculate the columns of the reduction.

-	Based on Algorithm 14.32 on pp.601 of HAC.

-*/

-_private_montgomery_reduce_comba :: proc(x, n: ^Int, rho: DIGIT, allocator := context.allocator) -> (err: Error) {

-	context.allocator = allocator

-	W: [_WARRAY]_WORD = ---

-

-	if x.used > _WARRAY { return .Invalid_Argument }

-

-	/*

-		Get old used count.

-	*/

-	old_used := x.used

-

-	/*

-		Grow `x` as required.

-	*/

-	internal_grow(x, n.used + 1)                                     or_return

-

-	/*

-		First we have to get the digits of the input into an array of double precision words W[...]

-		Copy the digits of `x` into W[0..`x.used` - 1]

-	*/

-	ix: int

-	for ix = 0; ix < x.used; ix += 1 {

-		W[ix] = _WORD(x.digit[ix])

-	}

-

-	/*

-		Zero the high words of W[a->used..m->used*2].

-	*/

-	zero_upper := (n.used * 2) + 1

-	if ix < zero_upper {

-		for ix = x.used; ix < zero_upper; ix += 1 {

-			W[ix] = {}

-		}

-	}

-

-	/*

-		Now we proceed to zero successive digits from the least significant upwards.

-	*/

-	for ix = 0; ix < n.used; ix += 1 {

-		/*

-			`mu = ai * m' mod b`

-

-			We avoid a double precision multiplication (which isn't required)

-			by casting the value down to a DIGIT.  Note this requires

-			that W[ix-1] have the carry cleared (see after the inner loop)

-		*/

-		mu := ((W[ix] & _WORD(_MASK)) * _WORD(rho)) & _WORD(_MASK)

-

-		/*

-			`a = a + mu * m * b**i`

-		

-			This is computed in place and on the fly.  The multiplication

-		 	by b**i is handled by offseting which columns the results

-		 	are added to.

-		

-			Note the comba method normally doesn't handle carries in the

-			inner loop In this case we fix the carry from the previous

-			column since the Montgomery reduction requires digits of the

-			result (so far) [see above] to work.

-

-			This is	handled by fixing up one carry after the inner loop.

-			The carry fixups are done in order so after these loops the

-			first m->used words of W[] have the carries fixed.

-		*/

-		for iy := 0; iy < n.used; iy += 1 {

-			W[ix + iy] += mu * _WORD(n.digit[iy])

-		}

-

-		/*

-			Now fix carry for next digit, W[ix+1].

-		*/

-		W[ix + 1] += (W[ix] >> _DIGIT_BITS)

-	}

-

-	/*

-		Now we have to propagate the carries and shift the words downward

-		[all those least significant digits we zeroed].

-	*/

-

-	for ; ix < n.used * 2; ix += 1 {

-		W[ix + 1] += (W[ix] >> _DIGIT_BITS)

-	}

-

-	/* copy out, A = A/b**n

-	 *

-	 * The result is A/b**n but instead of converting from an

-	 * array of mp_word to mp_digit than calling mp_rshd

-	 * we just copy them in the right order

-	 */

-

-	for ix = 0; ix < (n.used + 1); ix += 1 {

-		x.digit[ix] = DIGIT(W[n.used + ix] & _WORD(_MASK))

-	}

-

-	/*

-		Set the max used.

-	*/

-	x.used = n.used + 1

-

-	/*

-		Zero old_used digits, if the input a was larger than m->used+1 we'll have to clear the digits.

-	*/

-	internal_zero_unused(x, old_used)

-	internal_clamp(x)

-

-	/*

-		if A >= m then A = A - m

-	*/

-	if internal_gte_abs(x, n) {

-		return internal_sub(x, x, n)

-	}

-	return nil

-}

-

-/*

-	Computes xR**-1 == x (mod N) via Montgomery Reduction.

-	Assumes `x` and `n` not to be nil.

-*/

-_private_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.

-	*/

-	internal_clear_if_uninitialized(x, n) or_return

-

-	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)

-	_private_int_shr_leg(x, n.used)

-

-	/*

-		if x >= n then x = x - n

-	*/

-	if internal_gte_abs(x, n) {

-		return internal_sub(x, x, n)

-	}

-

-	return nil

-}

-

-/*

-	Shifts with subtractions when the result is greater than b.

-

-	The method is slightly modified to shift B unconditionally upto just under

-	the leading bit of b.  This saves alot of multiple precision shifting.

-

-	Assumes `a` and `b` not to be `nil`.

-*/

-_private_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.

-	*/

-	internal_clear_if_uninitialized(a, b) or_return

-

-	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)                                              or_return

-		bits = 1

-	}

-

-	/*

-		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_gte_abs(a, b) {

-			internal_sub(a, a, b)                                    or_return

-		}

-	}

-	return nil

-}

-

-/*

-	Sets up the Montgomery reduction stuff.

-*/

-_private_int_montgomery_setup :: proc(n: ^Int, allocator := context.allocator) -> (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

-	*/

-	internal_clear_if_uninitialized(n, allocator) or_return

-

-	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 _DIGIT_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 */

-	}

-

-	/*

-		rho = -1/m mod b

-	*/

-	rho = DIGIT(((_WORD(1) << _WORD(_DIGIT_BITS)) - _WORD(x)) & _WORD(_MASK))

-	return rho, nil

-}

-

-/*

-	Reduces `x` mod `m`, assumes 0 < x < m**2, mu is precomputed via reduce_setup.

-	From HAC pp.604 Algorithm 14.42

-

-	Assumes `x`, `m` and `mu` all not to be `nil` and have been initialized.

-*/

-_private_int_reduce :: proc(x, m, mu: ^Int, allocator := context.allocator) -> (err: Error) {

-	context.allocator = allocator

-

-	q := &Int{}

-	defer internal_destroy(q)

-	um := m.used

-

-	/*

-		q = x

-	*/

-	internal_copy(q, x)                                              or_return

-

-	/*

-		q1 = x / b**(k-1)

-	*/

-	_private_int_shr_leg(q, um - 1)

-

-	/*

-		According to HAC this optimization is ok.

-	*/

-	if DIGIT(um) > DIGIT(1) << (_DIGIT_BITS - 1) {

-		internal_mul(q, q, mu)                                       or_return

-	} else {

-		_private_int_mul_high(q, q, mu, um)                          or_return

-	}

-

-	/*

-		q3 = q2 / b**(k+1)

-	*/

-	_private_int_shr_leg(q, um + 1)

-

-	/*

-		x = x mod b**(k+1), quick (no division)

-	*/

-	internal_int_mod_bits(x, x, _DIGIT_BITS * (um + 1))              or_return

-

-	/*

-		q = q * m mod b**(k+1), quick (no division)

-	*/

-	_private_int_mul(q, q, m, um + 1)                                or_return

-

-	/*

-		x = x - q

-	*/

-	internal_sub(x, x, q)                                            or_return

-

-	/*

-		If x < 0, add b**(k+1) to it.

-	*/

-	if internal_is_negative(x) {

-		internal_set(q, 1)                                           or_return

-		_private_int_shl_leg(q, um + 1)                                or_return

-		internal_add(x, x, q)                                        or_return

-	}

-

-	/*

-		Back off if it's too big.

-	*/

-	for internal_gte(x, m) {

-		internal_sub(x, x, m)                                        or_return

-	}

-

-	return nil

-}

-

-/*

-	Reduces `a` modulo `n`, where `n` is of the form 2**p - d.

-*/

-_private_int_reduce_2k :: proc(a, n: ^Int, d: DIGIT, allocator := context.allocator) -> (err: Error) {

-	context.allocator = allocator

-

-	q := &Int{}

-	defer internal_destroy(q)

-

-	internal_zero(q)                                                 or_return

-

-	p := internal_count_bits(n)

-

-	for {

-		/*

-			q = a/2**p, a = a mod 2**p

-		*/

-		internal_shrmod(q, a, a, p)                                  or_return

-

-		if d != 1 {

-			/*

-				q = q * d

-			*/

-			internal_mul(q, q, d)                                    or_return

-		}

-

-		/*

-			a = a + q

-		*/

-		internal_add(a, a, q)                                        or_return

-		if internal_lt_abs(a, n)                                     { break }

-		internal_sub(a, a, n)                                        or_return

-	}

-

-	return nil

-}

-

-/*

-	Reduces `a` modulo `n` where `n` is of the form 2**p - d

-	This differs from reduce_2k since "d" can be larger than a single digit.

-*/

-_private_int_reduce_2k_l :: proc(a, n, d: ^Int, allocator := context.allocator) -> (err: Error) {

-	context.allocator = allocator

-

-	q := &Int{}

-	defer internal_destroy(q)

-

-	internal_zero(q)                                                 or_return

-

-	p := internal_count_bits(n)

-

-	for {

-		/*

-			q = a/2**p, a = a mod 2**p

-		*/

-		internal_shrmod(q, a, a, p)                                  or_return

-

-		/*

-			q = q * d

-		*/

-		internal_mul(q, q, d)                                        or_return

-

-		/*

-			a = a + q

-		*/

-		internal_add(a, a, q)                                        or_return

-		if internal_lt_abs(a, n)                                     { break }

-		internal_sub(a, a, n)                                        or_return

-	}

-

-	return nil

-}

-

-/*

-	Determines if `internal_int_reduce_2k` can be used.

-	Asssumes `a` not to be `nil` and to have been initialized.

-*/

-_private_int_reduce_is_2k :: proc(a: ^Int) -> (reducible: bool, err: Error) {

-	assert_if_nil(a)

-

-	if internal_is_zero(a) {

-		return false, nil

-	} else if a.used == 1 {

-		return true, nil

-	} else if a.used  > 1 {

-		iy := internal_count_bits(a)

-		iw := 1

-		iz := DIGIT(1)

-

-		/*

-			Test every bit from the second digit up, must be 1.

-		*/

-		for ix := _DIGIT_BITS; ix < iy; ix += 1 {

-			if a.digit[iw] & iz == 0 {

-				return false, nil

-			}

-

-			iz <<= 1

-			if iz > _DIGIT_MAX {

-				iw += 1

-				iz  = 1

-			}

-		}

-		return true, nil

-	} else {

-		return true, nil

-	}

-}

-

-/*

-	Determines if `internal_int_reduce_2k_l` can be used.

-	Asssumes `a` not to be `nil` and to have been initialized.

-*/

-_private_int_reduce_is_2k_l :: proc(a: ^Int) -> (reducible: bool, err: Error) {

-	assert_if_nil(a)

-

-	if internal_int_is_zero(a) {

-		return false, nil

-	} else if a.used == 1 {

-		return true, nil

-	} else if a.used  > 1 {

-		/*

-			If more than half of the digits are -1 we're sold.

-		*/

-		ix := 0

-		iy := 0

-

-		for ; ix < a.used; ix += 1 {

-			if a.digit[ix] == _DIGIT_MAX {

-				iy += 1

-			}

-		}

-		return iy >= (a.used / 2), nil

-	} else {

-		return false, nil

-	}

-}

-

-/*

-	Determines the setup value.

-	Assumes `a` is not `nil`.

-*/

-_private_int_reduce_2k_setup :: proc(a: ^Int, allocator := context.allocator) -> (d: DIGIT, err: Error) {

-	context.allocator = allocator

-

-	tmp := &Int{}

-	defer internal_destroy(tmp)

-	internal_zero(tmp)                                               or_return

-

-	internal_int_power_of_two(tmp, internal_count_bits(a))           or_return

-	internal_sub(tmp, tmp, a)                                        or_return

-

-	return tmp.digit[0], nil

-}

-

-/*

-	Determines the setup value.

-	Assumes `mu` and `P` are not `nil`.

-

-	d := (1 << a.bits) - a;

-*/

-_private_int_reduce_2k_setup_l :: proc(mu, P: ^Int, allocator := context.allocator) -> (err: Error) {

-	context.allocator = allocator

-

-	tmp := &Int{}

-	defer internal_destroy(tmp)

-	internal_zero(tmp)                                               or_return

-

-	internal_int_power_of_two(tmp, internal_count_bits(P))           or_return

-	internal_sub(mu, tmp, P)                                         or_return

-

-	return nil

-}

-

-/*

-	Pre-calculate the value required for Barrett reduction.

-	For a given modulus "P" it calulates the value required in "mu"

-	Assumes `mu` and `P` are not `nil`.

-*/

-_private_int_reduce_setup :: proc(mu, P: ^Int, allocator := context.allocator) -> (err: Error) {

-	context.allocator = allocator

-

-	internal_int_power_of_two(mu, P.used * 2 * _DIGIT_BITS)           or_return

-	return internal_int_div(mu, mu, P)

-}

-

-/*

-	Determines the setup value.

-	Assumes `a` to not be `nil` and to have been initialized.

-*/

-_private_int_dr_setup :: proc(a: ^Int) -> (d: DIGIT) {

-	/*

-		The casts are required if _DIGIT_BITS is one less than

-		the number of bits in a DIGIT [e.g. _DIGIT_BITS==31].

-	*/

-	return DIGIT((1 << _DIGIT_BITS) - a.digit[0])

-}

-

-/*

-	Determines if a number is a valid DR modulus.

-	Assumes `a` to not be `nil` and to have been initialized.

-*/

-_private_dr_is_modulus :: proc(a: ^Int) -> (res: bool) {

-	/*

-		Must be at least two digits.

-	*/

-	if a.used < 2 { return false }

-

-	/*

-		Must be of the form b**k - a [a <= b] so all but the first digit must be equal to -1 (mod b).

-	*/

-	for ix := 1; ix < a.used; ix += 1 {

-		if a.digit[ix] != _MASK {

-			return false

-		}

-	}

-	return true

-}

-

-/*

-	Reduce "x" in place modulo "n" using the Diminished Radix algorithm.

-	Based on algorithm from the paper

-

-		"Generating Efficient Primes for Discrete Log Cryptosystems"

-					Chae Hoon Lim, Pil Joong Lee,

-			POSTECH Information Research Laboratories

-

-	The modulus must be of a special format [see manual].

-	Has been modified to use algorithm 7.10 from the LTM book instead

-

-	Input x must be in the range 0 <= x <= (n-1)**2

-	Assumes `x` and `n` to not be `nil` and to have been initialized.

-*/

-_private_int_dr_reduce :: proc(x, n: ^Int, k: DIGIT, allocator := context.allocator) -> (err: Error) {

-	/*

-		m = digits in modulus.

-	*/

-	m := n.used

-

-	/*

-		Ensure that "x" has at least 2m digits.

-	*/

-	internal_grow(x, m + m)                                          or_return

-

-	/*

-		Top of loop, this is where the code resumes if another reduction pass is required.

-	*/

-	for {

-		i: int

-		mu := DIGIT(0)

-

-		/*

-			Compute (x mod B**m) + k * [x/B**m] inline and inplace.

-		*/

-		for i = 0; i < m; i += 1 {

-			r         := _WORD(x.digit[i + m]) * _WORD(k) + _WORD(x.digit[i] + mu)

-			x.digit[i] = DIGIT(r & _WORD(_MASK))

-			mu         = DIGIT(r >> _WORD(_DIGIT_BITS))

-		}

-

-		/*

-			Set final carry.

-		*/

-		x.digit[i] = mu

-

-		/*

-			Zero words above m.

-		*/

-		mem.zero_slice(x.digit[m + 1:][:x.used - m])

-

-		/*

-			Clamp, sub and return.

-		*/

-		internal_clamp(x)                                            or_return

-

-		/*

-			If x >= n then subtract and reduce again.

-			Each successive "recursion" makes the input smaller and smaller.

-		*/

-		if internal_lt_abs(x, n) { break }

-

-		internal_sub(x, x, n)                                        or_return

-	}

-	return nil

-}

-

-/*

-	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.

-*/

-_private_int_exponent_mod :: proc(res, G, X, P: ^Int, redmode: int, allocator := context.allocator) -> (err: Error) {

-	context.allocator = allocator

-

-	M := [_TAB_SIZE]Int{}

-	winsize: uint

-

-	/*

-		Use a pointer to the reduction algorithm.

-		This allows us to use one of many reduction algorithms without modding the guts of the code with if statements everywhere.

-	*/

-	redux: #type proc(x, m, mu: ^Int, allocator := context.allocator) -> (err: Error)

-

-	defer {

-		internal_destroy(&M[1])

-		for x := 1 << (winsize - 1); x < (1 << winsize); x += 1 {

-			internal_destroy(&M[x])

-		}

-	}

-

-	/*

-		Find window size.

-	*/

-	x := internal_count_bits(X)

-	switch {

-	case x <= 7:

-		winsize = 2

-	case x <= 36:

-		winsize = 3

-	case x <= 140:

-		winsize = 4

-	case x <= 450:

-		winsize = 5

-	case x <= 1303:

-		winsize = 6

-	case x <= 3529:

-		winsize = 7

-	case:

-		winsize = 8

-	}

-

-	winsize = min(_MAX_WIN_SIZE, winsize) if _MAX_WIN_SIZE > 0 else winsize

-

-	/*

-		Init M array.

-		Init first cell.

-	*/

-	internal_zero(&M[1])                                             or_return

-

-	/*

-		Now init the second half of the array.

-	*/

-	for x = 1 << (winsize - 1); x < (1 << winsize); x += 1 {

-		internal_zero(&M[x])                                         or_return

-	}

-

-	/*

-		Create `mu`, used for Barrett reduction.

-	*/

-	mu := &Int{}

-	defer internal_destroy(mu)

-	internal_zero(mu)                                                or_return

-

-	if redmode == 0 {

-		_private_int_reduce_setup(mu, P)                             or_return

-		redux = _private_int_reduce

-	} else {

-		_private_int_reduce_2k_setup_l(mu, P)                        or_return

-		redux = _private_int_reduce_2k_l

-	}

-

-	/*

-		Create M table.

-

-		The M table contains powers of the base, e.g. M[x] = G**x mod P.

-		The first half of the table is not computed, though, except for M[0] and M[1].

-	*/

-	internal_int_mod(&M[1], G, P)                                    or_return

-

-	/*

-		Compute the value at M[1<<(winsize-1)] by squaring M[1] (winsize-1) times.

-

-		TODO: This can probably be replaced by computing the power and using `pow` to raise to it

-		instead of repeated squaring.

-	*/

-	slot := 1 << (winsize - 1)

-	internal_copy(&M[slot], &M[1])                                   or_return

-

-	for x = 0; x < int(winsize - 1); x += 1 {

-		/*

-			Square it.

-		*/

-		internal_sqr(&M[slot], &M[slot])                             or_return

-

-		/*

-			Reduce modulo P

-		*/

-		redux(&M[slot], P, mu)                                       or_return

-	}

-

-	/*

-		Create upper table, that is M[x] = M[x-1] * M[1] (mod P)

-		for x = (2**(winsize - 1) + 1) to (2**winsize - 1)

-	*/

-	for x = slot + 1; x < (1 << winsize); x += 1 {

-		internal_mul(&M[x], &M[x - 1], &M[1])                        or_return

-		redux(&M[x], P, mu)                                          or_return

-	}

-

-	/*

-		Setup result.

-	*/

-	internal_one(res)                                                or_return

-

-	/*

-		Set initial mode and bit cnt.

-	*/

-	mode   := 0

-	bitcnt := 1

-	buf    := DIGIT(0)

-	digidx := X.used - 1

-	bitcpy := uint(0)

-	bitbuf := DIGIT(0)

-

-	for {

-		/*

-			Grab next digit as required.

-		*/

-		bitcnt -= 1

-		if bitcnt == 0 {

-			/*

-				If digidx == -1 we are out of digits.

-			*/

-			if digidx == -1 { break }

-

-			/*

-				Read next digit and reset the bitcnt.

-			*/

-			buf    = X.digit[digidx]

-			digidx -= 1

-			bitcnt = _DIGIT_BITS

-		}

-

-		/*

-			Grab the next msb from the exponent.

-		*/

-		y := buf >> (_DIGIT_BITS - 1) & 1

-		buf <<= 1

-

-		/*

-			If the bit is zero and mode == 0 then we ignore it.

-			These represent the leading zero bits before the first 1 bit

-			in the exponent.  Technically this opt is not required but it

-			does lower the # of trivial squaring/reductions used.

-		*/

-		if mode == 0 && y == 0 {

-			continue

-		}

-

-		/*

-			If the bit is zero and mode == 1 then we square.

-		*/

-		if mode == 1 && y == 0 {

-			internal_sqr(res, res)                                   or_return

-			redux(res, P, mu)                                        or_return

-			continue

-		}

-

-		/*

-			Else we add it to the window.

-		*/

-		bitcpy += 1

-		bitbuf |= (y << (winsize - bitcpy))

-		mode    = 2

-

-		if (bitcpy == winsize) {

-			/*

-				Window is filled so square as required and multiply.

-				Square first.

-			*/

-			for x = 0; x < int(winsize); x += 1 {

-				internal_sqr(res, res)                               or_return

-				redux(res, P, mu)                                    or_return

-			}

-

-			/*

-				Then multiply.

-			*/

-			internal_mul(res, res, &M[bitbuf])                       or_return

-			redux(res, P, mu)                                        or_return

-

-			/*

-				Empty window and reset.

-			*/

-			bitcpy = 0

-			bitbuf = 0

-			mode   = 1

-		}

-	}

-

-	/*

-		If bits remain then square/multiply.

-	*/

-	if mode == 2 && bitcpy > 0 {

-		/*

-			Square then multiply if the bit is set.

-		*/

-		for x = 0; x < int(bitcpy); x += 1 {

-			internal_sqr(res, res)                                   or_return

-			redux(res, P, mu)                                        or_return

-

-			bitbuf <<= 1

-			if ((bitbuf & (1 << winsize)) != 0) {

-				/*

-					Then multiply.

-				*/

-				internal_mul(res, res, &M[1])                        or_return

-				redux(res, P, mu)                                    or_return

-			}

-		}

-	}

-	return err

-}

-

-/*

-	Computes Y == G**X mod P, HAC pp.616, Algorithm 14.85

-

-	Uses a left-to-right `k`-ary sliding window to compute the modular exponentiation.

-	The value of `k` changes based on the size of the exponent.

-

-	Uses Montgomery or Diminished Radix reduction [whichever appropriate]

-

-	Assumes `res`, `G`, `X` and `P` to not be `nil` and for `G`, `X` and `P` to have been initialized.

-*/

-_private_int_exponent_mod_fast :: proc(res, G, X, P: ^Int, redmode: int, allocator := context.allocator) -> (err: Error) {

-	context.allocator = allocator

-

-	M := [_TAB_SIZE]Int{}

-	winsize: uint

-

-	/*

-		Use a pointer to the reduction algorithm.

-		This allows us to use one of many reduction algorithms without modding the guts of the code with if statements everywhere.

-	*/

-	redux: #type proc(x, n: ^Int, rho: DIGIT, allocator := context.allocator) -> (err: Error)

-

-	defer {

-		internal_destroy(&M[1])

-		for x := 1 << (winsize - 1); x < (1 << winsize); x += 1 {

-			internal_destroy(&M[x])

-		}

-	}

-

-	/*

-		Find window size.

-	*/

-	x := internal_count_bits(X)

-	switch {

-	case x <= 7:

-		winsize = 2

-	case x <= 36:

-		winsize = 3

-	case x <= 140:

-		winsize = 4

-	case x <= 450:

-		winsize = 5

-	case x <= 1303:

-		winsize = 6

-	case x <= 3529:

-		winsize = 7

-	case:

-		winsize = 8

-	}

-

-	winsize = min(_MAX_WIN_SIZE, winsize) if _MAX_WIN_SIZE > 0 else winsize

-

-	/*

-		Init M array

-		Init first cell.

-	*/

-	cap := internal_int_allocated_cap(P)

-	internal_grow(&M[1], cap)                                        or_return

-

-	/*

-		Now init the second half of the array.

-	*/

-	for x = 1 << (winsize - 1); x < (1 << winsize); x += 1 {

-		internal_grow(&M[x], cap)                                    or_return

-	}

-

-	/*

-		Determine and setup reduction code.

-	*/

-	rho: DIGIT

-

-	if redmode == 0 {

-		/*

-			Now setup Montgomery.

-		*/

-		rho = _private_int_montgomery_setup(P)                       or_return

-

-		/*

-			Automatically pick the comba one if available (saves quite a few calls/ifs).

-		*/

-		if ((P.used * 2) + 1) < _WARRAY && P.used < _MAX_COMBA {

-			redux = _private_montgomery_reduce_comba

-		} else {

-			/*

-				Use slower baseline Montgomery method.

-			*/

-			redux = _private_int_montgomery_reduce

-		}

-	} else if redmode == 1 {

-		/*

-			Setup DR reduction for moduli of the form B**k - b.

-		*/

-		rho = _private_int_dr_setup(P)

-		redux = _private_int_dr_reduce

-	} else {

-		/*

-			Setup DR reduction for moduli of the form 2**k - b.

-		*/

-		rho = _private_int_reduce_2k_setup(P)                        or_return

-		redux = _private_int_reduce_2k

-	}

-

-	/*

-		Setup result.

-	*/

-	internal_grow(res, cap)                                          or_return

-

-	/*

-		Create M table

-		The first half of the table is not computed, though, except for M[0] and M[1]

-	*/

-

-	if redmode == 0 {

-		/*

-			Now we need R mod m.

-		*/

-		_private_int_montgomery_calc_normalization(res, P)           or_return

-

-		/*

-			Now set M[1] to G * R mod m.

-		*/

-		internal_mulmod(&M[1], G, res, P)                            or_return

-	} else {

-		internal_one(res)                                            or_return

-		internal_mod(&M[1], G, P)                                    or_return

-	}

-

-	/*

-		Compute the value at M[1<<(winsize-1)] by squaring M[1] (winsize-1) times.

-	*/

-	slot := 1 << (winsize - 1)

-	internal_copy(&M[slot], &M[1])                                   or_return

-

-	for x = 0; x < int(winsize - 1); x += 1 {

-		internal_sqr(&M[slot], &M[slot])                             or_return

-		redux(&M[slot], P, rho)                                      or_return

-	}

-

-	/*

-		Create upper table.

-	*/

-	for x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x += 1 {

-		internal_mul(&M[x], &M[x - 1], &M[1])                        or_return

-		redux(&M[x], P, rho)                                         or_return

-	}

-

-	/*

-		Set initial mode and bit cnt.

-	*/

-	mode   := 0

-	bitcnt := 1

-	buf    := DIGIT(0)

-	digidx := X.used - 1

-	bitcpy := 0

-	bitbuf := DIGIT(0)

-

-	for {

-		/*

-			Grab next digit as required.

-		*/

-		bitcnt -= 1

-		if bitcnt == 0 {

-			/*

-				If digidx == -1 we are out of digits so break.

-			*/

-			if digidx == -1 { break }

-

-			/*

-				Read next digit and reset the bitcnt.

-			*/

-			buf    = X.digit[digidx]

-			digidx -= 1

-			bitcnt = _DIGIT_BITS

-		}

-

-		/*

-			Grab the next msb from the exponent.

-		*/

-		y := (buf >> (_DIGIT_BITS - 1)) & 1

-		buf <<= 1

-

-		/*

-			If the bit is zero and mode == 0 then we ignore it.

-			These represent the leading zero bits before the first 1 bit in the exponent.

-			Technically this opt is not required but it does lower the # of trivial squaring/reductions used.

-		*/

-		if mode == 0 && y == 0 { continue }

-

-		/*

-			If the bit is zero and mode == 1 then we square.

-		*/

-		if mode == 1 && y == 0 {

-			internal_sqr(res, res)                                   or_return

-			redux(res, P, rho)                                       or_return

-			continue

-		}

-

-		/*

-			Else we add it to the window.

-		*/

-		bitcpy += 1

-		bitbuf |= (y << (winsize - uint(bitcpy)))

-		mode    = 2

-

-		if bitcpy == int(winsize) {

-			/*

-				Window is filled so square as required and multiply

-				Square first.

-			*/

-			for x = 0; x < int(winsize); x += 1 {

-				internal_sqr(res, res)                               or_return

-				redux(res, P, rho)                                   or_return

-			}

-

-			/*

-				Then multiply.

-			*/

-			internal_mul(res, res, &M[bitbuf])                       or_return

-			redux(res, P, rho)                                       or_return

-

-			/*

-				Empty window and reset.

-			*/

-			bitcpy = 0

-			bitbuf = 0

-			mode   = 1

-		}

-	}

-

-	/*

-		If bits remain then square/multiply.

-	*/

-	if mode == 2 && bitcpy > 0 {

-		/*

-			Square then multiply if the bit is set.

-		*/

-		for x = 0; x < bitcpy; x += 1 {

-			internal_sqr(res, res)                                   or_return

-			redux(res, P, rho)                                       or_return

-

-			/*

-				Get next bit of the window.

-			*/

-			bitbuf <<= 1

-			if bitbuf & (1 << winsize) != 0 {

-				/*

-					Then multiply.

-				*/

-				internal_mul(res, res, &M[1])                        or_return

-				redux(res, P, rho)                                   or_return

-			}

-		}

-	}

-

-	if redmode == 0 {

-		/*

-			Fixup result if Montgomery reduction is used.

-			Recall that any value in a Montgomery system is actually multiplied by R mod n.

-			So we have to reduce one more time to cancel out the factor of R.

-		*/

-		redux(res, P, rho)                                           or_return

-	}

-

-	return nil

-}

-

-/*

-	hac 14.61, pp608

-*/

-_private_inverse_modulo :: proc(dest, a, b: ^Int, allocator := context.allocator) -> (err: Error) {

-	context.allocator = allocator

-	x, y, u, v, A, B, C, D := &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}

-	defer internal_destroy(x, y, u, v, A, B, C, D)

-

-	/*

-		`b` cannot be negative.

-	*/

-	if b.sign == .Negative || internal_is_zero(b) {

-		return .Invalid_Argument

-	}

-

-	/*

-		init temps.

-	*/

-	internal_init_multi(x, y, u, v, A, B, C, D) or_return

-

-	/*

-		`x` = `a` % `b`, `y` = `b`

-	*/

-	internal_mod(x, a, b) or_return

-	internal_copy(y, b) or_return

-

-	/*

-		2. [modified] if x,y are both even then return an error!

-	*/

-	if internal_is_even(x) && internal_is_even(y) {

-		return .Invalid_Argument

-	}

-

-	/*

-		3. u=x, v=y, A=1, B=0, C=0, D=1

-	*/

-	internal_copy(u, x) or_return

-	internal_copy(v, y) or_return

-	internal_one(A) or_return

-	internal_one(D) or_return

-

-	for {

-		/*

-			4.  while `u` is even do:

-		*/

-		for internal_is_even(u) {

-			/*

-				4.1 `u` = `u` / 2

-			*/

-			internal_int_shr1(u, u) or_return

-

-			/*

-				4.2 if `A` or `B` is odd then:

-			*/

-			if internal_is_odd(A) || internal_is_odd(B) {

-				/*

-					`A` = (`A`+`y`) / 2, `B` = (`B`-`x`) / 2

-				*/

-				internal_add(A, A, y) or_return

-				internal_add(B, B, x) or_return

-			}

-			/*

-				`A` = `A` / 2, `B` = `B` / 2

-			*/

-			internal_int_shr1(A, A) or_return

-			internal_int_shr1(B, B) or_return

-		}

-

-		/*

-			5.  while `v` is even do:

-		*/

-		for internal_is_even(v) {

-			/*

-				5.1 `v` = `v` / 2

-			*/

-			internal_int_shr1(v, v) or_return

-

-			/*

-				5.2 if `C` or `D` is odd then:

-			*/

-			if internal_is_odd(C) || internal_is_odd(D) {

-				/*

-					`C` = (`C`+`y`) / 2, `D` = (`D`-`x`) / 2

-				*/

-				internal_add(C, C, y) or_return

-				internal_add(D, D, x) or_return

-			}

-			/*

-				`C` = `C` / 2, `D` = `D` / 2

-			*/

-			internal_int_shr1(C, C) or_return

-			internal_int_shr1(D, D) or_return

-		}

-

-		/*

-			6.  if `u` >= `v` then:

-		*/

-		if internal_cmp(u, v) != -1 {

-			/*

-				`u` = `u` - `v`, `A` = `A` - `C`, `B` = `B` - `D`

-			*/

-			internal_sub(u, u, v) or_return

-			internal_sub(A, A, C) or_return

-			internal_sub(B, B, D) or_return

-		} else {

-			/* v - v - u, C = C - A, D = D - B */

-			internal_sub(v, v, u) or_return

-			internal_sub(C, C, A) or_return

-			internal_sub(D, D, B) or_return

-		}

-

-		/*

-			If not zero goto step 4

-		*/

-		if internal_is_zero(u) {

-			break

-		}

-	}

-

-	/*

-		Now `a` = `C`, `b` = `D`, `gcd` == `g`*`v`

-	*/

-

-	/*

-		If `v` != `1` then there is no inverse.

-	*/

-	if !internal_eq(v, 1) {

-		return .Invalid_Argument

-	}

-

-	/*

-		If its too low.

-	*/

-	if internal_is_negative(C) {

-		internal_add(C, C, b) or_return

-	}

-

-	/*

-		Too big.

-	*/

-	if internal_gte(C, 0) {

-		internal_sub(C, C, b) or_return

-	}

-

-	/*

-		`C` is now the inverse.

-	*/

-	swap(dest, C)

-

-	return

-}

-

-/*

-	Computes the modular inverse via binary extended Euclidean algorithm, that is `dest` = 1 / `a` mod `b`.

-

-	Based on slow invmod except this is optimized for the case where `b` is odd,

-	as per HAC Note 14.64 on pp. 610.

-*/

-_private_inverse_modulo_odd :: proc(dest, a, b: ^Int, allocator := context.allocator) -> (err: Error) {

-	context.allocator = allocator

-	x, y, u, v, B, D := &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}

-	defer internal_destroy(x, y, u, v, B, D)

-

-	sign: Sign

-

-	/*

-		2. [modified] `b` must be odd.

-	*/

-	if internal_is_even(b) { return .Invalid_Argument }

-

-	/*

-		Init all our temps.

-	*/

-	internal_init_multi(x, y, u, v, B, D) or_return

-

-	/*

-		`x` == modulus, `y` == value to invert.

-	*/

-	internal_copy(x, b) or_return

-

-	/*

-		We need `y` = `|a|`.

-	*/

-	internal_mod(y, a, b) or_return

-

-	/*

-		If one of `x`, `y` is zero return an error!

-	*/

-	if internal_is_zero(x) || internal_is_zero(y) { return .Invalid_Argument }

-

-	/*

-		3. `u` = `x`, `v` = `y`, `A` = 1, `B` = 0, `C` = 0, `D` = 1

-	*/

-	internal_copy(u, x) or_return

-	internal_copy(v, y) or_return

-

-	internal_one(D) or_return

-

-	for {

-		/*

-			4.  while `u` is even do.

-		*/

-		for internal_is_even(u) {

-			/*

-				4.1 `u` = `u` / 2

-			*/

-			internal_int_shr1(u, u) or_return

-

-			/*

-				4.2 if `B` is odd then:

-			*/

-			if internal_is_odd(B) {

-				/*

-					`B` = (`B` - `x`) / 2

-				*/

-				internal_sub(B, B, x) or_return

-			}

-

-			/*

-				`B` = `B` / 2

-			*/

-			internal_int_shr1(B, B) or_return

-		}

-

-		/*

-			5.  while `v` is even do:

-		*/

-		for internal_is_even(v) {

-			/*

-				5.1 `v` = `v` / 2

-			*/

-			internal_int_shr1(v, v) or_return

-

-			/*

-				5.2 if `D` is odd then:

-			*/

-			if internal_is_odd(D) {

-				/*

-					`D` = (`D` - `x`) / 2

-				*/

-				internal_sub(D, D, x) or_return

-			}

-			/*

-				`D` = `D` / 2

-			*/

-			internal_int_shr1(D, D) or_return

-		}

-

-		/*

-			6.  if `u` >= `v` then:

-		*/

-		if internal_cmp(u, v) != -1 {

-			/*

-				`u` = `u` - `v`, `B` = `B` - `D`

-			*/

-			internal_sub(u, u, v) or_return

-			internal_sub(B, B, D) or_return

-		} else {

-			/*

-				`v` - `v` - `u`, `D` = `D` - `B`

-			*/

-			internal_sub(v, v, u) or_return

-			internal_sub(D, D, B) or_return

-		}

-

-		/*

-			If not zero goto step 4.

-		*/

-		if internal_is_zero(u) { break }

-	}

-

-	/*

-		Now `a` = C, `b` = D, gcd == g*v

-	*/

-

-	/*

-		if `v` != 1 then there is no inverse

-	*/

-	if internal_cmp(v, 1) != 0 {

-		return .Invalid_Argument

-	}

-

-	/*

-		`b` is now the inverse.

-	*/

-	sign = a.sign

-	for internal_int_is_negative(D) {

-		internal_add(D, D, b) or_return

-	}

-

-	/*

-		Too big.

-	*/

-	for internal_gte_abs(D, b) {

-		internal_sub(D, D, b) or_return

-	}

-

-	swap(dest, D)

-	dest.sign = sign

-	return nil

-}

-

-

-/*

-	Returns the log2 of an `Int`.

-	Assumes `a` not to be `nil` and to have been initialized.

-	Also assumes `base` is a power of two.

-*/

-_private_log_power_of_two :: proc(a: ^Int, base: DIGIT) -> (log: int, err: Error) {

-	base := base

-	y: int

-	for y = 0; base & 1 == 0; {

-		y += 1

-		base >>= 1

-	}

-	log = internal_count_bits(a)

-	return (log - 1) / y, err

-}

-

-/*

-	Copies DIGITs from `src` to `dest`.

-	Assumes `src` and `dest` to not be `nil` and have been initialized.

-*/

-_private_copy_digits :: proc(dest, src: ^Int, digits: int, offset := int(0)) -> (err: Error) {

-	digits := digits

-	/*

-		If dest == src, do nothing

-	*/

-	if dest == src {

-		return nil

-	}

-

-	digits = min(digits, len(src.digit), len(dest.digit))

-	mem.copy_non_overlapping(&dest.digit[0], &src.digit[offset], size_of(DIGIT) * digits)

-	return nil

-}

-

-

-/*

-	Shift left by `digits` * _DIGIT_BITS bits.

-*/

-_private_int_shl_leg :: proc(quotient: ^Int, digits: int, allocator := context.allocator) -> (err: Error) {

-	context.allocator = allocator

-

-	if digits <= 0 { return nil }

-

-	/*

-		No need to shift a zero.

-	*/

-	if #force_inline internal_is_zero(quotient) {

-		return nil

-	}

-

-	/*

-		Resize `quotient` to accomodate extra digits.

-	*/

-	#force_inline internal_grow(quotient, quotient.used + digits) or_return

-

-	/*

-		Increment the used by the shift amount then copy upwards.

-	*/

-

-	/*

-		Much like `_private_int_shr_leg`, this is implemented using a sliding window,

-		except the window goes the other way around.

-	*/

-	#no_bounds_check for x := quotient.used; x > 0; x -= 1 {

-		quotient.digit[x+digits-1] = quotient.digit[x-1]

-	}

-

-	quotient.used += digits

-	mem.zero_slice(quotient.digit[:digits])

-	return nil

-}

-

-/*

-	Shift right by `digits` * _DIGIT_BITS bits.

-*/

-_private_int_shr_leg :: proc(quotient: ^Int, digits: int, allocator := context.allocator) -> (err: Error) {

-	context.allocator = allocator

-

-	if digits <= 0 { return nil }

-

-	/*

-		If digits > used simply zero and return.

-	*/

-	if digits > quotient.used { return internal_zero(quotient) }

-

-	/*

-		Much like `int_shl_digit`, this is implemented using a sliding window,

-		except the window goes the other way around.

-

-		b-2 | b-1 | b0 | b1 | b2 | ... | bb |   ---->

-					/\                   |      ---->

-					 \-------------------/      ---->

-	*/

-

-	#no_bounds_check for x := 0; x < (quotient.used - digits); x += 1 {

-		quotient.digit[x] = quotient.digit[x + digits]

-	}

-	quotient.used -= digits

-	internal_zero_unused(quotient)

-	return internal_clamp(quotient)

-}

-

-/*	

-	========================    End of private procedures    =======================

-

-	===============================  Private tables  ===============================

-

-	Tables used by `internal_*` and `_*`.

-*/

-

-_private_int_rem_128 := [?]DIGIT{

-	0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,

-	0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,

-	1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,

-	1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,

-	0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,

-	1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,

-	1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,

-	1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,

-}

-#assert(128 * size_of(DIGIT) == size_of(_private_int_rem_128))

-

-_private_int_rem_105 := [?]DIGIT{

-	0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1,

-	0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1,

-	0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1,

-	1, 0, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1,

-	0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1,

-	1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 0, 1, 1, 1, 1,

-	1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1,

-}

-#assert(105 * size_of(DIGIT) == size_of(_private_int_rem_105))

-

-_PRIME_TAB_SIZE :: 256

-_private_prime_table := [_PRIME_TAB_SIZE]DIGIT{

-	0x0002, 0x0003, 0x0005, 0x0007, 0x000B, 0x000D, 0x0011, 0x0013,

-	0x0017, 0x001D, 0x001F, 0x0025, 0x0029, 0x002B, 0x002F, 0x0035,

-	0x003B, 0x003D, 0x0043, 0x0047, 0x0049, 0x004F, 0x0053, 0x0059,

-	0x0061, 0x0065, 0x0067, 0x006B, 0x006D, 0x0071, 0x007F, 0x0083,

-	0x0089, 0x008B, 0x0095, 0x0097, 0x009D, 0x00A3, 0x00A7, 0x00AD,

-	0x00B3, 0x00B5, 0x00BF, 0x00C1, 0x00C5, 0x00C7, 0x00D3, 0x00DF,

-	0x00E3, 0x00E5, 0x00E9, 0x00EF, 0x00F1, 0x00FB, 0x0101, 0x0107,

-	0x010D, 0x010F, 0x0115, 0x0119, 0x011B, 0x0125, 0x0133, 0x0137,

-

-	0x0139, 0x013D, 0x014B, 0x0151, 0x015B, 0x015D, 0x0161, 0x0167,

-	0x016F, 0x0175, 0x017B, 0x017F, 0x0185, 0x018D, 0x0191, 0x0199,

-	0x01A3, 0x01A5, 0x01AF, 0x01B1, 0x01B7, 0x01BB, 0x01C1, 0x01C9,

-	0x01CD, 0x01CF, 0x01D3, 0x01DF, 0x01E7, 0x01EB, 0x01F3, 0x01F7,

-	0x01FD, 0x0209, 0x020B, 0x021D, 0x0223, 0x022D, 0x0233, 0x0239,

-	0x023B, 0x0241, 0x024B, 0x0251, 0x0257, 0x0259, 0x025F, 0x0265,

-	0x0269, 0x026B, 0x0277, 0x0281, 0x0283, 0x0287, 0x028D, 0x0293,

-	0x0295, 0x02A1, 0x02A5, 0x02AB, 0x02B3, 0x02BD, 0x02C5, 0x02CF,

-

-	0x02D7, 0x02DD, 0x02E3, 0x02E7, 0x02EF, 0x02F5, 0x02F9, 0x0301,

-	0x0305, 0x0313, 0x031D, 0x0329, 0x032B, 0x0335, 0x0337, 0x033B,

-	0x033D, 0x0347, 0x0355, 0x0359, 0x035B, 0x035F, 0x036D, 0x0371,

-	0x0373, 0x0377, 0x038B, 0x038F, 0x0397, 0x03A1, 0x03A9, 0x03AD,

-	0x03B3, 0x03B9, 0x03C7, 0x03CB, 0x03D1, 0x03D7, 0x03DF, 0x03E5,

-	0x03F1, 0x03F5, 0x03FB, 0x03FD, 0x0407, 0x0409, 0x040F, 0x0419,

-	0x041B, 0x0425, 0x0427, 0x042D, 0x043F, 0x0443, 0x0445, 0x0449,

-	0x044F, 0x0455, 0x045D, 0x0463, 0x0469, 0x047F, 0x0481, 0x048B,

-

-	0x0493, 0x049D, 0x04A3, 0x04A9, 0x04B1, 0x04BD, 0x04C1, 0x04C7,

-	0x04CD, 0x04CF, 0x04D5, 0x04E1, 0x04EB, 0x04FD, 0x04FF, 0x0503,

-	0x0509, 0x050B, 0x0511, 0x0515, 0x0517, 0x051B, 0x0527, 0x0529,

-	0x052F, 0x0551, 0x0557, 0x055D, 0x0565, 0x0577, 0x0581, 0x058F,

-	0x0593, 0x0595, 0x0599, 0x059F, 0x05A7, 0x05AB, 0x05AD, 0x05B3,

-	0x05BF, 0x05C9, 0x05CB, 0x05CF, 0x05D1, 0x05D5, 0x05DB, 0x05E7,

-	0x05F3, 0x05FB, 0x0607, 0x060D, 0x0611, 0x0617, 0x061F, 0x0623,

-	0x062B, 0x062F, 0x063D, 0x0641, 0x0647, 0x0649, 0x064D, 0x0653,

-}

-#assert(_PRIME_TAB_SIZE * size_of(DIGIT) == size_of(_private_prime_table))

-

-when MATH_BIG_FORCE_64_BIT || (!MATH_BIG_FORCE_32_BIT && size_of(rawptr) == 8) {

-	_factorial_table := [35]_WORD{

-/* f(00): */                                                     1,

-/* f(01): */                                                     1,

-/* f(02): */                                                     2,

-/* f(03): */                                                     6,

-/* f(04): */                                                    24,

-/* f(05): */                                                   120,

-/* f(06): */                                                   720,

-/* f(07): */                                                 5_040,

-/* f(08): */                                                40_320,

-/* f(09): */                                               362_880,

-/* f(10): */                                             3_628_800,

-/* f(11): */                                            39_916_800,

-/* f(12): */                                           479_001_600,

-/* f(13): */                                         6_227_020_800,

-/* f(14): */                                        87_178_291_200,

-/* f(15): */                                     1_307_674_368_000,

-/* f(16): */                                    20_922_789_888_000,

-/* f(17): */                                   355_687_428_096_000,

-/* f(18): */                                 6_402_373_705_728_000,

-/* f(19): */                               121_645_100_408_832_000,

-/* f(20): */                             2_432_902_008_176_640_000,

-/* f(21): */                            51_090_942_171_709_440_000,

-/* f(22): */                         1_124_000_727_777_607_680_000,

-/* f(23): */                        25_852_016_738_884_976_640_000,

-/* f(24): */                       620_448_401_733_239_439_360_000,

-/* f(25): */                    15_511_210_043_330_985_984_000_000,

-/* f(26): */                   403_291_461_126_605_635_584_000_000,

-/* f(27): */                10_888_869_450_418_352_160_768_000_000,

-/* f(28): */               304_888_344_611_713_860_501_504_000_000,

-/* f(29): */             8_841_761_993_739_701_954_543_616_000_000,

-/* f(30): */           265_252_859_812_191_058_636_308_480_000_000,

-/* f(31): */         8_222_838_654_177_922_817_725_562_880_000_000,

-/* f(32): */       263_130_836_933_693_530_167_218_012_160_000_000,

-/* f(33): */     8_683_317_618_811_886_495_518_194_401_280_000_000,

-/* f(34): */   295_232_799_039_604_140_847_618_609_643_520_000_000,

-	}

-} else {

-	_factorial_table := [21]_WORD{

-/* f(00): */                                                     1,

-/* f(01): */                                                     1,

-/* f(02): */                                                     2,

-/* f(03): */                                                     6,

-/* f(04): */                                                    24,

-/* f(05): */                                                   120,

-/* f(06): */                                                   720,

-/* f(07): */                                                 5_040,

-/* f(08): */                                                40_320,

-/* f(09): */                                               362_880,

-/* f(10): */                                             3_628_800,

-/* f(11): */                                            39_916_800,

-/* f(12): */                                           479_001_600,

-/* f(13): */                                         6_227_020_800,

-/* f(14): */                                        87_178_291_200,

-/* f(15): */                                     1_307_674_368_000,

-/* f(16): */                                    20_922_789_888_000,

-/* f(17): */                                   355_687_428_096_000,

-/* f(18): */                                 6_402_373_705_728_000,

-/* f(19): */                               121_645_100_408_832_000,

-/* f(20): */                             2_432_902_008_176_640_000,

-	}

-}

-

-/*

-	=========================  End of private tables  ========================

+/*
+	Copyright 2021 Jeroen van Rijn <nom@duclavier.com>.
+	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.
+
+	=============================    Private procedures    =============================
+
+	Private procedures used by the above low-level routines follow.
+
+	Don't call these yourself unless you really know what you're doing.
+	They include implementations that are optimimal for certain ranges of input only.
+
+	These aren't exported for the same reasons.
+*/
+
+
+package math_big
+
+import "base:intrinsics"
+import "core:mem"
+
+/*
+	Multiplies |a| * |b| and only computes upto digs digits of result.
+	HAC pp. 595, Algorithm 14.12  Modified so you can control how
+	many digits of output are created.
+*/
+_private_int_mul :: proc(dest, a, b: ^Int, digits: int, allocator := context.allocator) -> (err: Error) {
+	context.allocator = allocator
+
+	/*
+		Can we use the fast multiplier?
+	*/
+	if digits < _WARRAY && min(a.used, b.used) < _MAX_COMBA {
+		return #force_inline _private_int_mul_comba(dest, a, b, digits)
+	}
+
+	/*
+		Set up temporary output `Int`, which we'll swap for `dest` when done.
+	*/
+
+	t := &Int{}
+
+	internal_grow(t, max(digits, _DEFAULT_DIGIT_COUNT)) or_return
+	t.used = digits
+
+	/*
+		Compute the digits of the product directly.
+	*/
+	pa := a.used
+	for ix := 0; ix < pa; ix += 1 {
+		/*
+			Limit ourselves to `digits` DIGITs of output.
+		*/
+		pb    := min(b.used, digits - ix)
+		carry := _WORD(0)
+		iy    := 0
+
+		/*
+			Compute the column of the output and propagate the carry.
+		*/
+		#no_bounds_check for iy = 0; iy < pb; iy += 1 {
+			/*
+				Compute the column as a _WORD.
+			*/
+			column := _WORD(t.digit[ix + iy]) + _WORD(a.digit[ix]) * _WORD(b.digit[iy]) + carry
+
+			/*
+				The new column is the lower part of the result.
+			*/
+			t.digit[ix + iy] = DIGIT(column & _WORD(_MASK))
+
+			/*
+				Get the carry word from the result.
+			*/
+			carry = column >> _DIGIT_BITS
+		}
+		/*
+			Set carry if it is placed below digits
+		*/
+		if ix + iy < digits {
+			t.digit[ix + pb] = DIGIT(carry)
+		}
+	}
+
+	internal_swap(dest, t)
+	internal_destroy(t)
+	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 internal_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;
+	*/
+	_private_int_shl_leg(b1, 4 * B)             or_return
+	_private_int_shl_leg(S2, 3 * B)             or_return
+	internal_add(b1, b1, S2)                    or_return
+	_private_int_shl_leg(S1, 2 * B)             or_return
+	internal_add(b1, b1, S1)                    or_return
+	_private_int_shl_leg(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.
+
+	Let `B` represent the radix [e.g. 2**_DIGIT_BITS] and let `n` represent
+	half of the number of digits in the min(a,b)
+
+	`a` = `a1` * `B`**`n` + `a0`
+	`b` = `b`1 * `B`**`n` + `b0`
+
+	Then, a * b => 1b1 * B**2n + ((a1 + a0)(b1 + b0) - (a0b0 + a1b1)) * B + a0b0
+
+	Note that a1b1 and a0b0 are used twice and only need to be computed once.
+	So in total three half size (half # of digit) multiplications are performed,
+		a0b0, a1b1 and (a1+b1)(a0+b0)
+
+	Note that a multiplication of half the digits requires 1/4th the number of
+	single precision multiplications, so in total after one call 25% of the
+	single precision multiplications are saved.
+
+	Note also that the call to `internal_mul` can end up back in this function
+	if the a0, a1, b0, or b1 are above the threshold.
+
+	This is known as divide-and-conquer and leads to the famous O(N**lg(3)) or O(N**1.584)
+	work which is asymptopically lower than the standard O(N**2) that the
+	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_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{}
+	defer internal_destroy(x0, x1, y0, y1, t1, x0y0, x1y1)
+
+	/*
+		min # of digits, divided by two.
+	*/
+	B := min(a.used, b.used) >> 1
+
+	/*
+		Init all the temps.
+	*/
+	internal_grow(x0, B)          or_return
+	internal_grow(x1, a.used - B) or_return
+	internal_grow(y0, B)          or_return
+	internal_grow(y1, b.used - B) or_return
+	internal_grow(t1, B * 2)      or_return
+	internal_grow(x0y0, B * 2)    or_return
+	internal_grow(x1y1, B * 2)    or_return
+
+	/*
+		Now shift the digits.
+	*/
+	x0.used, y0.used = B, B
+	x1.used = a.used - B
+	y1.used = b.used - B
+
+	/*
+		We copy the digits directly instead of using higher level functions
+		since we also need to shift the digits.
+	*/
+	internal_copy_digits(x0, a, x0.used)
+	internal_copy_digits(y0, b, y0.used)
+	internal_copy_digits(x1, a, x1.used, B)
+	internal_copy_digits(y1, b, y1.used, B)
+
+	/*
+		Only need to clamp the lower words since by definition the
+		upper words x1/y1 must have a known number of digits.
+	*/
+	clamp(x0)
+	clamp(y0)
+
+	/*
+		Now calc the products x0y0 and x1y1,
+		after this x0 is no longer required, free temp [x0==t2]!
+	*/
+	internal_mul(x0y0, x0, y0)      or_return /* x0y0 = x0*y0 */
+	internal_mul(x1y1, x1, y1)      or_return /* x1y1 = x1*y1 */
+	internal_add(t1,   x1, x0)      or_return /* now calc x1+x0 and */
+	internal_add(x0,   y1, y0)      or_return /* t2 = y1 + y0 */
+	internal_mul(t1,   t1, x0)      or_return /* t1 = (x1 + x0) * (y1 + y0) */
+
+	/*
+		Add x0y0.
+	*/
+	internal_add(x0, x0y0, x1y1)    or_return /* t2 = x0y0 + x1y1 */
+	internal_sub(t1,   t1,   x0)    or_return /* t1 = (x1+x0)*(y1+y0) - (x1y1 + x0y0) */
+
+	/*
+		shift by B.
+	*/
+	_private_int_shl_leg(t1, B)       or_return /* t1 = (x0y0 + x1y1 - (x1-x0)*(y1-y0))<<B */
+	_private_int_shl_leg(x1y1, B * 2) or_return /* x1y1 = x1y1 << 2*B */
+
+	internal_add(t1, x0y0, t1)      or_return /* t1 = x0y0 + t1 */
+	internal_add(dest, t1, x1y1)    or_return /* t1 = x0y0 + t1 + x1y1 */
+
+	return nil
+}
+
+
+
+/*
+	Fast (comba) multiplier
+
+	This is the fast column-array [comba] multiplier.  It is
+	designed to compute the columns of the product first
+	then handle the carries afterwards.  This has the effect
+	of making the nested loops that compute the columns very
+	simple and schedulable on super-scalar processors.
+
+	This has been modified to produce a variable number of
+	digits of output so if say only a half-product is required
+	you don't have to compute the upper half (a feature
+	required for fast Barrett reduction).
+
+	Based on Algorithm 14.12 on pp.595 of HAC.
+*/
+_private_int_mul_comba :: proc(dest, a, b: ^Int, digits: int, allocator := context.allocator) -> (err: Error) {
+	context.allocator = allocator
+
+	/*
+		Set up array.
+	*/
+	W: [_WARRAY]DIGIT = ---
+
+	/*
+		Grow the destination as required.
+	*/
+	internal_grow(dest, digits) or_return
+
+	/*
+		Number of output digits to produce.
+	*/
+	pa := min(digits, a.used + b.used)
+
+	/*
+		Clear the carry
+	*/
+	_W := _WORD(0)
+
+	ix: int
+	for ix = 0; ix < pa; ix += 1 {
+		tx, ty, iy, iz: int
+
+		/*
+			Get offsets into the two bignums.
+		*/
+		ty = min(b.used - 1, ix)
+		tx = ix - ty
+
+		/*
+			This is the number of times the loop will iterate, essentially.
+			while (tx++ < a->used && ty-- >= 0) { ... }
+		*/
+		 
+		iy = min(a.used - tx, ty + 1)
+
+		/*
+			Execute loop.
+		*/
+		#no_bounds_check for iz = 0; iz < iy; iz += 1 {
+			_W += _WORD(a.digit[tx + iz]) * _WORD(b.digit[ty - iz])
+		}
+
+		/*
+			Store term.
+		*/
+		W[ix] = DIGIT(_W) & _MASK
+
+		/*
+			Make next carry.
+		*/
+		_W = _W >> _WORD(_DIGIT_BITS)
+	}
+
+	/*
+		Setup dest.
+	*/
+	old_used := dest.used
+	dest.used = pa
+
+	/*
+		Now extract the previous digit [below the carry].
+	*/
+	copy_slice(dest.digit[0:], W[:pa])	
+
+	/*
+		Clear unused digits [that existed in the old copy of dest].
+	*/
+	internal_zero_unused(dest, old_used)
+
+	/*
+		Adjust dest.used based on leading zeroes.
+	*/
+
+	return internal_clamp(dest)
+}
+
+/*
+	Multiplies |a| * |b| and does not compute the lower digs digits
+	[meant to get the higher part of the product]
+*/
+_private_int_mul_high :: proc(dest, a, b: ^Int, digits: int, allocator := context.allocator) -> (err: Error) {
+	context.allocator = allocator
+
+	/*
+		Can we use the fast multiplier?
+	*/
+	if a.used + b.used + 1 < _WARRAY && min(a.used, b.used) < _MAX_COMBA {
+		return _private_int_mul_high_comba(dest, a, b, digits)
+	}
+
+	internal_grow(dest, a.used + b.used + 1) or_return
+	dest.used = a.used + b.used + 1
+
+	pa := a.used
+	pb := b.used
+	for ix := 0; ix < pa; ix += 1 {
+		carry := DIGIT(0)
+
+		for iy := digits - ix; iy < pb; iy += 1 {
+			/*
+				Calculate the double precision result.
+			*/
+			r := _WORD(dest.digit[ix + iy]) + _WORD(a.digit[ix]) * _WORD(b.digit[iy]) + _WORD(carry)
+
+			/*
+				Get the lower part.
+			*/
+			dest.digit[ix + iy] = DIGIT(r & _WORD(_MASK))
+
+			/*
+				Carry the carry.
+			*/
+			carry = DIGIT(r >> _WORD(_DIGIT_BITS))
+		}
+		dest.digit[ix + pb] = carry
+	}
+	return internal_clamp(dest)
+}
+
+/*
+	This is a modified version of `_private_int_mul_comba` that only produces output digits *above* `digits`.
+	See the comments for `_private_int_mul_comba` to see how it works.
+
+	This is used in the Barrett reduction since for one of the multiplications
+	only the higher digits were needed.  This essentially halves the work.
+
+	Based on Algorithm 14.12 on pp.595 of HAC.
+*/
+_private_int_mul_high_comba :: proc(dest, a, b: ^Int, digits: int, allocator := context.allocator) -> (err: Error) {
+	context.allocator = allocator
+
+	W: [_WARRAY]DIGIT = ---
+	_W: _WORD = 0
+
+	/*
+		Number of output digits to produce. Grow the destination as required.
+	*/
+	pa := a.used + b.used
+	internal_grow(dest, pa) or_return
+
+	ix: int
+	for ix = digits; ix < pa; ix += 1 {
+		/*
+			Get offsets into the two bignums.
+		*/
+		ty := min(b.used - 1, ix)
+		tx := ix - ty
+
+		/*
+			This is the number of times the loop will iterrate, essentially it's
+			while (tx++ < a->used && ty-- >= 0) { ... }
+		*/
+		iy := min(a.used - tx, ty + 1)
+
+		/*
+			Execute loop.
+		*/
+		for iz := 0; iz < iy; iz += 1 {
+			_W += _WORD(a.digit[tx + iz]) * _WORD(b.digit[ty - iz])
+		}
+
+		/*
+			Store term.
+		*/
+		W[ix] = DIGIT(_W) & DIGIT(_MASK)
+
+		/*
+			Make next carry.
+		*/
+		_W = _W >> _WORD(_DIGIT_BITS)
+	}
+
+	/*
+		Setup dest
+	*/
+	old_used := dest.used
+	dest.used = pa
+
+	for ix = digits; ix < pa; ix += 1 {
+		/*
+			Now extract the previous digit [below the carry].
+		*/
+		dest.digit[ix] = W[ix]
+	}
+
+	/*
+		Zero remainder.
+	*/
+	internal_zero_unused(dest, old_used)
+
+	/*
+		Adjust dest.used based on leading zeroes.
+	*/
+	return internal_clamp(dest)
+}
+
+/*
+	Single-digit multiplication with the smaller number as the single-digit.
+*/
+_private_int_mul_balance :: proc(dest, a, b: ^Int, allocator := context.allocator) -> (err: Error) {
+	context.allocator = allocator
+	a, b := a, b
+
+	a0, tmp, r := &Int{}, &Int{}, &Int{}
+	defer internal_destroy(a0, tmp, r)
+
+	b_size   := min(a.used, b.used)
+	n_blocks := max(a.used, b.used) / b_size
+
+	internal_grow(a0, b_size + 2) or_return
+	internal_init_multi(tmp, r)   or_return
+
+	/*
+		Make sure that `a` is the larger one.
+	*/
+	if a.used < b.used {
+		a, b = b, a
+	}
+	assert(a.used >= b.used)
+
+	i, j := 0, 0
+	for ; i < n_blocks; i += 1 {
+		/*
+			Cut a slice off of `a`.
+		*/
+
+		a0.used = b_size
+		internal_copy_digits(a0, a, a0.used, j)
+		j += a0.used
+		internal_clamp(a0)
+
+		/*
+			Multiply with `b`.
+		*/
+		internal_mul(tmp, a0, b)                                     or_return
+
+		/*
+			Shift `tmp` to the correct position.
+		*/
+		_private_int_shl_leg(tmp, b_size * i)                          or_return
+
+		/*
+			Add to output. No carry needed.
+		*/
+		internal_add(r, r, tmp)                                      or_return
+	}
+
+	/*
+		The left-overs; there are always left-overs.
+	*/
+	if j < a.used {
+		a0.used = a.used - j
+		internal_copy_digits(a0, a, a0.used, j)
+		j += a0.used
+		internal_clamp(a0)
+
+		internal_mul(tmp, a0, b)                                     or_return
+		_private_int_shl_leg(tmp, b_size * i)                          or_return
+		internal_add(r, r, tmp)                                      or_return
+	}
+
+	internal_swap(dest, r)
+	return
+}
+
+/*
+	Low level squaring, b = a*a, HAC pp.596-597, Algorithm 14.16
+	Assumes `dest` and `src` to not be `nil`, and `src` to have been initialized.
+*/
+_private_int_sqr :: proc(dest, src: ^Int, allocator := context.allocator) -> (err: Error) {
+	context.allocator = allocator
+	pa := src.used
+
+	t := &Int{}; ix, iy: int
+	/*
+		Grow `t` to maximum needed size, or `_DEFAULT_DIGIT_COUNT`, whichever is bigger.
+	*/
+	internal_grow(t, max((2 * pa) + 1, _DEFAULT_DIGIT_COUNT)) or_return
+	t.used = (2 * pa) + 1
+
+	#no_bounds_check for ix = 0; ix < pa; ix += 1 {
+		carry := DIGIT(0)
+		/*
+			First calculate the digit at 2*ix; calculate double precision result.
+		*/
+		r := _WORD(t.digit[ix+ix]) + (_WORD(src.digit[ix]) * _WORD(src.digit[ix]))
+
+		/*
+			Store lower part in result.
+		*/
+		t.digit[ix+ix] = DIGIT(r & _WORD(_MASK))
+		/*
+			Get the carry.
+		*/
+		carry = DIGIT(r >> _DIGIT_BITS)
+
+		#no_bounds_check for iy = ix + 1; iy < pa; iy += 1 {
+			/*
+				First calculate the product.
+			*/
+			r = _WORD(src.digit[ix]) * _WORD(src.digit[iy])
+
+			/* Now calculate the double precision result. Nóte we use
+			 * addition instead of *2 since it's easier to optimize
+			 */
+			r = _WORD(t.digit[ix+iy]) + r + r + _WORD(carry)
+
+			/*
+				Store lower part.
+			*/
+			t.digit[ix+iy] = DIGIT(r & _WORD(_MASK))
+
+			/*
+				Get carry.
+			*/
+			carry = DIGIT(r >> _DIGIT_BITS)
+		}
+		/*
+			Propagate upwards.
+		*/
+		#no_bounds_check for carry != 0 {
+			r     = _WORD(t.digit[ix+iy]) + _WORD(carry)
+			t.digit[ix+iy] = DIGIT(r & _WORD(_MASK))
+			carry = DIGIT(r >> _WORD(_DIGIT_BITS))
+			iy += 1
+		}
+	}
+
+	err = internal_clamp(t)
+	internal_swap(dest, t)
+	internal_destroy(t)
+	return err
+}
+
+/*
+	The jist of squaring...
+	You do like mult except the offset of the tmpx [one that starts closer to zero] can't equal the offset of tmpy.
+	So basically you set up iy like before then you min it with (ty-tx) so that it never happens.
+	You double all those you add in the inner loop. After that loop you do the squares and add them in.
+
+	Assumes `dest` and `src` not to be `nil` and `src` to have been initialized.	
+*/
+_private_int_sqr_comba :: proc(dest, src: ^Int, allocator := context.allocator) -> (err: Error) {
+	context.allocator = allocator
+
+	W: [_WARRAY]DIGIT = ---
+
+	/*
+		Grow the destination as required.
+	*/
+	pa := uint(src.used) + uint(src.used)
+	internal_grow(dest, int(pa)) or_return
+
+	/*
+		Number of output digits to produce.
+	*/
+	W1 := _WORD(0)
+	_W  : _WORD = ---
+	ix := uint(0)
+
+	#no_bounds_check for ; ix < pa; ix += 1 {
+		/*
+			Clear counter.
+		*/
+		_W = {}
+
+		/*
+			Get offsets into the two bignums.
+		*/
+		ty := min(uint(src.used) - 1, ix)
+		tx := ix - ty
+
+		/*
+			This is the number of times the loop will iterate,
+			essentially while (tx++ < a->used && ty-- >= 0) { ... }
+		*/
+		iy := min(uint(src.used) - tx, ty + 1)
+
+		/*
+			Now for squaring, tx can never equal ty.
+			We halve the distance since they approach at a rate of 2x,
+			and we have to round because odd cases need to be executed.
+		*/
+		iy = min(iy, ((ty - tx) + 1) >> 1 )
+
+		/*
+			Execute loop.
+		*/
+		#no_bounds_check for iz := uint(0); iz < iy; iz += 1 {
+			_W += _WORD(src.digit[tx + iz]) * _WORD(src.digit[ty - iz])
+		}
+
+		/*
+			Double the inner product and add carry.
+		*/
+		_W = _W + _W + W1
+
+		/*
+			Even columns have the square term in them.
+		*/
+		if ix & 1 == 0 {
+			_W += _WORD(src.digit[ix >> 1]) * _WORD(src.digit[ix >> 1])
+		}
+
+		/*
+			Store it.
+		*/
+		W[ix] = DIGIT(_W & _WORD(_MASK))
+
+		/*
+			Make next carry.
+		*/
+		W1 = _W >> _DIGIT_BITS
+	}
+
+	/*
+		Setup dest.
+	*/
+	old_used := dest.used
+	dest.used = src.used + src.used
+
+	#no_bounds_check for ix = 0; ix < pa; ix += 1 {
+		dest.digit[ix] = W[ix] & _MASK
+	}
+
+	/*
+		Clear unused digits [that existed in the old copy of dest].
+	*/
+	internal_zero_unused(dest, old_used)
+
+	return internal_clamp(dest)
+}
+
+/*
+	Karatsuba squaring, computes `dest` = `src` * `src` using three half-size squarings.
+ 
+ 	See comments of `_private_int_mul_karatsuba` for details.
+ 	It is essentially the same algorithm but merely tuned to perform recursive squarings.
+*/
+_private_int_sqr_karatsuba :: proc(dest, src: ^Int, allocator := context.allocator) -> (err: Error) {
+	context.allocator = allocator
+
+	x0, x1, t1, t2, x0x0, x1x1 := &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}
+	defer internal_destroy(x0, x1, t1, t2, x0x0, x1x1)
+
+	/*
+		Min # of digits, divided by two.
+	*/
+	B := src.used >> 1
+
+	/*
+		Init temps.
+	*/
+	internal_grow(x0,   B) or_return
+	internal_grow(x1,   src.used - B) or_return
+	internal_grow(t1,   src.used * 2) or_return
+	internal_grow(t2,   src.used * 2) or_return
+	internal_grow(x0x0, B * 2       ) or_return
+	internal_grow(x1x1, (src.used - B) * 2) or_return
+
+	/*
+		Now shift the digits.
+	*/
+	x0.used = B
+	x1.used = src.used - B
+
+	#force_inline internal_copy_digits(x0, src, x0.used)
+	#force_inline mem.copy_non_overlapping(&x1.digit[0], &src.digit[B], size_of(DIGIT) * x1.used)
+	#force_inline internal_clamp(x0)
+
+	/*
+		Now calc the products x0*x0 and x1*x1.
+	*/
+	internal_sqr(x0x0, x0) or_return
+	internal_sqr(x1x1, x1) or_return
+
+	/*
+		Now calc (x1+x0)^2
+	*/
+	internal_add(t1, x0, x1) or_return
+	internal_sqr(t1, t1) or_return
+
+	/*
+		Add x0y0
+	*/
+	internal_add(t2, x0x0, x1x1) or_return
+	internal_sub(t1, t1, t2) or_return
+
+	/*
+		Shift by B.
+	*/
+	_private_int_shl_leg(t1, B) or_return
+	_private_int_shl_leg(x1x1, B * 2) or_return
+	internal_add(t1, t1, x0x0) or_return
+	internal_add(dest, t1, x1x1) or_return
+
+	return #force_inline internal_clamp(dest)
+}
+
+/*
+	Squaring using Toom-Cook 3-way algorithm.
+
+	Setup and interpolation from algorithm SQR_3 in Chung, Jaewook, and M. Anwar Hasan. "Asymmetric squaring formulae."
+	  18th IEEE Symposium on Computer Arithmetic (ARITH'07). IEEE, 2007.
+*/
+_private_int_sqr_toom :: proc(dest, src: ^Int, allocator := context.allocator) -> (err: Error) {
+	context.allocator = allocator
+
+	S0, a0, a1, a2 := &Int{}, &Int{}, &Int{}, &Int{}
+	defer internal_destroy(S0, a0, a1, a2)
+
+	/*
+		Init temps.
+	*/
+	internal_zero(S0) or_return
+
+	/*
+		B
+	*/
+	B := src.used / 3
+
+	/*
+		a = a2 * x^2 + a1 * x + a0;
+	*/
+	internal_grow(a0, B) or_return
+	internal_grow(a1, B) or_return
+	internal_grow(a2, src.used - (2 * B)) or_return
+
+	a0.used = B
+	a1.used = B
+	a2.used = src.used - 2 * B
+
+	#force_inline mem.copy_non_overlapping(&a0.digit[0], &src.digit[    0], size_of(DIGIT) * a0.used)
+	#force_inline mem.copy_non_overlapping(&a1.digit[0], &src.digit[    B], size_of(DIGIT) * a1.used)
+	#force_inline mem.copy_non_overlapping(&a2.digit[0], &src.digit[2 * B], size_of(DIGIT) * a2.used)
+
+	internal_clamp(a0)
+	internal_clamp(a1)
+	internal_clamp(a2)
+
+	/** S0 = a0^2;  */
+	internal_sqr(S0, a0) or_return
+
+	/** \\S1 = (a2 + a1 + a0)^2 */
+	/** \\S2 = (a2 - a1 + a0)^2  */
+	/** \\S1 = a0 + a2; */
+	/** a0 = a0 + a2; */
+	internal_add(a0, a0, a2) or_return
+	/** \\S2 = S1 - a1; */
+	/** b = a0 - a1; */
+	internal_sub(dest, a0, a1) or_return
+	/** \\S1 = S1 + a1; */
+	/** a0 = a0 + a1; */
+	internal_add(a0, a0, a1) or_return
+	/** \\S1 = S1^2;  */
+	/** a0 = a0^2; */
+	internal_sqr(a0, a0) or_return
+	/** \\S2 = S2^2;  */
+	/** b = b^2; */
+	internal_sqr(dest, dest) or_return
+	/** \\ S3 = 2 * a1 * a2  */
+	/** \\S3 = a1 * a2;  */
+	/** a1 = a1 * a2; */
+	internal_mul(a1, a1, a2) or_return
+	/** \\S3 = S3 << 1;  */
+	/** a1 = a1 << 1; */
+	internal_shl(a1, a1, 1) or_return
+	/** \\S4 = a2^2;  */
+	/** a2 = a2^2; */
+	internal_sqr(a2, a2) or_return
+	/** \\ tmp = (S1 + S2)/2  */
+	/** \\tmp = S1 + S2; */
+	/** b = a0 + b; */
+	internal_add(dest, a0, dest) or_return
+	/** \\tmp = tmp >> 1; */
+	/** b = b >> 1; */
+	internal_shr(dest, dest, 1) or_return
+	/** \\ S1 = S1 - tmp - S3  */
+	/** \\S1 = S1 - tmp; */
+	/** a0 = a0 - b; */
+	internal_sub(a0, a0, dest) or_return
+	/** \\S1 = S1 - S3;  */
+	/** a0 = a0 - a1; */
+	internal_sub(a0, a0, a1) or_return
+	/** \\S2 = tmp - S4 -S0  */
+	/** \\S2 = tmp - S4;  */
+	/** b = b - a2; */
+	internal_sub(dest, dest, a2) or_return
+	/** \\S2 = S2 - S0;  */
+	/** b = b - S0; */
+	internal_sub(dest, dest, S0) or_return
+	/** \\P = S4*x^4 + S3*x^3 + S2*x^2 + S1*x + S0; */
+	/** P = a2*x^4 + a1*x^3 + b*x^2 + a0*x + S0; */
+	_private_int_shl_leg(  a2, 4 * B) or_return
+	_private_int_shl_leg(  a1, 3 * B) or_return
+	_private_int_shl_leg(dest, 2 * B) or_return
+	_private_int_shl_leg(  a0, 1 * B) or_return
+
+	internal_add(a2, a2, a1) or_return
+	internal_add(dest, dest, a2) or_return
+	internal_add(dest, dest, a0) or_return
+	internal_add(dest, dest, S0) or_return
+	/** a^2 - P  */
+
+	return #force_inline internal_clamp(dest)
+}
+
+/*
+	Divide by three (based on routine from MPI and the GMP manual).
+*/
+_private_int_div_3 :: proc(quotient, numerator: ^Int, allocator := context.allocator) -> (remainder: DIGIT, err: Error) {
+	context.allocator = allocator
+
+	/*
+		b = 2^_DIGIT_BITS / 3
+	*/
+ 	b := _WORD(1) << _WORD(_DIGIT_BITS) / _WORD(3)
+
+	q := &Int{}
+	internal_grow(q, numerator.used) or_return
+	q.used = numerator.used
+	q.sign = numerator.sign
+
+	w, t: _WORD
+	#no_bounds_check for ix := numerator.used; ix >= 0; ix -= 1 {
+		w = (w << _WORD(_DIGIT_BITS)) | _WORD(numerator.digit[ix])
+		if w >= 3 {
+			/*
+				Multiply w by [1/3].
+			*/
+			t = (w * b) >> _WORD(_DIGIT_BITS)
+
+			/*
+				Now subtract 3 * [w/3] from w, to get the remainder.
+			*/
+			w -= t+t+t
+
+			/*
+				Fixup the remainder as required since the optimization is not exact.
+			*/
+			for w >= 3 {
+				t += 1
+				w -= 3
+			}
+		} else {
+			t = 0
+		}
+		q.digit[ix] = DIGIT(t)
+	}
+	remainder = DIGIT(w)
+
+	/*
+		[optional] store the quotient.
+	*/
+	if quotient != nil {
+		err = clamp(q)
+ 		internal_swap(q, quotient)
+ 	}
+	internal_destroy(q)
+	return remainder, nil
+}
+
+/*
+	Signed Integer Division
+
+	c*b + d == a [i.e. a/b, c=quotient, d=remainder], HAC pp.598 Algorithm 14.20
+
+	Note that the description in HAC is horribly incomplete.
+	For example, it doesn't consider the case where digits are removed from 'x' in
+	the inner loop.
+
+	It also doesn't consider the case that y has fewer than three digits, etc.
+	The overall algorithm is as described as 14.20 from HAC but fixed to treat these cases.
+*/
+_private_int_div_school :: proc(quotient, remainder, numerator, denominator: ^Int, allocator := context.allocator) -> (err: Error) {
+	context.allocator = allocator
+
+	error_if_immutable(quotient, remainder) or_return
+
+	q, x, y, t1, t2 := &Int{}, &Int{}, &Int{}, &Int{}, &Int{}
+	defer internal_destroy(q, x, y, t1, t2)
+
+	internal_grow(q, numerator.used + 2) or_return
+	q.used = numerator.used + 2
+
+	internal_init_multi(t1, t2) or_return
+	internal_copy(x, numerator) or_return
+	internal_copy(y, denominator) or_return
+
+	/*
+		Fix the sign.
+	*/
+	neg   := numerator.sign != denominator.sign
+	x.sign = .Zero_or_Positive
+	y.sign = .Zero_or_Positive
+
+	/*
+		Normalize both x and y, ensure that y >= b/2, [b == 2**MP_DIGIT_BIT]
+	*/
+	norm := internal_count_bits(y) % _DIGIT_BITS
+
+	if norm < _DIGIT_BITS - 1 {
+		norm = (_DIGIT_BITS - 1) - norm
+		internal_shl(x, x, norm) or_return
+		internal_shl(y, y, norm) or_return
+	} else {
+		norm = 0
+	}
+
+	/*
+		Note: HAC does 0 based, so if used==5 then it's 0,1,2,3,4, i.e. use 4
+	*/
+	n := x.used - 1
+	t := y.used - 1
+
+	/*
+		while (x >= y*b**n-t) do { q[n-t] += 1; x -= y*b**{n-t} }
+		y = y*b**{n-t}
+	*/
+
+	_private_int_shl_leg(y, n - t) or_return
+
+	gte := internal_gte(x, y)
+	for gte {
+		q.digit[n - t] += 1
+		internal_sub(x, x, y) or_return
+		gte = internal_gte(x, y)
+	}
+
+	/*
+		Reset y by shifting it back down.
+	*/
+	_private_int_shr_leg(y, n - t)
+
+	/*
+		Step 3. for i from n down to (t + 1).
+	*/
+	#no_bounds_check for i := n; i >= (t + 1); i -= 1 {
+		if i > x.used { continue }
+
+		/*
+			step 3.1 if xi == yt then set q{i-t-1} to b-1, otherwise set q{i-t-1} to (xi*b + x{i-1})/yt
+		*/
+		if x.digit[i] == y.digit[t] {
+			q.digit[(i - t) - 1] = 1 << (_DIGIT_BITS - 1)
+		} else {
+
+			tmp := _WORD(x.digit[i]) << _DIGIT_BITS
+			tmp |= _WORD(x.digit[i - 1])
+			tmp /= _WORD(y.digit[t])
+			if tmp > _WORD(_MASK) {
+				tmp = _WORD(_MASK)
+			}
+			q.digit[(i - t) - 1] = DIGIT(tmp & _WORD(_MASK))
+		}
+
+		/* while (q{i-t-1} * (yt * b + y{t-1})) >
+					xi * b**2 + xi-1 * b + xi-2
+
+			do q{i-t-1} -= 1;
+		*/
+
+		iter := 0
+
+		q.digit[(i - t) - 1] = (q.digit[(i - t) - 1] + 1) & _MASK
+		#no_bounds_check for {
+			q.digit[(i - t) - 1] = (q.digit[(i - t) - 1] - 1) & _MASK
+
+			/*
+				Find left hand.
+			*/
+			internal_zero(t1)
+			t1.digit[0] = ((t - 1) < 0) ? 0 : y.digit[t - 1]
+			t1.digit[1] = y.digit[t]
+			t1.used = 2
+			internal_mul(t1, t1, q.digit[(i - t) - 1]) or_return
+
+			/*
+				Find right hand.
+			*/
+			t2.digit[0] = ((i - 2) < 0) ? 0 : x.digit[i - 2]
+			t2.digit[1] = x.digit[i - 1] /* i >= 1 always holds */
+			t2.digit[2] = x.digit[i]
+			t2.used = 3
+
+			if internal_lte(t1, t2) {
+				break
+			}
+			iter += 1; if iter > 100 {
+				return .Max_Iterations_Reached
+			}
+		}
+
+		/*
+			Step 3.3 x = x - q{i-t-1} * y * b**{i-t-1}
+		*/
+		int_mul_digit(t1, y, q.digit[(i - t) - 1]) or_return
+		_private_int_shl_leg(t1, (i - t) - 1) or_return
+		internal_sub(x, x, t1) or_return
+
+		/*
+			if x < 0 then { x = x + y*b**{i-t-1}; q{i-t-1} -= 1; }
+		*/
+		if x.sign == .Negative {
+			internal_copy(t1, y) or_return
+			_private_int_shl_leg(t1, (i - t) - 1) or_return
+			internal_add(x, x, t1) or_return
+
+			q.digit[(i - t) - 1] = (q.digit[(i - t) - 1] - 1) & _MASK
+		}
+	}
+
+	/*
+		Now q is the quotient and x is the remainder, [which we have to normalize]
+		Get sign before writing to c.
+	*/
+	z, _ := is_zero(x)
+	x.sign = .Zero_or_Positive if z else numerator.sign
+
+	if quotient != nil {
+		internal_clamp(q)
+		internal_swap(q, quotient)
+		quotient.sign = .Negative if neg else .Zero_or_Positive
+	}
+
+	if remainder != nil {
+		internal_shr(x, x, norm) or_return
+		internal_swap(x, remainder)
+	}
+
+	return nil
+}
+
+/*
+	Direct implementation of algorithms 1.8 "RecursiveDivRem" and 1.9 "UnbalancedDivision" from:
+
+		Brent, Richard P., and Paul Zimmermann. "Modern computer arithmetic"
+		Vol. 18. Cambridge University Press, 2010
+		Available online at https://arxiv.org/pdf/1004.4710
+
+	pages 19ff. in the above online document.
+*/
+_private_div_recursion :: proc(quotient, remainder, a, b: ^Int, allocator := context.allocator) -> (err: Error) {
+	context.allocator = allocator
+
+	A1, A2, B1, B0, Q1, Q0, R1, R0, t := &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}
+	defer internal_destroy(A1, A2, B1, B0, Q1, Q0, R1, R0, t)
+
+	m := a.used - b.used
+	k := m / 2
+
+	if m < MUL_KARATSUBA_CUTOFF {
+		return _private_int_div_school(quotient, remainder, a, b)
+	}
+
+	internal_init_multi(A1, A2, B1, B0, Q1, Q0, R1, R0, t) or_return
+
+	/*
+		`B1` = `b` / `beta`^`k`, `B0` = `b` % `beta`^`k`
+	*/
+	internal_shrmod(B1, B0, b, k * _DIGIT_BITS) or_return
+
+	/*
+		(Q1, R1) =  RecursiveDivRem(A / beta^(2k), B1)
+	*/
+	internal_shrmod(A1, t, a, 2 * k * _DIGIT_BITS) or_return
+	_private_div_recursion(Q1, R1, A1, B1) or_return
+
+	/*
+		A1 = (R1 * beta^(2k)) + (A % beta^(2k)) - (Q1 * B0 * beta^k)
+	*/
+	_private_int_shl_leg(R1, 2 * k) or_return
+	internal_add(A1, R1, t) or_return
+	internal_mul(t, Q1, B0) or_return
+
+	/*
+		While A1 < 0 do Q1 = Q1 - 1, A1 = A1 + (beta^k * B)
+	*/
+	if internal_lt(A1, 0) {
+		internal_shl(t, b, k * _DIGIT_BITS) or_return
+
+		for {
+			internal_decr(Q1) or_return
+			internal_add(A1, A1, t) or_return
+			if internal_gte(A1, 0) { break }
+		}
+	}
+
+	/*
+		(Q0, R0) =  RecursiveDivRem(A1 / beta^(k), B1)
+	*/
+	internal_shrmod(A1, t, A1, k * _DIGIT_BITS) or_return
+	_private_div_recursion(Q0, R0, A1, B1) or_return
+
+	/*
+		A2 = (R0*beta^k) +  (A1 % beta^k) - (Q0*B0)
+	*/
+	_private_int_shl_leg(R0, k) or_return
+	internal_add(A2, R0, t) or_return
+	internal_mul(t, Q0, B0) or_return
+	internal_sub(A2, A2, t) or_return
+
+	/*
+		While A2 < 0 do Q0 = Q0 - 1, A2 = A2 + B.
+	*/
+	for internal_is_negative(A2) { // internal_lt(A2, 0) {
+		internal_decr(Q0) or_return
+		internal_add(A2, A2, b) or_return
+	}
+
+	/*
+		Return q = (Q1*beta^k) + Q0, r = A2.
+	*/
+	_private_int_shl_leg(Q1, k) or_return
+	internal_add(quotient, Q1, Q0) or_return
+
+	return internal_copy(remainder, A2)
+}
+
+_private_int_div_recursive :: proc(quotient, remainder, a, b: ^Int, allocator := context.allocator) -> (err: Error) {
+	context.allocator = allocator
+
+	A, B, Q, Q1, R, A_div, A_mod := &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}
+	defer internal_destroy(A, B, Q, Q1, R, A_div, A_mod)
+
+	internal_init_multi(A, B, Q, Q1, R, A_div, A_mod) or_return
+
+	/*
+		Most significant bit of a limb.
+		Assumes  _DIGIT_MAX < (sizeof(DIGIT) * sizeof(u8)).
+	*/
+	msb := (_DIGIT_MAX + DIGIT(1)) >> 1
+	sigma := 0
+	msb_b := b.digit[b.used - 1]
+	for msb_b < msb {
+		sigma += 1
+		msb_b <<= 1
+	}
+
+	/*
+		Use that sigma to normalize B.
+	*/
+	internal_shl(B, b, sigma) or_return
+	internal_shl(A, a, sigma) or_return
+
+	/*
+		Fix the sign.
+	*/
+	neg := a.sign != b.sign
+	A.sign = .Zero_or_Positive; B.sign = .Zero_or_Positive
+
+	/*
+		If the magnitude of "A" is not more more than twice that of "B" we can work
+		on them directly, otherwise we need to work at "A" in chunks.
+	*/
+	n := B.used
+	m := A.used - B.used
+
+	/*
+		Q = 0. We already ensured that when we called `internal_init_multi`.
+	*/
+	for m > n {
+		/*
+			(q, r) = RecursiveDivRem(A / (beta^(m-n)), B)
+		*/
+		j := (m - n) * _DIGIT_BITS
+		internal_shrmod(A_div, A_mod, A, j) or_return
+		_private_div_recursion(Q1, R, A_div, B) or_return
+
+		/*
+			Q = (Q*beta!(n)) + q
+		*/
+		internal_shl(Q, Q, n * _DIGIT_BITS) or_return
+		internal_add(Q, Q, Q1) or_return
+
+		/*
+			A = (r * beta^(m-n)) + (A % beta^(m-n))
+		*/
+		internal_shl(R, R, (m - n) * _DIGIT_BITS) or_return
+		internal_add(A, R, A_mod) or_return
+
+		/*
+			m = m - n
+		*/
+		m -= n
+	}
+
+	/*
+		(q, r) = RecursiveDivRem(A, B)
+	*/
+	_private_div_recursion(Q1, R, A, B) or_return
+
+	/*
+		Q = (Q * beta^m) + q, R = r
+	*/
+	internal_shl(Q, Q, m * _DIGIT_BITS) or_return
+	internal_add(Q, Q, Q1) or_return
+
+	/*
+		Get sign before writing to dest.
+	*/
+	R.sign = .Zero_or_Positive if internal_is_zero(Q) else a.sign
+
+	if quotient != nil {
+		swap(quotient, Q)
+		quotient.sign = .Negative if neg else .Zero_or_Positive
+	}
+	if remainder != nil {
+		/*
+			De-normalize the remainder.
+		*/
+		internal_shrmod(R, nil, R, sigma) or_return
+		swap(remainder, R)
+	}
+	return nil
+}
+
+/*
+	Slower bit-bang division... also smaller.
+*/
+@(deprecated="Use `_int_div_school`, it's 3.5x faster.")
+_private_int_div_small :: proc(quotient, remainder, numerator, denominator: ^Int) -> (err: Error) {
+
+	ta, tb, tq, q := &Int{}, &Int{}, &Int{}, &Int{}
+
+	defer internal_destroy(ta, tb, tq, q)
+
+	for {
+		internal_one(tq) or_return
+
+		num_bits, _ := count_bits(numerator)
+		den_bits, _ := count_bits(denominator)
+		n := num_bits - den_bits
+
+		abs(ta, numerator)   or_return
+		abs(tb, denominator) or_return
+		shl(tb, tb, n)       or_return
+		shl(tq, tq, n)       or_return
+
+		for n >= 0 {
+			if internal_gte(ta, tb) {
+				// ta -= tb
+				sub(ta, ta, tb) or_return
+				//  q += tq
+				add( q, q,  tq) or_return
+			}
+			shr1(tb, tb) or_return
+			shr1(tq, tq) or_return
+
+			n -= 1
+		}
+
+		/*
+			Now q == quotient and ta == remainder.
+		*/
+		neg := numerator.sign != denominator.sign
+		if quotient != nil {
+			swap(quotient, q)
+			z, _ := is_zero(quotient)
+			quotient.sign = .Negative if neg && !z else .Zero_or_Positive
+		}
+		if remainder != nil {
+			swap(remainder, ta)
+			z, _ := is_zero(numerator)
+			remainder.sign = .Zero_or_Positive if z else numerator.sign
+		}
+
+		break
+	}
+	return err
+}
+
+
+
+/*
+	Binary split factorial algo due to: http://www.luschny.de/math/factorial/binarysplitfact.html
+*/
+_private_int_factorial_binary_split :: proc(res: ^Int, n: int, allocator := context.allocator) -> (err: Error) {
+	context.allocator = allocator
+
+	inner, outer, start, stop, temp := &Int{}, &Int{}, &Int{}, &Int{}, &Int{}
+	defer internal_destroy(inner, outer, start, stop, temp)
+
+	internal_one(inner, false)                                       or_return
+	internal_one(outer, false)                                       or_return
+
+	bits_used := ilog2(n)
+
+	for i := bits_used; i >= 0; i -= 1 {
+		start := (n >> (uint(i) + 1)) + 1 | 1
+		stop  := (n >> uint(i)) + 1 | 1
+		_private_int_recursive_product(temp, start, stop, 0)         or_return
+		internal_mul(inner, inner, temp)                             or_return
+		internal_mul(outer, outer, inner)                            or_return
+	}
+	shift := n - intrinsics.count_ones(n)
+
+	return internal_shl(res, outer, int(shift))
+}
+
+/*
+	Recursive product used by binary split factorial algorithm.
+*/
+_private_int_recursive_product :: proc(res: ^Int, start, stop: int, level := int(0), allocator := context.allocator) -> (err: Error) {
+	context.allocator = allocator
+
+	t1, t2 := &Int{}, &Int{}
+	defer internal_destroy(t1, t2)
+
+	if level > FACTORIAL_BINARY_SPLIT_MAX_RECURSIONS {
+		return .Max_Iterations_Reached
+	}
+
+	num_factors := (stop - start) >> 1
+	if num_factors == 2 {
+		internal_set(t1, start, false)                               or_return
+		when true {
+			internal_grow(t2, t1.used + 1, false)                    or_return
+			internal_add(t2, t1, 2)                                  or_return
+		} else {
+			internal_add(t2, t1, 2)                                  or_return
+		}
+		return internal_mul(res, t1, t2)
+	}
+
+	if num_factors > 1 {
+		mid := (start + num_factors) | 1
+		_private_int_recursive_product(t1, start,  mid, level + 1)   or_return
+		_private_int_recursive_product(t2,   mid, stop, level + 1)   or_return
+		return internal_mul(res, t1, t2)
+	}
+
+	if num_factors == 1 {
+		return #force_inline internal_set(res, start, true)
+	}
+
+	return #force_inline internal_one(res, true)
+}
+
+/*
+	Internal function computing both GCD using the binary method,
+		and, if target isn't `nil`, also LCM.
+
+	Expects the `a` and `b` to have been initialized
+		and one or both of `res_gcd` or `res_lcm` not to be `nil`.
+
+	If both `a` and `b` are zero, return zero.
+	If either `a` or `b`, return the other one.
+
+	The `gcd` and `lcm` wrappers have already done this test,
+	but `gcd_lcm` wouldn't have, so we still need to perform it.
+
+	If neither result is wanted, we have nothing to do.
+*/
+_private_int_gcd_lcm :: proc(res_gcd, res_lcm, a, b: ^Int, allocator := context.allocator) -> (err: Error) {
+	context.allocator = allocator
+
+	if res_gcd == nil && res_lcm == nil {
+		return nil
+	}
+
+	/*
+		We need a temporary because `res_gcd` is allowed to be `nil`.
+	*/
+	if a.used == 0 && b.used == 0 {
+		/*
+			GCD(0, 0) and LCM(0, 0) are both 0.
+		*/
+		if res_gcd != nil {
+			internal_zero(res_gcd) or_return
+		}
+		if res_lcm != nil {
+			internal_zero(res_lcm) or_return
+		}
+		return nil
+	} else if a.used == 0 {
+		/*
+			We can early out with GCD = B and LCM = 0
+		*/
+		if res_gcd != nil {
+			internal_abs(res_gcd, b) or_return
+		}
+		if res_lcm != nil {
+			internal_zero(res_lcm) or_return
+		}
+		return nil
+	} else if b.used == 0 {
+		/*
+			We can early out with GCD = A and LCM = 0
+		*/
+		if res_gcd != nil {
+			internal_abs(res_gcd, a) or_return
+		}
+		if res_lcm != nil {
+			internal_zero(res_lcm) or_return
+		}
+		return nil
+	}
+
+	temp_gcd_res := &Int{}
+	defer internal_destroy(temp_gcd_res)
+
+	/*
+		If neither `a` or `b` was zero, we need to compute `gcd`.
+ 		Get copies of `a` and `b` we can modify.
+ 	*/
+	u, v := &Int{}, &Int{}
+	defer internal_destroy(u, v)
+	internal_copy(u, a) or_return
+	internal_copy(v, b) or_return
+
+ 	/*
+ 		Must be positive for the remainder of the algorithm.
+ 	*/
+	u.sign = .Zero_or_Positive; v.sign = .Zero_or_Positive
+
+ 	/*
+ 		B1.  Find the common power of two for `u` and `v`.
+ 	*/
+ 	u_lsb, _ := internal_count_lsb(u)
+ 	v_lsb, _ := internal_count_lsb(v)
+ 	k        := min(u_lsb, v_lsb)
+
+	if k > 0 {
+		/*
+			Divide the power of two out.
+		*/
+		internal_shr(u, u, k) or_return
+		internal_shr(v, v, k) or_return
+	}
+
+	/*
+		Divide any remaining factors of two out.
+	*/
+	if u_lsb != k {
+		internal_shr(u, u, u_lsb - k) or_return
+	}
+	if v_lsb != k {
+		internal_shr(v, v, v_lsb - k) or_return
+	}
+
+	for v.used != 0 {
+		/*
+			Make sure `v` is the largest.
+		*/
+		if internal_gt(u, v) {
+			/*
+				Swap `u` and `v` to make sure `v` is >= `u`.
+			*/
+			internal_swap(u, v)
+		}
+
+		/*
+			Subtract smallest from largest.
+		*/
+		internal_sub(v, v, u) or_return
+
+		/*
+			Divide out all factors of two.
+		*/
+		b, _ := internal_count_lsb(v)
+		internal_shr(v, v, b) or_return
+	}
+
+ 	/*
+ 		Multiply by 2**k which we divided out at the beginning.
+ 	*/
+ 	internal_shl(temp_gcd_res, u, k) or_return
+ 	temp_gcd_res.sign = .Zero_or_Positive
+
+	/*
+		We've computed `gcd`, either the long way, or because one of the inputs was zero.
+		If we don't want `lcm`, we're done.
+	*/
+	if res_lcm == nil {
+		internal_swap(temp_gcd_res, res_gcd)
+		return nil
+	}
+
+	/*
+		Computes least common multiple as `|a*b|/gcd(a,b)`
+		Divide the smallest by the GCD.
+	*/
+	if internal_lt_abs(a, b) {
+		/*
+			Store quotient in `t2` such that `t2 * b` is the LCM.
+		*/
+		internal_div(res_lcm, a, temp_gcd_res) or_return
+		err = internal_mul(res_lcm, res_lcm, b)
+	} else {
+		/*
+			Store quotient in `t2` such that `t2 * a` is the LCM.
+		*/
+		internal_div(res_lcm, b, temp_gcd_res) or_return
+		err = internal_mul(res_lcm, res_lcm, a)
+	}
+
+	if res_gcd != nil {
+		internal_swap(temp_gcd_res, res_gcd)
+	}
+
+	/*
+		Fix the sign to positive and return.
+	*/
+	res_lcm.sign = .Zero_or_Positive
+	return err
+}
+
+/*
+	Internal implementation of log.
+	Assumes `a` not to be `nil` and to have been initialized.
+*/
+_private_int_log :: proc(a: ^Int, base: DIGIT, allocator := context.allocator) -> (res: int, err: Error) {
+	bracket_low, bracket_high, bracket_mid, t, bi_base := &Int{}, &Int{}, &Int{}, &Int{}, &Int{}
+	defer internal_destroy(bracket_low, bracket_high, bracket_mid, t, bi_base)
+
+	ic := #force_inline internal_cmp(a, base)
+	if ic == -1 || ic == 0 {
+		return 1 if ic == 0 else 0, nil
+	}
+	defer if err != nil {
+		res = -1
+	}
+
+	internal_set(bi_base, base, true, allocator)       or_return
+	internal_clear(bracket_mid, false, allocator)      or_return
+	internal_clear(t, false, allocator)                or_return
+	internal_one(bracket_low, false, allocator)        or_return
+	internal_set(bracket_high, base, false, allocator) or_return
+
+	low := 0; high := 1
+
+	/*
+		A kind of Giant-step/baby-step algorithm.
+		Idea shamelessly stolen from https://programmingpraxis.com/2010/05/07/integer-logarithms/2/
+		The effect is asymptotic, hence needs benchmarks to test if the Giant-step should be skipped
+		for small n.
+	*/
+
+	for {
+		/*
+			Iterate until `a` is bracketed between low + high.
+		*/
+		if #force_inline internal_gte(bracket_high, a) { break }
+
+		low = high
+		#force_inline internal_copy(bracket_low, bracket_high) or_return
+		high <<= 1
+		#force_inline internal_sqr(bracket_high, bracket_high) or_return
+	}
+
+	for (high - low) > 1 {
+		mid := (high + low) >> 1
+
+		#force_inline internal_pow(t, bi_base, mid - low) or_return
+
+		#force_inline internal_mul(bracket_mid, bracket_low, t) or_return
+
+		mc := #force_inline internal_cmp(a, bracket_mid)
+		switch mc {
+		case -1:
+			high = mid
+			internal_swap(bracket_mid, bracket_high)
+		case  0:
+			return mid, nil
+		case  1:
+			low = mid
+			internal_swap(bracket_mid, bracket_low)
+		}
+	}
+
+	fc := #force_inline internal_cmp(bracket_high, a)
+	res = high if fc == 0 else low
+
+	return
+}
+
+/*
+	Computes xR**-1 == x (mod N) via Montgomery Reduction.
+	This is an optimized implementation of `internal_montgomery_reduce`
+	which uses the comba method to quickly calculate the columns of the reduction.
+	Based on Algorithm 14.32 on pp.601 of HAC.
+*/
+_private_montgomery_reduce_comba :: proc(x, n: ^Int, rho: DIGIT, allocator := context.allocator) -> (err: Error) {
+	context.allocator = allocator
+	W: [_WARRAY]_WORD = ---
+
+	if x.used > _WARRAY { return .Invalid_Argument }
+
+	/*
+		Get old used count.
+	*/
+	old_used := x.used
+
+	/*
+		Grow `x` as required.
+	*/
+	internal_grow(x, n.used + 1)                                     or_return
+
+	/*
+		First we have to get the digits of the input into an array of double precision words W[...]
+		Copy the digits of `x` into W[0..`x.used` - 1]
+	*/
+	ix: int
+	for ix = 0; ix < x.used; ix += 1 {
+		W[ix] = _WORD(x.digit[ix])
+	}
+
+	/*
+		Zero the high words of W[a->used..m->used*2].
+	*/
+	zero_upper := (n.used * 2) + 1
+	if ix < zero_upper {
+		for ix = x.used; ix < zero_upper; ix += 1 {
+			W[ix] = {}
+		}
+	}
+
+	/*
+		Now we proceed to zero successive digits from the least significant upwards.
+	*/
+	for ix = 0; ix < n.used; ix += 1 {
+		/*
+			`mu = ai * m' mod b`
+
+			We avoid a double precision multiplication (which isn't required)
+			by casting the value down to a DIGIT.  Note this requires
+			that W[ix-1] have the carry cleared (see after the inner loop)
+		*/
+		mu := ((W[ix] & _WORD(_MASK)) * _WORD(rho)) & _WORD(_MASK)
+
+		/*
+			`a = a + mu * m * b**i`
+		
+			This is computed in place and on the fly.  The multiplication
+		 	by b**i is handled by offseting which columns the results
+		 	are added to.
+		
+			Note the comba method normally doesn't handle carries in the
+			inner loop In this case we fix the carry from the previous
+			column since the Montgomery reduction requires digits of the
+			result (so far) [see above] to work.
+
+			This is	handled by fixing up one carry after the inner loop.
+			The carry fixups are done in order so after these loops the
+			first m->used words of W[] have the carries fixed.
+		*/
+		for iy := 0; iy < n.used; iy += 1 {
+			W[ix + iy] += mu * _WORD(n.digit[iy])
+		}
+
+		/*
+			Now fix carry for next digit, W[ix+1].
+		*/
+		W[ix + 1] += (W[ix] >> _DIGIT_BITS)
+	}
+
+	/*
+		Now we have to propagate the carries and shift the words downward
+		[all those least significant digits we zeroed].
+	*/
+
+	for ; ix < n.used * 2; ix += 1 {
+		W[ix + 1] += (W[ix] >> _DIGIT_BITS)
+	}
+
+	/* copy out, A = A/b**n
+	 *
+	 * The result is A/b**n but instead of converting from an
+	 * array of mp_word to mp_digit than calling mp_rshd
+	 * we just copy them in the right order
+	 */
+
+	for ix = 0; ix < (n.used + 1); ix += 1 {
+		x.digit[ix] = DIGIT(W[n.used + ix] & _WORD(_MASK))
+	}
+
+	/*
+		Set the max used.
+	*/
+	x.used = n.used + 1
+
+	/*
+		Zero old_used digits, if the input a was larger than m->used+1 we'll have to clear the digits.
+	*/
+	internal_zero_unused(x, old_used)
+	internal_clamp(x)
+
+	/*
+		if A >= m then A = A - m
+	*/
+	if internal_gte_abs(x, n) {
+		return internal_sub(x, x, n)
+	}
+	return nil
+}
+
+/*
+	Computes xR**-1 == x (mod N) via Montgomery Reduction.
+	Assumes `x` and `n` not to be nil.
+*/
+_private_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.
+	*/
+	internal_clear_if_uninitialized(x, n) or_return
+
+	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)
+	_private_int_shr_leg(x, n.used)
+
+	/*
+		if x >= n then x = x - n
+	*/
+	if internal_gte_abs(x, n) {
+		return internal_sub(x, x, n)
+	}
+
+	return nil
+}
+
+/*
+	Shifts with subtractions when the result is greater than b.
+
+	The method is slightly modified to shift B unconditionally upto just under
+	the leading bit of b.  This saves alot of multiple precision shifting.
+
+	Assumes `a` and `b` not to be `nil`.
+*/
+_private_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.
+	*/
+	internal_clear_if_uninitialized(a, b) or_return
+
+	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)                                              or_return
+		bits = 1
+	}
+
+	/*
+		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_gte_abs(a, b) {
+			internal_sub(a, a, b)                                    or_return
+		}
+	}
+	return nil
+}
+
+/*
+	Sets up the Montgomery reduction stuff.
+*/
+_private_int_montgomery_setup :: proc(n: ^Int, allocator := context.allocator) -> (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
+	*/
+	internal_clear_if_uninitialized(n, allocator) or_return
+
+	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 _DIGIT_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 */
+	}
+
+	/*
+		rho = -1/m mod b
+	*/
+	rho = DIGIT(((_WORD(1) << _WORD(_DIGIT_BITS)) - _WORD(x)) & _WORD(_MASK))
+	return rho, nil
+}
+
+/*
+	Reduces `x` mod `m`, assumes 0 < x < m**2, mu is precomputed via reduce_setup.
+	From HAC pp.604 Algorithm 14.42
+
+	Assumes `x`, `m` and `mu` all not to be `nil` and have been initialized.
+*/
+_private_int_reduce :: proc(x, m, mu: ^Int, allocator := context.allocator) -> (err: Error) {
+	context.allocator = allocator
+
+	q := &Int{}
+	defer internal_destroy(q)
+	um := m.used
+
+	/*
+		q = x
+	*/
+	internal_copy(q, x)                                              or_return
+
+	/*
+		q1 = x / b**(k-1)
+	*/
+	_private_int_shr_leg(q, um - 1)
+
+	/*
+		According to HAC this optimization is ok.
+	*/
+	if DIGIT(um) > DIGIT(1) << (_DIGIT_BITS - 1) {
+		internal_mul(q, q, mu)                                       or_return
+	} else {
+		_private_int_mul_high(q, q, mu, um)                          or_return
+	}
+
+	/*
+		q3 = q2 / b**(k+1)
+	*/
+	_private_int_shr_leg(q, um + 1)
+
+	/*
+		x = x mod b**(k+1), quick (no division)
+	*/
+	internal_int_mod_bits(x, x, _DIGIT_BITS * (um + 1))              or_return
+
+	/*
+		q = q * m mod b**(k+1), quick (no division)
+	*/
+	_private_int_mul(q, q, m, um + 1)                                or_return
+
+	/*
+		x = x - q
+	*/
+	internal_sub(x, x, q)                                            or_return
+
+	/*
+		If x < 0, add b**(k+1) to it.
+	*/
+	if internal_is_negative(x) {
+		internal_set(q, 1)                                           or_return
+		_private_int_shl_leg(q, um + 1)                                or_return
+		internal_add(x, x, q)                                        or_return
+	}
+
+	/*
+		Back off if it's too big.
+	*/
+	for internal_gte(x, m) {
+		internal_sub(x, x, m)                                        or_return
+	}
+
+	return nil
+}
+
+/*
+	Reduces `a` modulo `n`, where `n` is of the form 2**p - d.
+*/
+_private_int_reduce_2k :: proc(a, n: ^Int, d: DIGIT, allocator := context.allocator) -> (err: Error) {
+	context.allocator = allocator
+
+	q := &Int{}
+	defer internal_destroy(q)
+
+	internal_zero(q)                                                 or_return
+
+	p := internal_count_bits(n)
+
+	for {
+		/*
+			q = a/2**p, a = a mod 2**p
+		*/
+		internal_shrmod(q, a, a, p)                                  or_return
+
+		if d != 1 {
+			/*
+				q = q * d
+			*/
+			internal_mul(q, q, d)                                    or_return
+		}
+
+		/*
+			a = a + q
+		*/
+		internal_add(a, a, q)                                        or_return
+		if internal_lt_abs(a, n)                                     { break }
+		internal_sub(a, a, n)                                        or_return
+	}
+
+	return nil
+}
+
+/*
+	Reduces `a` modulo `n` where `n` is of the form 2**p - d
+	This differs from reduce_2k since "d" can be larger than a single digit.
+*/
+_private_int_reduce_2k_l :: proc(a, n, d: ^Int, allocator := context.allocator) -> (err: Error) {
+	context.allocator = allocator
+
+	q := &Int{}
+	defer internal_destroy(q)
+
+	internal_zero(q)                                                 or_return
+
+	p := internal_count_bits(n)
+
+	for {
+		/*
+			q = a/2**p, a = a mod 2**p
+		*/
+		internal_shrmod(q, a, a, p)                                  or_return
+
+		/*
+			q = q * d
+		*/
+		internal_mul(q, q, d)                                        or_return
+
+		/*
+			a = a + q
+		*/
+		internal_add(a, a, q)                                        or_return
+		if internal_lt_abs(a, n)                                     { break }
+		internal_sub(a, a, n)                                        or_return
+	}
+
+	return nil
+}
+
+/*
+	Determines if `internal_int_reduce_2k` can be used.
+	Asssumes `a` not to be `nil` and to have been initialized.
+*/
+_private_int_reduce_is_2k :: proc(a: ^Int) -> (reducible: bool, err: Error) {
+	assert_if_nil(a)
+
+	if internal_is_zero(a) {
+		return false, nil
+	} else if a.used == 1 {
+		return true, nil
+	} else if a.used  > 1 {
+		iy := internal_count_bits(a)
+		iw := 1
+		iz := DIGIT(1)
+
+		/*
+			Test every bit from the second digit up, must be 1.
+		*/
+		for ix := _DIGIT_BITS; ix < iy; ix += 1 {
+			if a.digit[iw] & iz == 0 {
+				return false, nil
+			}
+
+			iz <<= 1
+			if iz > _DIGIT_MAX {
+				iw += 1
+				iz  = 1
+			}
+		}
+		return true, nil
+	} else {
+		return true, nil
+	}
+}
+
+/*
+	Determines if `internal_int_reduce_2k_l` can be used.
+	Asssumes `a` not to be `nil` and to have been initialized.
+*/
+_private_int_reduce_is_2k_l :: proc(a: ^Int) -> (reducible: bool, err: Error) {
+	assert_if_nil(a)
+
+	if internal_int_is_zero(a) {
+		return false, nil
+	} else if a.used == 1 {
+		return true, nil
+	} else if a.used  > 1 {
+		/*
+			If more than half of the digits are -1 we're sold.
+		*/
+		ix := 0
+		iy := 0
+
+		for ; ix < a.used; ix += 1 {
+			if a.digit[ix] == _DIGIT_MAX {
+				iy += 1
+			}
+		}
+		return iy >= (a.used / 2), nil
+	} else {
+		return false, nil
+	}
+}
+
+/*
+	Determines the setup value.
+	Assumes `a` is not `nil`.
+*/
+_private_int_reduce_2k_setup :: proc(a: ^Int, allocator := context.allocator) -> (d: DIGIT, err: Error) {
+	context.allocator = allocator
+
+	tmp := &Int{}
+	defer internal_destroy(tmp)
+	internal_zero(tmp)                                               or_return
+
+	internal_int_power_of_two(tmp, internal_count_bits(a))           or_return
+	internal_sub(tmp, tmp, a)                                        or_return
+
+	return tmp.digit[0], nil
+}
+
+/*
+	Determines the setup value.
+	Assumes `mu` and `P` are not `nil`.
+
+	d := (1 << a.bits) - a;
+*/
+_private_int_reduce_2k_setup_l :: proc(mu, P: ^Int, allocator := context.allocator) -> (err: Error) {
+	context.allocator = allocator
+
+	tmp := &Int{}
+	defer internal_destroy(tmp)
+	internal_zero(tmp)                                               or_return
+
+	internal_int_power_of_two(tmp, internal_count_bits(P))           or_return
+	internal_sub(mu, tmp, P)                                         or_return
+
+	return nil
+}
+
+/*
+	Pre-calculate the value required for Barrett reduction.
+	For a given modulus "P" it calulates the value required in "mu"
+	Assumes `mu` and `P` are not `nil`.
+*/
+_private_int_reduce_setup :: proc(mu, P: ^Int, allocator := context.allocator) -> (err: Error) {
+	context.allocator = allocator
+
+	internal_int_power_of_two(mu, P.used * 2 * _DIGIT_BITS)           or_return
+	return internal_int_div(mu, mu, P)
+}
+
+/*
+	Determines the setup value.
+	Assumes `a` to not be `nil` and to have been initialized.
+*/
+_private_int_dr_setup :: proc(a: ^Int) -> (d: DIGIT) {
+	/*
+		The casts are required if _DIGIT_BITS is one less than
+		the number of bits in a DIGIT [e.g. _DIGIT_BITS==31].
+	*/
+	return DIGIT((1 << _DIGIT_BITS) - a.digit[0])
+}
+
+/*
+	Determines if a number is a valid DR modulus.
+	Assumes `a` to not be `nil` and to have been initialized.
+*/
+_private_dr_is_modulus :: proc(a: ^Int) -> (res: bool) {
+	/*
+		Must be at least two digits.
+	*/
+	if a.used < 2 { return false }
+
+	/*
+		Must be of the form b**k - a [a <= b] so all but the first digit must be equal to -1 (mod b).
+	*/
+	for ix := 1; ix < a.used; ix += 1 {
+		if a.digit[ix] != _MASK {
+			return false
+		}
+	}
+	return true
+}
+
+/*
+	Reduce "x" in place modulo "n" using the Diminished Radix algorithm.
+	Based on algorithm from the paper
+
+		"Generating Efficient Primes for Discrete Log Cryptosystems"
+					Chae Hoon Lim, Pil Joong Lee,
+			POSTECH Information Research Laboratories
+
+	The modulus must be of a special format [see manual].
+	Has been modified to use algorithm 7.10 from the LTM book instead
+
+	Input x must be in the range 0 <= x <= (n-1)**2
+	Assumes `x` and `n` to not be `nil` and to have been initialized.
+*/
+_private_int_dr_reduce :: proc(x, n: ^Int, k: DIGIT, allocator := context.allocator) -> (err: Error) {
+	/*
+		m = digits in modulus.
+	*/
+	m := n.used
+
+	/*
+		Ensure that "x" has at least 2m digits.
+	*/
+	internal_grow(x, m + m)                                          or_return
+
+	/*
+		Top of loop, this is where the code resumes if another reduction pass is required.
+	*/
+	for {
+		i: int
+		mu := DIGIT(0)
+
+		/*
+			Compute (x mod B**m) + k * [x/B**m] inline and inplace.
+		*/
+		for i = 0; i < m; i += 1 {
+			r         := _WORD(x.digit[i + m]) * _WORD(k) + _WORD(x.digit[i] + mu)
+			x.digit[i] = DIGIT(r & _WORD(_MASK))
+			mu         = DIGIT(r >> _WORD(_DIGIT_BITS))
+		}
+
+		/*
+			Set final carry.
+		*/
+		x.digit[i] = mu
+
+		/*
+			Zero words above m.
+		*/
+		mem.zero_slice(x.digit[m + 1:][:x.used - m])
+
+		/*
+			Clamp, sub and return.
+		*/
+		internal_clamp(x)                                            or_return
+
+		/*
+			If x >= n then subtract and reduce again.
+			Each successive "recursion" makes the input smaller and smaller.
+		*/
+		if internal_lt_abs(x, n) { break }
+
+		internal_sub(x, x, n)                                        or_return
+	}
+	return nil
+}
+
+/*
+	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.
+*/
+_private_int_exponent_mod :: proc(res, G, X, P: ^Int, redmode: int, allocator := context.allocator) -> (err: Error) {
+	context.allocator = allocator
+
+	M := [_TAB_SIZE]Int{}
+	winsize: uint
+
+	/*
+		Use a pointer to the reduction algorithm.
+		This allows us to use one of many reduction algorithms without modding the guts of the code with if statements everywhere.
+	*/
+	redux: #type proc(x, m, mu: ^Int, allocator := context.allocator) -> (err: Error)
+
+	defer {
+		internal_destroy(&M[1])
+		for x := 1 << (winsize - 1); x < (1 << winsize); x += 1 {
+			internal_destroy(&M[x])
+		}
+	}
+
+	/*
+		Find window size.
+	*/
+	x := internal_count_bits(X)
+	switch {
+	case x <= 7:
+		winsize = 2
+	case x <= 36:
+		winsize = 3
+	case x <= 140:
+		winsize = 4
+	case x <= 450:
+		winsize = 5
+	case x <= 1303:
+		winsize = 6
+	case x <= 3529:
+		winsize = 7
+	case:
+		winsize = 8
+	}
+
+	winsize = min(_MAX_WIN_SIZE, winsize) if _MAX_WIN_SIZE > 0 else winsize
+
+	/*
+		Init M array.
+		Init first cell.
+	*/
+	internal_zero(&M[1])                                             or_return
+
+	/*
+		Now init the second half of the array.
+	*/
+	for x = 1 << (winsize - 1); x < (1 << winsize); x += 1 {
+		internal_zero(&M[x])                                         or_return
+	}
+
+	/*
+		Create `mu`, used for Barrett reduction.
+	*/
+	mu := &Int{}
+	defer internal_destroy(mu)
+	internal_zero(mu)                                                or_return
+
+	if redmode == 0 {
+		_private_int_reduce_setup(mu, P)                             or_return
+		redux = _private_int_reduce
+	} else {
+		_private_int_reduce_2k_setup_l(mu, P)                        or_return
+		redux = _private_int_reduce_2k_l
+	}
+
+	/*
+		Create M table.
+
+		The M table contains powers of the base, e.g. M[x] = G**x mod P.
+		The first half of the table is not computed, though, except for M[0] and M[1].
+	*/
+	internal_int_mod(&M[1], G, P)                                    or_return
+
+	/*
+		Compute the value at M[1<<(winsize-1)] by squaring M[1] (winsize-1) times.
+
+		TODO: This can probably be replaced by computing the power and using `pow` to raise to it
+		instead of repeated squaring.
+	*/
+	slot := 1 << (winsize - 1)
+	internal_copy(&M[slot], &M[1])                                   or_return
+
+	for x = 0; x < int(winsize - 1); x += 1 {
+		/*
+			Square it.
+		*/
+		internal_sqr(&M[slot], &M[slot])                             or_return
+
+		/*
+			Reduce modulo P
+		*/
+		redux(&M[slot], P, mu)                                       or_return
+	}
+
+	/*
+		Create upper table, that is M[x] = M[x-1] * M[1] (mod P)
+		for x = (2**(winsize - 1) + 1) to (2**winsize - 1)
+	*/
+	for x = slot + 1; x < (1 << winsize); x += 1 {
+		internal_mul(&M[x], &M[x - 1], &M[1])                        or_return
+		redux(&M[x], P, mu)                                          or_return
+	}
+
+	/*
+		Setup result.
+	*/
+	internal_one(res)                                                or_return
+
+	/*
+		Set initial mode and bit cnt.
+	*/
+	mode   := 0
+	bitcnt := 1
+	buf    := DIGIT(0)
+	digidx := X.used - 1
+	bitcpy := uint(0)
+	bitbuf := DIGIT(0)
+
+	for {
+		/*
+			Grab next digit as required.
+		*/
+		bitcnt -= 1
+		if bitcnt == 0 {
+			/*
+				If digidx == -1 we are out of digits.
+			*/
+			if digidx == -1 { break }
+
+			/*
+				Read next digit and reset the bitcnt.
+			*/
+			buf    = X.digit[digidx]
+			digidx -= 1
+			bitcnt = _DIGIT_BITS
+		}
+
+		/*
+			Grab the next msb from the exponent.
+		*/
+		y := buf >> (_DIGIT_BITS - 1) & 1
+		buf <<= 1
+
+		/*
+			If the bit is zero and mode == 0 then we ignore it.
+			These represent the leading zero bits before the first 1 bit
+			in the exponent.  Technically this opt is not required but it
+			does lower the # of trivial squaring/reductions used.
+		*/
+		if mode == 0 && y == 0 {
+			continue
+		}
+
+		/*
+			If the bit is zero and mode == 1 then we square.
+		*/
+		if mode == 1 && y == 0 {
+			internal_sqr(res, res)                                   or_return
+			redux(res, P, mu)                                        or_return
+			continue
+		}
+
+		/*
+			Else we add it to the window.
+		*/
+		bitcpy += 1
+		bitbuf |= (y << (winsize - bitcpy))
+		mode    = 2
+
+		if (bitcpy == winsize) {
+			/*
+				Window is filled so square as required and multiply.
+				Square first.
+			*/
+			for x = 0; x < int(winsize); x += 1 {
+				internal_sqr(res, res)                               or_return
+				redux(res, P, mu)                                    or_return
+			}
+
+			/*
+				Then multiply.
+			*/
+			internal_mul(res, res, &M[bitbuf])                       or_return
+			redux(res, P, mu)                                        or_return
+
+			/*
+				Empty window and reset.
+			*/
+			bitcpy = 0
+			bitbuf = 0
+			mode   = 1
+		}
+	}
+
+	/*
+		If bits remain then square/multiply.
+	*/
+	if mode == 2 && bitcpy > 0 {
+		/*
+			Square then multiply if the bit is set.
+		*/
+		for x = 0; x < int(bitcpy); x += 1 {
+			internal_sqr(res, res)                                   or_return
+			redux(res, P, mu)                                        or_return
+
+			bitbuf <<= 1
+			if ((bitbuf & (1 << winsize)) != 0) {
+				/*
+					Then multiply.
+				*/
+				internal_mul(res, res, &M[1])                        or_return
+				redux(res, P, mu)                                    or_return
+			}
+		}
+	}
+	return err
+}
+
+/*
+	Computes Y == G**X mod P, HAC pp.616, Algorithm 14.85
+
+	Uses a left-to-right `k`-ary sliding window to compute the modular exponentiation.
+	The value of `k` changes based on the size of the exponent.
+
+	Uses Montgomery or Diminished Radix reduction [whichever appropriate]
+
+	Assumes `res`, `G`, `X` and `P` to not be `nil` and for `G`, `X` and `P` to have been initialized.
+*/
+_private_int_exponent_mod_fast :: proc(res, G, X, P: ^Int, redmode: int, allocator := context.allocator) -> (err: Error) {
+	context.allocator = allocator
+
+	M := [_TAB_SIZE]Int{}
+	winsize: uint
+
+	/*
+		Use a pointer to the reduction algorithm.
+		This allows us to use one of many reduction algorithms without modding the guts of the code with if statements everywhere.
+	*/
+	redux: #type proc(x, n: ^Int, rho: DIGIT, allocator := context.allocator) -> (err: Error)
+
+	defer {
+		internal_destroy(&M[1])
+		for x := 1 << (winsize - 1); x < (1 << winsize); x += 1 {
+			internal_destroy(&M[x])
+		}
+	}
+
+	/*
+		Find window size.
+	*/
+	x := internal_count_bits(X)
+	switch {
+	case x <= 7:
+		winsize = 2
+	case x <= 36:
+		winsize = 3
+	case x <= 140:
+		winsize = 4
+	case x <= 450:
+		winsize = 5
+	case x <= 1303:
+		winsize = 6
+	case x <= 3529:
+		winsize = 7
+	case:
+		winsize = 8
+	}
+
+	winsize = min(_MAX_WIN_SIZE, winsize) if _MAX_WIN_SIZE > 0 else winsize
+
+	/*
+		Init M array
+		Init first cell.
+	*/
+	cap := internal_int_allocated_cap(P)
+	internal_grow(&M[1], cap)                                        or_return
+
+	/*
+		Now init the second half of the array.
+	*/
+	for x = 1 << (winsize - 1); x < (1 << winsize); x += 1 {
+		internal_grow(&M[x], cap)                                    or_return
+	}
+
+	/*
+		Determine and setup reduction code.
+	*/
+	rho: DIGIT
+
+	if redmode == 0 {
+		/*
+			Now setup Montgomery.
+		*/
+		rho = _private_int_montgomery_setup(P)                       or_return
+
+		/*
+			Automatically pick the comba one if available (saves quite a few calls/ifs).
+		*/
+		if ((P.used * 2) + 1) < _WARRAY && P.used < _MAX_COMBA {
+			redux = _private_montgomery_reduce_comba
+		} else {
+			/*
+				Use slower baseline Montgomery method.
+			*/
+			redux = _private_int_montgomery_reduce
+		}
+	} else if redmode == 1 {
+		/*
+			Setup DR reduction for moduli of the form B**k - b.
+		*/
+		rho = _private_int_dr_setup(P)
+		redux = _private_int_dr_reduce
+	} else {
+		/*
+			Setup DR reduction for moduli of the form 2**k - b.
+		*/
+		rho = _private_int_reduce_2k_setup(P)                        or_return
+		redux = _private_int_reduce_2k
+	}
+
+	/*
+		Setup result.
+	*/
+	internal_grow(res, cap)                                          or_return
+
+	/*
+		Create M table
+		The first half of the table is not computed, though, except for M[0] and M[1]
+	*/
+
+	if redmode == 0 {
+		/*
+			Now we need R mod m.
+		*/
+		_private_int_montgomery_calc_normalization(res, P)           or_return
+
+		/*
+			Now set M[1] to G * R mod m.
+		*/
+		internal_mulmod(&M[1], G, res, P)                            or_return
+	} else {
+		internal_one(res)                                            or_return
+		internal_mod(&M[1], G, P)                                    or_return
+	}
+
+	/*
+		Compute the value at M[1<<(winsize-1)] by squaring M[1] (winsize-1) times.
+	*/
+	slot := 1 << (winsize - 1)
+	internal_copy(&M[slot], &M[1])                                   or_return
+
+	for x = 0; x < int(winsize - 1); x += 1 {
+		internal_sqr(&M[slot], &M[slot])                             or_return
+		redux(&M[slot], P, rho)                                      or_return
+	}
+
+	/*
+		Create upper table.
+	*/
+	for x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x += 1 {
+		internal_mul(&M[x], &M[x - 1], &M[1])                        or_return
+		redux(&M[x], P, rho)                                         or_return
+	}
+
+	/*
+		Set initial mode and bit cnt.
+	*/
+	mode   := 0
+	bitcnt := 1
+	buf    := DIGIT(0)
+	digidx := X.used - 1
+	bitcpy := 0
+	bitbuf := DIGIT(0)
+
+	for {
+		/*
+			Grab next digit as required.
+		*/
+		bitcnt -= 1
+		if bitcnt == 0 {
+			/*
+				If digidx == -1 we are out of digits so break.
+			*/
+			if digidx == -1 { break }
+
+			/*
+				Read next digit and reset the bitcnt.
+			*/
+			buf    = X.digit[digidx]
+			digidx -= 1
+			bitcnt = _DIGIT_BITS
+		}
+
+		/*
+			Grab the next msb from the exponent.
+		*/
+		y := (buf >> (_DIGIT_BITS - 1)) & 1
+		buf <<= 1
+
+		/*
+			If the bit is zero and mode == 0 then we ignore it.
+			These represent the leading zero bits before the first 1 bit in the exponent.
+			Technically this opt is not required but it does lower the # of trivial squaring/reductions used.
+		*/
+		if mode == 0 && y == 0 { continue }
+
+		/*
+			If the bit is zero and mode == 1 then we square.
+		*/
+		if mode == 1 && y == 0 {
+			internal_sqr(res, res)                                   or_return
+			redux(res, P, rho)                                       or_return
+			continue
+		}
+
+		/*
+			Else we add it to the window.
+		*/
+		bitcpy += 1
+		bitbuf |= (y << (winsize - uint(bitcpy)))
+		mode    = 2
+
+		if bitcpy == int(winsize) {
+			/*
+				Window is filled so square as required and multiply
+				Square first.
+			*/
+			for x = 0; x < int(winsize); x += 1 {
+				internal_sqr(res, res)                               or_return
+				redux(res, P, rho)                                   or_return
+			}
+
+			/*
+				Then multiply.
+			*/
+			internal_mul(res, res, &M[bitbuf])                       or_return
+			redux(res, P, rho)                                       or_return
+
+			/*
+				Empty window and reset.
+			*/
+			bitcpy = 0
+			bitbuf = 0
+			mode   = 1
+		}
+	}
+
+	/*
+		If bits remain then square/multiply.
+	*/
+	if mode == 2 && bitcpy > 0 {
+		/*
+			Square then multiply if the bit is set.
+		*/
+		for x = 0; x < bitcpy; x += 1 {
+			internal_sqr(res, res)                                   or_return
+			redux(res, P, rho)                                       or_return
+
+			/*
+				Get next bit of the window.
+			*/
+			bitbuf <<= 1
+			if bitbuf & (1 << winsize) != 0 {
+				/*
+					Then multiply.
+				*/
+				internal_mul(res, res, &M[1])                        or_return
+				redux(res, P, rho)                                   or_return
+			}
+		}
+	}
+
+	if redmode == 0 {
+		/*
+			Fixup result if Montgomery reduction is used.
+			Recall that any value in a Montgomery system is actually multiplied by R mod n.
+			So we have to reduce one more time to cancel out the factor of R.
+		*/
+		redux(res, P, rho)                                           or_return
+	}
+
+	return nil
+}
+
+/*
+	hac 14.61, pp608
+*/
+_private_inverse_modulo :: proc(dest, a, b: ^Int, allocator := context.allocator) -> (err: Error) {
+	context.allocator = allocator
+	x, y, u, v, A, B, C, D := &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}
+	defer internal_destroy(x, y, u, v, A, B, C, D)
+
+	// `b` cannot be negative.
+	if b.sign == .Negative || internal_is_zero(b) {
+		return .Invalid_Argument
+	}
+
+	// init temps.
+	internal_init_multi(x, y, u, v, A, B, C, D) or_return
+
+	// `x` = `a` % `b`, `y` = `b`
+	internal_mod(x, a, b) or_return
+	internal_copy(y, b) or_return
+
+	// 2. [modified] if x,y are both even then return an error!
+	if internal_is_even(x) && internal_is_even(y) {
+		return .Invalid_Argument
+	}
+
+	// 3. u=x, v=y, A=1, B=0, C=0, D=1
+	internal_copy(u, x) or_return
+	internal_copy(v, y) or_return
+	internal_one(A) or_return
+	internal_one(D) or_return
+
+	for {
+		// 4.  while `u` is even do:
+		for internal_is_even(u) {
+			// 4.1 `u` = `u` / 2
+			internal_int_shr1(u, u) or_return
+
+			// 4.2 if `A` or `B` is odd then:
+			if internal_is_odd(A) || internal_is_odd(B) {
+				// `A` = (`A`+`y`) / 2, `B` = (`B`-`x`) / 2
+				internal_add(A, A, y) or_return
+				internal_sub(B, B, x) or_return
+			}
+			// `A` = `A` / 2, `B` = `B` / 2
+			internal_int_shr1(A, A) or_return
+			internal_int_shr1(B, B) or_return
+		}
+
+		// 5.  while `v` is even do:
+		for internal_is_even(v) {
+			// 5.1 `v` = `v` / 2
+			internal_int_shr1(v, v) or_return
+
+			// 5.2 if `C` or `D` is odd then:
+			if internal_is_odd(C) || internal_is_odd(D) {
+				// `C` = (`C`+`y`) / 2, `D` = (`D`-`x`) / 2
+				internal_add(C, C, y) or_return
+				internal_sub(D, D, x) or_return
+			}
+			// `C` = `C` / 2, `D` = `D` / 2
+			internal_int_shr1(C, C) or_return
+			internal_int_shr1(D, D) or_return
+		}
+
+		// 6.  if `u` >= `v` then:
+		if internal_cmp(u, v) != -1 {
+			// `u` = `u` - `v`, `A` = `A` - `C`, `B` = `B` - `D`
+			internal_sub(u, u, v) or_return
+			internal_sub(A, A, C) or_return
+			internal_sub(B, B, D) or_return
+		} else {
+			// v - v - u, C = C - A, D = D - B
+			internal_sub(v, v, u) or_return
+			internal_sub(C, C, A) or_return
+			internal_sub(D, D, B) or_return
+		}
+
+		// If not zero goto step 4
+		if internal_is_zero(u) {
+			break
+		}
+	}
+
+	// Now `a` = `C`, `b` = `D`, `gcd` == `g`*`v`
+
+	// If `v` != `1` then there is no inverse.
+	if !internal_eq(v, 1) {
+		return .Invalid_Argument
+	}
+
+	// If its too low.
+	for internal_is_negative(C) {
+		internal_add(C, C, b) or_return
+	}
+
+	// Too big.
+	for internal_cmp_mag(C, b) > -1 {
+		internal_sub(C, C, b) or_return
+	}
+
+	// `C` is now the inverse.
+	swap(dest, C)
+	return
+}
+
+/*
+	Computes the modular inverse via binary extended Euclidean algorithm, that is `dest` = 1 / `a` mod `b`.
+
+	Based on slow invmod except this is optimized for the case where `b` is odd,
+	as per HAC Note 14.64 on pp. 610.
+*/
+_private_inverse_modulo_odd :: proc(dest, a, b: ^Int, allocator := context.allocator) -> (err: Error) {
+	context.allocator = allocator
+	x, y, u, v, B, D := &Int{}, &Int{}, &Int{}, &Int{}, &Int{}, &Int{}
+	defer internal_destroy(x, y, u, v, B, D)
+
+	sign: Sign
+
+	/*
+		2. [modified] `b` must be odd.
+	*/
+	if internal_is_even(b) { return .Invalid_Argument }
+
+	/*
+		Init all our temps.
+	*/
+	internal_init_multi(x, y, u, v, B, D) or_return
+
+	/*
+		`x` == modulus, `y` == value to invert.
+	*/
+	internal_copy(x, b) or_return
+
+	/*
+		We need `y` = `|a|`.
+	*/
+	internal_mod(y, a, b) or_return
+
+	/*
+		If one of `x`, `y` is zero return an error!
+	*/
+	if internal_is_zero(x) || internal_is_zero(y) { return .Invalid_Argument }
+
+	/*
+		3. `u` = `x`, `v` = `y`, `A` = 1, `B` = 0, `C` = 0, `D` = 1
+	*/
+	internal_copy(u, x) or_return
+	internal_copy(v, y) or_return
+
+	internal_one(D) or_return
+
+	for {
+		/*
+			4.  while `u` is even do.
+		*/
+		for internal_is_even(u) {
+			/*
+				4.1 `u` = `u` / 2
+			*/
+			internal_int_shr1(u, u) or_return
+
+			/*
+				4.2 if `B` is odd then:
+			*/
+			if internal_is_odd(B) {
+				/*
+					`B` = (`B` - `x`) / 2
+				*/
+				internal_sub(B, B, x) or_return
+			}
+
+			/*
+				`B` = `B` / 2
+			*/
+			internal_int_shr1(B, B) or_return
+		}
+
+		/*
+			5.  while `v` is even do:
+		*/
+		for internal_is_even(v) {
+			/*
+				5.1 `v` = `v` / 2
+			*/
+			internal_int_shr1(v, v) or_return
+
+			/*
+				5.2 if `D` is odd then:
+			*/
+			if internal_is_odd(D) {
+				/*
+					`D` = (`D` - `x`) / 2
+				*/
+				internal_sub(D, D, x) or_return
+			}
+			/*
+				`D` = `D` / 2
+			*/
+			internal_int_shr1(D, D) or_return
+		}
+
+		/*
+			6.  if `u` >= `v` then:
+		*/
+		if internal_cmp(u, v) != -1 {
+			/*
+				`u` = `u` - `v`, `B` = `B` - `D`
+			*/
+			internal_sub(u, u, v) or_return
+			internal_sub(B, B, D) or_return
+		} else {
+			/*
+				`v` - `v` - `u`, `D` = `D` - `B`
+			*/
+			internal_sub(v, v, u) or_return
+			internal_sub(D, D, B) or_return
+		}
+
+		/*
+			If not zero goto step 4.
+		*/
+		if internal_is_zero(u) { break }
+	}
+
+	/*
+		Now `a` = C, `b` = D, gcd == g*v
+	*/
+
+	/*
+		if `v` != 1 then there is no inverse
+	*/
+	if internal_cmp(v, 1) != 0 {
+		return .Invalid_Argument
+	}
+
+	/*
+		`b` is now the inverse.
+	*/
+	sign = a.sign
+	for internal_int_is_negative(D) {
+		internal_add(D, D, b) or_return
+	}
+
+	/*
+		Too big.
+	*/
+	for internal_gte_abs(D, b) {
+		internal_sub(D, D, b) or_return
+	}
+
+	swap(dest, D)
+	dest.sign = sign
+	return nil
+}
+
+
+/*
+	Returns the log2 of an `Int`.
+	Assumes `a` not to be `nil` and to have been initialized.
+	Also assumes `base` is a power of two.
+*/
+_private_log_power_of_two :: proc(a: ^Int, base: DIGIT) -> (log: int, err: Error) {
+	base := base
+	y: int
+	for y = 0; base & 1 == 0; {
+		y += 1
+		base >>= 1
+	}
+	log = internal_count_bits(a)
+	return (log - 1) / y, err
+}
+
+/*
+	Copies DIGITs from `src` to `dest`.
+	Assumes `src` and `dest` to not be `nil` and have been initialized.
+*/
+_private_copy_digits :: proc(dest, src: ^Int, digits: int, offset := int(0)) -> (err: Error) {
+	digits := digits
+	/*
+		If dest == src, do nothing
+	*/
+	if dest == src {
+		return nil
+	}
+
+	digits = min(digits, len(src.digit), len(dest.digit))
+	mem.copy_non_overlapping(&dest.digit[0], &src.digit[offset], size_of(DIGIT) * digits)
+	return nil
+}
+
+
+/*
+	Shift left by `digits` * _DIGIT_BITS bits.
+*/
+_private_int_shl_leg :: proc(quotient: ^Int, digits: int, allocator := context.allocator) -> (err: Error) {
+	context.allocator = allocator
+
+	if digits <= 0 { return nil }
+
+	/*
+		No need to shift a zero.
+	*/
+	if #force_inline internal_is_zero(quotient) {
+		return nil
+	}
+
+	/*
+		Resize `quotient` to accomodate extra digits.
+	*/
+	#force_inline internal_grow(quotient, quotient.used + digits) or_return
+
+	/*
+		Increment the used by the shift amount then copy upwards.
+	*/
+
+	/*
+		Much like `_private_int_shr_leg`, this is implemented using a sliding window,
+		except the window goes the other way around.
+	*/
+	#no_bounds_check for x := quotient.used; x > 0; x -= 1 {
+		quotient.digit[x+digits-1] = quotient.digit[x-1]
+	}
+
+	quotient.used += digits
+	mem.zero_slice(quotient.digit[:digits])
+	return nil
+}
+
+/*
+	Shift right by `digits` * _DIGIT_BITS bits.
+*/
+_private_int_shr_leg :: proc(quotient: ^Int, digits: int, allocator := context.allocator) -> (err: Error) {
+	context.allocator = allocator
+
+	if digits <= 0 { return nil }
+
+	/*
+		If digits > used simply zero and return.
+	*/
+	if digits > quotient.used { return internal_zero(quotient) }
+
+	/*
+		Much like `int_shl_digit`, this is implemented using a sliding window,
+		except the window goes the other way around.
+
+		b-2 | b-1 | b0 | b1 | b2 | ... | bb |   ---->
+					/\                   |      ---->
+					 \-------------------/      ---->
+	*/
+
+	#no_bounds_check for x := 0; x < (quotient.used - digits); x += 1 {
+		quotient.digit[x] = quotient.digit[x + digits]
+	}
+	quotient.used -= digits
+	internal_zero_unused(quotient)
+	return internal_clamp(quotient)
+}
+
+/*	
+	========================    End of private procedures    =======================
+
+	===============================  Private tables  ===============================
+
+	Tables used by `internal_*` and `_*`.
+*/
+
+_private_int_rem_128 := [?]DIGIT{
+	0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
+	0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
+	1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
+	1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
+	0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
+	1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
+	1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
+	1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1,
+}
+#assert(128 * size_of(DIGIT) == size_of(_private_int_rem_128))
+
+_private_int_rem_105 := [?]DIGIT{
+	0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1,
+	0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1,
+	0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1,
+	1, 0, 1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1,
+	0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1,
+	1, 1, 1, 1, 0, 1, 0, 1, 1, 0, 0, 1, 1, 1, 1,
+	1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1,
+}
+#assert(105 * size_of(DIGIT) == size_of(_private_int_rem_105))
+
+_PRIME_TAB_SIZE :: 256
+_private_prime_table := [_PRIME_TAB_SIZE]DIGIT{
+	0x0002, 0x0003, 0x0005, 0x0007, 0x000B, 0x000D, 0x0011, 0x0013,
+	0x0017, 0x001D, 0x001F, 0x0025, 0x0029, 0x002B, 0x002F, 0x0035,
+	0x003B, 0x003D, 0x0043, 0x0047, 0x0049, 0x004F, 0x0053, 0x0059,
+	0x0061, 0x0065, 0x0067, 0x006B, 0x006D, 0x0071, 0x007F, 0x0083,
+	0x0089, 0x008B, 0x0095, 0x0097, 0x009D, 0x00A3, 0x00A7, 0x00AD,
+	0x00B3, 0x00B5, 0x00BF, 0x00C1, 0x00C5, 0x00C7, 0x00D3, 0x00DF,
+	0x00E3, 0x00E5, 0x00E9, 0x00EF, 0x00F1, 0x00FB, 0x0101, 0x0107,
+	0x010D, 0x010F, 0x0115, 0x0119, 0x011B, 0x0125, 0x0133, 0x0137,
+
+	0x0139, 0x013D, 0x014B, 0x0151, 0x015B, 0x015D, 0x0161, 0x0167,
+	0x016F, 0x0175, 0x017B, 0x017F, 0x0185, 0x018D, 0x0191, 0x0199,
+	0x01A3, 0x01A5, 0x01AF, 0x01B1, 0x01B7, 0x01BB, 0x01C1, 0x01C9,
+	0x01CD, 0x01CF, 0x01D3, 0x01DF, 0x01E7, 0x01EB, 0x01F3, 0x01F7,
+	0x01FD, 0x0209, 0x020B, 0x021D, 0x0223, 0x022D, 0x0233, 0x0239,
+	0x023B, 0x0241, 0x024B, 0x0251, 0x0257, 0x0259, 0x025F, 0x0265,
+	0x0269, 0x026B, 0x0277, 0x0281, 0x0283, 0x0287, 0x028D, 0x0293,
+	0x0295, 0x02A1, 0x02A5, 0x02AB, 0x02B3, 0x02BD, 0x02C5, 0x02CF,
+
+	0x02D7, 0x02DD, 0x02E3, 0x02E7, 0x02EF, 0x02F5, 0x02F9, 0x0301,
+	0x0305, 0x0313, 0x031D, 0x0329, 0x032B, 0x0335, 0x0337, 0x033B,
+	0x033D, 0x0347, 0x0355, 0x0359, 0x035B, 0x035F, 0x036D, 0x0371,
+	0x0373, 0x0377, 0x038B, 0x038F, 0x0397, 0x03A1, 0x03A9, 0x03AD,
+	0x03B3, 0x03B9, 0x03C7, 0x03CB, 0x03D1, 0x03D7, 0x03DF, 0x03E5,
+	0x03F1, 0x03F5, 0x03FB, 0x03FD, 0x0407, 0x0409, 0x040F, 0x0419,
+	0x041B, 0x0425, 0x0427, 0x042D, 0x043F, 0x0443, 0x0445, 0x0449,
+	0x044F, 0x0455, 0x045D, 0x0463, 0x0469, 0x047F, 0x0481, 0x048B,
+
+	0x0493, 0x049D, 0x04A3, 0x04A9, 0x04B1, 0x04BD, 0x04C1, 0x04C7,
+	0x04CD, 0x04CF, 0x04D5, 0x04E1, 0x04EB, 0x04FD, 0x04FF, 0x0503,
+	0x0509, 0x050B, 0x0511, 0x0515, 0x0517, 0x051B, 0x0527, 0x0529,
+	0x052F, 0x0551, 0x0557, 0x055D, 0x0565, 0x0577, 0x0581, 0x058F,
+	0x0593, 0x0595, 0x0599, 0x059F, 0x05A7, 0x05AB, 0x05AD, 0x05B3,
+	0x05BF, 0x05C9, 0x05CB, 0x05CF, 0x05D1, 0x05D5, 0x05DB, 0x05E7,
+	0x05F3, 0x05FB, 0x0607, 0x060D, 0x0611, 0x0617, 0x061F, 0x0623,
+	0x062B, 0x062F, 0x063D, 0x0641, 0x0647, 0x0649, 0x064D, 0x0653,
+}
+#assert(_PRIME_TAB_SIZE * size_of(DIGIT) == size_of(_private_prime_table))
+
+when MATH_BIG_FORCE_64_BIT || (!MATH_BIG_FORCE_32_BIT && size_of(rawptr) == 8) {
+	_factorial_table := [35]_WORD{
+/* f(00): */                                                     1,
+/* f(01): */                                                     1,
+/* f(02): */                                                     2,
+/* f(03): */                                                     6,
+/* f(04): */                                                    24,
+/* f(05): */                                                   120,
+/* f(06): */                                                   720,
+/* f(07): */                                                 5_040,
+/* f(08): */                                                40_320,
+/* f(09): */                                               362_880,
+/* f(10): */                                             3_628_800,
+/* f(11): */                                            39_916_800,
+/* f(12): */                                           479_001_600,
+/* f(13): */                                         6_227_020_800,
+/* f(14): */                                        87_178_291_200,
+/* f(15): */                                     1_307_674_368_000,
+/* f(16): */                                    20_922_789_888_000,
+/* f(17): */                                   355_687_428_096_000,
+/* f(18): */                                 6_402_373_705_728_000,
+/* f(19): */                               121_645_100_408_832_000,
+/* f(20): */                             2_432_902_008_176_640_000,
+/* f(21): */                            51_090_942_171_709_440_000,
+/* f(22): */                         1_124_000_727_777_607_680_000,
+/* f(23): */                        25_852_016_738_884_976_640_000,
+/* f(24): */                       620_448_401_733_239_439_360_000,
+/* f(25): */                    15_511_210_043_330_985_984_000_000,
+/* f(26): */                   403_291_461_126_605_635_584_000_000,
+/* f(27): */                10_888_869_450_418_352_160_768_000_000,
+/* f(28): */               304_888_344_611_713_860_501_504_000_000,
+/* f(29): */             8_841_761_993_739_701_954_543_616_000_000,
+/* f(30): */           265_252_859_812_191_058_636_308_480_000_000,
+/* f(31): */         8_222_838_654_177_922_817_725_562_880_000_000,
+/* f(32): */       263_130_836_933_693_530_167_218_012_160_000_000,
+/* f(33): */     8_683_317_618_811_886_495_518_194_401_280_000_000,
+/* f(34): */   295_232_799_039_604_140_847_618_609_643_520_000_000,
+	}
+} else {
+	_factorial_table := [21]_WORD{
+/* f(00): */                                                     1,
+/* f(01): */                                                     1,
+/* f(02): */                                                     2,
+/* f(03): */                                                     6,
+/* f(04): */                                                    24,
+/* f(05): */                                                   120,
+/* f(06): */                                                   720,
+/* f(07): */                                                 5_040,
+/* f(08): */                                                40_320,
+/* f(09): */                                               362_880,
+/* f(10): */                                             3_628_800,
+/* f(11): */                                            39_916_800,
+/* f(12): */                                           479_001_600,
+/* f(13): */                                         6_227_020_800,
+/* f(14): */                                        87_178_291_200,
+/* f(15): */                                     1_307_674_368_000,
+/* f(16): */                                    20_922_789_888_000,
+/* f(17): */                                   355_687_428_096_000,
+/* f(18): */                                 6_402_373_705_728_000,
+/* f(19): */                               121_645_100_408_832_000,
+/* f(20): */                             2_432_902_008_176_640_000,
+	}
+}
+
+/*
+	=========================  End of private tables  ========================
 */
\ No newline at end of file