#### The Story

Around 2300 years ago, Archimedes developed a "BigInteger" scheme. As the Greeks used Roman numerals, there were limitations to the size of the numbers that were expressable. At the same time, there was another Greek dude who made a competing BigInteger package, and he was going around town telling everyone that his implementation was better than Archimedes'. So Archimedes sent him a letter challenging him to solve a problem in the BigInteger realm. In modern parlance, Archimedes likely sent the letter to troll him.

#### The Letter

If thou art diligent and wise, O stranger, compute the number of cattle of the Sun, who once upon a time grazed on the fields of the Thrinacian isle of Sicily, divided into four herds of different colours, one milk white, another a glossy black, a third yellow and the last dappled. In each herd were bulls, mighty in number according to these proportions: Understand, stranger, that the white bulls were equal to a half and a third of the black together with the whole of the yellow, while the black were equal to the fourth part of the dappled and a fifth, together with, once more, the whole of the yellow. Observe further that the remaining bulls, the dappled, were equal to a sixth part of the white and a seventh, together with all of the yellow. These were the proportions of the cows: The white were precisely equal to the third part and a fourth of the whole herd of the black; while the black were equal to the fourth part once more of the dappled and with it a fifth part, when all, including the bulls, went to pasture together. Now the dappled in four parts were equal in number to a fifth part and a sixth of the yellow herd. Finally the yellow were in number equal to a sixth part and a seventh of the white herd. If thou canst accurately tell, O stranger, the number of cattle of the Sun, giving separately the number of well-fed bulls and again the number of females according to each colour, thou wouldst not be called unskilled or ignorant of numbers, but not yet shalt thou be numbered among the wise.

But come, understand also all these conditions regarding the cattle of the Sun. When the white bulls mingled their number with the black, they stood firm, equal in depth and breadth, and the plains of Thrinacia, stretching far in all ways, were filled with their multitude. Again, when the yellow and the dappled bulls were gathered into one herd they stood in such a manner that their number, beginning from one, grew slowly greater till it completed a triangular figure, there being no bulls of other colours in their midst nor none of them lacking. If thou art able, O stranger, to find out all these things and gather them together in your mind, giving all the relations, thou shalt depart crowned with glory and knowing that thou hast been adjudged perfect in this species of wisdom.

#### Archimedes' Cattle Problem

There are black, white, yellow, and dappled bulls and cows. If we label the bulls by the upper-case letters B, W, Y, D and the cows by the lower-case letters b, w, y, d. The problem is given in two parts.

##### Part 1

The first paragraph reduces to solving a system of 7 equations in 8 variables.

 W = (1/2 + 1/3) B + Y rarr 6W = 5B + 6Y B = (1/4 + 1/5) D + Y rarr 20B = 9D + 20Y D = (1/6 + 1/7) W + Y rarr 42D = 13W + 42Y w = (1/3 + 1/4) (B + b) rarr 12w = 7B + 7b b = (1/4 + 1/5) (D + d) rarr 20b = 9D + 9d d = (1/5 + 1/6) (Y + y) rarr 30d = 11Y + 11y y = (1/6 + 1/7) (W + w) rarr 42y = 13W + 13w
Archimedes, paraphrased: "If you can solve this, you're not totally unskilled in math..."

This is an underdetermined system and there are infinitely many integer solutions. To find the solutions, we can reduce the matrix using the sympy package in python.

import sympy as sp

M = sp.Matrix([
[  6,  -5,   0,  -6,   0,   0,   0,   0],
[  0,  20,  -9, -20,   0,   0,   0,   0],
[-13,   0,  42, -42,   0,   0,   0,   0],
[  0,  -7,   0,   0,  12,  -7,   0,   0],
[  0,   0,  -9,   0,   0,  20,  -9,   0],
[  0,   0,   0, -11,   0,   0,  30, -11],
[-13,   0,   0,   0, -13,   0,   0,  42]
])

S = M.nullspace()
print(S)

The nullspace of M is a one-dimensional space spanned by the vector.

[Matrix([
[3455494 / 1813071],
[ 828946 /  604357],
[7358060 / 5439213],
[ 461043 /  604357],
[2402120 / 1813071],
[ 543694 /  604357],
[1171940 / 1813071],
[      1          ]])]

