ed/Odin
1
0
Fork 0
mirror of https://github.com/Ed94/Odin.git synced 2026-08-05 23:28:48 +00:00
Code Issues Packages Projects Releases Wiki Activity

big: Split up gcd + lcm.

Browse Source
This commit is contained in:
Jeroen van Rijn
2021-08-11 20:59:53 +02:00
parent 62dcccd7ef
commit c3db24f834
3 changed files with 185 additions and 197 deletions
Download Patch File Download Diff File Expand all files Collapse all files
core/math/big/basic.odin
+1 -195
View File
@@ -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 res_gcd == nil && res_lcm == nil { return nil; }
if err = clear_if_uninitialized(res_gcd, res_lcm, a, b); err != nil { return err; } if err = clear_if_uninitialized(res_gcd, res_lcm, a, b); err != nil { return err; }
az, _ := is_zero(a); bz, _ := is_zero(b); return #force_inline internal_int_gcd_lcm(res_gcd, res_lcm, a, 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);
} }
gcd_lcm :: proc { int_gcd_lcm, }; gcd_lcm :: proc { int_gcd_lcm, };
@@ -346,171 +317,6 @@ int_lcm :: proc(res, a, b: ^Int) -> (err: Error) {
} }
lcm :: proc { int_lcm, }; 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) remainder = numerator % (1 << bits)
*/ */
core/math/big/internal.odin
+178 -1
View File
@@ -680,7 +680,7 @@ internal_int_divmod :: proc(quotient, remainder, numerator, denominator: ^Int, a
// err = _int_div_recursive(quotient, remainder, numerator, denominator); // err = _int_div_recursive(quotient, remainder, numerator, denominator);
} else { } else {
when true { when true {
err = _private_int_div_school(quotient, remainder, numerator, denominator); err = #force_inline _private_int_div_school(quotient, remainder, numerator, denominator);
} else { } else {
/* /*
NOTE(Jeroen): We no longer need or use `_private_int_div_small`. 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; 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) { 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); 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 ======================= ======================== End of private procedures =======================
core/math/big/test.py
+6 -1
View File
@@ -446,7 +446,6 @@ TESTS = {
test_factorial: [ test_factorial: [
[ 6_000 ], # Regular factorial, see cutoff in common.odin. [ 6_000 ], # Regular factorial, see cutoff in common.odin.
[ 12_345 ], # Binary split factorial [ 12_345 ], # Binary split factorial
[ 100_000 ],
], ],
test_gcd: [ test_gcd: [
[ 23, 25, ], [ 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_passes = 0
total_failures = 0 total_failures = 0