diff --git a/core/math/big/basic.odin b/core/math/big/basic.odin index b6c49a5fd..6d1ab23e6 100644 --- a/core/math/big/basic.odin +++ b/core/math/big/basic.odin @@ -297,36 +297,7 @@ int_gcd_lcm :: proc(res_gcd, res_lcm, a, b: ^Int) -> (err: Error) { if res_gcd == nil && res_lcm == nil { return nil; } if err = clear_if_uninitialized(res_gcd, res_lcm, a, b); err != nil { return err; } - az, _ := is_zero(a); bz, _ := is_zero(b); - if az && bz { - if res_gcd != nil { - if err = zero(res_gcd); err != nil { return err; } - } - if res_lcm != nil { - if err = zero(res_lcm); err != nil { return err; } - } - return nil; - } - else if az { - if res_gcd != nil { - if err = abs(res_gcd, b); err != nil { return err; } - } - if res_lcm != nil { - if err = zero(res_lcm); err != nil { return err; } - } - return nil; - } - else if bz { - if res_gcd != nil { - if err = abs(res_gcd, a); err != nil { return err; } - } - if res_lcm != nil { - if err = zero(res_lcm); err != nil { return err; } - } - return nil; - } - - return #force_inline _int_gcd_lcm(res_gcd, res_lcm, a, b); + return #force_inline internal_int_gcd_lcm(res_gcd, res_lcm, a, b); } gcd_lcm :: proc { int_gcd_lcm, }; @@ -346,171 +317,6 @@ int_lcm :: proc(res, a, b: ^Int) -> (err: Error) { } lcm :: proc { int_lcm, }; -/* - Internal function computing both GCD using the binary method, - and, if target isn't `nil`, also LCM. - Expects the arguments to have been initialized. -*/ -_int_gcd_lcm :: proc(res_gcd, res_lcm, a, b: ^Int) -> (err: Error) { - /* - 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. - */ - if res_gcd == nil && res_lcm == nil { return nil; } - - /* - We need a temporary because `res_gcd` is allowed to be `nil`. - */ - az, _ := is_zero(a); bz, _ := is_zero(b); - if az && bz { - /* - GCD(0, 0) and LCM(0, 0) are both 0. - */ - if res_gcd != nil { - if err = zero(res_gcd); err != nil { return err; } - } - if res_lcm != nil { - if err = zero(res_lcm); err != nil { return err; } - } - return nil; - } else if az { - /* - We can early out with GCD = B and LCM = 0 - */ - if res_gcd != nil { - if err = abs(res_gcd, b); err != nil { return err; } - } - if res_lcm != nil { - if err = zero(res_lcm); err != nil { return err; } - } - return nil; - } else if bz { - /* - We can early out with GCD = A and LCM = 0 - */ - if res_gcd != nil { - if err = abs(res_gcd, a); err != nil { return err; } - } - if res_lcm != nil { - if err = zero(res_lcm); err != nil { return err; } - } - return nil; - } - - temp_gcd_res := &Int{}; - defer 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 destroy(u, v); - if err = copy(u, a); err != nil { return err; } - if err = copy(v, b); err != nil { return err; } - - /* - 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, _ := count_lsb(u); - v_lsb, _ := count_lsb(v); - k := min(u_lsb, v_lsb); - - if k > 0 { - /* - Divide the power of two out. - */ - if err = shr(u, u, k); err != nil { return err; } - if err = shr(v, v, k); err != nil { return err; } - } - - /* - Divide any remaining factors of two out. - */ - if u_lsb != k { - if err = shr(u, u, u_lsb - k); err != nil { return err; } - } - if v_lsb != k { - if err = shr(v, v, v_lsb - k); err != nil { return err; } - } - - for v.used != 0 { - /* - Make sure `v` is the largest. - */ - if c, _ := cmp_mag(u, v); c == 1 { - /* - Swap `u` and `v` to make sure `v` is >= `u`. - */ - swap(u, v); - } - - /* - Subtract smallest from largest. - */ - if err = sub(v, v, u); err != nil { return err; } - - /* - Divide out all factors of two. - */ - b, _ := count_lsb(v); - if err = shr(v, v, b); err != nil { return err; } - } - - /* - Multiply by 2**k which we divided out at the beginning. - */ - if err = shl(temp_gcd_res, u, k); err != nil { return err; } - 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 { - 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 c, _ := cmp_mag(a, b); c == -1 { - /* - Store quotient in `t2` such that `t2 * b` is the LCM. - */ - if err = div(res_lcm, a, temp_gcd_res); err != nil { return err; } - err = mul(res_lcm, res_lcm, b); - } else { - /* - Store quotient in `t2` such that `t2 * a` is the LCM. - */ - if err = div(res_lcm, a, temp_gcd_res); err != nil { return err; } - err = mul(res_lcm, res_lcm, b); - } - - if res_gcd != nil { - swap(temp_gcd_res, res_gcd); - } - - /* - Fix the sign to positive and return. - */ - res_lcm.sign = .Zero_or_Positive; - return err; -} - /* remainder = numerator % (1 << bits) */ diff --git a/core/math/big/internal.odin b/core/math/big/internal.odin index a3e2548db..995f37000 100644 --- a/core/math/big/internal.odin +++ b/core/math/big/internal.odin @@ -680,7 +680,7 @@ internal_int_divmod :: proc(quotient, remainder, numerator, denominator: ^Int, a // err = _int_div_recursive(quotient, remainder, numerator, denominator); } else { when true { - err = _private_int_div_school(quotient, remainder, numerator, denominator); + err = #force_inline _private_int_div_school(quotient, remainder, numerator, denominator); } else { /* NOTE(Jeroen): We no longer need or use `_private_int_div_small`. @@ -864,6 +864,18 @@ internal_int_factorial :: proc(res: ^Int, n: int) -> (err: Error) { return nil; } +/* + Returns GCD, LCM or both. + + Assumes `a` and `b` to have been initialized. + `res_gcd` and `res_lcm` can be nil or ^Int depending on which results are desired. +*/ +internal_int_gcd_lcm :: proc(res_gcd, res_lcm, a, b: ^Int) -> (err: Error) { + if res_gcd == nil && res_lcm == nil { return nil; } + + return #force_inline _private_int_gcd_lcm(res_gcd, res_lcm, a, b); +} + internal_int_zero_unused :: #force_inline proc(dest: ^Int, old_used := -1) { /* @@ -1466,6 +1478,171 @@ _private_int_recursive_product :: proc(res: ^Int, start, stop: int, level := int return #force_inline set(res, 1); } +/* + 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) -> (err: Error) { + 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 { + if err = zero(res_gcd); err != nil { return err; } + } + if res_lcm != nil { + if err = zero(res_lcm); err != nil { return err; } + } + return nil; + } else if a.used == 0 { + /* + We can early out with GCD = B and LCM = 0 + */ + if res_gcd != nil { + if err = abs(res_gcd, b); err != nil { return err; } + } + if res_lcm != nil { + if err = zero(res_lcm); err != nil { return err; } + } + return nil; + } else if b.used == 0 { + /* + We can early out with GCD = A and LCM = 0 + */ + if res_gcd != nil { + if err = abs(res_gcd, a); err != nil { return err; } + } + if res_lcm != nil { + if err = zero(res_lcm); err != nil { return err; } + } + return nil; + } + + temp_gcd_res := &Int{}; + defer 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 destroy(u, v); + if err = copy(u, a); err != nil { return err; } + if err = copy(v, b); err != nil { return err; } + + /* + 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, _ := count_lsb(u); + v_lsb, _ := count_lsb(v); + k := min(u_lsb, v_lsb); + + if k > 0 { + /* + Divide the power of two out. + */ + if err = shr(u, u, k); err != nil { return err; } + if err = shr(v, v, k); err != nil { return err; } + } + + /* + Divide any remaining factors of two out. + */ + if u_lsb != k { + if err = shr(u, u, u_lsb - k); err != nil { return err; } + } + if v_lsb != k { + if err = shr(v, v, v_lsb - k); err != nil { return err; } + } + + for v.used != 0 { + /* + Make sure `v` is the largest. + */ + if c, _ := cmp_mag(u, v); c == 1 { + /* + Swap `u` and `v` to make sure `v` is >= `u`. + */ + swap(u, v); + } + + /* + Subtract smallest from largest. + */ + if err = sub(v, v, u); err != nil { return err; } + + /* + Divide out all factors of two. + */ + b, _ := count_lsb(v); + if err = shr(v, v, b); err != nil { return err; } + } + + /* + Multiply by 2**k which we divided out at the beginning. + */ + if err = shl(temp_gcd_res, u, k); err != nil { return err; } + 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 { + 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 c, _ := cmp_mag(a, b); c == -1 { + /* + Store quotient in `t2` such that `t2 * b` is the LCM. + */ + if err = div(res_lcm, a, temp_gcd_res); err != nil { return err; } + err = mul(res_lcm, res_lcm, b); + } else { + /* + Store quotient in `t2` such that `t2 * a` is the LCM. + */ + if err = div(res_lcm, a, temp_gcd_res); err != nil { return err; } + err = mul(res_lcm, res_lcm, b); + } + + if res_gcd != nil { + swap(temp_gcd_res, res_gcd); + } + + /* + Fix the sign to positive and return. + */ + res_lcm.sign = .Zero_or_Positive; + return err; +} + /* ======================== End of private procedures ======================= diff --git a/core/math/big/test.py b/core/math/big/test.py index a06e8e119..6bdff91db 100644 --- a/core/math/big/test.py +++ b/core/math/big/test.py @@ -446,7 +446,6 @@ TESTS = { test_factorial: [ [ 6_000 ], # Regular factorial, see cutoff in common.odin. [ 12_345 ], # Binary split factorial - [ 100_000 ], ], test_gcd: [ [ 23, 25, ], @@ -464,6 +463,12 @@ TESTS = { ], } +if not FAST_TESTS: + TESTS[test_factorial].append( + # This one on its own takes around 800ms, so we exclude it for FAST_TESTS + [ 100_000 ], + ) + total_passes = 0 total_failures = 0