Partial fraction decomposition using SymPy / Python

You can do this using apart in sympy but apart will look for a rational factorisation by default so you have to tell it to work in Q(sqrt(3)):

In [37]: apart(1/(x**6+1))                                                                                                                     
Out[37]: 
        2                     
       x  - 2           1     
- ─────────────── + ──────────
    ⎛ 4    2    ⎞     ⎛ 2    ⎞
  3⋅⎝x  - x  + 1⎠   3⋅⎝x  + 1⎠

In [36]: apart(1/(x**6+1), extension=sqrt(3))                                                                                                  
Out[36]: 
       √3⋅x - 2            √3⋅x + 2           1     
- ───────────────── + ───────────────── + ──────────
    ⎛ 2           ⎞     ⎛ 2           ⎞     ⎛ 2    ⎞
  6⋅⎝x  - √3⋅x + 1⎠   6⋅⎝x  + √3⋅x + 1⎠   3⋅⎝x  + 1⎠

EDIT: The comments are asking for a way to find this in more general cases without needing to know that sqrt(3) generates the extension.

We can use .apart(full=True) to compute the full PFE over linear factors:

In [29]: e = apart(1/(x**6+1), full=True).doit()

In [30]: e
Out[30]: 
       √3   ⅈ           √3   ⅈ          √3   ⅈ           √3   ⅈ                            
     - ── - ─         - ── + ─          ── - ─           ── + ─                            
       2    2           2    2          2    2           2    2           ⅈ           ⅈ    
- ────────────── - ────────────── - ────────────── - ────────────── + ───────── - ─────────
    ⎛    √3   ⅈ⎞     ⎛    √3   ⅈ⎞     ⎛    √3   ⅈ⎞     ⎛    √3   ⅈ⎞   6⋅(x + ⅈ)   6⋅(x - ⅈ)
  6⋅⎜x + ── + ─⎟   6⋅⎜x + ── - ─⎟   6⋅⎜x - ── + ─⎟   6⋅⎜x - ── - ─⎟                        
    ⎝    2    2⎠     ⎝    2    2⎠     ⎝    2    2⎠     ⎝    2    2⎠ 

This is not what is wanted because you want to have real quadratic denominators rather than introducing complex numbers. The terms here come in complex conjugate pairs though so we can combine the pairs:

In [46]: terms1 = []

In [47]: for a in terms:
    ...:     if conjugate(a) not in terms1:
    ...:         terms1.append(a)
    ...: 

In [49]: terms_real = [(t+t.conjugate()) for t in terms1]

In [51]: Add(*(factor(cancel(t)) for t in terms_real))
Out[51]: 
      √3⋅x - 2            √3⋅x + 2            1     
- ───────────────── + ───────────────── + ──────────
    ⎛ 2           ⎞     ⎛ 2           ⎞     ⎛ 2    ⎞
  6⋅⎝x  - √3⋅x + 1⎠   6⋅⎝x  + √3⋅x + 1⎠   3⋅⎝x  + 1⎠

Note that in general there might not be any simple expressions for the roots (Abel-Ruffini) so this kind of expression for the partial fraction expansion will not succeed for all possible denominator polynomials. This is why .apart by default computes an expansion over irreducible denominators (something that can always succeed).