Crystal: integer overflow in gcd function

[jzakiya]

The Crystal gcd implementation has a bug. I haven't found its code yet to see what it is as I just discovered it because it was returning errors in my code.

Ex: All these results should be '1'

puts 614_889_782_588_491_410.gcd(53) # => 7
puts 614_889_782_588_491_410.gcd(59) # => 95
puts 614_889_782_588_491_410.gcd(61) # => 31

To fix function in my code I implemented Stein's gcd algorithm

https://en.wikipedia.org/wiki/Binary_GCD_algorithm
https://lemire.me/blog/2013/12/26/fastest-way-to-compute-the-greatest-common-divisor/

This is inside a module, inside of Struct Int so I can do: a.gcd b like in Ruby.

    def bbgcd(b)
      a, b = self.to_big_i, b.to_big_i
      return b if a == 0
      return a if b == 0

      k = typeof(a).zero
      while (a | b).even?
        a >>= 1
        b >>= 1
        k += 1
      end

      while a.even?; a >>= 1 end

      while b != 0
        while b.even?; b >>= 1 end
        a, b = b, a if a > b
        b -= a
      end

      a << k
    end

Here are the correct results using it:

puts 614_889_782_588_491_410.bbgcd(53) => 1
puts 614_889_782_588_491_410.bbgcd(59) => 1
puts 614_889_782_588_491_410.bbgcd(61) => 1

This is another instance where I was pulling my hair trying to figure out why the same Ruby code was producing incorrect results in Crystal, until I became resigned to the fact that the code was correct and this function was returning incorrect results. And sure enough it was for certain numbers. Again, I just verified this, and haven't had a chance to look at the Crystal implementation, but it may be better to use Stein's implementation. Whether to use it was also a discussion issue for Ruby.

[RX14]

This is the implementation: https://github.com/crystal-lang/crystal/blob/master/src/int.cr#L310

If you can write a PR which replaces the implementation (remember to run bin/crystal tool format) and adds some specs in spec/std/int_spec.cr it'd be great! You can even bind the cttz llvm intrinsics, similar to here except using cttz instead of ctpop. This allows translating the gcd function as-is from your second blog post.

[jzakiya]

Since gcd is always no bigger than the bigger number -- a.gcd(a) => a -- maybe you can compare each operand and use type of biggest for both. In my code example I set both to BigInts because I just wanted to make sure it always worked in my code.

[jzakiya]

Hey @RX14 all that C stuff is way above my paygrade, that's why I use Ruby/Crystal. :-)

[RX14]

The actual bug here is just an overflow, 614_889_782_588_491_410.gcd(53) performs 53i32 % 614_889_782_588_491_410i64, i.e. the exact same bug as #6584: https://carc.in/#/r/4x1h. This new algorithm is good though, since it's faster.

The simplest fix would be to always call self.class.new(other) in gcd, which would keep the same semantics as other operators, i.e. always keeping the left-hand type (here it actually can form a union).

[jzakiya]

The simplest fix would be to always call self.class.new(other) in gcd, which would keep the same semantics as other operators, i.e. always keeping the left-hand type (here it actually can form a union).

This is not sufficient. In doing self.class.new(other) if other is bigger type (say BigInt) than self (say UInt64) you can get the same problem. That's why inside my bbgcd function I set both operands to BigInt and eliminated having to account for all the type permutations (a real PITA) so I can use it universally for all possibilities (though there may be marginal performance gains having smaller types operate together).

In fact, I would recommend just doing that, and be done with it, otherwise you must use the larger type for both to always get correct results.

[RX14]

@jzakiya it's sufficient to make the behaviour the same as other operators - which will overflow if the rhs type is larger, as you already discovered in #6584.

Using bigint everywhere would destroy performance and it's not a feasible solution. In a release or two we will have exceptions on overflow, and this error would immediately have been caught as an exception and you could have adjusted your types. So overflow checks will solve this.

