What algorithm is Newton using in the "De analysi" to extract the square root of a polynomial?

Question

I'm reading a 1745 English translation of Newton's De analysi (apparently the most up-to-date there is, surprisingly). The Latin is here.

In this tract he shows how to use the integral power rule for rational exponents (together with a sum rule: the area under a sum of curves being the sum of the areas) to find quadratures, and he gives examples in three types of how to deal with the case of integrating $y$ when it is not a polynomial: examples using long division e.g. $y=a^2/(b+x)$, examples using the extraction of square roots e.g. $y=\sqrt{a^2+x^2}$, and examples using "the resolution of affected equations" when $y$ is defined implicitly via a polynomial $f(x,y)=0$.

For the second example, $y=\sqrt{a^2+x^2}$, he uses an algorithm to extract the root, but I don't recognize the algorithm. Here's an image:

Note that Newton is not applying his binomial series: he is doing something else.

I realize this is equivalent to the following procedure. Assume $y=\sum b_ix^i$ and we are given $y^2=\sum a_ix^i$; the goal is to solve for the $b_i$. Using the Cauchy product, squaring the first equation gives $y^2=\sum c_ix^i$ where $c_i=\sum_{k=0}^ib_kb_{i-k}$. Pattern matching the coefficients $c_i=a_i$ allows us to recursively solve for the $b_i$, because $c_i$ is the first term to use $b_i$. Thus

$$c_0=b_0^2$$ $$c_1=2b_0b_1$$ $$c_2=2b_0b_2+b_1^2$$ $$c_3=2b_0b_3+2b_1b_2$$ $$c_4=2b_0b_4+2b_1b_3+b_2^2$$ etc.

But Newton is clearly doing something else: drawing on some established algorithm with a visual representation. I don't think it's his algorithm. So:

Three Questions

• How does the algorithm work?
• Does the algorithm have a name?
• Who first used the algorithm?

EDITED

After puzzling it out a little, I at least see the answer to my first question. Here is the procedure applied to $\sqrt{a^2+x^2}$.

Step 1: Guess $a$. Store this as the latest estimate $S_0$.

Step 2: Square, giving $a^2$, then subtract from the radicand: $x^2$. Store this as the latest remainder $R_0$.

Step 3: Update $S_i$ to $S_{i+1}$ by adding a term $y_i$ to $S_i$ so that $2y_iS_0$ agrees with the lowest-degree term of $R_i$. In this case, we want to add $y_0$ so that the lowest-degree term of $2y_0a$ agrees with $x^2$. Hence $y_0=x^2/2a$. The estimate $S_0$ is thus updated to $S_1=a+x^2/2a$.

Step 4: Multiply $y_i$ by $(2S_i+y_i)$, then subtract from the last remainder and store the result as the latest remainder $R_{i+1}$. In this case, we multiply $x^2/2a$ by $(2a+x^2/2a)$, giving $x^2+x^4/4a^2$. Subtracting from the last remainder gives $-x^4/4a^2$.

Step 5: Repeat steps 3/4.

Just to illustrate for the next term, we want to add $y$ to that $2ya=-x^4/4a^2$. Thus we add $y=-x^4/8a^3$ and the estimate is updated to $a+x^2/2a-x^4/8a^3$. Now we multiply

$$-\frac{x^4}{8a^3}\left(2\left(a+\frac{x^2}{2a}\right)-\frac{x^4}{8a^3}\right)=-\frac{x^4}{4a^2}-\frac{x^6}{8a^4}+\frac{x^8}{64a^6}$$

Subtracting from the previous remainder $-x^4/4a^2$ yields the new remainder

$$\frac{x^6}{8a^4}-\frac{x^8}{64a^6}$$

and so on.

Answer 1

Partial answer with historical context

After googling "square root of a polynomial," I found that the algorithm is lucidly described in many standard 19th and early 20th century algebra textbooks for high schools and colleges: e.g. here, here, and here. (How poorly modern textbooks fare by comparison!)

Encouraged by this evidence, I decided to glance at George Chrystal's famous textbook Algebra (first published 1886). He gives a thorough discussion in section 17 ("Square root of an integral function of $x$") in Chapter XI ("Arithmetical theory of surds") in volume 1. He also leaves a historical footnote, dating the algorithm to 1557 in English. I reproduce this passage here for ease of reference:

When an integral function of $x$ is a complete square as regards $x$, its square root can be found by a method analogous to that explained in $\S12$, for finding the square root of a number. Although the method is of little interest, either theoretically or practically, we give a brief sketch of it here, because it illustrates at once the analogy and the fundamental difference between arithmetical and algebraical operations.*

I. We may restate Proposition I. of $\S 12$, understanding now $A$ and $B$ to mean integral functions of $x$.

II. If $F=p_ox^{2n}+p_1x^{2n-1}+\cdots+p_{2n}$, and if $$\sqrt{F}=(q_nx^n+q_1x^{n-1}+\cdots+q_{n-p+1}x^{n-p+1})+(q_{n-p}x^{n-p}+\cdots+q_0)=P+Q$$ say; and if we suppose the first $p$ terms, namely $P=q_nx^n+q_1x^{n-1}+\cdots+q_{n-p+1}x^{n-p+1}$, of $\sqrt{F}$ to be known, then the next $p$ terms will be the first $p$ terms in the integral part of $(F-P^2)/2P$.

For we have $$F=P^2+2PQ+Q^2$$ hence $$\frac{F-P^2}{2P}=Q+\frac{Q^2}{2P}$$ Now the degree of the integral part of $Q^2/2P$ is $2(n-p)-n=n-2p$. Hence $Q^2/2P$ will at most affect the term in $x^{n-2p}$. Hence $(F-P^2)/2P$ will be identical with $Q$ down to the term in $x^{n-2p+1}$ inclusive. In other words, the first $n-p-(n-2p)=p$ terms obtained by dividing $F-P^2$ by $2P$ will be the $p$ terms of the square root which follow $P$.

We may use this rule to obtain the whole of the terms one at a time, the highest being of course found by inspection as the square root of the highest term of the radicand; or we may obtain a certain number in this way, and then obtain the rest by division.$\dagger$

$*$ The method was probably obtained by analogy from the arithmetical process. It was first given by Recorde in The Whetstone of Witte (black letter, 1557), the earliest English work on algebra.

$\dagger$ Just as in division, we may, if we please, arrange the radicand according to ascending powers of $x$. The final result will be the same whichever arrangement be adopted, provided the radicand is a complete square. If this is not the case the operation may be prolonged indefinitely just as in continued division. We leave the learner to discover the meaning of the result obtained in such cases. The full discussion of the matter would require some reference to the theory of infinite series.