Add faster divison.

3 files changed, 177 insertions, 8 deletions

diff --git a/core/math/big/basic.odin b/core/math/big/basic.odin
index 69520a76c..28e339f26 100644
--- a/core/math/big/basic.odin
+++ b/core/math/big/basic.odin
@@ -681,10 +681,9 @@ int_divmod :: proc(quotient, remainder, numerator, denominator: ^Int) -> (err: E
 

 	if false && (denominator.used > 2 * _MUL_KARATSUBA_CUTOFF) && (denominator.used <= (numerator.used/3) * 2) {

 		// err = _int_div_recursive(quotient, remainder, numerator, denominator);

-	} else if false {

-		// err = _int_div_school(quotient, remainder, numerator, denominator);

 	} else {

-		err = _int_div_small(quotient, remainder, numerator, denominator);

+		err = _int_div_school(quotient, remainder, numerator, denominator);

+		// err = _int_div_small(quotient, remainder, numerator, denominator);

 	}

 

 	return err;

@@ -1312,6 +1311,176 @@ _int_div_3 :: proc(quotient, numerator: ^Int) -> (remainder: int, err: Error) {
 }

 

 /*

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

+*/

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

+	if err = error_if_immutable(quotient, remainder); err != nil { return err; }

+	if err = clear_if_uninitialized(quotient, numerator, denominator); err != nil { return err; }

+

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

+	defer destroy(q, x, y, t1, t2);

+

+	if err = grow(q, numerator.used + 2); err != nil { return err; }

+	q.used = numerator.used + 2;

+

+	if err = init_multi(t1, t2);   err != nil { return err; }

+	if err = copy(x, numerator);   err != nil { return err; }

+	if err = copy(y, denominator); err != nil { return err; }

+

+	/*

+		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, _ := count_bits(y);

+	norm %= _DIGIT_BITS;

+

+	if norm < _DIGIT_BITS - 1 {

+		norm = (_DIGIT_BITS - 1) - norm;

+		if err = shl(x, x, norm); err != nil { return err; }

+		if err = shl(y, y, norm); err != nil { return err; }

+	} 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}

+	*/

+

+	if err = shl_digit(y, n - t); err != nil { return err; }

+

+	c, _ := cmp(x, y);

+	for c != -1 {

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

+		if err = sub(x, x, y); err != nil { return err; }

+		c, _ = cmp(x, y);

+	}

+

+	/*

+		Reset y by shifting it back down.

+	*/

+	shr_digit(y, n - t);

+

+	/*

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

+	*/

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

+		for {

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

+

+			/*

+				Find left hand.

+			*/

+			zero(t1);

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

+			t1.digit[1] = y.digit[t];

+			t1.used = 2;

+			if err = mul(t1, t1, q.digit[(i - t) - 1]); err != nil { return err; }

+

+			/*

+				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 t1_t2, _ := cmp_mag(t1, t2); t1_t2 != 1 {

+

+				break;

+			}

+			iter += 1; if iter > 100 { return .Max_Iterations_Reached; }

+		}

+

+		/*

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

+		*/

+		if err = int_mul_digit(t1, y, q.digit[(i - t) - 1]); err != nil { return err; }

+		if err = shl_digit(t1, (i - t) - 1);       err != nil { return err; }

+		if err = sub(x, x, t1); err != nil { return err; }

+

+		/*

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

+		*/

+		if x.sign == .Negative {

+			if err = copy(t1, y); err != nil { return err; }

+			if err = shl_digit(t1, (i - t) - 1); err != nil { return err; }

+			if err = add(x, x, t1); err != nil { return err; }

+

+			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 {

+		clamp(q);

+		swap(q, quotient);

+		quotient.sign = .Negative if neg else .Zero_or_Positive;

+	}

+

+	if remainder != nil {

+		if err = shr(x, x, norm); err != nil { return err; }

+		swap(x, remainder);

+	}

+

+	return nil;

+}

+

+/*

 	Slower bit-bang division... also smaller.

 */

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

diff --git a/core/math/big/build.bat b/core/math/big/build.bat
index 2c1edfcec..d454fef4d 100644
--- a/core/math/big/build.bat
+++ b/core/math/big/build.bat
@@ -1,10 +1,10 @@
 @echo off

-odin run . -vet

+:odin run . -vet

 : -o:size -no-bounds-check

 :odin build . -build-mode:shared -show-timings -o:minimal -use-separate-modules

 :odin build . -build-mode:shared -show-timings -o:size -use-separate-modules -no-bounds-check

 :odin build . -build-mode:shared -show-timings -o:size -use-separate-modules

-:odin build . -build-mode:shared -show-timings -o:speed -use-separate-modules -no-bounds-check

+odin build . -build-mode:shared -show-timings -o:speed -use-separate-modules -no-bounds-check

 :odin build . -build-mode:shared -show-timings -o:speed -use-separate-modules

 

-:python test.py
\ No newline at end of file
+python test.py
\ No newline at end of file
diff --git a/core/math/big/helpers.odin b/core/math/big/helpers.odin
index a326b960b..f0eeea96b 100644
--- a/core/math/big/helpers.odin
+++ b/core/math/big/helpers.odin
@@ -340,7 +340,7 @@ zero  :: clear;
 	Set the `Int` to 1 and optionally shrink it to the minimum backing size.
 */
 int_one :: proc(a: ^Int, minimize := false, allocator := context.allocator) -> (err: Error) {
-	return copy(a, ONE, minimize, allocator);
+	return set(a, 1);
 }
 one :: proc { int_one, };
 
@@ -348,7 +348,7 @@ one :: proc { int_one, };
 	Set the `Int` to -1 and optionally shrink it to the minimum backing size.
 */
 int_minus_one :: proc(a: ^Int, minimize := false, allocator := context.allocator) -> (err: Error) {
-	return copy(a, MINUS_ONE, minimize, allocator);
+	return set(a, -1);
 }
 minus_one :: proc { int_minus_one, };