Square Root Process Similar to Long Division

[jeejohn]

I saw the thread "Third Order Convergence for Square Roots..."

I thought that I would offer up a process that allows calculations of square roots, similar to how long division is done:

Using Albert Chan's 12345 to calculate the square root:
100 + 10 + 1 + 0.1 + 0.00 + 0.008 + .0000 + .00005 = 111.10805
100 _________| 12345
(200) (*100) = -10000
200 __________2345
210 _________-2100
220 ___________245
221 __________-221
222 ____________24
222.1 _________-22.21
222.2 ___________1.79
222.20 __________0.0000 (222.20*.01=2.222 > 1.779)
222.20 __________1.79
222.208 ________-1.777664
222.216 _________0.012336
222.2160 ________0.000000 (222.216*.0001=.00222216 > .001336
222.2160 ________0.012336
222.21605 _______0.001111080525
222.21610 _______0.0012251975

etc.
Basically, we're keeping track of 2*X on the left, then we add a digit δ to get 2X+δ. Finally we multiple by δ to get 2Xδ+δ^2. We subtract this from the previous quotient to get a new quotient.

Note that for the final quotient = 12345-111.10805^2 = .0012251975 is the remainder left over from squaring 111.10805.

I used the process to find square roots of numbers manually to 25+ digits fairly quickly. I have thought of somehow programming this on the HP-48s, but I have not been able to get past the doubling of the one digit at each step in some quick ML way.

I'm sorry about how this looks, but this forum appears to take out all of the extra spaces I had placed to make it more readable. I used underlines to space things out some, but it doesn't quite look right.

Except for the first digit, the rest of the digits added are in bold (as δ), then multiplied by δ. In the next line the digit added is doubled (or add δ again?).

[Albert Chan]

(09-17-2021 11:00 PM)jeejohn Wrote:  I saw the thread "Third Order Convergence for Square Roots..."

I thought that I would offer up a process that allows calculations of square roots, similar to how long division is done

Nice. By the way, you can cut/paste table inside code block (like mine)

I just deleted a post about old-school square root from Third Order Convergence thread.
It really should belong here. Here goes.

For sqrt of value very close to 1, it is is even closer to 1. √(1+2ε) ≈ 1+ε
We can even *assume* result is 1, then adjust for it.

Example, √(0.9996)

Here is old school method for square root.
With nice divisor = 1, we do 3-digits at a time.
Numbers are scaled so that we work mostly with integer.

0.9996E6 = 1E3^2 - 400       // X = 1E3

Newtons correction (S-X*X)/2/X, for 2nd 3-digits.

-400E3/2 = -200E3 = -200X + 0

Negative digits OK, we will fix later. It just meant √(0.9996) ≈ 0.999800

Other digits followed below patterns:
1st column = quotient
2nd column = remainder, then corrected (shell-like cross multiply, except center)

Code:

1000
-200    0E3 - (-200)^2/2 = -20E3 = -20X + 0
-20     0E3 - (-20)*(-200) = -4E3 = -4X + 0
-4      0E3 - (-4)*(-200)-(-20)^2/2 = -1E3 = -X + 0
-1      0E3 - (-1)*(-200)-(-4)*(-20) = -280 = -X + 720
-1      720E3 - (-1)*(-200)-(-1)*(-20)-(-4)^2/2 = 719772 = 719X + 772
719     772E3 - 719*(-200)-(-1)*(-20)-(-1)*(-4) = 915776 = 915X + 776
915     776E3 - 915*(-200)-719*(-20)-(-1)*(-4)-(-1)^2/2 = 973375.5 = 973X + 375.5
973     375.5E3 - 973*(-200)-915*(-20)-719*(-4)-(-1)*(-1) = 591275 = 591X + 275
591     275E3 - ...

Normalize base-1000 digits:

√(.9996) = 0.999 799 979 995 998 999 719 915 973 591 ...

Ref: book "Dead Reckoning: Calculating without instruments", chapter on roots

[Albert Chan]

If divisor is not 1, it work the same way

12345 = 111^2 + 24       // X = 111

111       24E3/2 = 12E3 = 108X + 12
108       12E3 - 108^2/2 = 6168 = 55X + 63
055       63E3 - 55*108 = 57060 = 514X + 6

√(12345) ≈ 111. 108 055 514

[Albert Chan]

(09-17-2021 11:00 PM)jeejohn Wrote:  Basically, we're keeping track of 2*X on the left, then we add a digit δ to get 2X+δ. Finally we multiple by δ to get 2Xδ+δ^2. We subtract this from the previous quotient to get a new quotient.

Above is based on identity: (X+δ)^2 = X^2 + δ*(2X+δ)

N = 12345 = (X+δ)^2
X = 100
R = N-X^2 = 2345

R = δ*(2X+δ), but we cannot get δ without square root.

Newton's method assumed (2X+δ) ≈ 2X, then solve for δ
With this assumption, solved δ always over-estmate.
To compensate, we floored estimated result.

R/2/X = 11.725
δ = 11
R -= δ*(2X+δ) = 24

X = X+δ = 111
R/2/X ≈ 0.108108108
δ = 0.1081
R -= δ*(2X+δ) = -0.00988561

X = X+δ = 111.1081
R/2/X ≈ -4.448645058E-5
δ = -4.448646E-5             // floor of negative number, same sig. digits as X

√12345 ≈ 111.1081 + (-4.448646E-5) = 111.10805551354

[Albert Chan]

Newton's method converge fast, but required expensive mul/div.
HP has a digit-by-digit method for square-root, that is very efficient.

First, we halved R

R/2 = δ*(X+δ/2)

if δ=1, δ*(X+δ/2) = (X+0.5)
if δ=2, δ*(X+δ/2) = 2*(X+1) = (X+0.5) + (X+1.5)
if δ=3, δ*(X+δ/2) = 3*(X+1.5) = (X+0.5) + (X+1.5) + (X+2.5)
...

This way, we avoided division to get δ. When R/2 - δ*(X+δ/2) < 0, we gone too far.

√12345, digit-by-digit:

1.2345 - 1^2 = .2345
R/2 = .11725

11.725 - 10.5 = 1.225
1.225 - 11.5 = -10.275

122.5 - 110.5 = 12
12 - 111.5 = -99.5

1200 - 1110.5 = 89.5
89.5 - 1111.5 = -1022

8950 - 11110.5 = -2160.5

895000 - 111100.5 = 783899.5
783899.5 - 111101.5 = 672798
672798 - 111102.5 = 561695.5
561695.5 - 111103.5 = 450592
450592 - 111104.5 = 339487.5
339487.5 - 111105.5 = 228382
228382 - 111106.5 = 117275.5
117275.5 - 111107.5 = 6168
6168 - 111108.5 = -104940.5

616800 - 1111080.5 = -494280.5

61680000 - 11110800.5 = 50569199.5
50569199.5 - 11110801.5 = 39458398
39458398 - 11110802.5 = 28347595.5
28347595.5 - 11110803.5 = 17236792
17236792 - 11110804.5 = 6125987.5
6125987.5 - 11110805.5 = -4984818

√12345 ≈ 111.10805 (under-estimated due to positive remainder)

Except for initial setup, algorithm only use subtraction

Quote:Note that for the final quotient = 12345-111.10805^2 = .0012251975 is the remainder

We confirm by undo scaling and halving R
first digit 1.2345/12345 = 1e-4, other 7 digits 100^7 = 1e14

12345 - 111.10805^2 = 6125987.5 * 1e-10 * 2 = 0.0012251975

[jeejohn]

I had never taken the time to "reverse engineer" the HP square root algorithm. I knew it had to be something similar to what I posted as it calculated square roots only a little bit slower than doing divisions.

I take it that dividing R by 2 in DEC is still fast? In HEX it would just be a bit shift. I see that this allows each subtraction step to be incremented by 1. Nice.

Back in the days of the HP-48 newsgroup, I had posted a challenge over redoing the HP-48 ML multiplication in a different possibly faster way. It used one of the multipliers - say X and calculated 4X. Then using a bunch of GOTOs,
if digit=1 -> +X
digit=2 -> +X+X
digit=3 -> +4X-X (instead of +X+X+X)
digit=4 -> +4X (instead of +X+X+X+X)
etc. It was slightly faster, but not really worth the effort. All the tests and GOTOs steps added a lot of clock cycles.

For division, something similar can be done and also an idea was to add a second loop to follow a negative remainder to "Add Up" to a positive remainder.

Something similar could be done for the square root algorithm. When I did my manual square root calculations, I would keep a table of 2Xδ updated to grab quickly to add at the bottom and subtract to get the next remainder.

Interesting how the processor capabilities affects which algorithm is used to calculate irrational functions. HP finding out about the Cordic algorithm galvanized HP into creating the HP-9100 (the first desktop scientific calculator), which in turn, encourage Bill Hewlett to challenged his workers to create the HP-35.

[Thomas Klemm]

(09-18-2021 01:07 PM)Albert Chan Wrote:  Newton's method converge fast, but required expensive mul/div.

We can do as Newton did and use this equation instead:

\(
\begin{align}
\frac{1}{x^2} = d
\end{align}
\)

Using Newton's method we get the following recurrence:

\(
\begin{align}
x \leftarrow x + \frac{x(1 - d \cdot x^2)}{2}
\end{align}
\)

So we can get rid of the expensive division and only have to divide by 2.
To get \(\sqrt{d}\) we have to multiply the final result \(x\) by \(d\).

Program

This is the corresponding Python program:

Code:

def newton(d, x, n):
     for _ in range(n):
         u = d * x
         print(u)
         x += x * (1 - u * x) / 2

Examples

d = Decimal('.9996')
x = Decimal('1')
newton(d, x, 7)

0.9996
0.99979992
0.99979997999599359936
0.99979997999599899971991597354766255640239090891735007232
0.999799979995998999719915973591417139027263962386382736646308053824790220415186​0911590260475681681352
0.999799979995998999719915973591417139027263962386382736649180323587214856901603​0640762583498945948583
0.999799979995998999719915973591417139027263962386382736649180323587214856901603​0640762583498945948583

d = Decimal('12345')
x = Decimal('0.009')
newton(d, x, 7)

111.105
111.1080553875
111.1080555135405110305346994140625
111.1080555135405112450044387430752408689311675830734154751136313016974218865863​978862762451171875
111.1080555135405112450044387430752414899113774596977299764856731617825851115519​126275283305922167024
111.1080555135405112450044387430752414899113774596977299764856731617825903175167​629180920706522637384
111.1080555135405112450044387430752414899113774596977299764856731617825903175167​629180920706522637384

Depending on the use case only a few iterations might be needed.
Still the multiplications might get tedious after a while.

[Thomas Okken]

I learned the long-division-like square root algorithm from a book called Mathematics for Statistics, by W. L. Bashaw. The book was in English, a language I had only just started learning at the time, but the algorithm was explained with a worked-out example that made it easy to follow. I was a bit surprised that no one else seemed to know this algorithm and it wasn't taught in school, even though calculators were not standard equipment yet at the time.

Regarding the numerical approach, I knew the r -> (r + x/r)/2 iteration from applying Newton's method to r^2-x=0. It does require one division per step, but since it converges pretty quickly, that didn't seem like a big issue to me. I always assumed that the four-banger calculators back then used this algorithm, which would make sense because it's easy to implement, and it would be consistent with √ typically being slower than the arithmetic functions, but not very slow.