From 31c94bd7f83fcb1c2f32e6b168fd6e04cccbcfde Mon Sep 17 00:00:00 2001 From: Jeroen van Rijn Date: Sat, 24 Jul 2021 13:40:36 +0200 Subject: [PATCH] big: Finish `log`, fix `sqr`. --- core/math/big/basic.odin | 7 +-- core/math/big/example.odin | 10 +--- core/math/big/helpers.odin | 15 ++++- core/math/big/log.odin | 114 ++++++++++++++++++++++++++++++------- 4 files changed, 113 insertions(+), 33 deletions(-) diff --git a/core/math/big/basic.odin b/core/math/big/basic.odin index 18186b38c..53733fdb0 100644 --- a/core/math/big/basic.odin +++ b/core/math/big/basic.odin @@ -898,7 +898,7 @@ _int_sqr :: proc(dest, src: ^Int) -> (err: Error) { /* Grow `t` to maximum needed size, or `_DEFAULT_DIGIT_COUNT`, whichever is bigger. */ - if err = grow(t, min((2 * pa) + 1, _DEFAULT_DIGIT_COUNT)); err != .None { return err; } + if err = grow(t, max((2 * pa) + 1, _DEFAULT_DIGIT_COUNT)); err != .None { return err; } t.used = (2 * pa) + 1; for ix = 0; ix < pa; ix += 1 { @@ -906,13 +906,12 @@ _int_sqr :: proc(dest, src: ^Int) -> (err: Error) { /* First calculate the digit at 2*ix; calculate double precision result. */ - r := _WORD(t.digit[ix+ix]) + _WORD(src.digit[ix] * src.digit[ix]); + 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. */ @@ -924,7 +923,7 @@ _int_sqr :: proc(dest, src: ^Int) -> (err: Error) { */ r = _WORD(src.digit[ix]) * _WORD(src.digit[iy]); - /* Now calculate the double precision result. Nte we use + /* 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); diff --git a/core/math/big/example.odin b/core/math/big/example.odin index a6ba77667..c0de1d7c3 100644 --- a/core/math/big/example.odin +++ b/core/math/big/example.odin @@ -57,14 +57,8 @@ demo :: proc() { a, b, c := &Int{}, &Int{}, &Int{}; defer destroy(a, b, c); - for base in -3..=3 { - for power in -3..=3 { - err = pow(a, base, power); - fmt.printf("err: %v | pow(%v, %v) = ", err, base, power); print("", a, 10); - } - } - - + err = set(a, 5125); + print("a", a); } main :: proc() { diff --git a/core/math/big/helpers.odin b/core/math/big/helpers.odin index 20193e061..1d756e601 100644 --- a/core/math/big/helpers.odin +++ b/core/math/big/helpers.odin @@ -75,7 +75,7 @@ int_copy :: proc(dest, src: ^Int, allocator := context.allocator) -> (err: Error /* Copy everything over and zero high digits. */ - for v, i in src.digit[:src.used+1] { + for v, i in src.digit[:src.used] { dest.digit[i] = v; } dest.used = src.used; @@ -533,6 +533,19 @@ clear_if_uninitialized :: proc(dest: ^Int, minimize := false) -> (err: Error) { return .None; } +/* + Allocates several `Int`s at once. +*/ +_int_init_multi :: proc(integers: ..^Int) -> (err: Error) { + integers := integers; + for a in &integers { + if err = clear(a); err != .None { return err; } + } + return .None; +} + +_init_multi :: proc { _int_init_multi, }; + _copy_digits :: proc(dest, src: ^Int, digits: int) -> (err: Error) { digits := digits; if err = clear_if_uninitialized(src); err != .None { return err; } diff --git a/core/math/big/log.odin b/core/math/big/log.odin index 632973ec8..3bd1f4f25 100644 --- a/core/math/big/log.odin +++ b/core/math/big/log.odin @@ -13,36 +13,22 @@ int_log :: proc(a: ^Int, base: DIGIT) -> (res: int, err: Error) { if base < 2 || DIGIT(base) > _DIGIT_MAX { return -1, .Invalid_Argument; } - - if err = clear_if_uninitialized(a); err != .None { - return -1, err; - } - if n, _ := is_neg(a); n { - return -1, .Invalid_Argument; - } - if z, _ := is_zero(a); z { - return -1, .Invalid_Argument; - } + if err = clear_if_uninitialized(a); err != .None { return -1, err; } + if n, _ := is_neg(a); n { return -1, .Invalid_Argument; } + if z, _ := is_zero(a); z { return -1, .Invalid_Argument; } /* Fast path for bases that are a power of two. */ - if is_power_of_two(int(base)) { - return _log_power_of_two(a, base); - } + if is_power_of_two(int(base)) { return _log_power_of_two(a, base); } /* Fast path for `Int`s that fit within a single `DIGIT`. */ - if a.used == 1 { - return log(a.digit[0], DIGIT(base)); - } + if a.used == 1 { return log(a.digit[0], DIGIT(base)); } - // if (MP_HAS(S_MP_LOG)) { - // return s_mp_log(a, (mp_digit)base, c); - // } + return _int_log(a, base); - return -1, .Unimplemented; } log :: proc { int_log, int_log_digit, }; @@ -222,4 +208,92 @@ int_log_digit :: proc(a: DIGIT, base: DIGIT) -> (log: int, err: Error) { } else { return low, .None; } +} + + +/* + Internal implementation of log. +*/ +_int_log :: proc(a: ^Int, base: DIGIT) -> (res: int, err: Error) { + bracket_low, bracket_high, bracket_mid, t, bi_base := &Int{}, &Int{}, &Int{}, &Int{}, &Int{}; + + cnt := 0; + + ic, _ := cmp(a, base); + if ic == -1 || ic == 0 { + return 1 if ic == 0 else 0, .None; + } + + if err = set(bi_base, base); err != .None { return -1, err; } + if err = _init_multi(bracket_mid, t); err != .None { return -1, err; } + if err = one(bracket_low); err != .None { return -1, err; } + if err = set(bracket_high, base); err != .None { return -1, err; } + + 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 bc, _ := cmp(bracket_high, a); bc != -1 { + break; + } + + low = high; + if err = copy(bracket_low, bracket_high); err != .None { + destroy(bracket_low, bracket_high, bracket_mid, t, bi_base); + return -1, err; + } + high <<= 1; + if err = sqr(bracket_high, bracket_high); err != .None { + destroy(bracket_low, bracket_high, bracket_mid, t, bi_base); + return -1, err; + } + + cnt += 1; + if cnt == 7 { + destroy(bracket_low, bracket_high, bracket_mid, t, bi_base); + return -2, .Max_Iterations_Reached; + } + } + + for (high - low) > 1 { + mid := (high + low) >> 1; + + if err = pow(t, bi_base, mid - low); err != .None { + destroy(bracket_low, bracket_high, bracket_mid, t, bi_base); + return -1, err; + } + + if err = mul(bracket_mid, bracket_low, t); err != .None { + destroy(bracket_low, bracket_high, bracket_mid, t, bi_base); + return -1, err; + } + mc, _ := cmp(a, bracket_mid); + if mc == -1 { + high = mid; + swap(bracket_mid, bracket_high); + } + if mc == 1 { + low = mid; + swap(bracket_mid, bracket_low); + } + if mc == 0 { + destroy(bracket_low, bracket_high, bracket_mid, t, bi_base); + return mid, .None; + } + } + + fc, _ := cmp(bracket_high, a); + res = high if fc == 0 else low; + + destroy(bracket_low, bracket_high, bracket_mid, t, bi_base); + return; } \ No newline at end of file