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. :smiley: :+1:

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. :rage:)

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.