\\% Eta-products and Ramanujan functions.
\\ Author: Joerg Arndt
\\ License: GPL version 3 or later
\\ online at http://www.jjj.de/pari/
\\ version: 2011-August-15 (14:22)

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

\\read("lambert.gpi");  \\ for (commented out method for) Etax_xderiv()


Eta(q, S=-1)=
\\ Return (eta-product:) prod( n>=1, 1 - q^n ).
\\ Computation via the pentagonal number theorem.
\\ As series correct up to order floor(3/2 * S^2).
{
    local( et=1, qn=q, q3n2=q, q3=q^3, tt=1, ttq );
    S=chkS_(S);
\\    return( prod( n=1, S*S*3\2, 1-q^n ) );
\\    return( 1+sum(n=1, S, (-1)^n*(q^(n*(3*n-1)/2)+q^(n*(3*n+1)/2))) );
\\    return( 1+sum(n=1, S, (-1)^n*q^(n*(3*n-1)/2)*(1+q^n)) );
    for (n=1, S,
        tt *= q3n2;
        ttq = tt * (1+qn);
        if ( bitand(n, 1),  et-=ttq, et+=ttq );
        qn *= q;
        q3n2 *= q3;
    );
    return( et );
} /* ----- */

Eta_pow3(q, S=-1)=
\\ Return Eta(q)^3 = sum(n>=0, (-1)^n*(2*n+1)*q^(n*(n+1)\2) ).
{
    S=chkS_(S);
    return( sum(n=0, S, (-1)^n*(2*n+1)*q^(n*(n+1)\2) ) );
} /* ----- */


Eta_qderiv(q, S=-1)=
\\ Return derivative (with respect to q) of eta-product.
\\ As series correct up to order floor(3/2 * S^2).
{
    S=chkS_(S);
    return( sum(n=1, S,
        (-1)^n*(  (n*(3*n-1)/2) * q^(n*(3*n-1)/2-1) +
                  (n*(3*n+1)/2) * q^(n*(3*n+1)/2-1) ) )
    );
} /* ----- */

Eta_plus(q, S=-1)=
\\ Return prod( n>=1, 1 + q^n )
\\ As series correct up to order floor(3/2 * S^2).
{
    S=chkS_(S);
    return( Eta(q^2, S) / Eta(q, S) );
} /* ----- */

Eta_plus_qderiv(q, S=-1)=
\\ Return derivative (with respect to q) of eta-plus-product.
\\ As series correct up to order floor(3/2 * S^2).
{
    local(e1, e2, de1, de2);
    S=chkS_(S);
    e1 = Eta(q, S);
    e2 = Eta(q^2, S);
    de1 = Eta_qderiv(q, S);
    de2 = sum(n=1, S\2+1,
        (-1)^n*( (n*(3*n-1)) * q^(n*(3*n-1)-1) +
                 (n*(3*n+1)) * q^(n*(3*n+1)-1) )
    );
    return(  (de2*e1-e2*de1)/(e1^2) );
} /* ----- */


Etax_plus(q, x, S=-1)=
\\ Return q-exponential series Exp_q(x) = prod( n>=0, 1 + x*q^n )
\\ computed as sum(n=0, S, x^n*q^(n*(n-1)/2) / prod(k=1, n, 1-q^k) )
\\ As series in q correct up to order floor(S^2/2).
{
    local(xn=1, qn=1, qn2=1, pr=1, ret=1);
    S=chkS_(S);
\\    return( prod(n=0, S*S\2, 1+x*q^n) );
\\    return( sum(n=0, S, x^n*q^(n*(n-1)/2) / prod(k=1, n, 1-q^k) ) );
    for (n=1, S,
\\        ret += x^n*q^(n*(n-1)/2) / prod(k=1, n, 1-q^k);
        xn *= x;  \\ == x^n
        qn2 *= qn;  qn *= q;  \\ qn2 == q^(n*(n-1)/2)
        pr /= (1-qn);  \\ pr == 1/prod(k=1, n, 1-q^k);
        ret += ( xn*qn2 * pr );
    );
    return( ret );
} /* ----- */
\\
\\ Etax_plus(q, q) = Eta_plus(q).
\\ Etax_plus(q, x) = Etax(q, -x).
\\ 1 / Etax_plus(q, x) = sum(n>=0, (-x)^n/prod(k=1, n, 1-q^k ) ).
\\ Etax_plus(q, x) = Exp_q(q, x).
\\ Etax_plus(q, (1-q)*x) = Exp_Q(q, x).


