\\% Solving  X^2 - d*Y^2 == +-1
\\ Author: Joerg Arndt
\\ License: GPL version 3 or later
\\ online at http://www.jjj.de/pari/
\\ version: 2011-January-19 (12:49)

read("contfrac.gpi");

cf_solvepell(d)=
{ /* solve x^2-d*y^2 = +-1 using the current realprecision */
    local(a, qp, pv, qv );
    local(n, i, p, q, e);
    a = contfrac(sqrt(d));
    pq = cfab2pqvec(a);
    pv = pq[1];
    qv = pq[2];
    n = length(a);
    i = 0;  \\ P_i, Q_i is a solution
    for (k=1, n,
        p = pv[k];
        q = qv[k];
        e = p^2-d*q^2;
        if ( (+1==e) || (-1==e),  i=k; break() );
    );
    return( [i, e, p, q] );  \\ p^2-d*q^2==e, p=P_i, q=Q_i
} /* ----- */

solvepell(d)=
{ /* solve x^2-d*y^2 = +-1 */
    local(rp, wp, sp);
    rp = default(realprecision,,1); \\ save precision
    wp = rp;
    sp = [0];
    while ( 0==sp[1],
        sp = cf_solvepell(d);
\\        wp += 10; \\ linear
        wp = ceil(1.5*wp + 7);  \\ exponential
        default(realprecision, wp);
    );
    default(realprecision,rp); \\ restore precision
    return( sp );
} /* ----- */

chk_pell(d,v)=
{
/* Example
    d=53; s=solvepell(d) ==> [5, -1, 182, 25]
    chk_pell(d,s) ==> -1
    chk_pell(d,[7,7,7,7]) ==> 0
*/
    local(e, p, q);
    e=v[2]; p=v[3]; q=v[4];
    if ( e==p^2-d*q^2,  return(e), return(0) )
} /* ----- */

\\ ==== end of file ====