Add powmod to math library

[jzakiya]

Following on the discussions raised here:

#8612

#13244

I implemented an optimized universal powmod function that works with any integer type Int|UIint|BigInt.

As shown it has two main functions, powmodgmp and powmodint.

powmodint is the classical numerical algorithm to use for values that fit in UInt64 sizes.
This assumes a 64-bit cpu default architecture.

powmodgmp handles values that can exceed UInt64 to eliminate overflow cases.

They are combined in powmod, and selected appropriately as shown.
From extensive testing I found, if the base m is < UInt32::MAX then the
max power value at each stage fits within UInt64 so its fastest to set m to that.

For m > UInt32::MAX, doing math with UInt128 for values < UInt64::MAX in podmodint was
slower that just using BigInt with powmodgmp.

Another possible consideration is to return the type of m for the output value, since
by definition the result must be an integer from 0..m-1.

So in powmod can do for each branch: m.class.new(powmodxxx(b,e,m))

require "big"

@[Link("gmp")]
lib LibGMP
  fun mpz_powm = __gmpz_powm(rop : MPZ*, base : MPZ*, exp : MPZ*, mod : MPZ*)
end

def powmodgmp(b, e, m)
  y = BigInt.new
  LibGMP.mpz_powm(y, b.to_big_i, e.to_big_i, m.to_big_i)
  y
end

def powmodint(b, e, m)
  r = m.class.new(1)
  b = b % m;  b = m.class.new(b)
  while e > 0
    r = (r &* b) % m if e.odd?
    e >>= 1
    b = (b &* b) % m
  end
  r
end

def powmod(b, e, m)
  if m > UInt32::MAX
    powmodgmp(b, e, m) or  m.class.new(powmodgmp(b, e, m))
  else
    powmodint(b, e, m.to_u64) or m.class.new(powmodint(b ,e , m.to_u64))
  end
end

I assume this would go in the Math.cr library.

As per the referenced discussions, powmod is a frequently used function in various
fields of math, statistics, cryptography, data analysis, and more.

[jzakiya]

Oh I forgot to mention, Ruby has a fast (I think C implemented) native powmod function.
So adding this to Crystal would provide more out-of-the-box compatibility with it.

Also a shorter implementation of powmod could be

def powmod(b, e, m)
    return powmodgmp(b, e, m) if m > UInt32::MAX
    powmodint(b, e, m.to_u64)
end

though from the compiler's perspective I don't know if this is any faster.

[jzakiya]

I ran some old application code using powmod with Crystal 1.16.0-dev-1.

Converting the values with m.class.new(xxx) significantly slows the performance, especially when converting from BigInt from powmodgmp.

That why I didn't do that (I remembered after running the code again) as it's not only unnecessary, it's slower.