The Archimedes Cattle Problem
The Archimedes Cattle Problem is a well known mathematical problem attributed to Archimedes. The problem, which is described both here and here, is interesting because it suggests a remarkable early insight into mathematics that was only fully developed many centuries later in the middle ages. In particular, it suggests that, while those who posed the problem could not give a solution, they were confident enough in the mathematics involved to know that a solution did exist. The problem, the history involved its solution and methods for it solution are reviewed in this paper by H.W. Lenstra Jr.
I thought it would be interesting to produce a solution in Python.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 |
from fractions import Fraction as RF from functools import reduce from operator import mul from math import lcm from number_theory import factor, Qde # print full output if full = True full = False # Crout LU decomposition (which uses an implied unit # diagonal for the upper triangular matrix) # # S = source matrix, D = destination matrix # # if D is present return the LU decomposition with # S unchanged; otherwise the LU decomposition will # be computed in place and returned in S def crout(S, D=None): if D is None: D = S n = len(S) for k in range(n): for r in range(k, n): D[r][k] = S[r][k] - sum(D[r][j] * D[j][k] for j in range(k)) t = RF(1, D[k][k]) for c in range(k + 1, n): D[k][c] = t * (S[k][c] - sum(D[k][j] * D[j][c] for j in range(k))) # Solve the system L.U.x = v for the vector (v) def solve_crout(LU, v): n = len(LU) x = [0] * n for i in range(n): t = RF(1, LU[i][i]) x[i] = t * (v[i] - sum(LU[i][j] * x[j] for j in range(i))) for i in reversed(range(n)): x[i] -= sum(LU[i][j] * x[j] for j in range(i + 1, n)) return x # the fractions in the Archimedes Cattle Problem nd = ((5, 6), (9, 20), (13, 42), (7, 12), (9, 20), (11, 30)) f1, f2, f3, f4, f5, f6 = (RF(n, d) for n, d in nd) # First step: establish the equations set out on Wikipedia at: # https://en.wikipedia.org/wiki/Archimedes%27s_cattle_problem # W B D Y w b d y M = [ [ -1, f1, 0, 1, 0, 0, 0, 0 ], [ 0, -1, f2, 1, 0, 0, 0, 0 ], [ f3, 0, -1, 1, 0, 0, 0, 0 ], [ 0, f4, 0, 0, -1, f4, 0 , 0 ], [ 0, 0, f5, 0, 0, -1, f5, 0 ], [ 0, 0, 0, f6, 0, 0, -1, f6 ], [ f3, 0, 0, 0, f3, 0, 0, -1 ], [ 1, 1, 1, 1, 1, 1, 1, 1 ] ] # the right hand side R = [ 0, 0, 0, 0, 0, 0, 0, 1 ] # solve and scale the results to give integer values crout(M) S = solve_crout(M, R) d = lcm(*(x.denominator for x in S)) S = tuple(int(x * d) .denominator for x in S) if full: for n, s in zip('WBDYwbdy', S): print(f"{n} = {s:>10_} * k") print() # Part 2(1): B + W is a perfect square. All solutions are # multiples (k) of the above integer values for (B, W, ...) # and k must now be chosen such that k.(B + W) is a perfect # square. To ensure this, k itself must be the product of the # square free factors of B + W multiplied by any square (q) k_sqfp = reduce(mul, (f for f, x in factor(S[0] + S[1]) if x & 1)) # Part 2(2): k.(D + Y) is a triangular number t.(t + 1) / 2 # where k is derived above as k_sqfp.q^2, hence giving: # # (2t + 1)^2 - 8.(D + Y).bw_sqfp.q^2 = 1 # # This is a quadratic diophantine equation r^2 - C.q^2 = 1 with # C = 8.(D + Y).bw_sqfp; now use Amthor's strategy by expressing # C as a product of square and square-free components (with the # former as a square root) fctrs = factor(8 * (S[2] + S[3]) * k_sqfp) sqfc = reduce(mul, (f for f, x in fctrs if x & 1)) srtc = reduce(mul, (f ** t for f, x in fctrs if (t:= x // 2))) # we now have to solve the quadratic diophantine equation # r^2 - sqfc.(srtc.s)^2 where the second component of the # solution is a multiple of srtc p = Qde(sqfc, 1) for r, sm in p.solve(): # sm must be a multiple of srtc m, rem = divmod(sm, srtc) if m and not rem: # k = k_sqfp.m^2 - now compute the number of cattle in each # category and the overall total (only output the first and # last 20 digits) if full: for n, s in zip('WBDYwbdy', S): c = str(s * k_sqfp * m ** 2) print(f"{n} = {c[:20]} ..{len(c) - 40:}.. {c[-20:]}") c = str(sum(S) * k_sqfp * m ** 2) print(f"T = {c[:20]} ..{len(c) - 40}.. {c[-20:]}") break |
No comments yet