A Story of Computational Science: Colonel Titus’ Problem from the 17th Century

: Experimentation and the evaluation of algorithms have a long history in algebra. In this paper we follow a single test example over more than 250 years. In 1685, John Wallis published A treatise of algebra, both historical and practical , containing a solution of Colonel Titus’ problem that was proposed to him around 1650. The Colonel Titus problem consists of three algebraic quadratic equations in three unknowns, which Wallis transformed into the problem of ﬁnding the roots of a fourth-order (quartic) polynomial. When Joseph Raphson published his method in 1690, he demonstrated the method on 32 algebraic equations and one of the examples was this quartic equation. Edmund Halley later used the same polynomial as an example for his new methods in 1694. Although Wallis used the method of Vietè, which is a digit–by–digit method, the more efﬁcient methods of Halley and Raphson are clearly demonstrated in the works by Raphson and Halley. For more than 250 years the quartic equation has been used as an example in a wide range of solution methods for nonlinear equations. This paper provides an overview of the Colonel Titus problem and the equation ﬁrst derived by Wallis. The quartic equation has four positive roots and the equation has been found to be very useful for analyzing the number of roots and ﬁnding intervals for the individual roots, in the Cardan–Ferrari direct approach for solving quartic equations, and in Sturm’s method of determining the number of real roots of an algebraic equation. The quartic equation, together with two other algebraic equations, have likely been the ﬁrst set of test examples used to compare different iteration methods of solving algebraic equations.


Introduction
A problem brought to John Pell (1611-1685) in 1649, and discussed at the time with Silius Titus (1623-1704), was the following-to find numbers a, b, and c satisfying the equations a 2 + bc = 16, b 2 + ac = 17, and c 2 + ab = 22.
A solution with positive integers is easily seen to be a = 2, b = 3, and c = 4, but Pell decided to challenge himself by changing the final equation: a 2 + bc = 16, b 2 + ac = 17, and c 2 + ab = 18.
leading to a quartic equation, an equation that is not derived from a problem in geometry or trigonometry. A variant of the Colonel Titus problem is Question 113 in Ladies' Diary from 1725 shown in Figure 1 a 2 + bc = 920, b 2 + ac = 980, and c 2 + ab = 1000, (3) and was solved by John Turner in 1726. Turner only specifies the quartic equation to be solved and a solution a, b, c, of (3). Question 113 is also found in algebra textbooks in 1820 ( [4], p. 405) and in 1840 ( [5], p. 563).
A fourth variant of the problem was published in Question 209 in The Scientific Receptacle in 1796 and shown in Figure 2: a 2 + bc = 1 806 520, b 2 + ac = 2 225 275, and c 2 + ab = 5 567 720. (4) John Ryley solved the problem and introduced a new way to solve it by expressing a and b as a fraction of c [9]. Wallis ([2], pp. 225-256) eliminates the variables b and c in (2) and reduces the three equations to a fourth-order algebraic equation where x = 2a 2 . In the following we will use the term "Pell-Wallis equation" to refer to (5).
To determine a root x * , Wallis uses Viète's method and a is found through a = √ x * /2. To compute b, Wallis derives a cubic equation which follows from multiplying the first quadratic equation by a and the second by b and eliminates abc to obtain the cubic equation Having found a and b, the unknown c is found from the first quadratic a 2 + bc = 16.
One of the most classical problems in mathematics is the solution of systems of polynomial equations in several unknowns [10]. They arise in robotics, coding theory, optimization, mathematical biology, computer vision, game theory, statistics, machine learning, control theory, and numerous other areas [10]. Systems of quadratic polynomial equations appear in nearly every crypto-system [11] and in robotics [12].
For more than 250 years, the equation x 4 − 80x 3 + 1998x 2 − 14937x + 5000 = 0 has played an important role in the development of new methods, analyses of algebraic equations, and comparisons of methods for solving nonlinear equations.
In Section 2 we discuss four different approaches in solving Colones Titus' problem that have appeared in the literature and in Section 3 we discuss different techniques and methods using the Pell-Wallis Equation (5). For a modern treatment of numerical methods for roots of polynomials, see [13,14] and references therein. For solving systems of polynomial equations, see [10,11] and references therein.

