tribonacci_primality_test.sf

#!/usr/bin/ruby

# Author: Daniel "Trizen" Șuteu
# Date: 18 May 2019
# https://github.com/trizen

# A new primality test, using a Tribonacci-like sequence.

# Sequence definition:
#   T(0) = 0
#   T(1) = 0
#   T(2) = 9
#   T(n) = T(n-1) + 3*T(n-2) + 9*T(n-3)

# Closed form:
#   T(n) = (-9 sqrt(2) (-1 + i sqrt(2))^n + 2 (sqrt(2) + 4 i)×3^n + (7 sqrt(2) - 8 i) (-1 - i sqrt(2))^n)/(4 (sqrt(2) + 4 i))

# The sequence starts as:
#   0, 0, 9, 9, 36, 144, 333, 1089, 3384, 9648, 29601, 89001, 264636, 798048, ...

# When p is a prime > 5 congruent to {1,3} mod 8, then T(p) == 0 (mod p).
# When p is a prime > 5 congruent to {5,7} mod 8, then T(p) == 4 (mod p).

# The strong version also includes a base-2 strong Miller-Rabin test.

# When n ≡ 1 (mod 8), the first few counter-examples to the weak version are:
#   88561, 107185, 162401, 221761, 226801, 334153, 410041, 665281, 825265, 1569457, 1615681, 2727649, 3057601, 3363121, 4767841, 5148001, 5176153, 5632705, 6313681, 7207201...

# When n ≡ 3 (mod 8), several counter-examples to both versions include (with possible gaps):
#   80375707, 154287451, 326266051, 478614067, 4297753027, 4340265931, 11362973707, 11723962867, 13204140643, 19601470027, 24422578867, 34087042891, 37870128451, ...

# Several counter-examples to the strong version (with possible gaps):
#   29111881, 213035761, 1312332001, 1750412161, 8585937361, 9856290601, 19036199881, 26244842401, 35252113601, 39897898921, 58318467481, ...

# See also:
#   https://trizenx.blogspot.com/2020/01/primality-testing-algorithms.html

# Matrix for p == {1,3} (mod 8)
const A123 = Matrix(
                [1, 0, 0],
                [0, 1, 0],
                [0, 0, 1]
            )

# Matrix for p == {5,7} (mod 8)
const A57 = Matrix(
                [4, 2, 3],
                [1, 5, 3],
                [1, 2, 6],
            )

# Base matrix
const BASE = Matrix(
                [0, 3, 0],
                [0, 0, 3],
                [1, 1, 1],
            )

# Tribonacci-like sequence modulo m
func tribonaccimod(n, m) {
    BASE.powmod(n, m)
}

func is_tribonacci_pseudoprime(n) {

    return false if (n <= 1)
    return true  if (n <= 3)
    return false if (n.is_even)
    return true  if (n == 5)

    # Strong base-2 test (optional)
    n.is_strong_pseudoprime(2) || return false

    # The remainder modulo 8
    var r = n%8

    # For primes p congruent to {1,3} mod 8, we get 0.
    if (r ~~ [1,3]) {

        # Compute T_(n-1) (mod n)
        var A = tribonaccimod(n-1, n)
        return (A == A123)
    }

    # For primes p congruent to {5,7} mod 8, we get 4.
    if (r ~~ [5,7]) {

        # Compute T_(n+1) (mod n)
        var A = tribonaccimod(n+1, n)
        return (A == A57)
    }

    return false
}

var C1 = [561, 1105, 1729, 2465, 2821, 6601, 8911, 10585, 15841, 29341, 41041, 46657, 52633, 62745, 63973, 75361, 101101, 115921,
         126217, 162401, 172081, 188461, 252601, 278545, 294409, 314821, 334153, 340561, 399001, 410041, 449065, 488881, 512461]

var C2 = [5148001, 7519441, 8719921, 29111881, 178837201, 186393481, 221884001, 392099401, 399906001, 434932961, 438359041, 555465601, 558977761]
var C3 = [29111881, 213035761, 1312332001, 1750412161, 8585937361, 9856290601, 19036199881, 26244842401, 35252113601, 39897898921, 58318467481]

say C1.grep(is_tribonacci_pseudoprime)    #=> []
say C2.grep(is_tribonacci_pseudoprime)    #=> [29111881]
say C3.grep(is_tribonacci_pseudoprime)    #=> [29111881, 213035761, 1312332001, 1750412161, 8585937361, 9856290601, 19036199881, 26244842401, 35252113601, 39897898921, 58318467481]

say is_tribonacci_pseudoprime.first(25)   #=> [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]