+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Calculators (and very old HP Computers) (/forum-3.html)
+--- Forum: General Forum (/forum-4.html)
+--- Thread: Square Root Process Similar to Long Division (/thread-17484.html)
Square Root Process Similar to Long Division - jeejohn - 09-17-2021 11:00 PM
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?).
RE: Square Root Process Similar to Long Division - Albert Chan - 09-17-2021 11:58 PM
(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
RE: Square Root Process Similar to Long Division - Albert Chan - 09-18-2021 01:51 AM
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
RE: Square Root Process Similar to Long Division - Albert Chan - 09-18-2021 12:46 PM
(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
RE: Square Root Process Similar to Long Division - Albert Chan - 09-18-2021 01:07 PM
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
RE: Square Root Process Similar to Long Division - jeejohn - 09-18-2021 08:55 PM
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.
RE: Square Root Process Similar to Long Division - Thomas Klemm - 11-06-2022 04:33 PM
(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) / 2Examples
d = Decimal('.9996')
x = Decimal('1')
newton(d, x, 7)
0.9996
0.99979992
0.99979997999599359936
0.99979997999599899971991597354766255640239090891735007232
0.9997999799959989997199159735914171390272639623863827366463080538247902204151860911590260475681681352
0.9997999799959989997199159735914171390272639623863827366491803235872148569016030640762583498945948583
0.9997999799959989997199159735914171390272639623863827366491803235872148569016030640762583498945948583
d = Decimal('12345')
x = Decimal('0.009')
newton(d, x, 7)
111.105
111.1080553875
111.1080555135405110305346994140625
111.1080555135405112450044387430752408689311675830734154751136313016974218865863978862762451171875
111.1080555135405112450044387430752414899113774596977299764856731617825851115519126275283305922167024
111.1080555135405112450044387430752414899113774596977299764856731617825903175167629180920706522637384
111.1080555135405112450044387430752414899113774596977299764856731617825903175167629180920706522637384
Depending on the use case only a few iterations might be needed.
Still the multiplications might get tedious after a while.
RE: Square Root Process Similar to Long Division - Thomas Okken - 11-06-2022 07:25 PM
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.