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