Skip to content

Commit

Permalink
Improve precision of FSQRT - #42
Browse files Browse the repository at this point in the history
Test output after the change:

Running test of square root(x).
Testing if sqrt(X * X) == X for 20 Integers X.
Test for sqrt monotonicity.
sqrt has passed a test for Monotonicity.
Testing whether sqrt is rounded or chopped.
Square root is neither chopped nor correctly rounded.
Observed errors run from 0.E0 to 5.000000000000000000E-1  ulps.
  • Loading branch information
ikysil committed Sep 27, 2017
1 parent 07cba18 commit 9d8f991
Showing 1 changed file with 25 additions and 6 deletions.
31 changes: 25 additions & 6 deletions lib/~ik/float-ieee-binary/f-sqrt.4th
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
64 CONSTANT FSQRT-ITERATIONS-DEFAULT
32 CONSTANT FSQRT-ITERATIONS-DEFAULT


: FSQRT-APPROX (F r1 -- r2 )
\G r2 is the initial approximation of square root of r1.
Expand All @@ -8,31 +9,48 @@
'FPX FPE!
;


: FSQRT-NEWTON-STEP (F r1 r2 -- r3 )
\G Perform step in Newton approximation for a square root.
\G Calculate next approximation r3 for the square root of
\G r1 given the previous approximation r2.
?FPX0= IF FNIP EXIT THEN
?FPX0= IF FNIP EXIT THEN
FSWAP FOVER \ F: r2 r1 r2
F/ F+
FPX2/
FTWO F/
;

: FSQRT-NEWTON (S +n -- ) (F r1 r2 -- r3 )

: FSQRT-NEWTON-STEP2 (F r1 r2 -- r3 )
\G Perform step in Newton approximation for a square root.
\G Calculate next approximation r3 for the square root of
\G r1 given the previous approximation r2.
?FPX0= IF FNIP EXIT THEN
FSWAP FOVER FDUP F* \ F: r2 r1 r2**2
FSWAP F- \ F: r2 r2**2-r1
FOVER FTWO F* F/ \ F: r2 (r2**2-r1)/(2*r2)
F-
;


: FSQRT-NEWTON (S +n xt -- ) (F r1 r2 -- r3 )
\G Calculate square root of r1 given the initial approximation r2 and number of iterations +n.
\G xt calculates next approximation with stack effect (F r1 r2 -- r3 ).
1 ?FPSTACK-OVERFLOW
>R
BEGIN
DUP 0>
?FPX0= INVERT \ stop iterations if approximation eq. zero
AND
WHILE
FOVER FSWAP
\DEBUG S" FSQRT-NEWTON-A1: " CR TYPE CR F.DUMP CR
FSQRT-NEWTON-STEP
R@ EXECUTE
\DEBUG S" FSQRT-NEWTON-A2: " CR TYPE CR F.DUMP CR
1-
REPEAT
DROP
R> DROP
FNIP
\DEBUG S" FSQRT-NEWTON-B: " CR TYPE CR F.DUMP CR
;
Expand All @@ -45,5 +63,6 @@
?FPX0= IF EXIT THEN
?FPX0< IF FPX-NAN! EXIT THEN
?FPX-INF IF EXIT THEN
FSQRT-ITERATIONS-DEFAULT FDUP FSQRT-APPROX FSQRT-NEWTON
FDUP FSQRT-APPROX
FSQRT-ITERATIONS-DEFAULT ['] FSQRT-NEWTON-STEP2 FSQRT-NEWTON
;

0 comments on commit 9d8f991

Please sign in to comment.