fib-tzn_primality_test.sf

#!/usr/bin/ruby

# Daniel "Trizen" Șuteu
# Date: 01 August 2017
# https://github.com/trizen

# An efficient implementation of a new primilaty test, inspired from the AKS primality test.

# When n>2 is a (pseudo)prime:
#
#   (2 + sqrt(-1))^n - (sqrt(-1))^n - 2 = 0 (mod n)
#

# By breaking the formula into pieces, we get the following equivalent statements:
#
#   5^(n/2) * cos(n * atan(1/2)) = 2 (mod n)
#   5^(n/2) * sin(n * atan(1/2)) = { n-1   if n=3 (mod 4)
#                                      1   if n=1 (mod 4) } (mod n)
#

# Additionally, we have the following two identities:
#
#   cos(n * atan(1/2)) = (((2+i)/sqrt(5))^n + exp(-1 * log((2+i)/sqrt(5)) * n))/2
#   sin(n * atan(1/2)) = (((2+i)/sqrt(5))^n - exp(-1 * log((2+i)/sqrt(5)) * n))/(2i)
#

# For numbers of the form `2n+1`, the above formulas simplify to:
#
#   cos((2*n + 1) * atan(1/2)) = a(n)/(sqrt(5) * 5^n)
#   sin((2*n + 1) * atan(1/2)) = b(n)/(sqrt(5) * 5^n)
#
# where `a(n)` and `b(n)` are integers given by:
#
#   a(n) = real((2 + sqrt(-1))^n)
#   b(n) = imag((2 + sqrt(-1))^n)
#
# Defined recursively as:
#
#   a(1) = 2; a(2) = 3; a(n) = 4*a(n-1) - 5*a(n-2)
#   b(1) = 1; b(2) = 4; b(n) = 4*b(n-1) - 5*b(n-2)
#

# Currently, we use only the `b(n)` branch, as it is strong enough to reject most composites.

# Known counter-examples (for both branches) (in order):
#   [1105, 2465, 10585, 15841, 29341, 38081, 40501, 41041, 46657, 75361, ...]

# Known counter-examples (for `b(n)` branch) (in order):
#   [1105, 1729, 2465, 10585, 15841, 29341, 38081, 40501, 41041, 46657, ...]

###########################--CURRENT TEST--################################

func modulo_test(n, mod) {
    func (n) is cached {

        n.is_zero && return 1
        n.is_one  && return 4

        var k = (n >> 1)

        n.is_even
            ? (((__FUNC__(k) * __FUNC__(k    )) - (5 * __FUNC__(k - 1) * __FUNC__(k - 1))) % mod)
            : (((__FUNC__(k) * __FUNC__(k + 1)) - (5 * __FUNC__(k - 1) * __FUNC__(k    ))) % mod)

    }(n - 1)
}

func is_probably_prime(n) {

    var r = modulo_test(n, n)

    (n % 4 == 3) ? (r == n-1) : (r == 1)
}

###########################--OLD TEST--################################

func a(n) is cached {       # real((2 + sqrt(-1))^n)
    n == 1 && return 2
    n == 2 && return 3

    4*a(n-1) - 5*a(n-2)
}

func b(n) is cached {       # imag((2 + sqrt(-1))^n)
    n == 1 && return 1
    n == 2 && return 4

    4*b(n-1) - 5*b(n-2)
}

func is_probably_prime_old_1(n) {

    var _a_ = (a(n) % n) == 2 || return false
    var _b_ = (b(n) % n)

    (n % 4 == 3) ? (_b_ == n-1) : (_b_ == 1)
}

###########################--OLD TEST--################################

func is_probably_prime_old_2(n) {

    define r = 1/2
    local Number!PREC = ceil(log2(5) * (n*r) + 2*log(n))

    #var k = (isqrt(5**n) + r)
    var k = (5**(n >> 1) * sqrt(5))
    var p = (n * atan(r))

    var _a_ = (k * cos(p) -> round % n) == 2 || return false
    var _b_ = (k * sin(p) -> round % n)

    (n % 4 == 3) ? (_b_ == n-1) : (_b_ == 1)
}

######################################################################

say "=> Testing for counter-examples..."

each(3..100, { |n|
    if (is_probably_prime(n)) {

        if (!n.is_prime) {
            warn "Counter-example: #{n}"
        }
    }
    elsif (n.is_prime) {
        warn "Missed a prime: #{n}"
        Sys.sleep(1)
    }
})

#
## Run a few tests
#

say (is_probably_prime(6760517005636313)   ? 'prime' : 'error')    #=> prime
say (is_probably_prime(204524538079257577) ? 'prime' : 'error')    #=> prime
say (is_probably_prime(904935283655003749) ? 'prime' : 'error')    #=> prime

say "=> Testing a few large primes..."

say (is_probably_prime(90123127846128741241234935283655003749)                             ? 'prime' : 'error')    #=> prime
say (is_probably_prime(793534607085486631526003804503819188867498912352777)                ? 'prime' : 'error')    #=> prime
say (is_probably_prime(6297842947207644396587450668076662882608856575233692384596461)      ? 'prime' : 'error')    #=> prime
say (is_probably_prime(396090926269155174167385236415542573007935497117155349994523806173) ? 'prime' : 'error')    #=> prime

say "=> Testing a few large Mersenne primes...";

# Mersenne primes
say (is_probably_prime(2**127  - 1) ? 'prime' : 'error')   #=> prime
say (is_probably_prime(2**521  - 1) ? 'prime' : 'error')   #=> prime
say (is_probably_prime(2**1279 - 1) ? 'prime' : 'error')   #=> prime
say (is_probably_prime(2**3217 - 1) ? 'prime' : 'error')   #=> prime
say (is_probably_prime(2**4423 - 1) ? 'prime' : 'error')   #=> prime