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