#!/usr/bin/ruby

# Tests for the "Gauss" class.

var r = (-2 .. 2)

for a in (r), b in (r), c in (r), d in (r) {

    assert_eq([Complex(a,b) + Complex(c,d) -> reals], [cadd(a,b,c,d)])
    assert_eq([Complex(a,b) - Complex(c,d) -> reals], [csub(a,b,c,d)])
    assert_eq([Complex(a,b) * Complex(c,d) -> reals], [cmul(a,b,c,d)])

    if (c*c + d*d != 0) {
        assert_eq([Complex(a,b) / Complex(c,d) -> reals], [cdiv(a,b,c,d)].map{.float})
    }
}

for a in (r), b in (r) {

    var n = irand(0, 100)
    var m = irand(100, 1000)

    assert_eq([Complex(a,b)**n -> reals], [cpow(a,b,n)])
    assert_eq([Complex(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(Gauss(1,0) / Gauss(a, b), Gauss(a,b).inv)
    assert_eq(Mod(Gauss(a, b), m).inv * Gauss(a,b), Mod(Gauss(1,0), m))
    assert_eq(Mod(Gauss(a, b), m).inv * Gauss(a,b) -> lift, Gauss(1,0))
    assert_eq(Mod(Gauss(a, b), m)**(-1) * Gauss(a,b) -> lift, Gauss(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(Gauss(3,4) == Gauss(3,4))
assert(!(Gauss(3,4) == Gauss(3,3)))
assert(!(Gauss(3,3) == Gauss(3,4)))

assert(Gauss(3,3) != Gauss(3,4))
assert(Gauss(3,4) != Gauss(3,3))
assert(!(Gauss(3,4) != Gauss(3,4)))

assert_eq(Gauss(4,5) <=> Gauss(3,4), 1)
assert_eq(Gauss(3,4) <=> Gauss(3,4), 0)
assert_eq(Gauss(3,4) <=> Gauss(3,5), -1)

assert(Gauss(4,5) >  Gauss(3,4))
assert(Gauss(4,5) >  Gauss(4,4))
assert(Gauss(4,5) >= Gauss(3,4))
assert(Gauss(3,4) >= Gauss(3,4))
assert(Gauss(3,4) <= Gauss(3,4))

assert(!(Gauss(3,4) > Gauss(3,5)))
assert(!(Gauss(3,4) > Gauss(4,5)))
assert(!(Gauss(3,4) < Gauss(3,4)))
assert(!(Gauss(3,4) < Gauss(3,1)))
assert(!(Gauss(3,4) < Gauss(2,1)))
assert(!(Gauss(3,4) <= Gauss(3,3)))
assert(!(Gauss(2,1) >= Gauss(2,2)))

for a in (r), b in (r) {
    var y = irand(2, 100)
    var x = Gauss(a,b)
    assert_eq(x - floor(x/y)*y, x % y)
}

func gaussian_sum_2(n) {

    var i     = Gauss(0, 1)
    var total = Gauss(0)

    for k in (1..n) {
        total += (i**(k-1) / k)
    }

    total * n!
}

assert_eq(
    10.of(gaussian_sum_2),
    [Gauss(0, 0), Gauss(1, 0), Gauss(2, 1), Gauss(4, 3), Gauss(16, 6), Gauss(104, 30), Gauss(624, 300), Gauss(3648, 2100), Gauss(29184, 11760), Gauss(302976, 105840)]
)

assert_eq(powmod(Gauss(3,4), 1000, 1e6), Gauss(585313, 426784))
assert_eq([Mod(Gauss(3,4), 1e6)**1000 -> lift.reals], [585313, 426784])
assert_eq(Mod(Gauss(3,4), 1e6)**1000, Mod(Gauss(585313, 426784), 1e6))
assert(Mod(Gauss(3,4), 1e6)**1000 == Mod(Gauss(585313, 426784), 1e6))

assert_eq(Mod(43, 97).to_n, 43)
assert_eq(Gauss(3,4).to_n, 3+4i)
assert_eq(Gauss(3,4).to_c, 3+4i)

assert_eq(Gauss(42).invmod(2017), Gauss(1969, 0))
assert_eq(Gauss(3,4).invmod(2017), Gauss(1291, 968))
assert_eq(Gauss(91,23).invmod(2017), Gauss(590, 405))
assert_eq(Gauss(43, 99).invmod(2017), Gauss(1709,1272))
assert_eq(Gauss(43, 99).invmod(1234567), Gauss(1019551, 667302))

assert_eq(Mod(Gauss(42), 2017).inv, Mod(Gauss(1969, 0), 2017))
assert_eq(Mod(Gauss(3,4), 2017)**(-1), Mod(Gauss(1291, 968), 2017))
assert_eq(Mod(Gauss(91,23), 2017)**(-1), Mod(Gauss(590, 405), 2017))
assert_eq(Mod(Gauss(43, 99), 2017).inv, Mod(Gauss(1709,1272), 2017))
assert_eq(Mod(Gauss(43, 99), 1234567)**(-2), Mod(Gauss(1019551, 667302)**2, 1234567))
assert_eq(Mod(Gauss(43, 99), 1234567)**(-5), Mod(Gauss(1019551, 667302)**5, 1234567))

assert_eq(powmod(Gauss(43, 99), -4, 1234567), invmod(Gauss(43, 99)**4 % 1234567, 1234567))
assert_eq(powmod(Gauss(43, 99), -5, 1234567), invmod(Gauss(43, 99)**5 % Gauss(1234567), 1234567))
assert_eq(powmod(Gauss(43, 99), -5, 1234567), invmod(Gauss(43, 99)**5, 1234567))

assert_eq(powmod(Gauss(43, 99), -4, 1234567), invmod(powmod(Gauss(43, 99), 4, 1234567), 1234567))
assert_eq(powmod(Gauss(43, 99), -5, 1234567), invmod(powmod(Gauss(43, 99), 5, 1234567), 1234567))

assert_eq(Gauss(43, 97)**(-5), (Gauss(43,97)**5)**(-1))
assert_eq(Gauss(43, 97)**(-5), (Gauss(43,97)**5).inv)
assert_eq(Gauss(43, 97)**(-5), (Gauss(43,97).inv)**5)
assert_eq(Gauss(43, 97)**(-5), (Gauss(43,97)**(-1))**5)

assert_eq(Mod(Gauss(43, 97), 1234567)**(-5),  Mod(Gauss(43, 97), 1234567)**5 -> inv)
assert_eq(Mod(Gauss(43, 97), 1234567)**(-5), (Mod(Gauss(43, 97), 1234567)->inv)**5)

assert_eq(Mod(Gauss(43, 97), 1234567)**(-1234),  Mod(Gauss(43, 97), 1234567)**1234 -> inv)
assert_eq(Mod(Gauss(43, 97), 1234567)**(-1234), (Mod(Gauss(43, 97), 1234567)->inv)**1234)
assert_eq(Mod(Gauss(3,4), 1234567)**1234 -> lift, powmod(Gauss(3,4), 1234, 1234567))
assert_eq(Mod(Gauss(3,4), 1234567)**-1234 -> lift, powmod(Gauss(3,4), -1234, 1234567))

assert_eq(Gauss(3/5,11/4)**(-27), Gauss(3/5,11/4)**27 -> inv)
assert_eq(Gauss(3/5,11/4)**(-27), Gauss(3/5,11/4).inv**27)

assert_eq(Complex(Gauss(3,4), Gauss(17,19)), -16 + 21i)
assert_eq(Complex(Mod(13, 97), Mod(43, 97)), 13 + 43i)
assert_eq(Gauss(Mod(13, 97), Mod(43, 97)), Gauss(Mod(13, 97), Mod(43, 97)))
assert_eq(Mod(Gauss(3, 4), 97)*1234, Mod(Gauss(16, 86), 97))
assert_eq(Mod(Gauss(3/4, 5/6), 1234567)**10 * Mod(Gauss(3/4, 5/6), 1234567)**-10, Mod(Gauss(1,0), 1234567))

assert_eq((Mod(Gauss(43/3, 97/5), 127)**(-11) * Mod(Gauss(43/3, 97/5), 127))**+9, Mod(Gauss(52, 73), 127))
assert_eq((Mod(Gauss(43/3, 97/5), 127)**(-11) * Mod(Gauss(43/3, 97/5), 127))**-9, Mod(Gauss(81, 89), 127))

assert_eq((3 + Gauss(4, 5)), Gauss(7, 5))
assert_eq((3 - Gauss(4, 5)), Gauss(-1, -5))
assert_eq((3 * Gauss(4, 5)), Gauss(12, 15))
assert_eq((3 / Gauss(4, 5)), Gauss(12/41, -15/41))

var params = [
    [3, 4, 5, 6],
    [3, 4, 5, -2],
    [3,-11, 7, 23],
    [-9, -4, -1, -4],
    [0, -4, 1, 1],
    [0, -1, 13, 12],
    [5, 1, 7, 1],
    [1, 3, 0, 1],
]

params.each_2d {|a,b,c,d|

    var m = (2**64 + 1)

    for n in (-274176, 274176) {

        var x = powmod(Gauss(a*d, b*c), n, m)
        var y = powmod(b*d, -n, m)

        var r1 = (x * y)%m
        var r2 = powmod(Gauss(a/b, c/d), n, m)
        var r3 = Mod(Gauss(a/b, c/d), m)**n
        var r4 = Mod(Gauss(a/b, c/d), m)**(-n)

        say "Gauss(#{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, Gauss(1, 0))

        if (r1 != Gauss(1,0)) {
            assert_eq(gcd(r1.re-1, m), 274177)
            assert_eq(gcd(r2.re-1, m), 274177)
        }
    }
}

do {
    var a = Gauss(3,4)**10
    var b = Gauss(3,4).powmod(97, 1234)

    assert_eq([a.reals], [reals(Gauss(3,4)**10)])
    assert_eq([b.reals], [reals(Gauss(3,4).powmod(97, 1234))])
    assert_eq([reals(a+b)], [reals(Gauss(reals(a)) + Gauss(reals(b)))])
    assert_eq([reals(a-b)], [reals(Gauss(reals(a)) - Gauss(reals(b)))])

    assert_eq([a.reals], [reals(Quadratic(3,4,-1)**10)])
    assert_eq([b.reals], [reals(Quadratic(3,4,-1).powmod(97, 1234))])
    assert_eq([reals(a+b)], [reals(Quadratic(reals(a),-1) + Quadratic(reals(b),-1))])
    assert_eq([reals(a-b)], [reals(Quadratic(reals(a),-1) - Quadratic(reals(b),-1))])
}

assert_eq(  # OEIS: A006495
    with (Gauss(1, 2)) {|g| 20.of { g**_ -> real } },
    %n[1, 1, -3, -11, -7, 41, 117, 29, -527, -1199, 237, 6469, 11753, -8839, -76443, -108691, 164833, 873121, 922077, -2521451]
)

do {
    func lucasUVmod(P, Q, n, m) {

        var D = (P*P - 4*Q)

        var a = Quadratic(P, +1, D)/2
        var b = Quadratic(P, -1, D)/2

        var x = Mod(a, m)**n
        var y = Mod(b, m)**n

        ((x-y).lift.b, (x+y).lift.a)
    }

    var P = 43
    var Q = -11
    var n = 10000000
    var m = 1e9+1

    assert_eq([lucasUVmod(P, Q, n, m)], [::lucasUVmod(P, Q, n, m)])
}

do {
    func lucasUVmod(P, Q, n, m) {

        var D = (P*P - 4*Q)

        var a = Quadratic(Mod(P, m), +1, D)/2
        var b = Quadratic(Mod(P, m), -1, D)/2

        var x = a**n
        var y = b**n

        ((x-y).b.lift, (x+y).a.lift)
    }

    var P = 97
    var Q = 11
    var n = 10000
    var m = 1e9-1

    assert_eq([lucasUVmod(P, Q, n, m)], [::lucasUVmod(P, Q, n, m)])
}

do {
    assert_eq(43 % Gauss(13, 24), Gauss(43,0) % Gauss(13, 24))
    assert_eq(Gauss(Mod(Poly(3,4), 127), 923).lift, Gauss(Polynomial(3 => 4), 923))
}

assert_eq(Gauss(5040).factor.len, 14)
assert_eq(Gauss(0, 5040).factor.len, 13)

assert_eq(Gauss(5040).factor.prod, Gauss(5040))
assert_eq(Gauss(-5040).factor.prod, Gauss(-5040))
assert_eq(Gauss(0, 5040).factor.prod, Gauss(0, 5040))
assert_eq(Gauss(0, -5040).factor.prod, Gauss(0, -5040))

assert_eq(Gauss(38, 60).factor.len, 5)
assert_eq(Gauss(38, 60).factor.prod, Gauss(38, 60))

assert_eq(Gauss(5040).factor_exp.len, 6)
assert_eq(Gauss(0, 5040).factor_exp.len, 5)
assert_eq(Gauss(5040).factor_exp.prod_2d{|p,e| p**e }, Gauss(5040))
assert_eq(Gauss(0, 5040).factor_exp.prod_2d{|p,e| p**e }, Gauss(0, 5040))

assert_eq(Gauss(-38, -60).factor_exp.prod_2d{|p,e| p**e }, Gauss(-38, -60))
assert_eq(Gauss( 38, -60).factor_exp.prod_2d{|p,e| p**e }, Gauss( 38, -60))
assert_eq(Gauss(-38,  60).factor_exp.prod_2d{|p,e| p**e }, Gauss(-38,  60))

assert_eq(20.of{ Gauss(_).factor.prod.to_n }, 20.of { _ })

assert_eq(divisors(Gauss(0,-5040)).len, 864)
assert({|t| t.divisors.all { t.is_div(_) } }(Gauss(-38, 60)))
assert({|t| t.divisors.all { t.is_div(_) } }(Gauss(-41)))
assert({|t| t.divisors.all { t.is_div(_) } }(Gauss(0, 97)))

assert(Gauss(0, 100.primorial*503*863*43*97).factor.all{.is_prime})

say "** Test passed!"