Exercise 4.4.28* (Algorithm to compute $V_n$)

Question

Show that $V_{2n+1} = V_n V_{n+1} - a (-b)^n$. Explain how this formula and the duplication formula (4.13) can be used to compute the triple $(V_{2n}, V_{2n+1}, (-b)^{2n})$, if the triple $(V_n, V_{n+1}, (-b)^n)$ is known. Similarly, explain how the triple $(V_{2n+1}, V_{2n+2}, (-b)^{2n+1})$ can be computed in terms of the triple $(V_n, V_{n+1}, (-b)^n)$. Explain how this triple can be determined for general $n$ by using these two operations. (This method is not very much more efficient than the method described in the text, but it involves less work in the special case $b = -1$. By constructing congruential analogues of the identities in Problem 18 one may see that for purpose of constructing Lucas pseudoprime tests this does not involve any loss of generality.)

Richard Ganaye · Accepted Answer

ewcommand{\Z}{\mathbb{Z}}

ewcommand{\N}{\mathbb{N}}

\begin{proof} 
\begin{enumerate}
\item[(a)] By definition
\begin{align}
U_n = \frac{\lambda^n - \mu^n}{\lambda - \mu},\qquad V_n = \lambda^n + \mu^n,
\end{align}
where $\lambda, \mu$ are the roots of $x^2 - a x -b$, so that 
\begin{align}
\lambda + \mu = a,\qquad \lambda \mu = -b.
\end{align}
Then
\begin{align*}
V_{2n+1} & = \lambda^{2n+1}  + \mu^{2n+1}\
&= ( \lambda^{n}  + \mu^{n})( \lambda^{n+1}  + \mu^{n+1}) -(\lambda \mu)^n (\lambda + \mu)\
&= V_n V_{n+1} - a(-b)^n&(	ext{by (1)}).
\end{align*}
In conclusion,
\begin{align}
V_{2n+1} & = V_n V_{n+1} - a(-b)^n.
\end{align}
(Lucas defined $\lambda, \mu$ as roots of $x^2 - Px + Q$, which gives the more aesthetically pleasing formula $V_{2n+1} = V_n V_{n+1} - P Q^n$.)
\item[(b)] We recall the duplication formulas (4.13):
$$U_{2n} = U_n V_n,\qquad V_{2n} = V_n^2 - 2(-b)^n.$$
Suppose that we know the triple $(V_{n}, V_{n+1}, (-b)^n)$. Then
\begin{align}
V_{2n} \quad &= V_n^2 - 2(-b)^n,\
V_{2n+1} &= V_n V_{n+1} -a (-b)^n,\
(-b)^{2n} &=  [(-b)^n]^2.
\end{align}
This gives the triple $(V_{2n}, V_{2n+1},(-b)^{2n} )$ with a fixed number of multiplications and additions.

Similarly
\begin{align}
V_{2n+1} &= V_n V_{n+1} -a (-b)^n,\
V_{2n+2} &= V_{n+1}^2 + 2 b(-b)^n,\
(-b)^{2n+1} &= -b [(-b)^{n}]^2.
\end{align}

\item[(c)] We compute the triple $(V_k(a, b),  V_{k+1}(a, b), (-b)^n)$ recursively,. We use (4),(5),(6) with $n = k/2$ if $n$ is even, and (7),(8),(9)  with $n = (k-1)/ 2$ is $n$ is odd (see algorithm in the note 1), until we reach $n = 0$, for which $(V_0(a, b),  V_{1}(a, b), (-b)^0) = (2, a , 1)$. If we compute $V_k \pmod p)$ by reducing modulo $p$ at each step, the time is $T = O(\log(k))$.

\end{enumerate}
\end{proof}
Note 1: This algorithm gives in Python the following procedure

\begin{verbatim}
def triple(k, a, b):
    """return the tuple (V_k(a,b), V_(k+1)(a,b), (-b)^n)"""
    if k == 0: return (2, a, 1) 
    if k % 2 == 0:
        n = k // 2
        (V, W, P) = triple(n, a, b)
        return (V * V - 2 * P, V * W - a * P, P * P)
    else:
        n = (k - 1) // 2
        (V, W, P) = triple(n, a, b)
        return (V * W - a * P, W * W + 2 * b * P, -b * P * P)
        
