Partial fraction decomposition using SymPy / Python

Question

How do I find the constants A,B,C,D,K,S such that

1/(x**6+1) = (A*x+B)/(x**2+1) + (C*x+D)/(x**2-sqrt(3)*x+1) + (K*x+S)/(x**2+sqrt(3)*x+1)

is true for every real x.

I need some sympy code maybe, not sure. Or any other Python lib which could help here.

I tried by hand but it's not easy at all: after 1 hour of calculating, I found that I have probably made some mistake.

I tried partial fraction decomposition in SymPy but it does not go that far.

I tried Wolfram Alpha too, but it also does not decompose to that level of detail, it seems.

WA attempt

See the alternate forms which WA gives below.

Edit

I did a second try entirely by hand and I got these:

A = 0
B = 1/3
C = -1/(2*sqrt(3))
D = 1/3
K = 1/(2*sqrt(3))
S = 1/3

How can I verify if these are correct?

Edit 2

The main point of my question is: how to do this with some nice/reusable Python code?

Answer 1

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⎠