Archimedes's cattle problem

Question

Compute, O friend, the number of the cattle of the sun which once grazed upon the plains of Sicily, divided according to color into four herds, one milk-white, one black, one dappled and one yellow. The number of bulls is greater than the number of cows, and the relations between them are as follows:

White bulls \$= (\frac{1}{2} + \frac{1}{3})\$ black bulls + yellow bulls,
Black bulls \$= (\frac{1}{4} + \frac{1}{5})\$ dappled bulls + yellow bulls,
Dappled bulls \$= (\frac{1}{6} + \frac{1}{7})\$ white bulls + yellow bulls,
White cows \$= (\frac{1}{3} + \frac{1}{4})\$ black herd,
Black cows \$= (\frac{1}{4} + \frac{1}{5})\$ dappled herd,
Dappled cows \$= (\frac{1}{5} + \frac{1}{6})\$ yellow herd,
Yellow cows \$= (\frac{1}{6} + \frac{1}{7})\$ white herd.

If thou canst give, O friend, the number of each kind of bulls and cows, thou art no novice in numbers, yet can not be regarded as of high skill. Consider, however, the following additional relations between the bulls of the sun:

White bulls + black bulls = a square number,
Dappled bulls + yellow bulls = a triangular number.

If thou hast computed these also, O friend, and found the total number of cattle, then exult as a conqueror, for thou hast proved thyself most skilled in numbers.

- Archimedes

Some clarifications:

• black herd = black bulls + black cows, white herd = white bulls + white cows, etc.
• \$(\frac{1}{2} + \frac{1}{3})\$ black bulls + yellow bulls, means only the black bulls get the coefficient
• A square number is a number which can be represented as n * n where n is an integer
• A triangular number is a number which can be represented as 1 + 2 + ... + (n - 1) + n

Task

Write a program/function to output the size of the smallest herd that could satisfy both the first and second parts of the above problem.

Output

You must output this number. Standard Output rules for integers apply.

Scoring

This is code-golf so shortest bytes wins.

---

Hash of decimal representation of answer for easy checking:

• sha256 (with newline): 4fce7274734a18e9d27592e59480b50876da243de4caef4323f783339c456a7c
• sha256 (without newline): 8500b362f3d6b5e96c775ca16d1475b519ff18669b0f6f3d35e4fb23d2c18d01

---

Inspired by The Archimedes Number - Numberphile

Answer 1

Sledgehammer, 67 66 26 bytes

-40 bytes thanks to @GregMartin and I also no longer have any idea about how my answer works

Completes in less than a few seconds!

⡇⣄⠀⠇⣺⠇⢞⡞⣵⣍⠪⢺⡇⠜⢂⡒⢃⠦⠲⣎⠇⠣⡔⢻⡦⠔

Mathematica code:

Floor[Divide[25194541,184119152] * (NumberFieldFundamentalUnits@Sqrt[4729494])^4658]

(the reason for the Divide is that by default Mathematica represents x/y as x * y^-1, which is half a byte longer in Sledgehammer).

Answer 2

Wolfram Language (Mathematica), ⁸⁴ 83 bytes

224571490814418y^2/.{1}.FindInstance[x^2-410286423278424y^2==1&&x>1,{x,y},Integers]

Try it online!

^{-1 byte thanks to @J42161217.}

Gives the identical result to the existing Sledgehammer solution. Uses the Pell equation directly to find the required y, and substitutes into the formula for the desired value. One problem was that the OEIS didn't have the correct constant factor (which must be multiplied by 4456749).

This one should be easier to port to other languages, since the Pell equation can be brute-forced using just integers.

Deriving the formula

Start with the minimal solution of the linear equations, which is already present on MathWorld, where \$W,X,Y,Z\$ denote white, black, dappled, and yellow bulls and \$w,x,y,z\$ denote the cows, and \$S\$ is the grand total:

$$ W,X,Y,Z = 10366482, 7460514, 7358060, 4149387 \\ w,x,y,z = 7206360, 4893246, 3515820, 5439213 \\ W+X = 17826996, Y+Z = 11507447, S = 50389082 $$

Now we need to find the integer multiple of all the values such that \$(W+X)n\$ is a square and \$(Y+Z)n\$ is a triangular number (which can be alternatively written as "8 times that plus 1 is a square"):

$$ (W+X)n = x^2 \\ 8(Y+Z)n + 1 = y^2 $$

Note that \$W+X\$ is squarefree except for having the prime factor 2 twice, so \$x\$ must be a multiple of \$\frac{W+X}{2}\$. Substitute \$x=\frac{W+X}{2}x_1\$, then we get

$$ n = \frac{W+X}{4}x_1^2 $$

Substitute this into the equation for \$y\$ and rearrange a bit, then we get

$$ y^2 - 2(W+X)(Y+Z)x_1^2 = 1 $$

So the \$2(W+X)(Y+Z)\$ is where the number 410286423278424 comes from. Now assume we solved it; then we have the value for \$x_1\$. The final answer we want is \$S n\$, or

