用Python优雅地求解阿基米德分牛问题
文章目录
- 题目大意
- sympy求解
- 结果
题目大意
问 太阳神有一牛群，由白、黑、花、棕四种颜色的公、母牛组成，其间关系如下，求每种牛的个数。
- 公牛中，白牛多于棕牛，二者之差为黑牛的12+13\frac{1}{2}+\frac{1}{3}21+31；黑牛多于棕牛，二者之差为花牛的14+15\frac{1}{4}+\frac{1}{5}41+51；花牛多于棕牛，二者之差为白牛数的16+17\frac{1}{6}+\frac{1}{7}61+71
- 母牛中，白牛是全体黑牛的13+14\frac{1}{3}+\frac{1}{4}31+41；黑牛是全体花牛的14+15\frac{1}{4}+\frac{1}{5}41+51；花牛是全体棕牛的15+16\frac{1}{5}+\frac{1}{6}51+61；棕牛是全体白牛的16+17\frac{1}{6}+\frac{1}{7}61+71
如果用字母x0,x1,x2,x3x_0, x_1, x_2, x_3x0,x1,x2,x3分别表示白、黑、花、棕各色的公牛数；用y0,y1,y2,y3y_0, y_1, y_2, y_3y0,y1,y2,y3分别表示白、黑、花、棕各色母牛数，则得8 个未知数的如下7 个方程
x0−x3=(12+13)x1x1−x3=(14+15)x2x2−x3=(16+17)x0y0=(13+14)(x1+y1)y1=(14+15)(x2+y2)y2=(15+16)(x3+y3)y3=(16+17)(x0+y0)\begin{aligned} x_0-x_3=(\frac{1}{2}+\frac{1}{3})x_1\\ x_1-x_3=(\frac{1}{4}+\frac{1}{5})x_2\\ x_2-x_3=(\frac{1}{6}+\frac{1}{7})x_0\\ y_0=(\frac{1}{3}+\frac{1}{4})(x_1+y_1)\\ y_1=(\frac{1}{4}+\frac{1}{5})(x_2+y_2)\\ y_2=(\frac{1}{5}+\frac{1}{6})(x_3+y_3)\\ y_3=(\frac{1}{6}+\frac{1}{7})(x_0+y_0)\\ \end{aligned} x0−x3=(21+31)x1x1−x3=(41+51)x2x2−x3=(61+71)x0y0=(31+41)(x1+y1)y1=(41+51)(x2+y2)y2=(51+61)(x3+y3)y3=(61+71)(x0+y0)
这个题其实是毫无难度的，但非要用Python，那么难点主要如何优雅地表达这个过程，这里选用的是sympy符号计算。
所以第一步，先给定一些符号
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]
sympy求解
然后将阿基米德分牛问题转化为Python代码，其优雅之处在于，这些分数的构建遵循自然数递增的规律，故可通过循环来生成，非常便捷。
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) % 4fs.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
但是，8个未知数7个方程，显然没有唯一解，考虑到x3x_3x3貌似是最小的值，所以最后希望用x3x_3x3来表示其他数。
res = sympy.solve(fs, x[:3]+y)
结果
查看一下结果
for key in res:print(sympy.latex(key), "&=", sympy.latex(res[key]), r"\\")
x0=781x3336x1=89x356x2=55x342y0=2316515x31564752y1=1731719x31825544y2=1639880x32053737y3=806221x3684579\begin{aligned} x_{0} &= \frac{781 x_{3}}{336} \\ x_{1} &= \frac{89 x_{3}}{56} \\ x_{2} &= \frac{55 x_{3}}{42} \\ y_{0} &= \frac{2316515 x_{3}}{1564752} \\ y_{1} &= \frac{1731719 x_{3}}{1825544} \\ y_{2} &= \frac{1639880 x_{3}}{2053737} \\ y_{3} &= \frac{806221 x_{3}}{684579} \\ \end{aligned} x0x1x2y0y1y2y3=336781x3=5689x3=4255x3=15647522316515x3=18255441731719x3=20537371639880x3=684579806221x3
这道题到这里基本上就算解完了，但是牛至少得是个整数，所以接下来要做的是求解分母的最小公倍数。
在sympy中，对于一个分数r，r.p为分子，r.q为分母；lcm可求解其最小公倍数。
denominators = [(v/x3).q for v in res.values()]
x3Res = sympy.lcm(denominators)
# 32859792
然后让将x3的值加入fs，
fs.append(x3-x3Res)
res2 = sympy.solve(fs, x+y)
for key in res2:print(sympy.latex(key), "=", res2[key], r"\\")
结果如下
x0=76379457x1=52223598x2=43030680x3=32859792y0=48646815y1=31170942y2=26238080y3=38698608x_{0} = 76379457 \\ x_{1} = 52223598 \\ x_{2} = 43030680 \\ x_{3} = 32859792 \\ y_{0} = 48646815 \\ y_{1} = 31170942 \\ y_{2} = 26238080 \\ y_{3} = 38698608 \\ x0=76379457x1=52223598x2=43030680x3=32859792y0=48646815y1=31170942y2=26238080y3=38698608
这些牛加一起有349247972头，全世界大概有10万亿头，看来太阳神的牛还是比较多的。
