1 % Newton-Raphson Division
 2 % John Bryan 2019
 3 % Octave 5.1.0
 4 
 5 
 6 echo off;
 7 close all;
 8 clear all;
 9 format long;
10 
11 
12 function [d0,mul]=normalize(d)
13    % normalize to 0.5 <= d <= 1
14    % mul is denormalizing factor
15    mul = 1.0 % power of 2
16    d0=d;
17    while (d0 > 1)
18          d0 = d0/2.0;
19          mul = mul*2.0;
20    endwhile
21    while (d0 < 0.5)
22       d0 = d0*2.0;
23       mul = mul/2.0;
24    endwhile
25    %d0 now normalized
26    return;
27 endfunction
28 
29 
30 function x1=invD(d0)
31       % calculate 1/d0, 0.5 <= d0 <= 1
32       % initial approx. x0 = 48/17 - (32/17)*d0
33       x0 = 2.82353 - 1.88235*d0
34       P = 53;
35       iterations=ceil(log2((P+1)/log2(17)))
36       for i=[1:iterations]
37          x1 = x0 + x0 * (1-d0*x0);
38          x0 = x1;
39          printf('\nx[%d]=%.15f\n',i,x0);
40       endfor
41       return;
42 endfunction
43 
44 
45 function d2=denorm(X1,mul)
46     d2= X1/mul;
47     return;
48 endfunction
49 
50 
51 function [fpn,e]=fp(m)
52    fpn=m;
53    e=0;
54    if  (m<2) && (m>=1)
55       return
56    else
57       if (m>=2)
58          while(1)
59             fpn=fpn/2;
60             e=e+1;
61             if  (fpn<2)
62               return
63             endif
64          endwhile
65       endif
66       if (m<1)
67          while(1)
68             fpn=fpn*2;
69             e=e-1;
70             if  (fpn>=2)
71               break
72             endif
73           endwhile
74        endif
75    endif
76    return
77 endfunction
78 
79 
80 
81 n=[32.0,-15.27];
82 d=[27.0,34.34];
83 for i=[1:length(n)]
84    printf('\n\n');
85    [dfp,e]=fp(abs(d(i)))
86    d0=dfp/2
87    n0=abs(n(i)/(2^(e+1)))
88    X1=invD(d0)
89    signq=sign(n(i))*sign(d(i));
90    answer=signq*n0*X1
91    check=n(i)/d(i)
92    printf('\n\n');
93 endfor