123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869 |
- #!/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)
- }
|