[jzakiya]

Using bigint everywhere would destroy performance and it's not a feasible solution.

Has this been tested to be shown to be empirically true? My bbgcd function doesn't seem to have any measurable negative performance effects in my code, but more importantly, I know it will give correct results for all size numbers.

For me, and I hope you all hear me saying this as a user/programmer, the 3 most important things are:

1. Operational correctness of code
2. Time to write/test the code
3. Performance of the code

In my project of translating my primes-utils gem from Ruby to Crystal, I'm still stuck at verifying 1) while the time spent in 2) continues to increase. While this effort has been a learning experience, it hasn't been a very pleasurable one so far, especially coming from Ruby where Matz makes programmer happiness top priority.

I'm saying this to let you know, as just one user, if I can't create operationally provable correct code in a reasonable amount of time then I don't care how fast it runs to produce (unexpected) incorrect results.
And I'm beginning to really dislike languages that use number types (including Nim, et al) forcing me to be a manual compiler to juggle doing math to fit the language requirements. In fact, in Matz's EuRuKo keynote in August he specifically stated Ruby will never use static typing.

This whole exercise for me was to see what benefits Crystal would provide in comparison to Ruby. Yes, it can be faster (in code execution) but at a severe cost in code development and testing time. And with
Truffleruby already demonstrating significant performance improvements over MRI (while Ruby 2.6 is also faster than 2.5), and being able to provide true parallel programming capabilities (which Crystal currently lacks) , it will create even fewer incentives for me to use compiled/typed language, with all their headaches.

I'm not trying to be harsh. I'm just very frustrated right now, and don't know what else in my code may not give me expected results because of the implementation of the language. :(

[asterite]

I think we just need to restrict gcd to the same type of integers. The code right now is:

  def gcd(other : Int)
    self == 0 ? other.abs : (other % self).gcd(self)
  end

so the return type is typeof(self) | typeof(other) because of the other.abs. And that's not good. The return type should probably be the type of the "largest" integer type, but that's a bit inconsistent with the other operators. It gets tricky with signed vs. unsigned integers though... so probably restricting to the same types is good.

@jzakiya About using BigInt: it depends on libgmp. Not everyone has it installed, nor not everyone needs it for computing gcd of two numbers. In fact, another solution that I'm thinking right now is to completely remove gcd and lcm from the standard library. I personally never used them. They are probably much more useful for numeric/science/math stuff, and given that there are specialized algorithms for this it might be good to just delegate this to a shard.

Ideally BigInt and others should be implemented in Crystal, or libgmp should come bundled with a Crystal installation (and so libyaml and all libraries that are in the standard library). Then it would be fine to require BigInt by default because we know it'll always be available. But that's not the case right now.

Plus, I think BigInt is not needed for gcd, the result will always be equal or less than the lowest number.

Finally, Ruby is more than 20 years old. Of course it has a very polished API and it functions really well. I'm pretty sure in the past it wasn't like that (for example the rescue keyword used to be written resque).

For now, though, I think restricting the types is a good compromise. If anyone wants to do that, please open a PR.

[asterite]

Thinking it a bit more, I think we should remove gcd and lcm altogether, and just leave them in BigInt. For example:

100.lcm(100_000_001) # => 1410065508 ???
100_i64.lcm(100_000_001_i64) # => 10000000100

With lcm it's even worse because for two big primes the result is their multiplication, which will very likely overflow. I know we can just say "well, it throws overflow exception and that's it" but I don't think that's a correct solution here. I'm pretty sure that if one needs lcm and/or gcd in code, then one should probably be using BigInt to avoid overflow altogether. That is, for math/science stuff one should always use BigInt to be safe (or use Ruby/Python/Julia/Matlab/R which will always be safe by default).

And if we remove lcm because of the above reason, we should remove gcd, even though chances of overflowing are less likely.