These rational numbers have a common denominator of 5,439,213, so any integer solution to this system is just this vector scaled by a multiple of 5,439,213, and the total number of cattle is the sum of these numbers.

• W = 10355482 * k
• B = 7460514 * k
• D = 7358060 * k
• Y = 4149387 * k
• w = 7206360 * k
• b = 4893246 * k
• d = 3515820 * k
• y = 5439213 * k
• Total = 50389082 * k

In python:

def solution_vector(k):
return [5439213 * k * x for x in S[0]]

def total_cattle(k):
return sum(solution_vector(k))

The smallest solution has 50,389,082 cattle, but Archimedes is unimpressed. We are, by his accounting, not entirely unskilled in mathematics!

##### Part 2

In the second paragraph, Archimedes imposes two additional conditions. The first condition is that W + B must a square number. Using the computed solution vector, we have

W + B = 17826996 * k
W + B = 2^2 * 3 * 11 * 29 * 4657 * k

For W + B to be a square number, k must be equal to 3 * 11 * 29 * 4657 multiplied by a square number, ie, there must be some integer v such that

k = 4456749 * v^2

The second condition is that D + Y is a triangular number. Using the solution vector and value of k derived above, we have that

D + Y = 11507447 * k
D + Y = 410286423278424 * v^2

Additionally, D + Y must be a triangular number, so it must be of the form n (n + 1) / 2. Using the quadratic formula to solve for n and designating R = 410286423278424, we arrive at the equation

n = (-1 + sqrt(1 + R v^2)) / 2.

The inside of the square root must be a perfect square, say u^2, which means that we are left trying to find an integer solution to

u^2 - R v^2 = 1.

In modern mathematics, an equation of this type is call Pell's equation. It was misattributed to John Pell, who never actually worked on this problem, but the name stuck. Long before being misnamed, these equations were studied by the Ancient Greeks and solved by Indian mathematicians in the first millenium CE.

#### Continued Fractions

The Indian mathematician Brahmagupta (~600 CE) developed an identity for solving some instances of Pell's equation. Bhaskara II (~1150 CE) extended Brahmagupta's method to solve any instance of Pell's equation. But it wasn't until 1768 that Lagrange proved that Bhaskara's method will always terminate after finitely many steps and yield a solution. In Lagrange's work, he tied the Indian methods to continued fractions. This is the viewpoint we will adopt, although the procedure is essentially the same as Bhaskara II.

The continued fraction algorithm begins with a real number alpha, computes its integer part q, and its fractional part beta. This beta is between 0 and 1, so it's reciprocal is greater than 1. The reciprocal of beta becomes alpha in the next iteration and the procedure is continued indefinitely or until beta is equal to 0. The pseudo-code for the continued fraction algorithm is as follows.

beta = alpha - int(alpha)
while beta > 0 do
alpha = 1 / beta
beta = alpha - int(alpha)

This algorithm terminates after finitely many steps if and only if alpha is a rational number. In fact, for rational numbers, the continued fraction algorithm is essentially equivalent to the Euclidean algorithm.

The name continued fraction arises as follows. If we index the values of q as q_0, q_1, q_2, ..., then we write alpha = [q_0, q_1, q_2, ...], which means q_0 + 1 / (q_1 + 1 / (q_2 + 1 / ddots) ).

When alpha = sqrt(R) for a non-square integer R, it turns out that the sequence q_0, q_1, q_2, ... is periodic. For example, the continued fraction algorithm for alpha = sqrt(7) is the following.

alpha q beta
sqrt(7) 2 sqrt(7) - 2
(sqrt(7) + 2) / 3 1 (sqrt(7) - 1) / 3
(sqrt(7) + 1) / 2 1 (sqrt(7) - 1) / 2
(sqrt(7) + 1) / 3 1 (sqrt(7) - 2) / 3
sqrt(7) + 2 4 sqrt(7) - 2

At this point, we see that the final beta is equal to the initial beta, so the pattern will repeat indefinitely. We write sqrt(7) = [2, 1, 1, 1, 4, 1, 1, 1, 4, 1, 1, 1, 4, ...] which we abbreviate by sqrt(7) = [2; 1, 1, 1, 4]. The semicolon separates the initial part from the periodic part.

To compute the continued fraction in python, we represent the number (a sqrt(R) + b) / c by the tuple (a, b, c). Calculating m as the integer part of sqrt(R), we start with sqrt(R) = (1, 0, 1) and continue until alpha = (1, m, 1). This is enough, because the first beta is (1, -m, 1), which will also be the beta after alpha = (1, m, 1).

