> This is a reason I am sometimes confused if the series should even be addition of functions + an order term to make an object of the add class (rather than having its own class).

There is no need to be confused about this. Implementing series as an Add with an O term is a bad idea.

A Series should be a special object of its own and it should be like Poly rather than like Expr. The internal representation should be like a Poly i.e. with an actual ordered list of coefficients. The coefficients should be represented using domain elements in a canonical form, always expanded etc just like for Poly. There should be no need to waste time debating whether to call expand, cancel etc because the domain elements are always canonical.

The ring_series module shows how to do this but is too low-level (gh-26100). It needs to be wrapped in something higher-level like Poly. Also python-flint should be used when possible:

In [7]: from flint import fmpq_series

In [8]: x = fmpq_series([0, 1])

In [9]: x
Out[9]: x + O(x^10)

In [10]: 1/(x + 1)
Out[10]: 1 + (-1)*x + x^2 + (-1)*x^3 + x^4 + (-1)*x^5 + x^6 + (-1)*x^7 + x^8 + (-1)*x^9 + O(x^10)

In [11]: x.exp()
Out[11]: 1 + x + 1/2*x^2 + 1/6*x^3 + 1/24*x^4 + 1/120*x^5 + 1/720*x^6 + 1/5040*x^7 + 1/40320*x^8 + 1/362880*x^9 + O(x^10)

It is much faster:

In [1]: from flint import ctx

In [2]: ctx.cap = 1000

In [3]: from flint import fmpq_series

In [4]: x = fmpq_series([0, 1])

In [5]: %time r = x.exp()
CPU times: user 5.04 ms, sys: 79 μs, total: 5.12 ms
Wall time: 6.18 ms

In [6]: len(r)
Out[6]: 1000

Compare:

In [12]: %time ok = exp(x).series(x, n=1000)
CPU times: user 19.2 s, sys: 23.3 ms, total: 19.3 s
Wall time: 19.3 s

That's 3000x slower.

It should be possible to have multivariate series as well (gh-25483).

Originally posted by @oscarbenjamin in #26941 (comment)