$$ S n = \frac{W+X}{4}Sx_1^2 $$

where \$\frac{W+X}{4}S\$ evaluates to 224571490814418.

Answer 3

JavaScript (Node.js),  373 ... 304  296 bytes

Returns a BigInt of 206545 digits.

Most BigInt literals in the code are stored as strings in base 119. This saves 8 bytes but leads to many unprintable characters. Below is a sanitized version without this compression scheme.

_=>[...1e9+[9542]].map(i=>M.push(m=([a,b,c]=m,[d,e,f]=M[i-2]||m,[v=a*d+b*f,a*e+b*d,u=c*d+a*f])),m=[0x4EDF512CD794532694B80D70C2648ADB08931n,0x119A739D926824D28537A6B609D64903072098n,0x16163EA6FB9A566AD5B17C9614A6476F10864n],M=[[0x104941B82B6E51BED5n,0n,0x48F880D7EEB3F6CAAn]])&&v*u*v*u*48222351474n/4657n

Try it online!

^{Note: Interestingly, it takes 6 times longer to turn the number into a string than to compute it (~2400ms and ~400ms on TIO respectively).}

How?

Preamble

Node is definitely not the right tool for the job. We need to use an algorithm relying on integers exclusively and we don't have any advanced math built-in at our disposal.

It would probably be possible to solve the Pell equation used by Bubbler, but given the magnitude of the number we're looking for, I doubt that a size-optimized version of such an algorithm would return anything in a reasonable amount of time.

In order to get something that actually works, I've decided to use the formula described in the next paragraph.

Formula

This is an implementation of the method described in A simple solution to Archimedes' cattle problem by Antti Nygrén, whose final formula is:

$$\begin{pmatrix}v\\u\end{pmatrix}=\begin{pmatrix}a_1&a_2\\a_3&a_1\end{pmatrix}^{1164}\times\begin{pmatrix}r_1\\r_2\end{pmatrix}$$

$$T=6299\times 957\times u^2v^2+\frac{21054639\times 957\times u^2v^2}{4657}$$

with:

$$\begin{align}&a_1=109931986732829734979866232821433543901088049\\ &a_2=392567302329690546856394748066206816187916440\\ &a_3=30784636507697855142356992218944109072681060\\ &r_1=300426607914281713365\\ &r_2=84129507677858393258\end{align}$$

Implementation

We start with m = [a1, a2, a3] and M = [[r1, 0, r2]].

We iterate through the list [...1e9 + [9542]] which expands to:

['1','0','0','0','0','0','0','0','0','0','9','5','4','2']

For each value \$i\$ in the above list:

• we build the vector [a,b,c] from m

• we build the vector [d,e,f] from either M[i-2] if it exists, or m otherwise

• we update m to [v,v',u] where:

  $$\begin{pmatrix}v&v'\\u&\color{grey}{u'}\end{pmatrix}=\begin{pmatrix}a&b\\c&a\end{pmatrix}\times\begin{pmatrix}d&e\\f&d\end{pmatrix}=\begin{pmatrix}ad+bf&ae+bd\\cd+af&\color{grey}{ce+ad}\end{pmatrix}$$

  ^{Note: we don't need \$u'\$, so it's not computed at all}

• we push the new instance of m into M

During the first 10 iterations, M[i-2] is undefined and m is used instead. So we essentially just square the original matrix 10 times. In other words, we compute:

$$\begin{pmatrix}a_1&a_2\\a_3&a_1\end{pmatrix}^{1024}$$

During the next 3 iterations, we multiply by previous occurrences of m: M[9-2], M[5-2] and M[4-2], which correspond to the original matrix to the power of \$128\$, \$8\$ and \$4\$ respectively in order to reach:

$$\begin{pmatrix}a_1&a_2\\a_3&a_1\end{pmatrix}^{1024+128+8+4}=\begin{pmatrix}a_1&a_2\\a_3&a_1\end{pmatrix}^{1164}$$

The last iteration is a multiplication by M[0] where we've stored [r1, 0, r2] at the very beginning of the process, leading to:

$$\begin{pmatrix}v&\color{grey}{v'}\\u&\color{grey}{u'}\end{pmatrix}=\begin{pmatrix}a_1&a_2\\a_3&a_1\end{pmatrix}^{1164}\times\begin{pmatrix}r_1&\color{grey}{0}\\r_2&\color{grey}{r_1}\end{pmatrix}$$

At this point, \$v\$ and \$u\$ are set correctly and we just need to compute:

$$\frac{v^2u^2\times 48222351474}{4657}$$

Answer 4

Python 2, 166 bytes

u=0x104941b82b6e51bed5
v=0x48f880d7eeb3f6caa
k=2*u*v
a=0x4edf512cd794532694b80d70c2648adb08931
exec"u,v=a*u+7766*k*v,a*v+609*k*u;"*1164
print u*u*v*v*48222351474/4657

Try it online!

Based on Arnauld's formula and method. I start with an initial u,v, and update them by a 2*2 matrix operation 1164 times. I looked for constants sharing large common factors and extracted those.