We also need the greatest common divisor to simplify the fraction. The rest of the algorithm is just keeping track of the arithmetic of rationalizing the denominator after taking the reciprocal.

import math

def gcd(a, b):
if b == 0:
return abs(a)
else:
return gcd(b, a % b)

def continued_fraction(R):
m = int(math.sqrt(R))
cf = []
alpha = (1, 0, 1)
while alpha != (1, m, 1):
a, b, c = alpha
q = int((a*m + b) / c)
x = q*c - b
d = R * a**2 - x**2
D = gcd(c, d)
alpha = (a*c/D, x*c/D, d/D)
cf.append(q)
return cf

We compute the continued fraction for sqrt(R) = sqrt(410286423278424) and find that its period to equal 203253. The first 10 terms are [20255528; 4, 1, 1, 1, 4, 1, 2, 3, 4, ...]. In python:

R = 410286423278424
cf = continued_fraction(R)

print('first ten terms: {}'.format(cf[:10]))
print('period: {}'.format(len(cf) -1))

#### Convergents

The convergents of a continued fraction are the rational numbers obtained by stopping after finitely many q's. For example, the first five convergents of sqrt(7) are:

• 1st convergent = 2
• 2nd convergent = 3 = 2 + 1/1
• 3rd convergent = 3/2 = 2 + 1 / (1 + 1/1)
• 4th convergent = 8/3 = 2 + 1 / (1 + 1 / (1 + 1/1))
• 5th convergent = 37/14 = 2 + 1 / (1 + 1 / (1 + 1 / (1 + 1/4)))

Two important facts:

1. the convergents form a sequence of increasingly accurate rational approximations to the real number in question.
2. the convergent right before the end of the first period of the continued fration of sqrt(R) gives a solution to x^2 - R y^2 = pm 1.

In our example for sqrt(7), if we write each convergent as the fraction a/b, we can plug in x = a and y = b to x^2 - 7 y^2:

• 1st convergent = 2/1 gives 2^2 - 7(1^2) = -3
• 2nd convergent = 3/1 gives 3^2 - 7(1^2) = 2
• 3rd convergent = 3/2 gives 3^2 - 7(2^2) = -19
• 4th convergent = 8/3 gives 8^2 - 7(3^2) = 1

In short, the continued fraction convergents of sqrt(R) can be used to find a solution to Pell's equation.

To calculate the convergents efficiently, we use another property of sqrt(R), namely that the numerators and denominators of a convergent both satisfy the following recurrence relation:

x_(n+1) = q_n * x_n + x_(n-1)
def convergent(cf, n):
a0, a1 = 1, cf[0]
b0, b1 = 0, 1
for i in range(1, n):
a0, a1 = a1, cf[i] * a1 + a0
b0, b1 = b1, cf[i] * b1 + b0
return a1, b1

To find a solution to our Pell's equation, we compute the 2nd to last convergent, and check its value in x^2 - Ry^2. By the way, for the solution found, x has 103273 digits and y has 103266 digits.

a, b = convergent(cf, len(cf))

assert( a**2 - R * b**2 == 1 )

#### The Solution

We have found a solution u = a and v = b to u^2 - R v^2 = 1, where a, b are the extremely large numbers just calculated. We plug this b into the equation k = 4456749 * v^2 to get k. Finally, we plug this k into 50389082 * k to get the total number of cattle. The total number of digits is in the solution 206,545.

#### History

• 1880: Amthor was the first to find the solution. He used a few techniques from algebraic number theory to find an exact formula for the solution. The computation of the actual number was impossible for him to perform by hand. However, he used logarithms to calculate the correct number of digits and calculated, and he attempted to compute the first 4 digits. He computed 7766 as the first four digits, the first three of which were correct.
• 1965: Mathematicians at the University of Waterloo used an IBM 7040 to compute the full solution. It took 7.5 hours.
• 1981: A Cray-1 computer was used to compute the solution. It took 10 minutes.
• This year: you can run the code contained here to find the full solution. On my laptop, it took 3.74 seconds.
%%time
R = 410286423278424
cf = continued_fraction(R)
a, b = convergent(cf, len(cf))
k = 4456749 * b**2
total = total_cattle(k)