\\% Basic manipulations of vectors.
\\ Author: Ralf Stephan, (reformatted by Joerg Arndt)
\\ License: GPL version 3 or later
\\ online at http://www.jjj.de/pari/
\\ version: 2011-April-08 (19:18)


\\ gp/pari code for Doron's  challenge at http://mathoverflow.net/questions/71092/how-many-integer-partitions-of-a-googol-10100-into-at-most-60-parts/71132#71132
\\Needs http://www.jjj.de/pari/fastrec.gpi

ggf(v)=
{ /* Guess generating function for supplied vector of integers */
/* Example: ggf(vector(10,j,fibonacci(j)^2)) ==> (-x + 1)/(x^3 - 2*x^2 - 2*x + 1) */
    local(l,m,p,q,B);
    l = length(v);
    B = floor(l/2);
    if ( B<4,  return(0) );
    m = matrix(B, B, x, y, v[x-y+B+1]);
    q = qflll(m, 4)[1];
    if ( length(q)==0,  return(0) );
\\    p = sum(k=1, B, x^(k-1)*q[k,1]);
    p = sum(k=1, B, 'x^(k-1)*q[k,1]);
\\    q = Pol( Pol( vector(l, n, v[l-n+1]) )*p + O(x^(B+1)) );
    q = Pol( Pol( vector(l, n, v[l-n+1]) )*p + O('x^(B+1)) );
    if ( polcoeff(p,0)<0,  q=-q; p=-p);
    q = q/p;
\\    p = Ser(q+O(x^(l+1)));
    p = Ser(q+O('x^(l+1)));
    for (m=1, l, if( polcoeff(p,m-1)!=v[m],  return(0)) );
    return(q);
} /* ------- */


\\ aux:
strsgn(a)=if(a>=0, Str("+" a), Str(a));
\\prtpol(c)=
\\{ /* print coefficients with ascending degrees */
\\    print1("(");
\\    for (d=(poldegree(c)-...)
\\    print1(")");
\\} /* ----- */
\\
\\prtgf(v)=
\\{ /* Print Generating function */
\\    local(gf, t);
\\    gf = ggf(v);
\\    if ( gf==0, return(0) );
\\    t = numerator(gf);
\\    print1
\\    return(1);
\\} /* ----- */


prtrec(v)=
{ /* Print recurrence relation */
    local(c,d,c0,t);
    c=ggf(v);
    if ( c==0, return(0) );
    c = denominator(c);
\\    print(c);
\\    c = polrecip(c);
\\    print(c);
    d = poldegree(c);
    c0 = polcoeff(c,0);
    print1("a(n) =");
    for (k=1, d,
        t=polcoeff(c,k)/c0;
        if (t, print1(" "); print1( strsgn( -t ), "*a(n-",k,")" ) );
    );
    print();
    return(1);
} /* ----- */

prtrec2(c)=
{ /* Print recurrence relation */
    local(d,c0,t,r,t2);
    \\c=ggf(v);
    if ( c==0, return(0) );
    c = denominator(c);
\\    print(c);
\\    c = polrecip(c);
\\    print(c);
    d = poldegree(c);
    r=[];\\listcreate(d);
    c0 = polcoeff(c,0);
    print1("a(n) =");
    forstep (k=1, d,1,
        t=polcoeff(c,k)/c0;
        t2=-polcoeff(c,k)/c0;
	r=concat(r,t2);
        if (t, print1(" "); print1( strsgn( -t ), "*aa(n-",k,")" ) );
    );
    print();
    return(r);
} /* ----- */

k4=60;
d=prod(n=1,k4,(1-x^n));
{aa(n)=return(polcoeff(1/d+O(x^(n+1)),n ));};
zz=prtrec2(1/d);
a0=[];\\listcreate(poldegree(d));
for(n=0,poldegree(d)-1,a0=concat(a0,aa(n)); );
{getd()=d;}
{geta0()=a0;}
{getzz()=zz;}

X=frec(a0,zz,10^100)[1];

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