#!/usr/bin/ruby

# Tests for Lucas sequences U_n(P, Q) and V_n(P, Q).

func is(Number a, b) {

    if (a.is_nan) {
        return __FUNC__(a.to_s, b.to_s)
    }

    assert_eq(a, b)
}

func is(a, b){
    assert_eq(a, b)
}

is(lucas(-1),  NaN)
is(lucas(0),     2)
is(lucas(1),     1)
is(lucas(15), 1364)
is(lucas(15), 1364)
is(lucas(-4),  NaN)

is(lucasU(4, 4, 13), 53248)
is(lucasV(4, 4, 13), 16384)

is(lucasV(-4, 4,  24), 33554432)
is(lucasV(4,  -4, 24), 25783171895263232)
is(lucasV(4,  -4, 25), 124492166541082624)
is(lucasV(-4, 4,  25), -67108864)
is(lucasV(-4, 4,  26), 134217728)
is(lucasV(-4, -4, 16), 87275208704)
is(lucasV(-4, -4, 48), 664771952980691802270648617664512)

is(lucasU(4,  -4, 24), 4557863921909760)
is(lucasU(-4, 4,  24), -201326592)
is(lucasU(4,  -4, 25), 22007313791451136)
is(lucasU(-4, 4,  25), 419430400)
is(lucasU(-4, 4,  26), -872415232)
is(lucasU(-4, -4, 16), -15428222976)
is(lucasU(-4, -4, 48), -117516188973817974394087349944320)

is(lucasUmod(-4, -4, 48,    700237709),  424976705)
is(lucasUmod(4,  -4, 4171,  619442071),  93411036)
is(lucasUmod(4,  -4, 54223, 1238884142), 793124222)

is(lucasUmod(5, 3, 1234, 789347),  134939)
is(lucasUmod(5, 3, 1234, -789347), 134939)   # maybe this should be: 134939 - 789347 = -654408

is({ lucasUmod(3, -2, _, 24) }.map(0 .. 9), %n[0, 1, 3, 11, 15, 19, 15, 11, 15, 19])
is({ lucasVmod(3, -2, _, 24) }.map(0 .. 9), %n[2, 3, 13, 21, 17, 21, 1, 21, 17, 21])

is({ [lucasUVmod(3, -2, _, 24)] }.map(0 .. 9),
    %n[0, 1, 3, 11, 15, 19, 15, 11, 15, 19] ~Z
    %n[2, 3, 13, 21, 17, 21, 1, 21, 17, 21])

is({ lucasUmod(4, 4, _, 25) }.map(0 .. 9), %n[0, 1, 4, 12, 7, 5, 17, 23, 24, 4])
is({ lucasVmod(4, 4, _, 25) }.map(0 .. 9), %n[2, 4, 8, 16, 7, 14, 3, 6, 12, 24])

is({ [lucasUVmod(4, 4, _, 25)] }.map(0 .. 9),
    %n[0, 1, 4, 12, 7, 5, 17, 23, 24, 4] ~Z
    %n[2, 4, 8, 16, 7, 14, 3, 6, 12, 24])

is({ lucasUmod(-5, -3, _, 48) }.map(0 .. 9), %n[0, 1, 43, 28, 37, 43, 40, 25, 43, 4])
is({ lucasVmod(-5, -3, _, 48) }.map(0 .. 9), %n[2, 43, 31, 22, 31, 7, 10, 19, 31, 46])

is({ [lucasUVmod(-5, -3, _, 48)] }.map(0 .. 9),
    %n[0, 1, 43, 28, 37, 43, 40, 25, 43, 4] ~Z
    %n[2, 43, 31, 22, 31, 7, 10, 19, 31, 46])

is({ lucasUmod( 5, 3, _, 789347) }.map(0 .. 14), %n[0, 1, 5, 22, 95, 409, 1760, 7573, 32585, 140206, 603275, 227716, 118102, 696709, 761198])
is({ lucasUmod(-5, 3, _, 789347) }.map(0 .. 14), %n[0, 1, 789342, 22, 789252, 409, 787587, 7573, 756762, 140206, 186072, 227716, 671245, 696709, 28149])

is({ lucasUmod( 5, 3, _, 789347) }.map(0 .. 14), { lucasU( 5, 3, _) % 789347 }.map(0 .. 14))
is({ lucasUmod(-5, 3, _, 789347) }.map(0 .. 14), { lucasU(-5, 3, _) % 789347 }.map(0 .. 14))

is({ lucasVmod( 5, 3, _, 789347) }.map(0 .. 14), { lucasV( 5, 3, _) % 789347 }.map(0 .. 14))
is({ lucasVmod(-5, 3, _, 789347) }.map(0 .. 14), { lucasV(-5, 3, _) % 789347 }.map(0 .. 14))

is(fibmod(0, 100), 0)
is(fibmod(1, 100), 1)
is(lucasmod(0, 100), 2)
is(lucasmod(1, 100), 1)

is(lucasUmod(1, -1, 0, 100), 0)
is(lucasUmod(1, -1, 1, 100), 1)

is(lucasVmod(1, -1, 0, 100), 2)
is(lucasVmod(1, -1, 1, 100), 1)

is(lucasUmod(4, 4, 0, 100), 0)
is(lucasUmod(4, 4, 1, 100), 1)

is(lucasUmod(4,  4, 50, 1000), lucasU(4,  4, 50) % 1000)
is(lucasUmod(4, -4, 50, 1000), lucasU(4, -4, 50) % 1000)

is(lucasVmod(4, 4, 0, 100), 2)
is(lucasVmod(4, 4, 1, 100), 4)