Colonel Titus' Problem
Using the notation in Wallis algebra book from 1685 [2], Ch. LX-LXI, the general Colonel Titus problem is as follows. For given positive real numbers l, m, and n, find a, b, and c such that We review several solution techniques, mainly using what could be described as high-school algebra [15]. An elegant solution is given in Solutions of the principal questions of Dr. Hutton's course of mathematics by Thomas Stephens Davies, and we follow his solution technique.
From (6) and (7) we have Equating the two expressions for c, we have a cubic equation in b b 3 − mb + a(l − a 2 ) = 0.
From (8) and the two expression for c above, we have Multiply the quadratic equation by b and the cubic equation by l − 2a 2 and subtract the two expressions to eliminate the cubic term. We now have two quadratic equations in b To eliminate b 2 , multiply the first quadratic equation by n and the second by l − 2a 2 and subtract the two resulting quadratic equations. The result is a linear equation in b. Solve Substitute the value for b in Equate the two expressions in nb − ma = (l − a 2 )(l − 2a 2 )/b and simplify Multiply the equation by 2 and let x = 2a 2 and we have the equation For each real root, x * of (9) a, b, and c, can be computed in the following way; a in 2a 2 , and c in a 2 + bc = l. For l = 16, m = 17, and n = 18 we have the equation Different techniques for solving Colonel Titus' problem have been suggested in the literature by philomaths and mathematicians, school teachers of mathematics, and professors of mathematics. The different solution techniques can mainly be divided in two groups; the first group is based on elimination and the second group on first reformulating the problem and then performing an elimination.
The first solution to Colonel Titus' problem was published by J. Wallis [2] and this was an elimination of the unknowns that results in the quartic Equation (5). To find the four positive roots of (5) Wallis used a digit-by-digit computation method. The solution of Colonel Titus' problem by Wallis was republished by Francis Maseres (1731-1824) in 1800, including numerous details ( [16], pp. 187-239). However, Maseres did not use a digit-by-digit method to find the roots, but rather the Newton-Raphson method. Similar solutions using explicit elimination are found in [5,[17][18][19], all leading to the same quartic Equation (5). J. Kirkby [20] in 1735 and A. Cayley [21] in 1860 used a general elimination theory, leading to the same quartic equation.
The method of introducing two new variables expressing the unknowns as a fraction of one of the other variable was studied by J. Ryley [9] in 1796, and made popular by William Frend [22] in 1800. Variations of this technique are found in [23][24][25]. Ivory expressed two of the unknowns as a difference of the third [26,27]. All these reformulations lead to quartic equations that are different from the quartic Equation (5). These quartic equations never reached the same popularity as (5).
Using iterative methods to solve the three equations simultaneously was suggested in the Diarian Repository [6] in 1774 and by Whitley [28] in 1824.

Ladies' Diary 1725 Question 113
We find a variation of Colonel Titus' problem in the journal Ladies' Diary from 1725 in Question 113 shown in Figure 1 where l = 920, m = 980, and m = 1000.
In Ladies' Diary in 1726 John Turner (active in Ladies' Diary from 1726 to 1750 ( [17], p. 423)) gives a solution of the problem and states the equation (using the notation in Wallis) Let x = 2u 2 and multiplying the equation by two gives (9). Turner gives the solution of Question 113 in Ladies' Diary to be 19.5991, 22.7788, and 23.5276. There are three minor typographical errors in the solution by Turner ( [29], p. 7): These three typographical errors are repeated in Diarian Miscellany [7] and Diarian Repository [6] and one error is pointed out in the Errata of [8].
For l = 920, m = 980, and n = 1000, the Equation (9) has four positive roots approximately equal to 1937.6, 1881.6, 768.0, and 12.7, and the only root that gives reasonable ages is 768.0, and the ages are approximately a = 19.5965, b = 22.7799, and c = 23.5286.
Leybourn also presents an additional solution of Colones Titus' problem provided by Mark Noble (a mathematician at Royal Military College (Sandhurst)) in the appendix in the fourth volume [17], pp. 255-259. The contribution is signed N and in the preface of Leybourn's first volume ( [8], Preface page X) it is signed "this is Mark Noble".This is an elimination technique and it leads to the same quartic equation as in Wallis The problem is find positive numbers (using the notation in Wallis) a, b, and c so that a 2 + bc = 1, 806, 520, b 2 + ac = 2, 225, 275, and c 2 + ab = 5, 567, 720 with a solution published in a later issue in the same volume ( [9], p. 95).
and derived the equation From a root of (11), all other quantities can be determined. However, Ryley does not compute any root of (11) or values of a, b, and c for the given l, m, and n.

