modular_lucas_sequence_V.sf

#!/usr/bin/ruby

# Algorithm due to Aleksey Koval for efficiently computing the modular Lucas V sequence.

# See also:
#   https://en.wikipedia.org/wiki/Lucas_sequence
#   https://file.scirp.org/pdf/IJCNS20101200011_90376712.pdf

func lucasVmod(P, Q, n, m) {

    var(V1 = 2, V2 = P)
    var(Q1 = 1, Q2 = 1)

    var s = n.valuation(2)

    for bit in ((n>>(s+1)).as_bin.chars) {

        Q1 = mulmod(Q1, Q2, m)

        if (bit) {
            Q2 = mulmod(Q1, Q, m)
            V1 = mulsubmulmod(V2, V1, P, Q1, m)
            V2 = mulsubmulmod(V2, V2, 2, Q2, m)
        }
        else {
            Q2 = Q1
            V2 = mulsubmulmod(V2, V1, P, Q1, m)
            V1 = mulsubmulmod(V1, V1, 2, Q2, m)
        }
    }

    Q1 = mulmod(Q1, Q2, m)
    Q2 = mulmod(Q1, Q, m)
    V1 = mulsubmulmod(V2, V1, P, Q1, m)
    Q1 = mulmod(Q1, Q2, m)

    s.times {
        V1 = mulsubmulmod(V1, V1, 2, Q1, m)
        Q1 = mulmod(Q1, Q1, m)
    }

    return V1
}

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

say {|n| lucasVmod(1, -1, n, n) }.map(1..30)        # Lucas(n) modulo n
say 25.by {|n| lucasVmod(1, -1, n, n) == 1 }        # First 25 prime numbers

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

for k in (1..20) {      # run some tests

    var p = 1e30.random_prime

    assert_eq(lucasVmod(1, -1, p, p), 1)

    var n = 1e30.irand
    var m = 1e30.irand

    var P = (100.irand * [-1,1].rand)
    var Q = (100.irand * [-1,1].rand)

    var V1 = lucasVmod(P, Q, n, m)
    var V2 = ::LucasVmod(P, Q, n, m)

    assert_eq(V1, V2)
}