Etax(q, x, S=-1)=
\\ Return inverse q-exponential series 1/exp_q(x) = prod( n>=0, 1 - x*q^n )
\\ As series in q correct up to order floor(S^2/2).
{
    S=chkS_(S);
\\    return( prod(n=0, S*S\2, 1-x*q^n) );
    return( Etax_plus(q, -x, S) );
} /* ----- */
\\
\\ Write E(q, x) for Etax(q, x) and P(q, x) for Etax_plus(q, x), then
\\
\\ E(q, +x) = P(q, -x).
\\ E(q, q) = Eta(q).
\\
\\ E(q,x) = sum(n=0, S, (-x)^n * q^(n*(n-1)/2) / prod(k=1, n, 1 - q^k) )
\\
\\ 1 / E(q, x)  = sum(n>=0, x^n / prod(k=1, n, 1-q^k ) ) =
\\   = exp_q(q, x) =
\\   = exp( sum(n>=1, ( 1/n *  x^n ) / (1 - q^n) ) ).
\\
\\ 1 / E(q, (1-q)*x)  = 1 / prod(n>=0, 1 - x*(1-q)*q^n )
\\   = exp_Q(q, x) =
\\   = sum(n>=0, x^n/prod(k=1, n, (1 - q^k)/(1 - q) ) ) =
\\   = exp( sum(n>=1, ( 1/n * x^n * (1 - q)^n ) / (1 - q^n) ) ).
\\
\\ E(+q, +x) * E(+q, -x) = P(+q^2, -x^2) = E(+q^2, +x^2)
\\ E(+q, +I*x) * E(+q, -I*x) = P(+q^2, +x^2)
\\
\\ E(+q, +x^2) = E(+q, +x) * E(+q, -x) * E(+q^2, +q*x^2)
\\ Now replace x by I*x to obtain
\\ P(+q, +x^2) = E(+q, +I*x) * E(+q, -I*x) * E(+q^2, -q*x^2) =
\\   = P(+q^2, +x^2) * E(+q^2, -q*x^2) =
\\   = P(+q^2, +x^2) * P(+q^2, +q*x^2) =
\\   = E(+q^2, -x^2) * E(+q^2, -q*x^2)
\\
\\ Multiplicative splitting:
\\ E(q, x) = prod( s=0, M-1, E(q^M, q^s * x) )  for M>=1
\\ For example,
\\ E(q, x) = E(q^2, x) * E(q^2, q * x)
\\ E(q, x) = E(q^3, x) * E(q^3, q * x) * E(q^3, q^2 * x)
\\ The identities also hold for P(q, x).



Etax_xderiv(q, x, S=-1)=
\\ Derivative of Etax(q, x) with respect to x.
\\ Computation as sum(n>=1, -n * (-x)^(n-1) * q^(n*(n-1)/2) / prod(k=1, n, 1 - q^k) ).
\\ Etax_xderiv(q, x) = Etax(q, x)^2 * sum(n>=1, -n * x^(n-1) / prod(k=1, n, 1 - q^k) )
{
     \\ Alternative method using generalized Lambert series:
\\    local(L, E);
\\    L=lambert_G0(q, x, S);
\\    E=Etax(q, x, S);
\\    return( (-1/x)*E*L );

    local(xn=-1, qn=1, qn2=1, pr=1, ret=0);
    S=chkS_(S);
\\    return( sum(n=1, S, -n * (-x)^(n-1) * q^(n*(n-1)/2) / prod(k=1, n, 1-q^k) ) );
    for (n=1, S,
\\        ret += ( -n * (-x)^(n-1) * q^(n*(n-1)/2) ) / prod(k=1, n, 1-q^k);
        qn2 *= qn;  qn *= q;  \\ qn2 == q^(n*(n-1)/2)
        pr /= (1-qn);  \\ pr == 1/prod(k=1, n, 1-q^k);
        ret += ( n * xn*qn2 * pr );
        xn *= -x;  \\ == x^(n-1)
    );
    return( ret );
} /* ----- */

