# Lucas sequence implementation

Could you suggest any improvements in this Lucas sequence's implementation:

from mpmath.libmp import bitcount as _bitlength

def _int_tuple(*i):
return tuple(int(_) for _ in i)

def _lucas_sequence(n, P, Q, k):
"""Return the modular Lucas sequence (U_k, V_k, Q_k).
Given a Lucas sequence defined by P, Q, returns the kth values for
U and V, along with Q^k, all modulo n.
"""
D = P*P - 4*Q
if n < 2:
raise ValueError("n must be >= 2")
if k < 0:
raise ValueError("k must be >= 0")
if D == 0:
raise ValueError("D must not be zero")

if k == 0:
return _int_tuple(0, 2, Q)
U = 1
V = P
Qk = Q
b = _bitlength(k)
if Q == 1:
# For strong tests
while b > 1:
U = (U*V) % n
V = (V*V - 2) % n
b -= 1
if (k >> (b - 1)) & 1:
t = U*D
U = U*P + V
if U & 1:
U += n
U >>= 1
V = V*P + t
if V & 1:
V += n
V >>= 1
elif P == 1 and Q == -1:
# For Selfridge parameters
while b > 1:
U = (U*V) % n
if Qk == 1:
V = (V*V - 2) % n
else:
V = (V*V + 2) % n
Qk = 1
b -= 1
if (k >> (b-1)) & 1:
t = U*D
U = U + V
if U & 1:
U += n
U >>= 1
V = V + t
if V & 1:
V += n
V >>= 1
Qk = -1
else:
# The general case with any P and Q
while b > 1:
U = (U*V) % n
V = (V*V - 2*Qk) % n
Qk *= Qk
b -= 1
if (k >> (b - 1)) & 1:
t = U*D
U = U*P + V
if U & 1:
U += n
U >>= 1
V = V*P + t
if V & 1:
V += n
V >>= 1
Qk *= Q
Qk %= n
U %= n
V %= n
return _int_tuple(U, V, Qk)

• The indentation of the last line is plain wrong. Can you check that this is the only one? Also, can you add the formal definition of this sequence? – Mathias Ettinger Sep 11 '16 at 13:37
• Yes, my bad. And, you can read more about Lucas Sequence here : Lucas Sequence – Mike Sep 12 '16 at 13:01
• I'm not familiar with python, so forgive my ignorance, but where is the else to the if (k >> (b - 1)) & 1: statement? If the next bit is 0, there is still an operation you need to perform. Can you confirm that the code outputs correct results in its current state? – Myridium Sep 27 '16 at 12:16

## 3 Answers

    if D == 0:
raise ValueError("D must not be zero")


D is not a parameter, so it seems a bit strange to me to use its name in the exception error message.

Also, given that the reference you give defines Lucas sequences only for D > 0, it seems a bit strange that you're happy to have D < 0. It might be worth adding a comment to say that you're supporting negative determinants.

        # For strong tests


I deduce from this that you're using the sequences for a Lucas or Lucas-Lehmer primality test, but that comment seems a bit out of place.

Having three cases also seems a bit overkill. Have you profiled the code to see how much of a difference the special cases make? You seem to mainly save a few multiplications by 1, and I'm not convinced that the maintenance effort is a worthwhile tradeoff if it's only a 2% speedup in a function which probably won't be the bottleneck anyway.

b = _bitlength(k)
...
while b > 1:
...
b -= 1
if (k >> (b - 1)) & 1:
...


I had to test this code to be convinced that you hadn't got the bit-shifting backwards (i.e. that you really did want to work from the most-significant bit down to the least-significant). It does seem correct, but it's definitely too clever to go without a comment explaining why it's correct.

                t = U*D
U = U*P + V
if U & 1:
U += n
U >>= 1
V = V*P + t
if V & 1:
V += n
V >>= 1


One of the great features of Python which I wish more languages had is parallel assignment. I would ditch the temporary variable and expose the parallels between U and V as

                U, V = U*P + V, V*P + U*D
if U & 1:
U += n
if V & 1:
V += n
U, V = U >> 1, V >> 1


    U %= n
V %= n
return _int_tuple(U, V, Qk)


In my opinion it's not putting too much into a line to make that (following Graipher's suggestion)

    return (U % n, V % n, Qk)

• I can't answer for the OP code, but in a C implementation, the P=1,Q=-1 case comes out about 30% faster than the generic case, and is hit by about 50% of inputs to Lucas/Strong Lucas test. Generic: 3 mulmod + 4/odd bit. P=1,Q=-1: 2 mulmod + 1/odd bit. Q=1: 2 mulmod + 3/odd bit. (OP's code matches this also). The Q=1 case is used for the extra strong Lucas test. – DanaJ Sep 23 '16 at 15:47

Your _int_tuple function does not seem to serve any real purpose. U and V are always integers (you initialize them to 1 and do only arithmetic on them afterwards). And what exactly your program will do with a Q which is e.g. a string is hard to follow, but I think it will run into the else block (Q is neither 1 nor -1) and then try to do Qk = Q, followed by Qk *= Qk which fails for almost everything but number types.

I also can't see any parts of your code which would produce a float (though I might have overlooked one).

So I would get rid of that function and return just a tuple and see what it does.

In general it looks good to me. I'd clarify the comments a bit more. Q=1 is typically used for the extra strong Lucas test, P=1,Q=-1 is used by about half the cases for the standard or strong Lucas test (the rest have a different Q). Some info on the various tests can be seen on my Pseudoprime Statistics page.

What happens if n is even? For example, Lucas_5(6,1) = (1189,6726). Mod 1000 should yield (189,726). But the routine gives us (689,726). This is an artifact of the halving method used. My solution was to use a different, slower, method in that case, since it's rarely used (never by a standard primality test, but there are more uses that just those).

For D=0, we can return a solution. See the Wikipedia page.

Check the cases where P and/or Q and/or D exceeds n. It looks from a quick test that the returned Q_k values are out of range for k={0,1}. Doing a check and mod up front may improve performance.

Consider k.bit_length() instead of using libmp. Up to you, but I like fewer dependencies. I don't know all the tradeoffs between the two however.

For the Q=1 case (used by the extra strong Lucas test) there is yet another method that can be faster. If D can be inverted mod n, then we can calculate V_k with 2 mulmods per bit, then compute U_k using the inverse. I hate to add more cases, and of course you'd need to benchmark to see if it worked better for your infrastructure and use.