sidef-scripts/Math/quadratic_integers.sf

#!/usr/bin/ruby

# Simple implementation of quadratic integers.

# See also:
#   https://en.wikipedia.org/wiki/Quadratic_integer

class QuadraticInteger(a, b, w = 2) {   # represents: a + b*sqrt(w)

    method to_s {
        "QuadraticInteger(#{a}, #{b}, #{w})"
    }

    method ==(QuadraticInteger c) {
        (a == c.a) && (b == c.b) && (w == c.w)
    }

    method conjugate {
        QuadraticInteger(a, -b, w)
    }

    method norm {
        a*a - w*b*b
    }

    method add (Number c) {
        QuadraticInteger(a+c, b, w)
    }

    method add (QuadraticInteger z) {
        var (c,d) = (z.a, z.b)
        QuadraticInteger(a+c, b+d, w)
    }

    __CLASS__.alias_method(:add, '+')

    method sub (Number c) {
        QuadraticInteger(a-c, b, w)
    }

    method sub (QuadraticInteger z) {
        var (c,d) = (z.a, z.b)
        QuadraticInteger(a-c, b-d, w)
    }

    __CLASS__.alias_method(:sub, '-')

    method mul (Number c) {
        QuadraticInteger(a*c, b*c, w)
    }

    method mul (QuadraticInteger z) {
        var (c,d) = (z.a, z.b)
        QuadraticInteger(a*c + b*d*w, a*d + b*c, w)
    }

    __CLASS__.alias_method(:mul, '*')

    method mod (Number m) {
        QuadraticInteger(a % m, b % m, w)
    }

    __CLASS__.alias_method(:mod, '%')

    method pow(Number n) {
        var x = self
        var c = QuadraticInteger(1, 0, w)

        for bit in (n.digits(2)) {
            c *= x if bit
            x *= x
        }

        return c
    }

    __CLASS__.alias_method(:pow, '**')

    method powmod(Number n, Number m) {

        var x = self
        var c = QuadraticInteger(1, 0, w)

        for bit in (n.digits(2)) {
            (c *= x) %= m if bit        #=
            (x *= x) %= m               #=
        }

        return c
    }
}

#---------------------------------------------------------------------

# Determine if a given number is probably a prime number.
func is_quadratic_pseudoprime (n, r=2) {

    return false if (n <= 1)
    return true  if (n <= 3)

    return true if (r <= 0)

    var x = QuadraticInteger(r, 1, r+2).powmod(n, n)

    x.a == r || return false

    var y = QuadraticInteger(r, -1, r+2).powmod(n, n)

    y.a == r || return false

    (x.b + y.b == n) && __FUNC__(n, r-1)
}

say is_quadratic_pseudoprime(43)    #=> true
say is_quadratic_pseudoprime(97)    #=> true

with (QuadraticInteger(1, 1, 2)) {|q|
    say 15.of { q.pow(_).a }        #=> A001333
    say 15.of { q.pow(_).b }        #=> A000129
}

with (QuadraticInteger(1, 1, 3)) {|q|
    say 15.of { q.pow(_).a }        #=> A026150
    say 15.of { q.pow(_).b }        #=> A002605
}

var n = (274177-1)
var m = (2**64 + 1)

with (QuadraticInteger(3, 4, 2)) {|q|
    var r = q.powmod(n, m)
    say gcd(r.a-1, m)       #=> 2741177
    say gcd(r.b, m)         #=> 2741177
}

#---------------------------------------------------------------------

func is_gaussian_quadratic_pseudoprime(n) {

    return false if (n <= 1)
    return true  if (n <= 3)

    static x = QuadraticInteger(1, -1, -2)

    given (n%8) {
        when ([1,3]) {
            var t = x.powmod(n-1, n)
            (t.a==1 && t.b==0)
        }
        when ([5, 7]) {
            var t = x.powmod(n+1, n)
            (t.a==3 && t.b==0)
        }
        else {
            false
        }
    }
}

assert([88561,107185,162401,221761,226801,334153,410041,665281,825265,1569457,1615681,2727649].all(is_gaussian_quadratic_pseudoprime))
assert([80375707,154287451,267559627,326266051,478614067,573183451,643767931,2433943891,4297753027].all(is_gaussian_quadratic_pseudoprime))

for n in (1..1e2) {
    if (is_gaussian_quadratic_pseudoprime(n)) {
        if (n.is_composite) {
            say "Counter-example: #{n}"
        }
    }
    elsif (n.is_prime) {
        say "Missed prime: #{n}"
    }
}