1 % Square root using Newton's Method
 2 % Octave 5.1.0
 3 % John Bryan 2019
 4 
 5 
 6 echo off;
 7 close all;
 8 clear all;
 9 format long;
10 
11 
12 function [r0,mul]=normalize(r)
13    % normalize to 0.25 <= r <= 1
14    % mul is denormalizing factor
15    mul = 1.0; % power of 2
16    r0=r;
17    while (r0 > 1)
18          r0 = r0/4.0;
19          mul = mul*2.0;
20    endwhile
21    while (r0 < 0.25)
22       r0 = r0*4.0;
23       mul = mul/2.0;
24    endwhile
25    %r0 now normalized
26    return;
27 endfunction
28 
29 
30 function x1=sqrt0(r0)
31       % initial approx x0: (2/3)r0 + 0.3506
32       x0 = 0.6666666666666667*r0+0.3506;
33       printf('\nx[0]=%.14f\n',x0);
34       i=1;
35       while(1)
36          x1 = x0/2.0 + r0/(2.0*x0);
37          printf('\nx[%d]=%.14f\n',i,x1);
38          if (abs(x1 - x0) < 1.0e-15) return; endif
39          x0 = x1;
40          i=i+1;
41       endwhile
42       return;
43 endfunction
44 
45 
46 function r2=denorm(sqrtr0,mul)
47     r2= sqrtr0*mul;
48     return;
49 endfunction
50 
51 
52 [r0,mul]=normalize(13.0)
53 x1=sqrt0(r0);
54 r2=denorm(x1,mul);
55 printf('\napproximate answer of sqrt(13.0) = %.14f\n',r2);
56 printf('\ncorrect answer of sqrt(13.0) = %.14f\n\n',sqrt(13.00000000000000));