def V(k, a, b):
    """return (V_k(a,b)"""
    v , _ , _ = triple(k, a, b)
    return v
\end{verbatim}
To test this program, we verify that $V_p(1, 1) \equiv a = 1 \pmod p$ for some prime $p$.
\begin{verbatim}
p = 11213
V(p, 1, 1)

239006726.....34499771   (2344 digits)

V(p, 1, 1) % p
1
\end{verbatim}
We give the modular version:
\begin{verbatim}
def triple_mod(k, a, b , p):
    """return the tuple (V_k(a,b), V_(k+1)(a,b), (-b)^n) mod p"""
    a, b = a % p, b % p
    if k == 0: return (2, a, 1) 
    if k % 2 == 0:
        n = k // 2
        (V, W, P) = triple_mod(n, a, b, p)
        return ((V * V - 2 * P) % p, (V * W - a * P) % p, (P * P) % p)
    else:
        n = (k - 1) // 2
        (V, W, P) = triple_mod(n, a, b, p)
        return ((V * W - a * P) % p, (W * W + 2 * b * P) % p, (-b * P * P) % p)
        
def V_mod(k, a, b, p):
    """return (V_k(a,b) mod p"""
    v , _ , _ = triple_mod(k, a, b, p)
    return v
    
p = 1234567891234567891234567891
# p is prime!
V_mod(p, 1, 1, p)
1
\end{verbatim}

\bigskip
Note 2: A more natural and generalizable method to compute $V_n$ (or $V_n \mod p$) is to write the relation between the $V_n$ in matrix form. Since $V_0 = 2,\ V_1 = a$ and $V_{n+2} = a V_{n+1} + b V_n$, then
\begin{align*}
\begin{pmatrix} V_{n+1} \ V_{n+2} \end{pmatrix} &= \begin{pmatrix} 0 & 1 \ b & a  \end{pmatrix} \begin{pmatrix} V_{n} \ V_{n+1} \end{pmatrix},
\end{align*} 
so
\begin{align*}
\begin{pmatrix} V_{n} \ V_{n+1} \end{pmatrix} = \begin{pmatrix} 0 & 1 \ b & a  \end{pmatrix}^n \begin{pmatrix} 2 \ a \end{pmatrix}.
\end{align*} 
By using fast exponentiation for matrices, we obtain the result $V_n \mod p$ in a time $T = O( \log(n))$.

With Sagemath:
\begin{verbatim}
def V(n, a, b, p):
    """return (V_n(a,b) mod p"""
    a , b = Mod(a, p), Mod(b, p)
    A = matrix([[0, 1], [b, a]])
    B = A^n * matrix([[2], [a]])
    return B[0][0]
    
p = 1234567891234567891234567891
V(p, 1, 1, p)

1
%timeit V(p, 1, 1, p)
1000 loops, best of 3: 981 micro-s per loop

%timeit V_mod(p, 1, 1, p)
The slowest run took 60.92 times longer than the fastest. 
This could mean that an intermediate result is being cached.
1000 loops, best of 3: 244 micro-s per loop

\end{verbatim}
We see in this example that the execution speeds are comparable, nevertheless with an advantage for the first function V\_mod.

Note 3: For Problem 27, the second method becomes indispensable. This gives
\begin{verbatim}
sage: def u(n, p):
....:     """return u_n mod p, see Problem 4.4.27* from Niven"""
....:     R = Integers(p)
....:     A = matrix(R, [[0, 1,  0], [0, 0 , 1], [1, 1, 0]])
....:     B = A^n * matrix(R, [[3], [0], [2]])
....:     return B[0][0]
....: 
\end{verbatim}
To verify the affirmation in the statement of Problem 27:
\begin{verbatim}
sage:  for n in range(2, 1000000):
....:      if u(n,n) == 0 and not is_pseudoprime(n):
....:          print(n)
....:         
    271441
    904631
\end{verbatim}
So the least composite number with the property $p \mid u_p$ is $271,441 = 521^2$.

Exercise 4.4.28* (Algorithm to compute $V_n$)

Answers

Comments

Exercise 4.4.28* (Algorithm to compute $V_n$)

Answers

Comments

Add answer