0

I want to invert a power series with PARI/GP:
if $y=a[1]*x+a[2]*x^2+...+O(x^N)$ is encoded by a given array, say a=vector(N-1), I want PARI to produce the vector b such that $x=b[1]*y+b[2]*y^2+...+O(y^N)$.
This can be done using the Bell polynomials. The Pari handbook has

Bell(k,n=-1)=
{
my(var(i)=eval(Str("X",i)));
my(x,v,dv);
v=vector(k,i,if(i==1,’E,var(i-1)));
dv=vector(k,i,if(i==1,’X*var(1)*’E,var(i)));
x=diffop(’E,v,dv,k)/’E;
if(n<0,subst(x,’X,1),polcoeff(x,n,’X))
}

producing e.g. 

gp > Bell(3)
%3 = X1^3 + 3*X2*X1 + X3

But I have no idea how to use those, i.e. to attribute values to X1, X2, ... in a subsequent formula in Pari (in fact, I hardly understand how eval() and subst() work in the above formula!). That should be trivial for someone who knows... Help please!

Adrian Mole
  • 49,934
  • 160
  • 51
  • 83
Wolfgang
  • 103
  • 2

1 Answers1

1

PARI has a built in function to do power series reversion, so there is no need to get into a low level implementation which is certainly very sub-optimal and complicated.

Example use:

serreverse( x / (1 + x)^2 + O(x^20))

A second example (based on an example in the link given in the question)

serreverse(atan(x + O(x^20)))
tan(x + O(x^20))

If your series is in a vector v you would first convert it to a power series and then call serreverse. For example:

serreverse(Ser(v))

To convert back to a vector later use the Vec function. When converting from a power series to a vector PARI will strip any leading zeros, which can be annoying. To prevent this there is a second argument which you typically give as -n or -(n+1) where n is the number of terms in the power series.

Andrew
  • 873
  • 4
  • 10