First Reformulation and then Elimination
J. Ryley was the first to express two of the unknowns as a fraction of the third. W. Frend (1757-1841) ( [22], pp. 240-246) in 1800 provided a different derivation and introduced x and y so that b = x a and c = y a, and derived the equation A minor improvement of Frend's solution, avoiding a square root, was given by John Hellins (c. 1749-1827) in the introduction of the same volume in which Frend's solution was found [16], pp. lxxi-lxxii. By interchanging the variables, Equation (13) can be derived from (11). For l = 16, m = 17, and n = 18 we obtain the equation The equation has four real solutions (or roots), of which one is negative. Maseres [16] (pp. 246-275) finds the three positive roots to be approximately 1.027179787, 1.17565, and 9.3388519 using Newton-Raphson iteration. Maseres regards the root 1.027179787 as "impossible" since y is negative. For a given root, y and the unknowns a, b, and c are easily found. Maseres ends the tract with a comment that Mr. Frend's solution has the advantage that it saves the trouble of those very tedious and perplexing algebraic multiplications and divisions necessary in Dr. Wallis's solution [16], p. 275. A similar solution to Frend's was given by Tebay [25] in 1845. A third variation is to express a and b in terms of c [23]. James Ivory (1765-1842) [26] wrote that the solution provided by Wallis to the problem (2) was remarkably operose and inellegant and a solution of the same problem by Frend [22] is preferable to Wallis's solution. Ivory expressed two of the unknowns as a difference of the third and the analysis was printed in 1804 in [26], but with no numerical solution. Ivory restricts his analysis to the specific choice l = 16, m = 17, and n = 18. Ivory's analysis was mailed to Baron Maseres [27] p. 360 and Maseres added many details and a numerical solution based on the Newton-Raphson method [32].
Whitley [24], p. 121 wrote in 1824 that Mr. Ivory's solution was an elegant specimen of analysis and Davis [18], p. 274 in 1840 called it an exceedingly elegant investigation. Cockle [3] speculated that the analysis can be carried over to the case where m = (n + l)/2. It can be shown that the derivation by Ivory can be extended to the general case of l, m, and n. Maseres [32], pp. 371-395 computed the two positive roots of the quartic equation derived by Ivory and these correspond to the positive a, b, and c values provided by Wallis.

Simultaneous Solution of the Three Unknowns
In the Repository Solution section in the Diarian Repository [6], pp. 190-191 an iterative approach was suggested. First, find an approximate solution (in this case 23, 22.5, and 21.1); then, find a correction (x, y, and z) that solves the (linear) equations where the second order terms are eliminated. . . . then via the solution of the resulting equations, x, y, and z will be determined to a sufficient degree of exactness; if not, the operation must be again repeated with the last found values. . . [6], pp. 190-191. This is Newton's method but no actual computations of a, b, and c are shown, except for finding the starting point for the iteration.
J. H. Swales, the editor of the Liverpool Apollonius asked its readers in 1823 to find a simpler solution than those given by Ivory [26], p. 156 and Frend [22], p. 240. Three The method proposed by Whitley [28], pp. 127-128 in 1824 is the fixed-point iteration with the starting point given by a 0 = b 0 = c 0 = 3, where l = 16, m = 17, and n = 18. Table 1 which are three homogeneous equations of second order in three unknowns [21]. However, Cayley did not solve the homogeneous equations. Schumacher solved this problem [33] in 1911.

