############################################################
# Manipulations with power series 
#
# Following "Handbook of Mathematical Functions"
#    Abramowitz and Stegun (A & S)
#    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::anything) 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::anything)
   solve({coeffs(P(s[3]) - 
                 P(series(series_in,x=0,5)),x)},unknowns);
end:

############################################################
# Use 'series_op' procedure to reproduce various A & S
# formulae.
# 
# Note:  Due to an apparent bug in 'series' and 'taylor'
# in Maple V.5, the computations of AS_3_6_18 and AS_3_6_19 
# must be repeated---the first time through, 'series_op'
# fails, the second time it works. 
############################################################
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_18 := series_op( sqrt(s[1]) ):
AS_3_6_19 := series_op( 1 / 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]) ):