#!/usr/bin/ruby
# Tests for the "Quadratic" class.
# 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 = Quadratic(r, 1, r+2).powmod(n, n)
    x.a == r || return false
    var y = Quadratic(r, -1, r+2).powmod(n, n)
    y.a == r || return false
    (x.b + y.b == n) && __FUNC__(n, r-1)
}
assert(is_quadratic_pseudoprime(43))
assert(is_quadratic_pseudoprime(97))
with (Quadratic(1, 1, 2)) {|q|
    assert_eq(
        15.of { q.pow(_).a }        #=> A001333
        %n[1, 1, 3, 7, 17, 41, 99, 239, 577, 1393, 3363, 8119, 19601, 47321, 114243]
    )
    assert_eq(
        15.of { q.pow(_).b }        #=> A000129
        %n[0, 1, 2, 5, 12, 29, 70, 169, 408, 985, 2378, 5741, 13860, 33461, 80782]
    )
}
with (Quadratic(1, 1, 3)) {|q|
    assert_eq(
        15.of { q.pow(_).a }        #=> A026150
        %n[1, 1, 4, 10, 28, 76, 208, 568, 1552, 4240, 11584, 31648, 86464, 236224, 645376]
    )
    assert_eq(
        15.of { q.pow(_).b }        #=> A002605
        %n[0, 1, 2, 6, 16, 44, 120, 328, 896, 2448, 6688, 18272, 49920, 136384, 372608]
    )
}
var n = (274177-1)
var m = (2**64 + 1)
with (Quadratic(3, 4, 2)) {|q|
    var r = q.powmod(n, m)
    assert_eq(gcd(r.a-1, m), 274177)
    assert_eq(gcd(r.b, m), 274177)
}
do {
    var a = Quadratic(5, 8, 10)
    var b = Quadratic(3, 9, 10)
    assert(a > b)
    assert(!(a < b))
    assert(b < a)
    assert(!(b > a))
    assert(a == a)
    assert(!(a != a))
    assert(b == b)
    assert(a != b)
    assert(b != a)
    assert(!(b == a))
    assert_eq(a, a)
    assert_eq(b, b)
    assert_ne(b, a)
    assert_ne(a, b)
    assert(a > 4)
    assert(a >= 5)
    assert(a < 6)
    assert_eq(a.a, 5)
    assert_eq(a.b, 8)
    assert_eq(a.w, 10)
    assert_eq(a+5, Quadratic(a.a+5, a.b, a.w))
    assert_eq(a-5, Quadratic(a.a-5, a.b, a.w))
    assert_eq(a*5, Quadratic(a.a*5, a.b*5, a.w))
    assert_eq(a/5, Quadratic(a.a/5, a.b/5, a.w))
    assert_eq(a+b, Quadratic(a.a + b.a, a.b + b.b, a.w))
    assert_eq(a*b, Quadratic(a.a*b.a + a.b*b.b*a.w, a.b*b.a + a.a*b.b, a.w))
    assert_eq(a.inv, a**(-1))
    assert_eq(b.inv, b**(-1))
    assert_eq(a.invmod(43), a.powmod(-1, 43))
    assert_eq(b.invmod(97), b.powmod(-1, 97))
    assert_eq(b.invmod(146).sqr.mod(146), b.powmod(-2, 146))
    assert_eq(b.sqr.invmod(146), b.powmod(-2, 146))
    assert_eq(a+b -> to_n.round(-30), a.to_n+b.to_n -> round(-30))
    assert_eq(a-b -> to_n.round(-30), a.to_n-b.to_n -> round(-30))
    assert_eq(a*b -> to_n.round(-30), a.to_n*b.to_n -> round(-30))
    assert_eq(a/b -> to_n.round(-30), a.to_n/b.to_n -> round(-30))
}
func Gaussian(a,b=0) {
    Quadratic(a, b, -1)
}
var r = (-2 .. 2)
for a in (r), b in (r), c in (r), d in (r) {
    assert_eq([Gaussian(a,b) + Gaussian(c,d) -> reals], [cadd(a,b,c,d)])
    assert_eq([Gaussian(a,b) - Gaussian(c,d) -> reals], [csub(a,b,c,d)])
    assert_eq([Gaussian(a,b) * Gaussian(c,d) -> reals], [cmul(a,b,c,d)])
    if (c*c + d*d != 0) {
        assert_eq([Gaussian(a,b) / Gaussian(c,d) -> reals], [cdiv(a,b,c,d)])
    }
}
for a in (r), b in (r) {
    var n = irand(0, 100)
    var m = irand(100, 1000)
    assert_eq([Gaussian(a,b)**n -> reals], [cpow(a,b,n)])
    assert_eq([Gaussian(a,b)**n -> reals].map { .mod(m) } , [cpowmod(a,b,n,m)])
}
func gaussian_sum(n) {
    var total = [0, 0]
    for k in (1..n) {
        total = [cadd(total..., cdiv(cpow(0, 1, k-1), k))]
    }
    [cmul(total..., n!)]
}
var arr = 10.of(gaussian_sum)
assert_eq(arr.map{.head}, %n[0, 1, 2, 4, 16, 104, 624, 3648, 29184, 302976])
assert_eq(arr.map{.tail}, %n[0, 0, 1, 3, 6, 30, 300, 2100, 11760, 105840])
do {
    var m = 10001
    var a = 43
    var b = 97
    assert_eq(Gaussian(1,0) / Gaussian(a, b), Gaussian(a,b).inv)
    assert_eq(Mod(Gaussian(a, b), m).inv * Gaussian(a,b), Mod(Gaussian(1,0), m))
    assert_eq(Mod(Gaussian(a, b), m).inv * Gaussian(a,b) -> lift, Gaussian(1,0))
    assert_eq(Mod(Gaussian(a, b), m)**(-1) * Gaussian(a,b) -> lift, Gaussian(1,0))
    assert_eq([cdiv(1, 0, a, b)], [a / (a*a + b*b), -b / (a*a + b*b)])
    assert_eq([a * invmod(a*a + b*b, m), -b * invmod(a*a + b*b, m)].map{.mod(m)}, [complex_invmod(a, b, m)])
    assert_eq([cmod(cmul(a, b, complex_invmod(a, b, m)), m)], [1, 0])
}
assert(Gaussian(3,4) == Gaussian(3,4))
assert(!(Gaussian(3,4) == Gaussian(3,3)))
assert(!(Gaussian(3,3) == Gaussian(3,4)))
assert(Gaussian(3,3) != Gaussian(3,4))
assert(Gaussian(3,4) != Gaussian(3,3))
assert(!(Gaussian(3,4) != Gaussian(3,4)))
assert_eq(Gaussian(4,5) <=> Gaussian(3,4), 1)
assert_eq(Gaussian(3,4) <=> Gaussian(3,4), 0)
assert_eq(Gaussian(3,4) <=> Gaussian(3,5), -1)
assert(Gaussian(4,5) >  Gaussian(3,4))
assert(Gaussian(4,5) >  Gaussian(4,4))
assert(Gaussian(4,5) >= Gaussian(3,4))
assert(Gaussian(3,4) >= Gaussian(3,4))
assert(Gaussian(3,4) <= Gaussian(3,4))
assert(!(Gaussian(3,4) > Gaussian(3,5)))
assert(!(Gaussian(3,4) > Gaussian(4,5)))
assert(!(Gaussian(3,4) < Gaussian(3,4)))
assert(!(Gaussian(3,4) < Gaussian(3,1)))
assert(!(Gaussian(3,4) < Gaussian(2,1)))
assert(!(Gaussian(3,4) <= Gaussian(3,3)))
assert(!(Gaussian(2,1) >= Gaussian(2,2)))
for a in (r), b in (r) {
    var y = irand(2, 100)
    var x = Gaussian(a,b)
    assert_eq(x - floor(x/y)*y, x % y)
}
func gaussian_sum_2(n) {
    var i     = Gaussian(0, 1)
    var total = Gaussian(0)
    for k in (1..n) {
        total += (i**(k-1) / k)
    }
    total * n!
}
assert_eq(
    10.of(gaussian_sum_2),
    [Gaussian(0, 0), Gaussian(1, 0), Gaussian(2, 1), Gaussian(4, 3), Gaussian(16, 6), Gaussian(104, 30), Gaussian(624, 300), Gaussian(3648, 2100), Gaussian(29184, 11760), Gaussian(302976, 105840)]
)
assert_eq(powmod(Gaussian(3,4), 1000, 1e6), Gaussian(585313, 426784))
assert_eq([Mod(Gaussian(3,4), 1e6)**1000 -> lift.reals], [585313, 426784])
assert_eq(Mod(Gaussian(3,4), 1e6)**1000, Mod(Gaussian(585313, 426784), 1e6))
assert(Mod(Gaussian(3,4), 1e6)**1000 == Mod(Gaussian(585313, 426784), 1e6))
assert_eq(Mod(43, 97).to_n, 43)
assert_eq(Gaussian(3,4).to_n, Gaussian(3,4).to_n)
assert_eq(Gaussian(3,4).to_c, 3+4i)
assert_eq(Gaussian(42).invmod(2017), Gaussian(1969, 0))
assert_eq(Gaussian(3,4).invmod(2017), Gaussian(1291, 968))
assert_eq(Gaussian(91,23).invmod(2017), Gaussian(590, 405))
assert_eq(Gaussian(43, 99).invmod(2017), Gaussian(1709,1272))
assert_eq(Gaussian(43, 99).invmod(1234567), Gaussian(1019551, 667302))
assert_eq(Mod(Gaussian(42), 2017).inv, Mod(Gaussian(1969, 0), 2017))
assert_eq(Mod(Gaussian(3,4), 2017)**(-1), Mod(Gaussian(1291, 968), 2017))
assert_eq(Mod(Gaussian(91,23), 2017)**(-1), Mod(Gaussian(590, 405), 2017))
assert_eq(Mod(Gaussian(43, 99), 2017).inv, Mod(Gaussian(1709,1272), 2017))
assert_eq(Mod(Gaussian(43, 99), 1234567)**(-2), Mod(Gaussian(1019551, 667302)**2, 1234567))
assert_eq(Mod(Gaussian(43, 99), 1234567)**(-5), Mod(Gaussian(1019551, 667302)**5, 1234567))
assert_eq(powmod(Gaussian(43, 99), -4, 1234567), invmod(Gaussian(43, 99)**4 % 1234567, 1234567))
assert_eq(powmod(Gaussian(43, 99), -5, 1234567), invmod(Gaussian(43, 99)**5 % Gaussian(1234567), 1234567))
assert_eq(powmod(Gaussian(43, 99), -5, 1234567), invmod(Gaussian(43, 99)**5, 1234567))
assert_eq(powmod(Gaussian(43, 99), -4, 1234567), invmod(powmod(Gaussian(43, 99), 4, 1234567), 1234567))
assert_eq(powmod(Gaussian(43, 99), -5, 1234567), invmod(powmod(Gaussian(43, 99), 5, 1234567), 1234567))
assert_eq(Quadratic(Gauss(3,4), Gauss(5,6), 2).powmod(43, 97), Quadratic(Gauss(52, 44), Gauss(92, 8), 2))
assert_eq(Gaussian(43, 97)**(-5), (Gaussian(43,97)**5)**(-1))
assert_eq(Gaussian(43, 97)**(-5), (Gaussian(43,97)**5).inv)
assert_eq(Gaussian(43, 97)**(-5), (Gaussian(43,97).inv)**5)
assert_eq(Gaussian(43, 97)**(-5), (Gaussian(43,97)**(-1))**5)
assert_eq(Mod(Gaussian(43, 97), 1234567)**(-5),  Mod(Gaussian(43, 97), 1234567)**5 -> inv)
assert_eq(Mod(Gaussian(43, 97), 1234567)**(-5), (Mod(Gaussian(43, 97), 1234567)->inv)**5)
assert_eq(Mod(Gaussian(43, 97), 1234567)**(-1234),  Mod(Gaussian(43, 97), 1234567)**1234 -> inv)
assert_eq(Mod(Gaussian(43, 97), 1234567)**(-1234), (Mod(Gaussian(43, 97), 1234567)->inv)**1234)
assert_eq(Mod(Gaussian(3,4), 1234567)**1234 -> lift, powmod(Gaussian(3,4), 1234, 1234567))
assert_eq(Mod(Gaussian(3,4), 1234567)**-1234 -> lift, powmod(Gaussian(3,4), -1234, 1234567))
assert_eq(Gaussian(3/5,11/4)**(-27), Gaussian(3/5,11/4)**27 -> inv)
assert_eq(Gaussian(3/5,11/4)**(-27), Gaussian(3/5,11/4).inv**27)
assert_eq(Gaussian(Gaussian(3,4), Gaussian(17,19)).to_n.to_n, -16 + 21i)
assert_eq(Gaussian(Mod(13, 97), Mod(43, 97)).to_n.to_n, 13+43i)
assert_eq(Gaussian(Mod(13, 97), Mod(43, 97)), Gaussian(Mod(13, 97), Mod(43, 97)))
assert_eq(Mod(Gaussian(3, 4), 97)*1234, Mod(Gaussian(16, 86), 97))
assert_eq(Mod(Gaussian(3/4, 5/6), 1234567)**10 * Mod(Gaussian(3/4, 5/6), 1234567)**-10, Mod(Gaussian(1,0), 1234567))
assert_eq((Mod(Gaussian(43/3, 97/5), 127)**(-11) * Mod(Gaussian(43/3, 97/5), 127))**+9, Mod(Gaussian(52, 73), 127))
assert_eq((Mod(Gaussian(43/3, 97/5), 127)**(-11) * Mod(Gaussian(43/3, 97/5), 127))**-9, Mod(Gaussian(81, 89), 127))
assert_eq((3 + Gaussian(4, 5)), Gaussian(7, 5))
assert_eq((3 - Gaussian(4, 5)), Gaussian(-1, -5))
assert_eq((3 * Gaussian(4, 5)), Gaussian(12, 15))
assert_eq((3 / Gaussian(4, 5)), Gaussian(12/41, -15/41))
var params = [
    %n[3, 4, 5, 6],
    %n[3, 4, 5, -2],
    %n[3,-11, 7, 23],
    %n[-9, -4, -1, -4],
    %n[0, -4, 1, 1],
    %n[0, -1, 13, 12],
    %n[5, 1, 7, 1],
    %n[1, 3, 0, 1],
]
params.each_2d {|a,b,c,d|
    var m = (2**64 + 1)
    for n in (-274176, 274176) {
        var x = powmod(Gaussian(a*d, b*c), n, m)
        var y = powmod(b*d, -n, m)
        var r1 = (x * y)%m
        var r2 = powmod(Gaussian(a/b, c/d), n, m)
        var r3 = Mod(Gaussian(a/b, c/d), m)**n
        var r4 = Mod(Gaussian(a/b, c/d), m)**(-n)
        say "Gaussian(#{a}/#{b}, #{c}/#{d})^#{n} == #{r1} (mod m)"
        assert_eq(r1, r2)
        assert_eq(r3.lift, r1)
        assert(r1 == r2)
        assert(r1 == r3.lift)
        assert_eq(r3 * r4 -> lift, Gaussian(1, 0))
        if (r1 != Gaussian(1,0)) {
            assert_eq(gcd(r1.re-1, m), 274177)
            assert_eq(gcd(r2.re-1, m), 274177)
        }
    }
}
func is_gaussian_quadratic_pseudoprime(n) {
    return false if (n <= 1)
    return true  if (n <= 3)
    static x = Quadratic(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(%n[88561,107185,162401,221761,226801,334153,410041,665281,825265,1569457,1615681,2727649].all(is_gaussian_quadratic_pseudoprime))
assert(%n[80375707,154287451,267559627,326266051,478614067,573183451,643767931,2433943891,4297753027].all(is_gaussian_quadratic_pseudoprime))
func is_quadratic_pseudoprime (n, r=2) {
    return false if (n <= 1)
    return true  if (n <= 3)
    return true if (r <= 0)
    var x = Quadratic(r, 1, r+2).powmod(n, n)
    x.a == r || return false
    var y = Quadratic(r, -1, r+2).powmod(n, n)
    y.a == r || return false
    (x.b + y.b == n) && __FUNC__(n, r-1)
}
assert(is_quadratic_pseudoprime(43))
assert(is_quadratic_pseudoprime(97))
with (Quadratic(1, 1, 2)) {|q|
    assert_eq(15.of { q.pow(_).a }, %n[1, 1, 3, 7, 17, 41, 99, 239, 577, 1393, 3363, 8119, 19601, 47321, 114243])      # A001333
    assert_eq(15.of { q.pow(_).b }, %n[0, 1, 2, 5, 12, 29, 70, 169, 408, 985, 2378, 5741, 13860, 33461, 80782])        # A000129
}
with (Quadratic(1, 1, 3)) {|q|
    assert_eq(15.of { q.pow(_).a }, %n[1, 1, 4, 10, 28, 76, 208, 568, 1552, 4240, 11584, 31648, 86464, 236224, 645376]) # A026150
    assert_eq(15.of { q.pow(_).b }, %n[0, 1, 2, 6, 16, 44, 120, 328, 896, 2448, 6688, 18272, 49920, 136384, 372608])    # A002605
}
var n = (274177-1)
var m = (2**64 + 1)
with (Quadratic(3, 4, 2)) {|q|
    var r = q.powmod(n, m)
    assert_eq(gcd(r.a-1, m), 274177)
    assert_eq(gcd(r.b, m), 274177)
}
do {
    func f(n) {
        var t = ((((Quadratic(Gauss(0, 1), 1, Gauss(-1, 4))**(n+1))) - ((-Quadratic(Gauss(0, -1), 1, Gauss(-1, 4)))**(n+1))) / 2**(n+1) / Quadratic(0, 1, Gauss(-1, 4)))
        t.norm.norm.isqrt
    }
    # OEIS: A105309
    assert_eq(20.of(f), %n[1, 1, 2, 5, 9, 20, 41, 85, 178, 369, 769, 1600, 3329, 6929, 14418, 30005, 62441, 129940, 270409, 562725])
    assert_eq(f(1000), 953007415142879751205548648717471976378723433989590065332724291393966564557241346033182551724015717399625822029716392984323783678994296947484979558415164662370712797688946091924852312992157861694515487367535997242571046016868798763497148464793595996730200474727493070930998306358829398004034629488647010035263698492841)
}
do {
    var phi   = Quadratic(1, 1, 5)/2
    var ihp   = 1-phi
    var sqrt5 = Quadratic(0, 1, 5)
    assert_eq(ihp, -1/phi)
    assert_eq(ihp, Quadratic(1, -1, 5)/2)
    func f(n) {     # n-th Fibonacci number
        (phi**n - ihp**n) / sqrt5
    }
    assert_eq(20.of(f), 20.of { .fib })
    assert_eq(f( 999),  999.fib)
    assert_eq(f(1000), 1000.fib)
}
assert_eq(  # OEIS: A006495
    with (Quadratic(1, 2, -1)) { |q| 20.of { q**_ -> a } },
    %n[1, 1, -3, -11, -7, 41, 117, 29, -527, -1199, 237, 6469, 11753, -8839, -76443, -108691, 164833, 873121, 922077, -2521451]
)
do {  # OEIS: A066408
    #var ω = Quadratic(-1, 1, -3)/2
    var ω = (-1/2 + sqrtQ(-3)/2)
    assert(%n[2, 5, 7, 11, 17, 19, 79, 163, 193, 239, 317, 353, 659, 709, 1049, 1103].all {|n|
        (1 - ω)**n - 1 -> norm.is_prime
    })
}
do {
    var ω = Quadratic(-1, 1, -3)
    with (1 - 1*ω) {|q|
        assert_eq(  # OEIS: A213421
            20.of {|n| q**n -> a },
            %n[1, 2, 1, -10, -47, -118, -143, 254, 2017, 6290, 11041, 134, -76751, -307942, -694511, -622450, 2371777, 13844258, 38774593, 58188566],
        )
        assert_eq(  # OEIS: A168175
            20.of {|n| q**(n+1) -> b.neg },
            %n[1, 4, 9, 8, -31, -180, -503, -752, 513, 7316, 25673, 51480, 26209, -255524, -1205559, -3033568, -3695359, 6453540, 51681673, 161551912],
        )
    }
}
do {
    var ω = Quadratic(-1, 1, Gauss(1, Quadratic(0, 1, 1)))
    with (2 + 1*ω) {|q|
        assert_eq(  # OEIS: A138766
            20.of {|n| q**n -> real.real },
            %n[1, 1, 2, 4, 7, 11, 14, 8, -31, -167, -558, -1572, -4025, -9645, -21922, -47536, -98431, -193935, -360094, -617100]
        )
    }
}
do {
    func a(n) {  # OEIS: A109516
        var s = Quadratic(0, 1, (n-1)*(n+3))
        ((n + s - 1)**n - (n - s - 1)**n) / (2**n * s) -> a
    }
    assert_eq(
        a.map(2..15),
        %n[1, 6, 45, 464, 6000, 93528, 1707111, 35721216, 843160671, 22165100000, 642268811184, 20339749638144, 698946255836933, 25903663544572800]
    )
    assert_eq(a(100), 980535248623862021263601221274544258060956062594742050353606992495209528021598427086799187928951344994564615926105529217544063038419696049119411195562750041858776047253609373010165706844861440666701)
}
do {
    func a(n) {  # OEIS: A097691
        var q = Quadratic(0, 1, 4 - n*n)
        var i = Gauss(0, 1)
        ((2**(-n) * (q + i*n)**n) - (2**n * (-q - i*n)**(-n))) / q -> a.norm.isqrt
    }
    assert_eq(
        a.map(3..20),
        %n[8, 56, 551, 6930, 105937, 1905632, 39424240, 922080050, 24057287759, 692686638072, 21817946138353, 746243766783074, 27543862067299424, 1091228270370045824, 46187969968474139807, 2080128468827570457762, 99318726126650358502921, 5011361251329169946919800]
    )
    assert_eq(a(100), 990246417437806351190277754425603307809612897254870471251275960965725797918060784184076057065314451820651920841907304582494680760348836023974638203217830928705082772363290340579154282010420824995000)
}
do {
    var a = Quadratic(3,4,5)
    var b = Quadratic(7,11,43)
    assert_eq(a + b, Quadratic(a.a + b, a.b, a.w))
    assert_eq(a - b, Quadratic(a.a - b, a.b, a.w))
    assert_eq(a * b, Quadratic(a.a * b, a.b * b, a.w))
    assert_eq(a / b, Quadratic(a.a / b, a.b / b, a.w))
    assert_eq(a ^ b, Quadratic(a.a ^ b, a.b, a.w))
    assert_eq(a | b, Quadratic(a.a | b, a.b, a.w))
    assert_eq(a & b, Quadratic(a.a & b, a.b, a.w))
    assert_eq(a <=> b, -1)
    assert_eq(b <=> a, +1)
    assert_eq(a.to_n <=> b.to_n, -1)
}
do {
    assert_eq(43 % Quadratic(97, 120, 19), Quadratic(43, 0, 19) % Quadratic(97, 120, 19))
    assert_eq(Str(Quadratic(Quadratic(3,5),4).to_n), "12")
}
say "** Test passed!"