Computing Lucas sequence for large n

Question

I've been trying to write a test function for Fibonacci pseudo-primes with large $n$. Fibonacci pseudoprimes are composite numbers such that $V_n(P,Q) \equiv P \mod n$ for $P>0$ and $Q =\pm 1$, with $V_n$ Lucas sequence.

As such I need to compute Lucas sequence for large $n$. They are defined by:

$U_n(P,Q)=\frac{a^n-b^n}{a-b}$ et $V_n(P,Q)=a^n+b^n$ with $a,b=(P\pm \sqrt{P^2-4Q})/2$

Direct computation is not possible because I need to work with integers and $\sqrt{P^2-4Q}$ may not be one.

There are also binomial formulations of the form $U_n(P,Q)=2^{1-n}\sum_{k=0}^{\lfloor (n-1)/2 \rfloor}\binom{n}{2k+1}P^{n-2k-1}(P^2-4Q)^k$. But computing $U_n$ that way seems to be $O(n^2)$ which is not acceptable for me.

Are there any ways to compute $U_n$ and $V_n$ in $O(1)$ or $O(\log n)$ which involve only integers and not floating point ?

Thanks in advance

Answer 1

Following Emil remark, I used other relations to obtain a $O(\log n)$ algorithm:

$\begin{cases}U_{2n}=U_nV_n\\ V_{2n}=V_n^2-2Q^n\\ U_{2n+1} = U_{n+1}V_n - Q^n\\ V_{2n+1} = V_{n+1}V_n -PQ^n\end{cases}$

This concluded in the following python code:

def lucas_sequence(p, q, n):
    if n == 0:
        return (0, 2)
    if n == 1:
        return (1, p)
    if n == 2:
        return (p, p**2 - 2* q)
    ns2 = n // 2
    un, vn = lucas_sequence(p, q, ns2)
    if n % 2 == 0:
        return (un*vn, vn**2 - 2 * (q**ns2))
    else:
        unp1, vnp1 = lucas_sequence(p, q, ns2 + 1)
        return (unp1*vn - q**ns2, vnp1*vn - p * (q**ns2))