Erroneous Solution
The achievements of Adrien Quentin Buée (1748-1825), also called Abbé Buée, are important in relation to the conceptual development of the negative numbers and for the graphical representation of the complex numbers. In [34], he considers Colonel Titus' problem and makes an attempt to solve it using geometry and complex numbers. He claims that the solution must be a = 3.25x, b = 4.25x, and c = 5.25x, where x is the area of a circle in the geometric construction. However, he does not find any correct solution to the problem.

The Pell-Wallis Equation
In the late 17th and early 18th centuries, there were numerous collections of algebraic equations [35]. Most practical algebraic equations were derived from geometric or trigonometric problems. An algebra book by John Ward from 1695 contains ten geometric problems with corresponding algebraic equations [36] and The Young Mathematician's guide from 1707 contains more than 20 practical problems from geometry and trigonometry, leading to algebraic equations [37]. However, The Pell-Wallis equation is derived from a different type of problem. The equation has been in use for 270 years, from the first time it appeared in print in 1685 to the most recent reference to the equation in a paper from 1955.

Digit-by-Digit Methods
The root finding method used by Wallis in 1685 was a digit-by-digit computation method [2]. The method used by Wallis was based on Vietè's method but it deviated from Vietè's method in the divisor used to compute the next digit [38]. In this method, the roots are computed with a very high degree of accuracy. With Horner's technique to compute shifted polynomials, the digit-by-digit approach became more efficient using Holdred's and Horner's divisor [38]. The Pell-Wallis equation is used as an example in Holdred [39], pp. 55-56 and Nicholson [40]

Bracketing Methods
In Vietè's method, the first digit of a root must be specified. This will normally lead to the determination of the intervals of the roots. Intervals of the real roots may also provide a starting point for linear interpolation. Cardano's golden rule and regula falsi are methods in which a root is bracketed. Application of the Newton-Raphson method and the Halley method, which are iterative methods, requires a starting point sufficiently close to a root/solution and this point is often determined to be in an interval including the root.

Linear Interpolation
The first use of the Pell-Wallis equation and interpolation occurred in 1732. Graaf [55], pp. 33-35 considered (5) and scaled the variable x ← x/10 in the interval of 0 to 3.6 and plotted the graph (x, f (x)), where f (x) is the left-hand side of (5). Based on the graph, an interval where a solution exists was identified, and then linear interpolation. This is a variation of regula falsi [56] and Cardano's regula aurea [57], Chapter 30 methods, since both end points of the interval are changed in de Graaf's approach.
The method of John Davidson, a teacher in mathematics in Burntisland, involves a bracketing approach and linear interpolation [58], p. 114, [59], p. 38, as shown in his textbooks from 1814 and 1852. This is Cardano's golden rule [57], Chapter 30.

The Newton-Raphson Method
Wallis published his algebra book in 1685 [2] and it contained the first printed version of Newton's method. When Raphson presented his method in 1690 it was regarded as a different method. It was not until the mid-18th century that it became clear that the two methods generated the same sequence of iterations [35]. In Volume III of Scriptores logarithmici from 1796, Francis Maseres used the Newton-Raphson method. First an approximation 0.3507 to the smallest root is found by using a series expansion and then two iterations are performed [61], pp. 718-725. Maseres writes ". . . and this I take to be the very best method that can be employed to find the value of x to this degree of exactness".
Lockhart [62] in 1839 argues that the numbers of digits required to compute an approximate solution using the Newton-Raphson method is not worse than Horner's digit-by-digit method, as presented by De Morgan [42] in 1839.