Etax_plus_xderiv(q, x, S=-1)=
\\ Derivative of Etax_plus(q, x) with respect to x.
\\ Etax_plus_xderiv(q, x) = -Etax_xderiv(q, -x, S) ).
{
\\    return( -Etax_xderiv(q, -x, S) );

    local(xn=1, qn=1, qn2=1, pr=1, ret=0);
    S=chkS_(S);
\\    return( sum(n=1, S, n * x^(n-1) * q^(n*(n-1)/2) / prod(k=1, n, 1-q^k) ) );
    for (n=1, S,
\\        ret += n * x^(n-1) * q^(n*(n-1)/2) / prod(k=1, n, 1-q^k);
        qn2 *= qn;  qn *= q;  \\ qn2 == q^(n*(n-1)/2)
        pr /= (1-qn);  \\ pr == 1/prod(k=1, n, 1-q^k);
        ret += ( n * xn*qn2 * pr );
        xn *= x;  \\ == x^(n-1)
    );
    return( ret );
} /* ----- */


Eta_msect(q, s, m, S=-1)=
\\ Multisection of eta-product:
\\ return series prod( n>=1,  1-q^n )
\\ but only terms q^k such that k == s (mod m)
{
    local(et_, qq);
    S=chkS_(S);
    \\ eta(q) == 1 + sum(n=1, S, (-1)^n * q^(n*(3*n-1)/2) * (1+q^n) )
    et_ = if(s==0, 1, 0);
    for (n=1, S,
        qq = (n*(3*n-1)/2);
        if ( (qq%m)==s, et_ += (-1)^n*q^qq);
        qq = (n*(3*n+1)/2);
        if ( (qq%m)==s, et_ += (-1)^n*q^qq);
    );
    return( et_ );
} /* ----- */

\\ ellj(q^2):
\\ellj_q2(q, S)=
\\{
\\    local(kp2, kp4);
\\    kp2 = ellkprime(q, S)^2;
\\    kp4 = kp2^2;
\\    return( 256 * (1-kp2+kp4)^3 / (kp4*(1-kp2)^2) );
\\} /* ----- */

ellj_q(q, S)=
\\ Same as ellj(q), but this one does not switch to ellj(tau) for complex arguments
{

\\    local(e1, e2);
\\    e1=Eta(q, S);  e2=Eta(q^2, S);
\\    ((e1^24 + 256 * q *e2^24) / (e1 * e2)^8)^3 / (q * e1^24)
\\    (256*e2^24*q + e1^24)^3/(e2^24*e1^48*q)
    local(P1, P2);
    P1=Eta(q, S)^24;  P2=Eta(q^2, S)^24;
    return( (P1 + 256*q*P2)^3 / (q*P1^2*P2) );
} /* ----- */

\\ Ramanujan functions:
ram_ff(a, b, S=-1)={S=chkS_(S); sum(n=-S, S, a^(n*(n+1)/2) * b^(n*(n-1)/2));}
\\ == (-a, a*b) * (-b, a*b) * (a*b, a*b)  where (a, b) = Etax(b, a) = Etax_plus(b, -a)
\\ == ram_ff(b, a)
\\ ram_ff(q*w, q/w) == sum(n>=0, q^(n^2)*( w^n + 1/w^n )
\\ Theta3(q) = ram_ff(+q, +q)
\\ Theta4(q) = ram_ff(+q, -q)
\\ Theta4(q^4) = ram_ff(-q, -q)
ram_f(q, S=-1)=Eta(-q, S);  \\ ram_f(-q) == ram_ff(-q, -q^2)
ram_phi(q, S=-1)=Theta3(q, S);
ram_chi(q, S=-1)=Eta(q^2, S)^2/(Eta(q, S)*Eta(q^4, S)); \\ == (1+q)*(1+q^3)*(1+q^5)*(1+q^7)*...
ram_psi(q, S=-1)={S=chkS_(S); sum(n=0, S, q^(n*(n+1)/2));} \\ == Theta2_q4(sqrt(q))/2


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