4

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?

northerner
  • 449
  • Please do not delete after having received an answer. – quid Aug 29 '19 at 09:57
  • @quid I did that because there was a valuable answer but it got deleted for some reason. As it is now my preference is this gets deleted as I don't think it would be useful to anyone. – northerner Aug 29 '19 at 10:30
  • I don't oversee all the details but it seems you had pointed out some error in that answer. Maybe that's why but the situation is not exactly clear to me. – quid Aug 29 '19 at 10:32

2 Answers2

0

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.

John Omielan
  • 47,976
  • 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 Aug 19 '19 at 08:54
  • 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 Aug 19 '19 at 08:54
  • @northerner I'm quite certain you're correct. As you can see for the Uk1 and Vk1 calculations, it uses "// 2" at the end. I believe the code is in Python, where the double slash forces integer division (i.e., it always rounds down, basically implementing the floor function). For example, the answer to Stack Overflow's Python division describes this in some detail. – John Omielan Aug 20 '19 at 04:46
  • @northerner With $V_{2k}$, note it's also equal to $V_k^2 - 2Q^k$, which is an integer as both $V_k$ and $Q^k$ should always be integers. Thus, the next part of the equation, i.e., $\frac{V_k^2 + DU_k^2}{2}$ must also be an integer. As there's no mention of any kind of truncation or other adjustment, I'm quite certain it's because the numerator part must be even and, thus, can't be odd. – John Omielan Aug 20 '19 at 05:01
  • I think you're mistaking. I tested it out and some input $n$ cause the last term ($v_{n+1}$) to have an odd numerator. This seems to happen often when n is prime and $v_{n+1-1}=0$ because $\frac{v_k^2+Du^2_k}{2}$ becomes $\frac{Du^2_{n+1}}{2}$ where D is always odd. – northerner Aug 26 '19 at 08:49
  • Take for example n=7 where D=5. The numerator of $v_8$ is odd. – northerner Aug 26 '19 at 08:50
  • I (for now) undeleted this answer, since it's visibility might help to resolved some confusion related to the thread. (See comments on main.) – quid Aug 29 '19 at 10:35
  • Thank you for the detailed reply this clarifies the situation quite a bit. If you don't mind, could we leave this up for a bit. – quid Aug 29 '19 at 16:43
0

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.

Misha Lavrov
  • 142,276