Here is a rule set that will do the job:
[ a*(b + c) := a*b + a*c,
opt(a) O(x^n) + opt(b) O(x^m) := O(x^n) :: n <= m
:: constant(a) :: constant(b),
opt(a) O(x^n) + opt(b) x^m := O(x^n) :: n <= m
:: constant(a) :: constant(b),
a O(x^n) := O(x^n) :: constant(a),
x^opt(m) O(x^n) := O(x^(n+m)),
O(x^n) O(x^m) := O(x^(n+m)) ]
If we really want the + and * keys to operate naturally
on power series, we should put these rules in EvalRules. For
testing purposes, it is better to put them in a different variable,
say, O, first.
The first rule just expands products of sums so that the rest of the
rules can assume they have an expanded-out polynomial to work with.
Note that this rule does not mention `O' at all, so it will
apply to any product-of-sum it encounters--this rule may surprise
you if you put it into EvalRules!
In the second rule, the sum of two O's is changed to the smaller O. The optional constant coefficients are there mostly so that `O(x^2) - O(x^3)' and `O(x^3) - O(x^2)' are handled as well as `O(x^2) + O(x^3)'.
The third rule absorbs higher powers of `x' into O's.
The fourth rule says that a constant times a negligible quantity is still negligible. (This rule will also match `O(x^3) / 4', with `a = 1/4'.)
The fifth rule rewrites, for example, `x^2 O(x^3)' to `O(x^5)'. (It is easy to see that if one of these forms is negligible, the other is, too.) Notice the `x^opt(m)' to pick up terms like `x O(x^3)'. Optional powers will match `x' as `x^1' but not 1 as `x^0'. This turns out to be exactly what we want here.
The sixth rule is the corresponding rule for products of two O's.
Another way to solve this problem would be to create a new "data type"
that represents truncated power series. We might represent these as
function calls `series(coefs, x)' where coefs is
a vector of coefficients for x^0, x^1, x^2, and so
on. Rules would exist for sums and products of such series
objects, and as an optional convenience could also know how to combine a
series object with a normal polynomial. (With this, and with a
rule that rewrites `O(x^n)' to the equivalent series form,
you could still enter power series in exactly the same notation as
before.) Operations on such objects would probably be more efficient,
although the objects would be a bit harder to read.
Some other symbolic math programs provide a power series data type
similar to this. Mathematica, for example, has an object that looks
like `PowerSeries[x, x0, coefs, nmin,
nmax, den]', where x0 is the point about which the
power series is taken (we've been assuming this was always zero),
and nmin, nmax, and den allow pseudo-power-series
with fractional or negative powers. Also, the PowerSeries
objects have a special display format that makes them look like
`2 x^2 + O(x^4)' when they are printed out. (See section Compositions,
for a way to do this in Calc, although for something as involved as
this it would probably be better to write the formatting routine
in Lisp.)