As I said somewhere else, now I'd like Crystal to have integers automatically convert themselves to BigInt on overflow, so just having one Number type that can hold any number type, but I think it's a bit too late for that, plus it will usually kill most of performance gains we have, I think.

So my suggestion is that for this use case one should, at least for now, use another language.

[straight-shoota]

@asterite s/use another language/use BigInt/ ;)

[asterite]

Sure, if you don't mind having to write BigInt.new(...) or .to_big_i all over the place ^_^

[gmarcais]

Not to get into a language comparison, but Julia is in the same situation as Crystal here. The default type is Int64 and it does not promote fixed sized integer to large integer automatically, as the simple example of factorial implementation below shows:

> julia
               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.0.0 (2018-08-08)
 _/ |\__'_|_|_|\__'_|  |  
|__/                   |

julia> typeof(1)
Int64

julia> f(x) = x == 1 ? 1 : x * f(x-1)
f (generic function with 1 method)

julia> f(3)
6

julia> f(5)
120

julia> f(100)
0

[asterite]

@gmarcais Interesting! I thought Julia behaved more correctly...

[jzakiya]

julia> f(x) = (x | 1) == 1 ? 1 : x * f(x-1)
f (generic function with 1 method)

julia> f(1)
1
julia> f(0)
1
julia> f(BigInt(100))
93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000

[jzakiya]

Ruby

2.5.1 :075 > def f(x); x = (x | 1) == 1 ? 1 : x * f(x-1) end
 => :f 
2.5.1 :076 > f(100)
 => 93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000 
2.5.1 :077 > f(1)
 => 1 
2.5.1 :078 > f(0)
 => 1

[jzakiya]

In file fact.cr

require "big"

def f(x); x = (x | 1) == 1 ? 1 : x * f(x-1) end

puts "f(100) = #{ f(100) }"

print "f(BigInt.new(100)) = ", f(BigInt.new("100"))

puts

puts "f(1) = #{ f(1) }"

puts "f(0) = #{ f(0) }"

puts

def f1(x); x = x.to_big_i;  x = (x | 1) == 1 ? 1 : x * f1(x-1) end

puts "f1(100) = #{ f1(100) }"

puts "f1(200) = #{ f1(200) }"

puts "f1(1) = #{ f1(1) }"

puts "f1(0) = #{ f1(0) }"

Using 0.26.1

  Crystal Projects crystal run fact.cr --release
f(100) = 0
f(BigInt.new(100)) = 93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000
f(1) = 1
f(0) = 1

f1(100) = 93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000
f1(200) = 788657867364790503552363213932185062295135977687173263294742533244359449963403342920304284011984623904177212138919638830257642790242637105061926624952829931113462857270763317237396988943922445621451664240254033291864131227428294853277524242407573903240321257405579568660226031904170324062351700858796178922222789623703897374720000000000000000000000000000000000000000000000000
f1(1) = 1
f1(0) = 1

[yxhuvud]

Ok, so I did some benchmarking (on my 5 year old stationary computer. If someone has a more modern computer it would be nice to have a verification that the results are similar):
https://gist.github.com/yxhuvud/1ea800431d62f0d273e7dd2e06f8b65e#file-gcd

Results:

(master *>): crystal run gcd_bench.cr --release -D preview_overflow
Verifying binary_gcd
no overflows
Int8
      current   2.51M (398.56ns) (± 0.59%)    0.0B/op   2.21× slower
BigInt/LibGMP 448.73k (  2.23µs) (± 6.96%)  1.69kB/op  12.35× slower
   binary_gcd   5.54M (180.39ns) (± 0.78%)    0.0B/op        fastest

UInt8
      current   1.88M (532.22ns) (± 0.40%)    0.0B/op   3.61× slower
BigInt/LibGMP 483.50k (  2.07µs) (± 1.95%)  1.69kB/op  14.02× slower
   binary_gcd   6.78M (147.48ns) (± 1.66%)    0.0B/op        fastest

