Skip to content

Commit 679cbf1

Browse files
committed
math.big: restore gdc_euclid, use it for smaller numbers, fix bench_euclid.v .
1 parent d6db4f9 commit 679cbf1

File tree

4 files changed

+65
-23
lines changed

4 files changed

+65
-23
lines changed

vlib/math/big/big_test.v

Lines changed: 5 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -655,11 +655,14 @@ fn test_big_mod_pow() {
655655
}
656656
}
657657

658-
fn test_gcd() {
658+
fn test_gcd_vs_gcd_binary_vs_gcd_euclid() {
659659
for t in gcd_test_data {
660660
a := t.a.parse()
661661
b := t.b.parse()
662-
assert a.gcd(b) == t.expected.parse()
662+
expected := t.expected.parse()
663+
assert a.gcd(b) == expected
664+
assert a.gcd_binary(b) == expected
665+
assert a.gcd_euclid(b) == expected
663666
}
664667
}
665668

vlib/math/big/integer.v

Lines changed: 45 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -1021,6 +1021,19 @@ fn bi_max(a Integer, b Integer) Integer {
10211021

10221022
// gcd returns the greatest common divisor of the two integers `a` and `b`.
10231023
pub fn (a Integer) gcd(b Integer) Integer {
1024+
// The cutoff is determined empirically, using vlib/v/tests/bench/math_big_gcd/bench_euclid.v .
1025+
if b.digits.len < 8 {
1026+
return a.gcd_euclid(b)
1027+
}
1028+
return a.gcd_binary(b)
1029+
}
1030+
1031+
// gcd_binary returns the greatest common divisor of the two integers `a` and `b`.
1032+
// Note that gcd_binary is faster than gcd_euclid, for large integers (over 8 bytes long).
1033+
// Inspired by the 2013-christmas-special by D. Lemire & R. Corderoy https://en.algorithmica.org/hpc/analyzing-performance/gcd/
1034+
// For more information, refer to the Wikipedia article: https://en.wikipedia.org/wiki/Binary_GCD_algorithm
1035+
// Discussion and further information: https://lemire.me/blog/2013/12/26/fastest-way-to-compute-the-greatest-common-divisor/
1036+
pub fn (a Integer) gcd_binary(b Integer) Integer {
10241037
if a.signum == 0 {
10251038
return b.abs()
10261039
}
@@ -1031,24 +1044,42 @@ pub fn (a Integer) gcd(b Integer) Integer {
10311044
return one_int
10321045
}
10331046

1034-
return gcd_binary(a.abs(), b.abs())
1035-
}
1036-
1037-
// Inspired by the 2013-christmas-special by D. Lemire & R. Corderoy https://en.algorithmica.org/hpc/analyzing-performance/gcd/
1038-
// For more information, refer to the Wikipedia article: https://en.wikipedia.org/wiki/Binary_GCD_algorithm
1039-
// Discussion and further information: https://lemire.me/blog/2013/12/26/fastest-way-to-compute-the-greatest-common-divisor/
1040-
fn gcd_binary(x Integer, y Integer) Integer {
1041-
mut a, az := x.rsh_to_set_bit()
1042-
mut b, bz := y.rsh_to_set_bit()
1047+
mut aa, az := a.abs().rsh_to_set_bit()
1048+
mut bb, bz := b.abs().rsh_to_set_bit()
10431049
shift := umin(az, bz)
10441050

1045-
for a.signum != 0 {
1046-
diff := b - a
1047-
b = bi_min(a, b)
1048-
a, _ = diff.abs().rsh_to_set_bit()
1051+
for aa.signum != 0 {
1052+
diff := bb - aa
1053+
bb = bi_min(aa, bb)
1054+
aa, _ = diff.abs().rsh_to_set_bit()
10491055
}
1056+
return bb.left_shift(shift)
1057+
}
10501058

1051-
return b.left_shift(shift)
1059+
// gcd_euclid returns the greatest common divisor of the two integers `a` and `b`.
1060+
// Note that gcd_euclid is faster than gcd_binary, for very-small-integers up to 8-byte/u64.
1061+
pub fn (a Integer) gcd_euclid(b Integer) Integer {
1062+
if a.signum == 0 {
1063+
return b.abs()
1064+
}
1065+
if b.signum == 0 {
1066+
return a.abs()
1067+
}
1068+
if a.signum < 0 {
1069+
return a.neg().gcd_euclid(b)
1070+
}
1071+
if b.signum < 0 {
1072+
return a.gcd_euclid(b.neg())
1073+
}
1074+
mut x := a
1075+
mut y := b
1076+
mut r := x % y
1077+
for r.signum != 0 {
1078+
x = y
1079+
y = r
1080+
r = x % y
1081+
}
1082+
return y
10521083
}
10531084

10541085
// mod_inverse calculates the multiplicative inverse of the integer `a` in the ring `ℤ/nℤ`.

vlib/v/tests/bench/math_big_gcd/bench_euclid.v

Lines changed: 13 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -33,7 +33,8 @@ fn main() {
3333
mut clocks := Clocks(map[string]benchmark.Benchmark{})
3434

3535
for algo in [
36-
'euclid',
36+
'gcd',
37+
'gcd_euclid',
3738
'gcd_binary',
3839
//'u32binary',
3940
//'u64binary'
@@ -105,10 +106,10 @@ fn main() {
105106
println(clocks['euclid'].total_message('both algorithms '))
106107

107108
msg := [
108-
'Seems to me as if euclid in big.Integer.gcd() performs better on ',
109+
'Seems to me, that big.Integer.gcd_euclid in big.Integer.gcd() performs better on ',
109110
'very-small-integers up to 8-byte/u64. The tests #-1..5 show this.',
110-
'The gcd_binary-algo seems to be perform better, the larger the numbers/buffers get.',
111-
'On my machine, i see consistent gains between ~10-30-percent with :',
111+
'The big.Integer.gcd_binary algo seems to be perform better, the larger the numbers/buffers get.',
112+
'On my machine, I see consistent gains between ~10-30-percent with :',
112113
'\n',
113114
'v -prod -cg -gc boehm bench_euclid.v',
114115
'\n',
@@ -177,13 +178,20 @@ fn run_benchmark(data []DataI, heap bool, mut clocks Clocks) bool {
177178
for cycles < rounds {
178179
for set in testdata {
179180
match algo {
180-
'euclid' {
181+
'gcd' {
181182
if set.r != set.aa.gcd(set.bb) {
182183
eprintln('${algo} failed ?')
183184
clock.fail()
184185
break
185186
}
186187
}
188+
'gcd_euclid' {
189+
if set.r != set.aa.gcd_euclid(set.bb) {
190+
eprintln('${algo} failed ?')
191+
clock.fail()
192+
break
193+
}
194+
}
187195
'gcd_binary' {
188196
if set.r != set.aa.gcd_binary(set.bb) {
189197
eprintln('${algo} failed ?')

vlib/v/tests/bench/math_big_gcd/prime/maker.v

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -120,7 +120,7 @@ pub fn random_list(cfg []string) []string {
120120
}
121121
mut num := ''
122122
for prime_count != 0 {
123-
num = primes[prime_size][rand.int_in_range(0, primes[prime_size].len - 1)]
123+
num = primes[prime_size][rand.int_in_range(0, primes[prime_size].len - 1) or { 0 }]
124124
if num in p_list {
125125
continue
126126
}
@@ -140,7 +140,7 @@ pub fn random_list(cfg []string) []string {
140140
return p_list
141141
}
142142

143-
pub fn random_set(cfg PrimeCfg) ?[]PrimeSet {
143+
pub fn random_set(cfg PrimeCfg) ![]PrimeSet {
144144
p_lists := [
145145
cfg.r.split('.'),
146146
cfg.a.split('.'),

0 commit comments

Comments
 (0)