\\% Manipulations of continued fractions.
\\ Author: Joerg Arndt
\\ License: GPL version 3 or later
\\ online at http://www.jjj.de/pari/
\\ version: 2011-January-19 (12:51)


cfab2pq(a,b,m=-2)=
{
/* Example:
    a=vector(5,j,1+(j>1)) ==> [1, 2, 2, 2, 2]
    b=vector(5,j,1)       ==> [1, 1, 1, 1, 1]
    cfab2pq(a,b) ==> [41, 29]  (41/29 approx sqrt(2))
 */
    local(n, p, p1, p2, q, q1, q2, i);
    n = min(length(a), length(b))-1;
    if ( m==-2,  m = n ); \\ default: m=n
    if ( m<0,  return( [1, 0] ) ); \\ infinity
    p1 = 1;
    q1 = 0;
    p = a[1];
    q = 1;  \\ b[1] unused
    for (k=1, m,
        i = k+1;
        p2 = p1;  p1 = p;
        q2 = q1;  q1 = q;
        p = a[i]*p1 + b[i]*p2;
        q = a[i]*q1 + b[i]*q2;
    );
    return( [p,q] );
} /* ------ */

cfab2pqvec(a,b=0)=
{
/* Example:
    a=vector(5,j,1+(j>1)) ==> [1, 2, 2, 2, 2]
    b=vector(5,j,1)       ==> [1, 1, 1, 1, 1]
    cfab2pqvec(a,b) ==> [[1, 3, 7, 17, 41], [1, 2, 5, 12, 29]]
 */
    local(n, p, p1, p2, q, q1, q2, i, pv, qv);
    n = length(a)-1;
    if ( n<=0,  return( [[], []] ) ); \\ infinity
    if ( b==0,  b=vector(n+1,j,1); );
    pv = vector(n+1);
    qv = vector(n+1);
    p1 = 1;
    q1 = 0;
    p = a[1];
    q = 1;  \\ b[1] unused
    pv[1] = p;  qv[1] = q;
    for (k=1, n,
        i = k+1;
        p2 = p1;  p1 = p;
        q2 = q1;  q1 = q;
        p = a[i]*p1 + b[i]*p2;
        q = a[i]*q1 + b[i]*q2;
        pv[k+1] = p;  qv[k+1] = q;
    );
    return( [pv,qv] );
} /* ------ */


cfab2r(a,b, m=-2)=
{
/* Example:
    a=vector(5,j,1+(j>1)) ==> [1, 2, 2, 2, 2]
    b=vector(5,j,1)       ==> [1, 1, 1, 1, 1]
    cfab2r ==> 41/29
 */
    local(n, r);
    n = min(length(a), length(b));
    if ( m==-2, m = n-1 ); \\ default: m=n
    if ( m>=n,  m=n-1 );
    if ( m<0, return( 0 ) ); \\ infinity
    r = 0;
    m += 1;
    forstep (k=m, 2, -1,  r = b[k]/(a[k]+r); );
    r += a[1];  \\ b[1] unused
    return( r );
} /* ------ */


cfab2rvec(a,b=0)=
{
    local(t,r);
    t = cfab2pqvec(a,b);
    r = vector( length(t[1]), j, t[1][j]/t[2][j] );
    return( r );
} /* ------ */

cf2simple(A,B)=
{
    local(c);  c=1;
    for (j=2, #A-1,
        c = 1/(B[j]);
\\        print([j-1,c]);
        B[j] *= c;  \\ B[j]==1
        B[j+1] *= c;
        A[j] *= c;
    );
    \\ note last term of B[] != 1 in general
    return([A,B]);
} /* ----- */

cf2simpleB(A,B)=
{
    local(c);  c=1;
    for (j=2, #B-1,  \\ leave a0 as it is
        c = 1/(A[j]);
\\        print([j-1,c]);
        A[j] *= c;  \\ A[j]==1
        B[j] *= c;
        B[j+1] *= c;
    );
    \\ note first and last term of A[] != 1 in general
    return([A,B]);
} /* ----- */


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