The partial fraction decomposition of a univariate rational function:

\[f(x) = \frac{p(x)}{q(x)}\]

where \(p\) and \(q\) are co-prime and \(\deg(p) < \deg(q)\), is an expression of the form:

\[\sum_{i=1}^k \sum_{j=1}^{n_i} \frac{a_{ij}(x)}{q_i^j(x)}\]

where \(q_i\) for \(i=1 \ldots k\) are factors (e.g. over rationals or Gaussian rationals) of \(q\):

\[q(x) = \prod_{i=1}^k q_i^{n_i}\]

If \(p\) and \(q\) aren’t co-prime, we can use `cancel()` to remove common
factors and if \(\deg(p) >= \deg(q)\), then `div()` can be used to extract
the polynomial part of \(f\) and reduce the degree of \(p\).

Suppose we would like to compute partial fraction decomposition of:

```
>>> f = 1/(x**2*(x**2 + 1))
>>> f
1
───────────
2 ⎛ 2 ⎞
x ⋅⎝x + 1⎠
```

This can be achieved with SymPy’s built-in function `apart()`:

```
>>> apart(f)
1 1
- ────── + ──
2 2
x + 1 x
```

We can use `together()` to verify this result:

```
>>> together(_)
1
───────────
2 ⎛ 2 ⎞
x ⋅⎝x + 1⎠
```

Now we would like to compute this decomposition step-by-step. The rational
function \(f\) is already in factored form and has two factors \(x^2\) and
\(x^2 + 1\). If \(f\) was in expanded from, we could use `factor()` to
obtain the desired factorization:

```
>>> numer(f)/expand(denom(f))
1
───────
4 2
x + x
>>> factor(_)
1
───────────
2 ⎛ 2 ⎞
x ⋅⎝x + 1⎠
```

Based on the definition, the partial fraction expansion of \(f\) will be of the following form:

\[\frac{A}{x} + \frac{B}{x^2} + \frac{C x + D}{x^2 + 1}\]

Let’s do this with SymPy. We will use undetermined coefficients method to solve this problem. Let’s start by defining some symbols:

```
>>> var('A:D')
(A, B, C, D)
```

We use here the lexicographic syntax of `var()`. Next we can define three
rational functions:

```
>>> p1 = A/x
>>> p2 = B/x**2
>>> p3 = (C*x + D)/(x**2 + 1)
>>> p1, p2, p3
⎛A B C⋅x + D⎞
⎜─, ──, ───────⎟
⎜x 2 2 ⎟
⎝ x x + 1⎠
```

Let’s add them together to get the desired form:

```
>>> h = sum(_)
>>> h
A B C⋅x + D
─ + ── + ───────
x 2 2
x x + 1
```

The next step is to rewrite this expression as rational function in \(x\):

```
>>> together(h)
⎛ 2 ⎞ ⎛ 2 ⎞ 2
A⋅x⋅⎝x + 1⎠ + B⋅⎝x + 1⎠ + x ⋅(C⋅x + D)
────────────────────────────────────────
2 ⎛ 2 ⎞
x ⋅⎝x + 1⎠
>>> factor(_, x)
3 2
A⋅x + B + x ⋅(A + C) + x ⋅(B + D)
─────────────────────────────────
2 ⎛ 2 ⎞
x ⋅⎝x + 1⎠
```

Let’s now visually compare the last expression with \(f\):

```
>>> Eq(_, f)
3 2
a⋅x + b + x ⋅(a + c) + x ⋅(b + d) 1
───────────────────────────────── = ───────────
2 ⎛ 2 ⎞ 2 ⎛ 2 ⎞
x ⋅⎝x + 1⎠ x ⋅⎝x + 1⎠
```

Our task boils down to finding \(A\), \(B\), \(C\) and \(D\). We notice that denominators are equal so we will proceed only with numerators:

```
>>> eq = Eq(numer(_.lhs), numer(_.rhs))
>>> eq
3 2
a⋅x + b + x ⋅(a + c) + x ⋅(b + d) = 1
```

To solve this equation, we use `solve_undetermined_coeffs()`:

```
>>> solve_undetermined_coeffs(eq, [A, B, C, D], x)
{A: 0, B: 1, C: 0, D: -1}
```

This gave us values for our parameters, which now can be put into the initial expression:

```
>>> h.subs(_)
1 1
- ────── + ──
2 2
x + 1 x
```

This result is identical to the result we got from `apart(f)`. Suppose
however, we would like to see how undetermined coefficients method works.
First we have to extract coefficients of \(x\) of both sides of the equation:

```
>>> lhs, rhs = Poly(eq.lhs, x), Poly(eq.rhs, x)
>>> lhs
Poly((A + C)*x**3 + (B + D)*x**2 + A*x + B, x, domain='ZZ[A,B,C,D]')
>>> rhs
Poly(1, x, domain='ZZ')
```

Now we can use `Poly.nth()` to obtain coefficients of \(x\):

```
>>> [ Eq(lhs.nth(i), rhs.nth(i)) for i in xrange(4) ]
[b = 1, a = 0, b + d = 0, a + c = 0]
```

Solving this system of linear equations gives the same solution set as previously:

```
>>> solve(_)
{a: 0, b: 1, c: 0, d: -1}
>>> f.subs(_)
1 1
- ────── + ──
2 2
x + 1 x
```

There are several other ways we can approach undetermined coefficients
method. For example we could use `collect()` for this:

```
>>> collect(eq.lhs - eq.rhs, x, evaluate=False)
⎧ 2 3 ⎫
⎨1: B - 1, x: A, x : B + D, x : A + C⎬
⎩ ⎭
>>> solve(_.values())
{A: 0, B: 1, C: 0, D: -1}
```

Notice that even though the expressions were not `Eq()`‘s, this still
worked. This is because SymPy assumes by default that expressions are
identically equal to 0, so `solve(Eq(expr, 0))` is the same as
`solve(expr)`.

This approach is even simpler than using `Poly.nth()`. Finally we use a
little trick with `Symbol` and visually present solution to partial
fraction decomposition of \(f\):

```
>>> Eq(Symbol('apart')(f), f.subs(_))
⎛ 1 ⎞ 1 1
apart⎜───────────⎟ = - ────── + ──
⎜ 2 ⎛ 2 ⎞⎟ 2 2
⎝x ⋅⎝x + 1⎠⎠ x + 1 x
```