up vote 3 down vote favorite

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)
share|improve this question
3  
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
1  
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
up vote 2 down vote accepted 
    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)
share|improve this answer
    
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
up vote 2 down vote

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.

share|improve this answer
up vote 2 down vote

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.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.