# Archimedes's cattle problem

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 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

• In each formula, there aren't parens around, say, black bulls + yellow bulls, meaning only the black bulls get the coefficient?
– xnor
Jun 25 '20 at 5:22
• Helios should have had very small cattle or a quite big idea of the size of Sicily... ;-) Jun 27 '20 at 15:00
• @Rmano Very very very very very small cattle - see this Feb 2 '21 at 9:58

# Sledgehammer, 6766 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).

• OEIS A096151 has some information: the value is 50389082*t(1)^2, where (s(1), t(1)) is the smallest nontrivial solution to s^2 - D*t^2 = 1, with D=410286423278424 (or smallest solution t divisible by 9314 for squarefree D=4729494). So while it is a primary solution to the Pell's equation with large D, it is also >2000th solution for a smaller D. How that turns into the exact formula is still a mystery. Jun 25 '20 at 6:24
• Found a shorter expression in Mathematica that gives the same result: 224571490814418*y^2 /. FindInstance[x^2-410286423278424*y^2==1 && x>1, {x, y}, Integers][[1]]. I don't know how it'll compress well, but it's definitely easier to port to other languages. Jun 25 '20 at 8:30
• @Bubbler After minor Sledgehammer-specific golfing, that seems to be 41 bytes. (besides, I think FindInstance is harder to port to other languages than simple mathematical operations with arbitrary precision!) Jun 25 '20 at 11:36
• I don't know Sledgehammer, but I'm guessing it will save many bytes to replace 109931986732829734979866232821433543901088049 + 50549485234315033074477819735540408986340 Sqrt[4729494] by the mathematically equivalent NumberFieldFundamentalUnits[Sqrt[4729494]]. (prepend with Last if you don't like the fact that this returns the answer as a single-element list) Jun 25 '20 at 17:31
• Jun 26 '20 at 1:35

# Wolfram Language (Mathematica), 84 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.

• 83 bytes Jun 25 '20 at 12:17

# 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}$$

• I might be missing something, but would it be shorter to just update [u,v] by the matrix 1164 times?
– xnor
Jun 27 '20 at 2:22
• @xnor Well, your much more clever implementation is the proof that it can be made significantly shorter (at least 111 bytes shorter). I think I've been worrying too much from the beginning about the poor performance of Node with large integers on that one. It's a bit too late now, so I'll leave my answer as-is. :-) Jun 27 '20 at 10:05
• surely the shorter form would be worth posting another solution? Dec 21 '20 at 7:03
• @ASCII-only Well, the shorter form is really just a port of xnor's answer. Dec 21 '20 at 9:49

# 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.