Dead paper link in comments

[brubsby]

Dead link in a comment found here: https://github.com/aleaxit/gmpy/blob/205c4fdbd28fb64d2159822103bba9b4d7ed09ec/src/gmpy_mpz_lucas.c#L40C41-L40C84

link: http://joye.site88.net/papers/JQ96lucas.pdf

I assume it is this '96 paper: https://citeseerx.ist.psu.edu/document?doi=95e1b5b6cc356ecc6c860be9f8cca55502e37fce

It also seems that an improvement was made to this algorithm in 2010, if anyone were willing to port it to gmpy: https://www.researchgate.net/publication/220099375_On_Lucas_Sequences_Computation

I was going to port whatever algorithm that I had been using in gmpy2 to pari/gp, so that I could more easily create SNFS sieving polynomials for generalized lucas numbers, (which makes use of pari/gp's lindep function) but it turns out the original paper included a pari/gp implementation!

However, the 30 year old link in the paper is, surprisingly, dead. But I was able to find what I believe to be a copy in the homework directory of a number theory file server from my alma mater:

\\lucas2.gp - Fast computation of Lucas sequences.
\\You can get a copy of this file  via anonymous ftp at 
\\math.math.ucl.ac.be in the directory /pub/joye/pari.

\\lucas2(k,P,Q) returns the k-th terms of the Lucas sequences
\\{U_k(P,Q)} and {V_k(P,Q)} with parameters P and Q. 

\\Usage:
        print("lucas2(k,P,Q)=k-th term of U_k(P,Q) and V_k(P,Q)");

{
        lucas2(k,P,Q,     n,s,j,Ql,Qh,Uh,Vl,Vh) =
        n = floor(log(k)/log(2)) + 1;
        s = 0;
        while( bittest(k,s) == 0, s = s + 1 );
        Uh = 1; Vl = 2; Vh = P;
        Ql = 1; Qh = 1;
        forstep( j = n-1, s+1, -1,
                Ql = Ql * Qh;
                if( bittest(k,j),
                        Qh = Ql * Q;
                        Uh = Uh * Vh;
                        Vl = Vh * Vl - P * Ql;
                        Vh = Vh * Vh - 2 * Qh,

                        Uh = Uh * Vl - Ql;
                        Vh = Vh * Vl - P * Ql;
                        Vl = Vl * Vl - 2 * Ql;
                        Qh = Ql;
                );
        );
        Ql = Ql * Qh; Qh = Ql * Q;
        Uh = Uh * Vl - Ql;
        Vl = Vh * Vl - P * Ql;
        Ql = Ql * Qh;
        for( j = 1, s,
                Uh = Uh * Vl;
                Vl = Vl * Vl - 2 * Ql;
                Ql = Ql * Ql;
        );
        [Uh,Vl]
}

Small world!