Ferrari-Cardano Approach
The linear shift x − 20 in the Pell-Wallis equation makes the term x 3 vanish and the depressed quartic equation is Taking two slightly different approaches, Francis Maseres first finds the depressed quartic (14) and then, with reference to Ferrari, finds the resolvent cubic and, with reference to Descartes [68], p. 142, the resolvent cubic (in e 2 ) is The four roots of the Pell-Wallis equation can then found [68], pp. 134-182. Maseres points out that the use of linear interpolation and one iteration with Newton-Raphson will require fewer arithmetic operations than the use of Ferrari-Cardano approach [68] p. 178. William Rutherford (1798Rutherford ( -1871 found the depressed quartic (14) and then derived the resolvent cubic equation (in u 2 ) From the resolvent cubic (17), using Horner's method, Rutherford found one root and the four roots of the Pell-Wallis equation [69], pp. 17-18. Orson Pratt (1811-1881) [70], pp. 130-131 used the depressed quartic (14) and derived the resolvent depressed cubic y 3 − 1,401,372y − 633,074,427 = 0.
A root of the depressed cubic is found using a digit-by-digit approach with a modified divisor in Vietè's method, to eleven decimal places, and then the roots of the Pell-Wallis equation are given.
Christian Heinrich Schnuse (1800-1878) considered the Pell-Wallis equation and derived the depressed quartic (14) and the resolvent cubic (17). Using a digit-by-digit approach, he found the same root of (17) as Rutherford in 1849 [46], pp. 358-359.  [71], pp. 42-44. In Gräffe's method a sequence of polynomials is generated and the method is a "root-squaring" process and approximations to the roots can be computed from the coefficients of of the generated polynomials. The method works well for the Pell-Wallis equation since the roots are real, positive, and separated. The method is suitable for computation by hand, whereas computer implementations usually exhibit overflow after only a few steps. After a few steps, the estimates of the roots are good and suitable for a correction by means of Newton-Raphson iterations. The two smallest roots are correct with four decimal digits after four steps in Gräffe's method. Given that the two smallest roots have been accurately computed, the remaining two roots can be computed [72], pp. 74-75. Encke in 1839, Merino in 1879, and Rey Pastor [72] (1888-1962) in 1924 found it more convenient to work with the log of the coefficients of the polynomials.  (14) and scaled the variable √ 402x and obtained the equation

Miscellaneous Methods and Comments
By graphical inspection, the roots of (18) are located in intervals of length 0.01.For a point in the interval, a first correction method is a Newton-Raphson iteration, then a correction based on the next term in the Taylor expansion. The four roots are computed using two or three corrections. • Silvestre François Lacroix (1765-1843) [78], p. 261 discussed the Pell-Wallis equation as a problem of scaling the coefficients and found that the two roots are between 0 and 10 and 10 and 20.

An Early Comparison of Four Algorithms on Three Examples
One of the first comparisons of the use of several algorithms on different problems is found in [16]. The methods used were Newton-Raphson, Halley's two methods, and regula falsi or linear interpolation. The latter method is called the the differential method in [16] or the method of double position. Maseres [16] p. 109 provides a reference to A Course of Mathematics in Two Volumes, Composed for the Use of the Royal Military Academy by Charles Hutton for the equivalence between the differential method and the method of double position.

Concluding Remarks
We have shown that the three quadratic equations in three unknowns forming Colonel Titus' problem can be reduced to a single quartic equation using standard high school algebra. The different derivations of a quartic equation have been suggested by philomaths and mathematicians, school teachers of mathematics, and professors of mathematics over a period ranging from the mid-17th to the early 20th century. We find systems of quadratic equations in modern crypto-systems or robotics. Today, solutions can easily be obtained through the use of computer algebra systems implemented in Maple, Mathematica, or Wolfram. The modern theory related to solution methods, such as the use of a Gröbner basis, has not yet been explored in relation to Colonel Titus' problem.
We have seen that the quartic equation, the Pell-Wallis Equation (5), derived from Colonel Titus' problem, has been used for more than 250 years as a test example to develop methods to solve algebraic equations, techniques to determine the number of roots, or intervals of the roots, as well as in numerous textbooks. As a well-known equation, it has been included in the early numerical comparisons of root finding methods.
The references in this paper do not form a complete list of the use of this equation and Colonel Titus' problem.