From 47397a6a488ad9cae625e531e5236513b81660a0 Mon Sep 17 00:00:00 2001 From: Jeroen van Rijn Date: Wed, 4 Aug 2021 00:59:15 +0200 Subject: [PATCH] Add faster divison. --- core/math/big/basic.odin | 175 ++++++++++++++++++++++++++++++++++++- core/math/big/build.bat | 6 +- core/math/big/helpers.odin | 4 +- 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; @@ -1311,6 +1310,176 @@ _int_div_3 :: proc(quotient, numerator: ^Int) -> (remainder: int, err: Error) { 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. +*/ +_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. */ 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, };