\\% Conversion from power series to infinite products
\\ Author: Joerg Arndt
\\ License: GPL version 3 or later
\\ online at http://www.jjj.de/pari/
\\ version: 2011-August-07 (08:22)

read("ser2lambert.gpi");


ser2prod(t)=
{
/* Let t=[1,a1,a2,a3, ...], n=length(v), where t(x)=1+sum_{k=1}^{n}{a_k*x^k};
 * Return  p=[p1,p2,p3,...] so that (up to order n)
 * t(x)=\prod_{j=1}^{n}{(1-x^j)^{p_j}}
 */

    local(v);  \\ version by Michael Somos
    v = Ser(t);
    v = v'/v;
    v = vector(#t-1, j, polcoeff(v, j-1));
    v = ser2lambert(v);
    v = vector(#v, j, -v[j]/j);
    return( v );

\\    local(a1, a, q, v);
\\    a1 = t[2];
\\    if ( a1==0,  \\ must have nonzero linear term
\\        for (k=2, #t, t[k]+=t[k-1]; );  \\ divide by (1-x)
\\    );
\\    q = Vec( 'x*deriv(Ser(t),'x)/Ser(t) );
\\    v = ser2lambert(q);
\\    v = vector(#v, j, -v[j]/j);
\\    if ( a1==0,  v[1]+=1; );  \\ correct with zero linear term
\\    return( v );
} /* ----- */

prod2ser(p, w='y)=
{
    local(n, v);
    n = length(p);
    v = prod(k=1,#p, (1-(w+O(w^n))^k)^(p[k]));
    return( v );
} /* ----- */

ser2etaprod(v)=
{
/* Let t=[1,a1,a2,a3, ...], n=length(v), where t(x)=1+sum_{k=1}^{n}{a_k*x^k};
 * Return  p=[p1,p2,p3,...] so that (up to order n)
 * t(x)=\prod_{j=1}^{n}{eta(x^j)^{p_j}}
 * where eta(x) = prod(k>0, (1-x^k))
 */
    local(n, t);
    v = ser2prod(v);
    n = length(v);
    for (k=1, n,
        t = v[k];
        forstep (j=k+k, n, k, v[j]-=t; );
    );
    return( v );
} /* ----- */

etaprod2ser(p, w='y)=
{
    local(n, v);
    n = length(p);
    v = prod(k=1,#p, eta( (w+O(w^n))^k )^p[k] );
    return( v );
} /* ----- */

ser2prodplus(t)=
{
/* Let t=[1,a1,a2,a3, ...], n=length(v), where t(x)=1+sum_{k=1}^{n}{a_k*x^k};
 * Return  p=[p1,p2,p3,...] so that (up to order n)
 * t(x)=\prod_{j=1}^{n}{(1+x^j)^{p_j}}
 */
    local(v);  \\ version by Michael Somos
    v = Ser(t);
    v = v'/v;
    v = vector(#t-1, j, polcoeff(v, j-1));
    v = ser2lambertplus(v);
    v = vector(#v, j, v[j]/j);
    return( v );

\\    local(a1, q, v);
\\    a1 = t[2];
\\    if ( a1==0,   \\ must have nonzero linear term
\\        for (k=2, #t, t[k]-=t[k-1]; );  \\ divide by (1+x)
\\    );
\\    q = Vec( 'x*deriv(Ser(t),'x)/Ser(t) );
\\    v = ser2lambertplus(q);
\\    v = vector(#v, j, v[j]/j);
\\    if ( a1==0,  v[1]+=1; );  \\ correct with zero linear term
\\    return( v );
} /* ----- */

prodplus2ser(p, w='y)=
{
    local(n, v);
    n = length(p);
    v = prod(k=1,#p, (1+(w+O(w^n))^k)^(p[k]));
    return( v );
} /* ----- */


ser2etaprodplus(v)=
{
/* Let t=[1,a1,a2,a3, ...], n=length(v), where t(x)=1+sum_{k=1}^{n}{a_k*x^k};
 * Return  p=[p1,p2,p3,...] so that (up to order n)
 * t(x)=\prod_{j=1}^{n}{eta_+(x^j)^{p_j}}
 * where eta_+(x) = prod(k>0, (1+x^k))
 */
    local(n, t);
    v = ser2prodplus(v);
    n = length(v);
    for (k=1, n,
        t = v[k];
        forstep (j=k+k, n, k, v[j]-=t; );
    );
    return( v );
} /* ----- */

etaprodplus2ser(p, w='y)=
{
    local(n, v);
    n = length(p);
    v = prod(k=1,#p, ( eta( (w+O(w^n))^(2*k) )/eta( (w+O(w^n))^k ) )^p[k] );
    return( v );
} /* ----- */


aux_printvpos(v, es, red=1)=
{
    local(p, ct=0);
\\    print1("( ");
\\    print(" :: v,es,red= ",[v,es,red]);
    for (k=1, #v,
        p = v[k];
        p = eval("Vec(p)[1]");  \\ trick to cast type to int
\\        print(" :: p= ",p);
\\        print(" :: p= ",type(p));

        if ( p > 0,
            if ( ct!=0,  print1("*") );
            ct+=1;
            print1(es,k/red);
            if ( p!=1,
                print1("^");
                if ( type(p)==type(1),
                    print1(p);
                , /*else */
                    print1("(",p,")");
                );
            );
        );
    );
\\    print1(" )");
    return( ct );  \\ == 0 if nothing printed
} /* ----- */


printetaprod(v, es="E", red=1)=
{
    local( f );
    f = aux_printvpos(v,es,red);
    if ( 0==f, print1("1") );
    print1(" / ");
    print1("( ");
    f = aux_printvpos(-v,es,red);
    if ( 0==f, print1("1") );
    print1(" ) ");
    print();
} /* ---- */

aux_printvpos_tex(v, es)=
{
    local(p);
    print1("{ 1");
    for (k=1, #v,
        p = v[k];
        if ( p>0, print1("\\,",es,"_{",k,"}");
            if ( p!=1, print1("^{",p,"}"); );
        );
    );
    print1(" }");
} /* ----- */


printetaprod_tex(v, es="E")=
{
    print1("\\frac");
    aux_printvpos_tex(v,es);
    aux_printvpos_tex(-v,es);
    print();
} /* ---- */

printprod(p, texq=0, pm=-1)=
{
    local(dq, op);
    op = if ( pm<0, "-", "+");
    if ( 0==texq,
        print1("  1");
        for (j=1, #p, if ( p[j], print1(" * (1", op, "y^",j,")^{",p[j],"}") ) );
    , /* else */
        dq = 0; \\ have terms in denom?
        for (j=1, #p, if ( p[j]<0, dq=1 ) );
        if ( dq, print1("\\frac{") );
        for (j=1, #p, if ( p[j]>0, print1("\\, (1", op, "y^",j,")^{",p[j],"}" ) ) );
        if ( dq,
            print1("}\n  {");
            for (j=1, #p, if ( p[j]<0, print1("\\, (1", op, "y^",j,")^{",p[j],"}" ) ) );
            print1("}");
        );
    );
    print();
    return();
} /* ----- */


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