1 | initial version |

One more possible approach is to use a Groebner basis. Unfortunately, after finding such a basis, I'm still getting non-real solutions. But, again, maybe you can modify this approach to find what you want.

```
R = PolynomialRing(QQ, 4, 'abcd', order='lp')
a,b,c,d=R.gens()
I=(a*b*c*d - a*b - a*c - b*c - a*d - b*d + 2*a + 2*b - 448, -a*b*c - a*b*d - a*c*d - b*c*d + 3*a + 3*b + 2*c + 2*d + 452, a*b + a*c + b*c + a*d + b*d + c*d - 159, -a - b - c - d + 21)*R
B = I.groebner_basis()
print B
```

Copyright Sage, 2010. Some rights reserved under creative commons license. Content on this site is licensed under a Creative Commons Attribution Share Alike 3.0 license.