# Lucas sequence - how to know which terms to calculate when testing if a number is prime?

To check if a number is a Lucas pseudoprime, a Lucas sequence is computed. The Lucas sequence is based on a recurrence relation, but there's a method involving inspecting the bits of the number being checked for primatliy to know which terms in the sequence to compute. How exactly does inspecting the bits work?

Wikipedia gives the example

For example, if n+1 = 44 (= 101100 in binary), then, taking the bits one at a time from left to right, we obtain the sequence of indices to compute: $$1_2$$ = 1, $$10_2$$ = 2, $$100_2$$ = 4, $$101_2$$ = 5, $$1010_2$$ = 10, $$1011_2$$ = 11, $$10110_2$$ = 22, $$101100_2$$ = 44. Therefore, we compute U1, U2, U4, U5, U10, U11, U22, and U44

It says "taking the bits one at a time" but the example clearly is not doing that. Though the end result makes sense because the index each term is either double or one more from the previous, and we have formals to easily compute that. However in example code I've seen it appears one bit at a time is inspected. So how exactly are the bits used to decide which terms to compute? Are the terms known and advanced or determined while checking each individual bit?

The bits are being checked one at a time from left to right, just not in quite the way you're expecting. In particular, if the bit is $$0$$, just use that bit. However, if it's $$1$$, use $$0$$ and then $$1$$. One way to look at this is that it's checking all of the possible bit values from $$0$$ up to the actual bit value, inclusive.

With the Wikipedia article's example of $$44 = 101100_2$$, the first bit is $$1$$. Since $$0$$ is not checked, it starts at $$1_2 = 1$$. The next bit is $$0$$, so it just checks $$10_2 = 2$$. The next bit after this is $$1$$, so it checks first with $$0$$, i.e., $$100_2 = 4$$, and then $$1$$, i.e., $$101_2 = 5$$. After this, the next bit is $$1$$ again, so it checks $$1010_2 = 10$$ and $$1011_2 = 11$$. The second last bit is $$0$$, so it just checks $$10110_2 = 22$$. Finally, the last bit is $$0$$, so it just checks $$101100_2 = 44$$.

Regarding the example code you've linked to, I took a look but found it somewhat hard to follow regarding exactly what it's doing. Although it does inspect one bit at a time, as you state, it seems to handle the $$1$$ bits by basically repeating certain calculations compared to the $$0$$ bits. First, note the following comment part of the LucasUVthTerm function:

    '''
BIT MANIPULATION OF P TO CREATE LUCAS CHAIN.
Lucas chain help us to know when to double and when to add one to 0btain next terms
example:let nth term be 24 and 112, below are its lucas chains
24 <- 12 <- 6 <- 3 <- 2 <- 1 <- 0
112 <- 61 <- 60 <- 30 <- 15 <- 14 <- 7 <- 6 <- 3 <- 2 <- 1 <- 0
We need to read this sequence from right to left.
here is how we do it without pre calculation of chain
1. we double the value add one when nth bit in the bitstring of p is 1,
so one carry two processes
2. we only double when nth bit in bitstring is zero
Note our terms start from zero index
'''


In the first example of the Lucas chain starting at $$24 = 11000_2$$, reading from right to left, you get $$0 = 0_2$$, $$1 = 1_2$$, $$2 = 10_2$$, $$3 = 11_2$$, $$6 = 110_2$$, $$12 = 1100_2$$ and $$24 = 11000_2$$. Apart from starting at $$0$$, this matches what Wikipedia states and I described above. However, their second example has a typo as the initial number should be $$122$$. With this correction, the rest of the values follow the expected pattern. Since $$122 = 1111010_2$$, the values are $$0 = 0_2$$, $$1 = 1_2$$, $$2 = 10_2$$, $$3 = 11_2$$, $$6 = 110_2$$, $$7 = 111_2$$, $$14 = 1110_2$$, $$15 = 1111_2$$, $$30 = 11110_2$$, $$60 = 111100_2$$, $$61 = 111101_2$$ and $$122 = 1111010_2$$.

Also, their "1." point's (for if the bit is $$1$$) second line says

so one carry two processes

Once again, this implies the process is repeated twice.

Later, in the code itself, there's

        # when nth bit is !, double and add one
if bit == "1":
#doubling from k to 2k
Uk = V * U
Vk = V**2 - 2 * Qk
k = 2 * k

#adding one from 2k to 2k + 1
Uk1 = (P * Uk + Vk) * (p + 1) // 2
Vk1 = (D * Uk + P * Vk)  * (p + 1) // 2
#cahnging alues of our Lucas parameters
U, V = Uk1 % p, Vk1 % p
Qk = (Qk**2) % p
Qk = (Qk * Q) % p
k = k + 1
else:
# doubling only
Uk = V * U
Vk = V**2 - 2 * Qk
U, V = Uk % p, Vk % p
Qk = (Qk**2) % p
k = 2 * k


As you can see, the $$1$$ bit conditional code section involves basically repeating certain calculations compared to the $$0$$ bit part.

• When dividing by 2 what should be done if the numerator is an odd number? In the example code I think it always rounds down. – northerner 3 hours ago
• On Wikipedia it states "If $P ⋅ U 2 k + V 2 k {\displaystyle P\cdot U_{2k}+V_{2k}} {\displaystyle P\cdot U_{2k}+V_{2k}}$ is odd, replace it with $P ⋅ U 2 k + V 2 k + n {\displaystyle P\cdot U_{2k}+V_{2k}+n} {\displaystyle P\cdot U_{2k}+V_{2k}+n}$; this is even so it can now be divided by 2. The numerator of $V 2 k + 1 {\displaystyle V_{2k+1}} {\displaystyle V_{2k+1}}$ is handled in the same way." but what about $V_{2k}$ having an odd numerator? – northerner 3 hours ago

Maybe the recursive solution to this problem is more straightforward. We have $$(U_k, V_k) = \begin{cases} (U_{k/2} \cdot V_{k/2}, \frac12(V_{k/2}^2 + D U_{k/2}^2)) & \text{ if k is even,} \\ (\frac12(P \cdot U_{k-1} + V_{k-1}), \frac12 (D \cdot U_{k-1} + P \cdot V_{k-1})) & \text{ if k is odd.} \end{cases}$$ This can be directly turned into code in the form of a recursive function.

Your ultimate goal is to computer $$(U_{44}, V_{44})$$; this recurrence will tell you that you need to know $$(U_{22}, V_{22})$$ first, and for that you need to know $$(U_{11}, V_{11})$$, and so forth.

This is, in disguise, exactly the logarithmic algorithm for computing the $$k^{\text{th}}$$ power of some kind of object $$\alpha$$: $$\alpha^k = \begin{cases} (\alpha^{k/2})^2 & \text{ if k is even,} \\ \alpha^{k-1} \cdot \alpha & \text{ if k is odd.} \end{cases}$$ We turn that into an algorithm for computing $$(U_k, V_k)$$ by using the relation $$\left(\frac{P + \sqrt D}{2}\right)^k = \frac{V_k + U_k \sqrt D}{2}.$$ Equivalently, this is a matrix power: $$\begin{bmatrix}U_k \\ V_k\end{bmatrix} = \begin{bmatrix}P/2 & 1/2 \\ D/2 & P/2\end{bmatrix}^k \begin{bmatrix}1 \\ P\end{bmatrix}.$$ Either way, we're taking the $$k^{\text{th}}$$ power of some fancy object, and what multiplying or squaring it does is exactly the thing in the original algorithm.