Post Reply 
(12C) Square Root
10-02-2019, 10:23 AM (This post was last modified: 10-03-2019 04:10 AM by Gamo.)
Post: #1
(12C) Square Root
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
Find all posts by this user
Quote this message in a reply
10-02-2019, 02:36 PM (This post was last modified: 10-15-2020 11:45 PM by Albert Chan.)
Post: #2
RE: (12C) Square Root
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
Find all posts by this user
Quote this message in a reply
10-03-2019, 02:01 AM (This post was last modified: 10-03-2019 04:13 AM by Gamo.)
Post: #3
RE: (12C) Square Root
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. Smile
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
Find all posts by this user
Quote this message in a reply
09-28-2020, 05:18 PM
Post: #4
RE: (12C) Square Root
(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
Find all posts by this user
Quote this message in a reply
09-28-2020, 07:14 PM
Post: #5
RE: (12C) Square Root
(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
Find all posts by this user
Quote this message in a reply
09-28-2020, 09:45 PM
Post: #6
RE: (12C) Square Root
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
Find all posts by this user
Quote this message in a reply
06-08-2021, 03:36 PM (This post was last modified: 06-09-2021 11:24 AM by Albert Chan.)
Post: #7
RE: (12C) Square Root




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
Find all posts by this user
Quote this message in a reply
Post Reply 




User(s) browsing this thread: 1 Guest(s)