Int16
      current 981.41k (  1.02µs) (± 0.55%)    0.0B/op   2.82× slower
BigInt/LibGMP 412.77k (  2.42µs) (± 2.14%)  1.69kB/op   6.72× slower
   binary_gcd   2.77M (360.72ns) (± 0.51%)    0.0B/op        fastest

UInt16
      current   1.50M (665.17ns) (± 0.39%)    0.0B/op   1.16× slower
BigInt/LibGMP 428.11k (  2.34µs) (± 1.45%)  1.69kB/op   4.06× slower
   binary_gcd   1.74M (575.06ns) (± 0.54%)    0.0B/op        fastest

Int32
      current 587.63k (  1.70µs) (± 0.69%)    0.0B/op   2.56× slower
BigInt/LibGMP 368.09k (  2.72µs) (± 1.47%)  1.69kB/op   4.09× slower
   binary_gcd   1.50M (664.99ns) (± 1.03%)    0.0B/op        fastest

UInt32
      current 698.56k (  1.43µs) (± 0.51%)    0.0B/op   1.92× slower
BigInt/LibGMP 357.24k (  2.80µs) (± 2.33%)  1.69kB/op   3.75× slower
   binary_gcd   1.34M (746.54ns) (± 1.37%)    0.0B/op        fastest

Int64
      current 163.39k (  6.12µs) (± 0.41%)    0.0B/op   4.64× slower
BigInt/LibGMP 270.68k (  3.69µs) (± 1.63%)  1.69kB/op   2.80× slower
   binary_gcd 757.45k (  1.32µs) (± 0.64%)    0.0B/op        fastest

UInt64
      current 188.22k (  5.31µs) (± 0.48%)    0.0B/op   1.97× slower
BigInt/LibGMP 273.87k (  3.65µs) (± 1.20%)  1.69kB/op   1.35× slower
   binary_gcd 370.63k (  2.70µs) (± 0.98%)    0.0B/op        fastest

I also tried lemire's variant, but a: it was not conclusively faster on anything but Int8. b: misbehaved in the benchmarks somehow and got unreasonably fast so I'm fairly certain something that shouldn't be optimized away was removed.

Binary gcd seems the way to go. Next step would be to figure out the desired return type - I'm thinking for GCD, to take the greater size of self and the argument, and for LCM match whatever regular multiplication ends up as.

EDIT: Oh, and no 128 bit tests as there are no working random numbers there.

[RX14]

If it's just plain faster can someone PR it?

[yxhuvud]

@RX14 I'm currently working on getting the types right (mainly by casting stuff in ways that make certain the result can't be a union type), as well as adding actual tests for that. But it is quite a lot of cases to go through as I'm trying to keep down the amount of overflows that one would get in the case of different integer sizes and signedness.

[asterite]

I'm thinking... do we need gcd in the standard library? I never used it in my life (in any language). I'm sure it has its use cases in math-related programs, but crystal's std isn't math-oriented (for instance we don't have a factorial function).

[yxhuvud]

I've used it (one cute task in an AOC that reduced to a oneliner using gcd), but you may be correct that gcd (and lcm which depend on it) may not be used often enough to motivate its inclusion.

It is certainly true that they are needed less seldom than, for example, priority queues, which are used a lot more often but which isn't included.

And even if we would want it in stdlib, it could be that it would rather belong in a Math module rather than polluting the Integer namespace.

[yxhuvud]

Anyhow, I fixed up my implementation of a fix in #7996. Go forward with that one or go through with the removal?

[vlazar]

Fixed? Original example now shows 1 in all cases.

Crystal 0.31.1 on MacOS

$ crystal --version
Crystal 0.31.1 (2019-10-02)

LLVM: 8.0.1
Default target: x86_64-apple-macosx

[asterite]

This seems to work now but I don't understand why.

The result is still incorrectly a union, though.

7996 will eventually fix this.