利用Python求解阿基米德分牛问题

import sympy
x0,x1,x2,x3 = sympy.symbols("x0,x1,x2,x3")
y0,y1,y2,y3 = sympy.symbols("y0,y1,y2,y3")
x = [x0,x1,x2,x3]
y = [y0,y1,y2,y3]

frac = lambda x : sympy.Rational(1,x)
fs = []
for i in range(3):
    fs.append(x[i]-x[3]-(frac(2*i+2)+frac(2*i+3))*x[i+1])

for i in range(4):
    ind = (i + 1) % 4
    fs.append(y[i]-(frac(i+3)+frac(i+4))*(x[ind]+y[ind]))

>>> for f in fs: print(f)
...
x0 - 5*x1/6 - x3
x1 - 9*x2/20 - x3
x2 - 55*x3/42
-7*x1/12 + y0 - 7*y1/12
-9*x2/20 + y1 - 9*y2/20
-11*x3/30 + y2 - 11*y3/30
-13*x0/42 - 13*y0/42 + y3

res = sympy.solve(fs, x[:3]+y)

for key in res:
    print(sympy.latex(key), "&=", sympy.latex(res[key]), r"\\")

denominators = [(v/x3).q for v in res.values()]
x3Res = sympy.lcm(denominators)
# 32859792

fs.append(x3-x3Res)
res2 = sympy.solve(fs, x+y)
for key in res2:
    print(sympy.latex(key), "=", res2[key], r"\\")