\\% Fast evaluation of (generalized) Lambert series and sum(n>=1,1/fibonacci(n))
\\ Author: Joerg Arndt
\\ License: GPL version 3 or later
\\ online at http://www.jjj.de/pari/
\\ version: 2012-January-28 (16:14)

read("aux-chks.gpi");  \\ for chkS_()

lambert(q, S=-1)=
{
/*
 Return L(q) = sum(n>=1, q^n/(1-q^n) (Lambert series)
 computed as sum(n=1, S, q^(n^2)*(1+q^n)/(1-q^n) )
 As series correct up to order S^2.
 We also have  L(q) = sum(n>=1, sigma(n,0)*q^n ), cf. A000005
 L(0.5) = 1.6066951524152917...
  = sum(n>=1, 1/(2^n-1)) (Erdos-Borwein constant)
 Vec(lambert(q/(1-q),S))  cf. A101509
 Vec(lambert(q/(1+q),S))  cf. A001659
*/
    S=chkS_(S);
    return( sum(n=1,S, q^(n^2)*(1+q^n)/(1-q^n) ) );
} /* ----- */

lambert_G(q, x, S=-1)=
{
/*
 Return G(q,x) = sum(n>=1, x*q^n/(1-x*q^n) (generalized Lambert series)
 computed as sum(n=1,S, q^(n^2)*x^n*( 1/(1-q^n) + x*q^n/(1-x*q^n) ) )
 As series in q correct up to order S^2.
 We also have  G(q,x) = sum(n>=1, x^n*q^n/(1-q^n) )
 G(q,x) = sum(n>=1, sumdiv(n,d,x^d)*q^n )
 G(q,x) = -x/(1-x) +  sum(n>=1, x^n/(1-q^n) )
 G(q,1) = L(q) = lambert(q)
 G(q^2,1/q) = sum(n>=1, sum_odd_divs(n)*q^n ), cf. A001227
 G(q^3,1/q) = divs 3k+2, cf. A001822
 G(q^3,1/q^2) = divs 3k+1, cf. A001817
 G(q,2) cf. A055895
 G(2*q,1/2) cf. A074854 a(n)=Sum(d|n)(2^(n-d)).
 G(q,3)/3 cf. A034730 Dirichlet convolution of b_n=1 with c_n=3^(n-1).
 G(q,1/(1-q)) cf. A097939 Sum of smallest parts of all compositions of n.
*/
    S=chkS_(S);
    return( sum(n=1,S, q^(n^2)*x^n*( 1/(1-q^n) + x*q^n/(1-x*q^n) ) ) );
} /* ----- */

lambert_G0(q, x, S=-1)=
{
/*
 Return G0(q,x) = sum(n>=0, x*q^n/(1-x*q^n) (generalized Lambert series)
 As series in q correct up to order S^2.
 We have G0(q,x) = -x*deriv(E,x)/E = -x*deriv(log(E),'x)  where  E = prod(n>=0, 1-x*q^n)
 and G0(q,x) = sum(n>=1, x^n/(1-q^n) )
*/
    S=chkS_(S);

 \\ alternative methods:
\\    return( -x * Etax_xderiv(q,x,S) / Etax(q,x,S) );
\\    return( sum(n=0,S, q^(n^2+n)*x^(n+1)*(1 - x*q^(2*n+1)) / ((1 - x*q^n)*(1 - q*q^n))) );

    return( x/(1-x) + lambert_G(q, x, S) );
} /* ----- */


lambert_qxt(q, x, t, S=-1)=
\\ Return sum(n>=0, x*t^n / (1-x*q^n) )   (3-parameter generalized Lambert series)
\\ As series in q correct up to order S^2.
\\ lambert_G0(q, x) is the special case lambert_qxt(q, x, q).
{
    S=chkS_(S);
\\    return( sum(n=0, N, x*t^n / (1-x*q^n) ) );
    return( sum(n=0, S, (1-x*t*q^(2*n))*x*(x*t)^n*q^(n^2) / ( (1-x*q^n)*(1-t*q^n)) ) );
} /* ----- */



inv_rec_sum(m1=1, m2=1, S=-1)=
{
/*
 Return sum(n>=1, 1/f(n) ) where f(0)=0, f(1)=1, f(n) = m1*f(n-1) + m2*f(n-1)
 computed using generalized Lambert series.
 Must have m1^2+4*m2 > 0
 Examples:
 inv_rec_sum(1, 1) == 3.35988566624318..., cf. A079586
 inv_rec_sum(3, 1) == 1.53537050883625..., cf. A153386
*/
    local(al,be);
    S=chkS_(S);
    \\ Note: the m2 here is -q in the following paper:
    \\  A. F. Horadam: "Elliptic functions and Lambert series in the summation
    \\   of reciprocals in certain recurrence-generated sequences",
    \\ The Fibonacci Quarterly, vol.26, no.2, pp.98-114, (May-1988)
    \\ The following numerical examples are for m1=m2=1:
    al=1/2*(m1+sqrt(m1^2+4*m2));  \\ == +1.6180339887498...
    be=1/2*(m1-sqrt(m1^2+4*m2));  \\ == -0.6180339887498...
\\    print([al,be, m1^2-4*m2]);
\\    print([be/al, 1/al]);  \\ == [-0.381966..., 0.618033...]
\\    print([-be^2/m2, -be/m2]);  \\ same

    return( (al-be)*( 1/(al-1) + lambert_G(be/al, 1/al, S) ) ); \\ == 3.35988566624318...
    \\ same as:
\\    return( (al-be)*( 1/(al-1) + lambert_G(-be^2/m2, -be/m2, S) ) );
    \\ another way to compute the same is
\\    return( (al-be)*( 1/(al-1) + lambert_G((be/al)^2, 1/al, S) + lambert_G( (be/al)^2, 1/be, S) ) );
} /* ----- */



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