is(lucasVmod(4, 4, 0, 100), lucasV(4, 4, 0) % 100)
is(lucasVmod(4, 4, 1, 100), lucasV(4, 4, 1) % 100)

is(lucasVmod( 4, 4, 50, 1000), 248)
is(lucasVmod(-4, 4, 50, 1000), 248)

is(lucasVmod( 4, 4, 50, 1000), lucasV(4, 4, 50) % 1000)
is(lucasVmod(-4, 4, 50, 1000), lucasV(4, 4, 50) % 1000)

is(lucasVmod( 4, 4, 50, 1001), lucasV( 4, 4, 50) % 1001)
is(lucasVmod(-4, 4, 50, 1001), lucasV(-4, 4, 50) % 1001)

is(lucasUmod( 4, 4, 50, 1001), lucasU( 4, 4, 50) % 1001)
is(lucasUmod(-4, 4, 50, 1001), lucasU(-4, 4, 50) % 1001)

is(lucasVmod( 4, -4, 50, 1000), 624)
is(lucasVmod(-4, -4, 50, 1000), 624)

is(lucasVmod( 4, -4, 50, 1000), lucasV( 4, -4, 50) % 1000)
is(lucasVmod(-4, -4, 50, 1000), lucasV(-4, -4, 50) % 1000)

is(lucasVmod(3, -4, 50, 1000), lucasV(3, -4, 50) % 1000)
is(lucasUmod(3, -4, 50, 1000), lucasU(3, -4, 50) % 1000)

is(fibmod(1234, 987654), fibonacci(1234) % 987654)
is(lucasmod(1234, 987654), lucas(1234) % 987654)

is(lucasUmod(1, -1, -43, 42), NaN)
is(lucasUmod(1, -1, -2, 43), NaN)

is(lucasVmod(1, -1, -43, 42), NaN)
is(lucasVmod(1, -1, -2, 43), NaN)

is([lucasUVmod(1, -1, -41, 97)].join(' '), [NaN, NaN].join(' '))
is([lucasUVmod(1, -1, -12, 42)].join(' '), [NaN, NaN].join(' '))

is(fibmod(-1, 1234), NaN)
is(lucasmod(-1, 1234), NaN)

is(fibmod(42, 0), NaN)
is(lucasmod(42, 0), NaN)

is(lucasUmod(2**95 + 1, 2**75 - 1, 1234, 98765), 35621)
is(lucasVmod(2**95 + 1, 2**75 - 1, 1234, 98765), 3022)

is(lucasUmod(2**95 + 1, 2**75 - 1, 98765, 1234), 477)
is(lucasVmod(2**95 + 1, 2**75 - 1, 98765, 1234), 29)

is([lucasUVmod(2**95 + 1, 2**75 - 1, 1234, 98765)], [35621, 3022])
is([lucasUVmod(2**95 + 1, 2**75 - 1, 98765, 1234)], [477, 29])

func pisano_period_pp(P, Q, p, k=1) {
    (p - kronecker(P*P - 4*Q, p)).divisors.first_by {|d| lucasUmod(P, Q, d, p) == 0 } * p**(k-1)
}

func pisano_period(n, P=1, Q=-1) {

    return 0 if (n <= 0)
    return 1 if (n == 1)

    var d = n.factor_map {|p,k| pisano_period_pp(P, Q, p, k) }.lcm

    3.times {|k|
        var t = d<<k
        if ((lucasUmod(P, Q, t, n) == 0) && (lucasUmod(P, Q, t+1, n) == 1)) {
            return t
        }
    }

    die "Error for n = #{n} with P = #{P} and Q = #{Q}"
}

{|n| pisano_period(n, 1, -1) }.map(1..30)       # Periods of Fibonacci numbers
{|n| pisano_period(n, 2, -1) }.map(1..30)       # Periods of Pell numbers
{|n| pisano_period(n, 3, -1) }.map(1..30)       # Periods of 3-Fibonacci numbers

for k in (0..10) {
    var m = 43*97
    assert_eq(20.of { .fibmod(k, m) }, 20.of { .fib(k) % m })
}

assert_eq(fibmod(43, 3, 2**128 + 1), fib(43, 3) % (2**128 + 1))
assert_eq(fibmod(2**128 + 1, 3, 2**128 - 1), 220193395512101899249036009241856381864)
assert_eq(fibmod(2**128 + 1, 3, 2**127 - 1), 21277636584740672870826766395962370082)
assert_eq(fibmod(2**127 - 1, 3, 2**128 + 1), 85447148131207702798335464791538973804)
assert_eq(fibmod(2**127 - 1, 4, 2**128 + 1), 279405613631561682761560761315872506950)
assert_eq(fibmod(2**128 + 1, 2, 2**128 - 1), fibmod(2**128 + 1, 2**128 - 1))

assert(%n[913, 150267335403, 430558874533, 14760229232131, 936916995253453].all { .is_bfw_psp })
assert(%n[7056721, 79397009999, 443372888629441, 582920080863121, 2491924062668039, 14522256850701599, 39671149333495681, 242208715337316001, 729921147126771599, 842526563598720001, 1881405190466524799, 2380296518909971201, 3188618003602886401, 33711266676317630401, 54764632857801026161, 55470688965343048319, 72631455338727028799, 122762671289519184001, 361266866679292635601, 734097107648270852639].all{.is_chebyshev_psp})

assert(is_strong_fibonacci_pseudoprime(443372888629441))
assert(is_strong_fibonacci_pseudoprime(39671149333495681))

say "** Test passed!"