(12C) Square Root - Printable Version +- HP Forums (https://www.hpmuseum.org/forum) +-- Forum: HP Software Libraries (/forum-10.html) +--- Forum: General Software Library (/forum-13.html) +--- Thread: (12C) Square Root (/thread-13750.html) (12C) Square Root - Gamo - 10-02-2019 10:23 AM For case study purpose here is the very interesting algorithm to compute the square root of a number by using this equation. B = [A+(n÷A)] ÷ 2 This equation used Successive Approximation Algorithms that continue to better approximate some desired result by providing the initial guess then this program will converge on to the correct answer if your initial guess is bad then more iteration will be required to achieve a giving accuracy. With this program it doesn't matter where you start your initial guess because this algorithm also used the Change in the Approximation: ABS [(B-A) ÷ A] and Tolerance of 0.00000001 to guarantee that the result will be possible. This equation solves for B, which is the new improved approximation for the square root of n until we are happy with this approximation, we continue repeating the process by coping the new B value into A and then re-executing the same equation to get a new B value. ----------------------------------- Procedure: FIX 8 n [ENTER] A [R/S] display Answer [X<>Y] number of iterations ----------------------------------- Example: √10 used 3 as a guess 10 [ENTER] 3 [R/S] 3.16227766 [X<>Y] 3.00000000 Answer: √10 = 3.16227766 and program took 3 iterations to get this answer. ------------------------------------ Program: Code:  01 STO 2  // A 02 R↓ 03 STO 1  // n 04  0 05 STO 0  // Counter  06 RCL 1   07 RCL 2   08  ÷    09 LSTx   10  + 11  2 12  ÷  13 STO 3  // Store B 14 RCL 2 15  - 16 LSTx 17  ÷ 18  0 19 X<>Y 20 X≤Y 21 CHS 22 EEX 23 CHS 24  8 25 X<>Y 26 X≤Y 27 GTO 33 28 RCL 3  // B 29 STO 2  // Store A 30  1 31 STO+0 32 GTO 06 33 RCL 0 34 RCL 3 35 GTO 00 Gamo RE: (12C) Square Root - Albert Chan - 10-02-2019 02:36 PM It might be better not asking the user to supply a possibly bad guess. Instead, user is limited to enter value between 0.1 to 10 For example, N=12345, user entered n=1.2345, √N = 100 √n A naive guess = 0.5*(1+n) might create big relative errors. Worst case at the edge, n=0.1 or 10 max_rel_err = |1 - 0.5*11/√10| ≈ 0.7393 ≈ 74% Each iteration reduce max_rel_err to around ½(prev_rel_err)², we got: Ref: http://www.azillionmonkeys.com/qed/sqroot.html 0.7393 → 0.2732 → 0.03733 → 0.0006968 → 2.428e-7 → 2.947e-14 Thus, a maximum of 5 iterations to reach 10 digits accuracy. A better guess = k * (1+n), such that rel_err(n=1) = - rel_err(n=10). In other words, |rel_error| for n = 0.1 to 10 have 3 peaks, with W shape. XCas> simplify(solve(1-k*2/1 = k*11/sqrt(10) - 1, k)) ﻿ ﻿ ﻿ ﻿ → k = (22√10 - 40)/81 ≈ 0.365 guess = 0.365 * (1 + n) reduced max_rel_err to 27%, thus required at most 4 iterations. Ref: Numerical Analysis by Ridgway Scott, section 1.3.2, the best start RE: (12C) Square Root - Gamo - 10-03-2019 02:01 AM Thanks Albert Chan Very interesting article. I have streamline this program a bit, now no counter, with an educated guess, program will take on average less than 5 iterations. This program update included Pause to observe each iterations until converge to an answer. Program: Code:  01 STO 2 02 R↓ 03 STO 1 04 RCL 1 05 RCL 2 06  ÷ 07 LSTx 08  + 09  2 10  ÷ 11 PSE 12 STO 3 13 RCL 2 14  - 15 LSTx 16  ÷ 17  0 18 X<>Y 19 X≤Y 20 CHS 21 EEX 22 CHS 23  8 24 X<>Y 25 X≤Y 26 GTO 30 27 RCL 3 28 STO 2 29 GTO 04 30 RCL 3 Hints: Program line 17 to 20 is the [ABS] function with this routine we avoid using n [ENTER] [x] [√x] because it is kind of funny to use [√x] in program while this program is looking for the Square Root. This [ABS] function routine is also work good on HP-10C I have post earlier here https://www.hpmuseum.org/forum/thread-12138.html Gamo RE: (12C) Square Root - Albert Chan - 09-28-2020 05:18 PM (10-02-2019 02:36 PM)Albert Chan Wrote:  guess = 0.365 * (1 + n) reduced max_rel_err to 27%, thus required at most 4 iterations. Turns out above is not the optimal. It should be guess = 0.37913 * (1+n), 0.1 ≤ n ≤ 10 For increasing function, Newton's method preferred a guess on the right. I learned about this from "Fast Inverse Square Root", by Chris Lomon One application of Newton's method automatically gives a guess, to the right of actual root. Example, sqrt(n = 0.7) ≈ 0.83666 x = 0.37913*(1+n) ≈ 0.64452 (guess to the left of root) x = (x + n/x) / 2 ﻿ ﻿ ﻿ ﻿ ≈ 0.86530 (guess, now to the right) XCas> nt(x,n) := (x + n/x)/2 XCas> relerr(k,n) := 1 - nt(k*(1+n),n) / sqrt(n) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ // note: relative error ≤ 0 XCas> solve(relerr(k,1) = relerr(k,10) and k>0, k) ﻿ ﻿ ﻿ → $$[\sqrt{\sqrt{10} / 22}]$$ XCas> k0 := approx(ans()[0])﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → 0.379130444101 Relative error (worst case) at n=0.1, or 1., or 10. With guess = k0*(1+n), how many iterations to ensure full convergence ? Try n=1 (note: √1=1) XCas> nt(k0*2﻿, 1) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿→ 1.03853409762 XCas> nt(ans(), 1) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → 1.00071489067 XCas> nt(ans(), 1) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → 1.00000025535 XCas> nt(ans(), 1) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → 1.0 4 Newton iterations will converged to full 12 digits P.S. above shows Newton's method on square-root, convergence rate ≈ ε²/2 RE: (12C) Square Root - Albert Chan - 09-28-2020 07:14 PM (10-02-2019 02:36 PM)Albert Chan Wrote:  Ref: Numerical Analysis by Ridgway Scott, section 1.3.2, the best start Scott's best start = (1+n) * a, where a = -8 + 6√2 ≈ 0.485281374239 was a little off. Minimized relative errors of guess and actual root is not enough. Taking account the effect of Newton's method, reusing previous post defined relerr() XCas> solve(relerr(k,1) = relerr(k,2) and k > 0, k) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → $$[\sqrt{\sqrt{2} / 6}]$$ XCas> a := approx(ans()[0]) ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → 0.485491771707 I tested this in lua, with added 3 Newton's iterations. Note: Newton's iterations treated in pairs, to optimize away 1 multiply. Code: function fastsqrt(x)     local g, e = frexp(x)     if e%2 ~= 0 then g, e = g+g, e-1 end     x = (1+g) * 0.97098 -- factor to miminize relative error     x = x*0.25 + g/x    -- abs(1- x/g) < 4.337e-4     x = x + g/x         -- abs(1-2x/g) < 9.400e-8     return ldexp(x*0.25 + g/x, e*0.5) end Expected maximum relative error ≈ ε²/2 = (9.4e-8)²/2, ≈ 4.4e-15 And, it should occur when x is pow-of-2, or close to it. lua> x = 1 lua> for i=1,11 do print(x, fastsqrt(x)); x=x*100 end 1﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ 1.0000000000000044 100 ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ 10 10000 ﻿ ﻿ ﻿ ﻿ 100 1e+006 ﻿ ﻿ ﻿ 1000.0000000000041 1e+008 ﻿ ﻿ ﻿ 10000 1e+010 ﻿ ﻿ ﻿ 100000 1e+012 ﻿ ﻿ ﻿ 1000000.0000000033 1e+014 ﻿ ﻿ ﻿ 1e+007 1e+016 ﻿ ﻿ ﻿ 1e+008 1e+018 ﻿ ﻿ ﻿ 1000000000.0000021 1e+020 ﻿ ﻿ ﻿ 1e+010 RE: (12C) Square Root - SlideRule - 09-28-2020 09:45 PM Report No. NADC-91067-50, the Square Root CORDIC, NAVAL AIR SYSTEMS COMMAND, NAVAL AIR DEVELOPMENT CENTER, Mission Avionics Technology Department, AD-A242 318, JUL 1991, may be of minor interest: "     The CORDIC (Coordinate Rotation Digital Computer) algorithm, computes certain functions such as the sine, cosine, and √(x² + y²) using only additions and bit shifting operations.      We have implemented an integer math CORDIC algorithm on a high speed RISC processor. During the course of this work, we identified a convergence problem with the √(x² + y²) CORDIC.  A solution to this problem is presented along with an overview of this algorithm. " BEST! SlideRule RE: (12C) Square Root - Albert Chan - 06-08-2021 03:36 PM Inspired by ∫(1/(a - cos(x)), x=0..pi) = pi/√(a²-1), I find another way to get 1/√(a²-1) For x = 0 to pi, cos(x) = 1 to -1, we can approximate integral without cos(x) term: ∫(1/(a - cos(x)), x=0..pi) ≈ pi/a ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → 1/√(a²-1) ≈ 1/a ∫(1/(a - cos(x)), x=0..pi) = ∫(1/(a - cos(x)) + 1/(a + cos(x)), x=0 .. pi/2); ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ // fold the integral = 2*a * ∫(1/ (a² - cos(x)²), x = 0..pi/2) = 4*a * ∫(1/ (2*a² - (1 + cos(2*x))), x = 0..pi/2); ﻿ ﻿ ﻿ ﻿ ﻿ // let y = 2*x = 2*a * ∫(1/ (a2-cos(y)), y = 0..pi); ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ // let a2=2*a²-1 Again, dropping cos(y), we have: ∫(1/(a - cos(x)), x=0..pi) ≈ 2*a/a2 * pi ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → 1/√(a²-1) ≈ (1/a)·(1+1/a2) Rinse and repeat, let ak = 2*ak-1² - 1, we have 1/√(a²-1) = (1/a)·(1+1/a2)·(1+1/a3) ... Convergence is quadratic, example, with a=2 1/√3 = (1/2)·(1+1/7)·(1+1/97)·(1+1/18817)·(1+1/708158977) ... Or, we flip both side: √3 = 2·(1-1/8)·(1-1/98)·(1-1/18818)·(1-1/708158978) ... function seq(a) -- sequence converging to sqrt(a*a-1), a > 1 local t = a return function() a=2*a*a; t=t-t/a; a=a-1; return t end end lua> g = seq(2) lua> for i=1,4 do print(i, g()) end 1 1.75 2 1.7321428571428572 3 1.7320508100147276 4 1.7320508075688774 This is equivalent to Newton's method for √N, N=a*a-1, guess x=a: lua> N, x = 3, 2 lua> for i=1,4 do x=(x+N/x)/2; print(i, x) end 1 1.75 2 1.7321428571428572 3 1.7320508100147274 4 1.7320508075688772 RE: (12C) Square Root - depor - 12-21-2023 11:50 PM I will calculate quickly on a regular calculator, so that there are more numbers on the display. There is a very simple calculation algorithm that was used in Soviet Iskra (Spark) computers back in the 1970s. All that is required is the arithmetic operation of subtraction and multiplication by 10 (in essence, this is the addition of the inverse of the number and the shift registr of the left). We divide the number into groups of 2 digits before and after the decimal point. If one digit remains, it counts as a group. Next, we use the amazing property of odd numbers — their sum is equal to the square of the number of items. For example 1=1²; 1+3=4=2²; 1+3+5=9=3²; 1+3+5+7=16=4². We proceed as follows - from the first group on the left, we subtract odd numbers, starting with 1. Example. 10-1=9-3=6-5=1 The number of successful subtractions of 3 is the first digit of the root. We add the numbers of the next group. 100 We subtract odd numbers, start with the number obtained in this way — the number that gave a negative result remains, multiply by 10 and subtract 9. 7*10-9=61 100-61=39-63=-24 The second digit of the root is 1. The root is 3.1 New group 3900 The new number is 621 3900-621=3279-623=2656-625=2031-627=1404-629=775-631=144-633=-489 The third digit of the root is 6. The root is 3.16 New group 14400 The new number is 6321 14400-6321=8079-6323=1756 1756-6325=-4569 The root is 3.162 New group 175600 The new number is 63241 175600-63241=112359-63243=49116 The root is 3.1622 Group 4911600 The number is 632421 4911600-632421=4279179-632423=3646756-632425=3014331-632427=2381904-632429=1749475-11170444-632431=484611 The root is 3.1627 Group 48461100 The number is 6324321 And so on. The number of digits after the decimal point is determined by the size of the number with which the device works. For a pen and pencil, this is a lot, and for a calculator, it is 4-6 characters, because only up to 13 digits of the mantissa. RE: (12C) Square Root - Dave Hicks - 12-23-2023 01:48 AM Interesting algorithm! I think you have a typo at: “The root is 3.1627”. It should be 3.16227. Just a typo, the algorithm correctly produced the 2nd “2” in the previous iteration. RE: (12C) Square Root - Thomas Klemm - 12-23-2023 04:26 AM (12-23-2023 01:48 AM)Dave Hicks Wrote:  Interesting algorithm! This is the good old digit by digit method. Here's a program for the HP-42S: Code: 00 { 50-Byte Prgm } 01▸LBL "√" 02 1 03 X<>Y 04 0 05 X<>Y 06▸LBL 00 07 RCL- ST Z 08 X<0? 09 GTO 01 10 1 11 STO+ ST Z 12 STO+ ST X 13 STO+ ST T 14 R↓ 15 GTO 00 16▸LBL 01 17 RCL+ ST Z 18 10 19 STO× ST T 20 STO× ST Z 21 X↑2 22 × 23 9 24 STO- ST T 25 R↓ 26 STOP 27 GTO 00 28 END Example 10 XEQ "√" y: 30 x: 100 R/S y: 310 x: 3900 R/S y: 3160 x: 14400 R/S y: 31620 x: 175600 R/S y: 316220 x: 4911600 R/S y: 3162270 x: 48447100 But we can do better using David Cochran's trick, multiplying the number by 5. Therefore we only have to keep track of the subtrahend. We just ignore the trailing 5 in the result. Code: 00 { 45-Byte Prgm } 01▸LBL "√" 02 5 03 × 04 LASTX 05 X<>Y 06▸LBL 00 07 RCL- ST Y 08 X<0? 09 GTO 01 10 10 11 STO+ ST Z 12 R↓ 13 GTO 00 14▸LBL 01 15 RCL+ ST Y 16 10 17 STO× ST Z 18 X↑2 19 × 20 45 21 STO- ST Z 22 R↓ 23 STOP 24 GTO 00 25 END Example 10 XEQ "√" y: 305 x: 500 R/S y: 3105 x: 19500 R/S y: 31605 x: 72000 R/S y: 316205 x: 878000 R/S y: 3162205 x: 24558000 R/S y: 31622705 x: 242235500 This method is used in most HP calculators.