# Floating point exception error, is Math.log2 the problem

This code does a Mersenne Prime search.
It fails with an error such as: `10690 floating point exception`

The only thing I can think is causing the problem is the function: `Math.log2(...)`
There doesn’t seem to be a `BigFloat` equivalent.

Here’s the code - `mpsearch2.cr`

``````# Compile as: \$ crystal build mpsearch2.cr --release
require "big"

module Primes
module Utils

# Returns true if +self+ is a prime number, else returns false.
def prime? (k = 10)
primes = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47}
return primes.includes? self if self <= primes.last
return false if 614_889_782_588_491_410.to_big_i.gcd(self.to_big_i) != 1
# wits = [range, [wit_prms]] or nil
wits = WITNESS_RANGES.find { |range, wits| range > self }
witnesses = wits && wits[1] || k.times.map{ 2 + rand(self - 4) }
miller_rabin_test(witnesses)
end

# Returns true if +self+ passes Miller-Rabin Test on witnesses +b+
private def miller_rabin_test(witnesses)
neg_one_mod = n = d = self - 1
d >>= d.trailing_zeros_count
witnesses.each do |b|
next if (b % self).zero?
y = powmod(b, d, self)
s = d
until y == 1 || y == neg_one_mod || s == n
y = (y &* y) % self        # y = (y**2) mod self
s <<= 1
end
return false unless y == neg_one_mod || s.odd?
end
true
end

private def powmod(b, e, m)
r, b = typeof(m).new(1), b.to_big_i
while e > 0
r = (r &* b) % m if e.odd?
e >>= 1
b = (b &* b) % m
end
r
end

# Best known deterministic witnnesses for given range and set of bases
# https://miller-rabin.appspot.com/
# https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test
private WITNESS_RANGES = {
341_531u64 => {9345883071009581737u64},
1_050_535_501u64 => {336781006125u64, 9639812373923155u64},
350_269_456_337u64 => {4230279247111683200u64, 14694767155120705706u64, 16641139526367750375u64},
55_245_642_489_451u64 => {2u64, 141889084524735u64, 1199124725622454117u64, 11096072698276303650u64},
7_999_252_175_582_851u64 => {2u64, 4130806001517u64, 149795463772692060u64, 186635894390467037u64, 3967304179347715805u64},
585_226_005_592_931_977u64 => {2u64, 123635709730000u64, 9233062284813009u64, 43835965440333360u64, 761179012939631437u64, 1263739024124850375u64},
18_446_744_073_709_551_615u64 => {2, 325, 9375, 28178, 450775, 9780504, 1795265022},
318_665_857_834_031_151_167_461u128   => {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37},
3_317_044_064_679_887_385_961_981u128 => {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41}
}
end
end

struct Int; include Primes::Utils end

def mp1(k) 1 + Math.log2(15*k + 1).to_big_i end
def mp1_candidates(k)
x = mp1(k);
pp [Math.log2(30*k + 2).to_big_i] if x.prime? && (2.to_big_i**x - 1).prime?
end

def mp2(k) 1 + Math.log2(15*k + 4).to_big_i end
def mp2_candidates(k)
x = mp2(k)
pp [Math.log2(30*k + 8).to_big_i] if x.prime? && (2.to_big_i**x - 1).prime?
end

k = 0.to_big_i

loop do
if k.even?
mp2_candidates(k)
else
mp1_candidates(k)
end
k = k.even? ? k*4 + 1 : k*4
end
``````

Here’s typical output (the next should be [1279])

``````➜  crystal-projects time ./mpsearch2
[3]
[5]
[7]
[13]
[17]
[19]
[31]
[61]
[89]
[107]
[127]
[521]
[607]
[2]    10690 floating point exception  ./mpsearch2
./mpsearch2  0.06s user 0.00s system 99% cpu 0.057 total
``````

`Math.log2` calls `#to_f64` on its argument, so when a `BigInt` is larger than `Float64::MAX`, it becomes infinity and then GMP intentionally raises SIGFPE on `Float64::INFINITY.to_big_i`. One possible workaround is:

``````def big_i_log2(k : BigInt)
normalized, exp = Math.frexp(k.to_big_f)
(Math.log2(normalized) + exp).to_f64
end
``````

This gives `[1279]`, `[2203]`, `[2281]`, `[3217]`, and so on.

`Math.log2(BigFloat)` would only be available in MPFR: https://www.mpfr.org/mpfr-current/mpfr.html#Transcendental-Functions

2 Likes

Thank you, again, @HertzDevil.

Crystal would be a world class language for Number Theory, Cryptography, Data Analytics, etc, if it could do arbitrary size math for all the transcendental, trigonometric, algebraic, etc, functions,
provided in the MPFR library.

It would be so much easier to do these algorithms in Crystal, than having to use Rust, C|C++, etc.
(Don’t get me started about not being able to do true parallel thread algorithms yet too. )

3 Likes

Similarly, how would you wrap MPFR for this:

``````      #  r = b^e mod m
def powmod(b, e, m)
r, b = typeof(m).new(1), b.to_big_i
while e > 0
r = (r &* b) % m if e.odd?
e >>= 1
b = (b &* b) % m
end
r
end
``````

I’ve had it on my todo list to implement a crystal native implementation of BigInt (and I guess BigFloat at some point too), because I have the same thoughts too. Using crystal for a bunch of scientific or mathematical calculations just feels nicer than any other language (except Ruby, for obvious reasons).

2 Likes

At the risk of sounding like a one-note song, I just want to say: I agree!

Life hasn’t let me pursue this for a while, but I’d like to know who else out there is interested in doing scientific computing in crystal.

I wrote this group theory library as part of my PhD CrystalSymm · GitLab and have always wanted to try to make a library that derives point groups and space groups instead of hard-coding them.

My current work in trapped-ion quantum computing has taught me that many people rely on python because of the scientific community. I would love to make it realistic for them to do their work in Crystal and see if we can’t improve the developer experience.

Some day!

6 Likes

I think it could be useful to set up a working group on scientific|numeric language development.

At minimum, the devs would need to:

1. setup a process to teach people how use git.
2. setup course|documentation on how to wrap|test|benchmark libraries (MPFR, GMP, etc).
3. setup process to create necessary documentation, tutorials, videos, etc.

I have no real knowledge how to do these things, but I would be willing to learn, to contribute.
I also have a number of algorithms already done in Crystal I could contribute|document.

So if there’s a standard process for people to contribute I’m sure many will learn it to do so.
This process could apply to other specific fields of software development, and ease the burden on the devs, who would primarily just need to review|approve the submitted work.

Again, these are my thoughts off the top of my head, but I would really like to feel confident I know what to do to contribute toward these goals in a tangible way, even if it’s just writing documentation.

2 Likes

GMP supports powmod, we just haven’t really agreed upon the API for that

Yeh, I remember seeing that PR now.

FWIW, my ideal preference would be: `a.powmod(e, m)`

@jzakiya that makes sense to me! My mind tends to be more on the “killer framework” side of building a community, but that won’t work without the guts/bindings!

I don’t know anything about forming working groups and can’t commit time right now, but I’m definitely interested in supporting and maybe contributing in the future.