modular_lucas_sequences_U_V.sf

#!/usr/bin/ruby

# Algorithm due to M. Joye and J.-J. Quisquater for computing the Lucas U and V sequences.

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

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

    var(U1 = 1)
    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)
            U1 = mulmod(U1, V2, m)
            V1 = mulsubmulmod(V2, V1, P, Q1, m)
            V2 = mulsubmulmod(V2, V2, 2, Q2, m)
        }
        else {
            Q2 = Q1
            U1 = mulsubmod(U1, V1, Q1, m)
            V2 = mulsubmulmod(V2, V1, P, Q1, m)
            V1 = mulsubmulmod(V1, V1, 2, Q2, m)
        }
    }

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

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

    return (U1, V1)
}

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

var p = 1e6.irand
var q = 1e6.irand
var n = 1e20.irand
var m = 4123412312

var (U, V) = lucasUVmod(p, q, n, m)

say "lucasU(P=#{p}, Q=#{q}, N=#{n}) = #{U} (mod #{m})"
say "lucasV(P=#{p}, Q=#{q}, N=#{n}) = #{V} (mod #{m})"

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

    var p = 1e30.random_prime

    assert_eq([lucasUVmod(1, -1, p, p)][1], 1)
    assert_eq([lucasUVmod(1, -1, p - kronecker(p, 5), p)][0], 0)

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

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

    var (U1, V1) = lucasUVmod(P, Q, n, m)
    var (U2, V2) = (::LucasUmod(P, Q, n, m), ::LucasVmod(P, Q, n, m))
    var (U3, V3) = ::LucasUVmod(P, Q, n, m)

    assert_eq(U1, U2)
    assert_eq(U1, U3)
    assert_eq(V1, V2)
    assert_eq(V1, V3)
}

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