############################################################
# Manipulations with power series 
#
# Following "Handbook of Mathematical Functions"
#    Abramowitz and Stegun
#    Section 3.6, Infinite Series
############################################################

s[1] := 1 + a[1]*x + a[2]*x^2 + a[3]*x^3 + a[4]*x^4;
s[2] := 1 + b[1]*x + b[2]*x^2 + b[3]*x^3 + b[4]*x^4;
s[3] := 1 + c[1]*x + c[2]*x^2 + c[3]*x^3 + c[4]*x^4;

unknowns := {c[1],c[2],c[3],c[4]};

############################################################
# Define a 'shorthand' procedure for converting an
# expression (in the current case a series) to a polynomial
############################################################
P := proc(x) convert(x,polynom) end;

############################################################
# Solves for coefficients c_i in terms of a_i and b_i.
# 'series_in' should be an expression in s[1] and s[2].
# Note the use of global variables (s[3],unknowns).
############################################################
series_op := proc(series_in)
   solve({coeffs(P(s[3]) - 
                 P(series(series_in,x=0,5)),x)},unknowns);
end:

AS_3_6_16 := series_op( 1 / s[1] );
AS_3_6_17 := series_op( 1 / s[1]^2 );
AS_3_6_18 := series_op( sqrt(s[1]) );
AS_3_6_19 := series_op( 1 / sqrt(s[1]) );
AS_3_6_20 := series_op( s[1]^n );
AS_3_6_21 := series_op( s[1] * s[2] );
AS_3_6_22 := series_op( s[1] / s[2] );
AS_3_6_23 := series_op( exp(s[1] - 1) );
AS_3_6_24 := series_op( 1 + ln(s[1]) );