Integer square root
Sometimes you just need to have the square root of a number, but you need it as an INTEGER, not as a REAL. And you don't want to use the detour over Math.Sqrt. There's an easy way around this: we just calculate the root ourself. No big deal at all. Just some basic math required.
There's a simple method to get the job done. It's called the "Newton Raphson approximation". In the graph on the right it is shown. In this case for the function f(x) = x^2 - n. If we find the interception of the function with the X axis, we have the square root of 'n'.
The method is based on the tangent of the function at any point 'x'. Below is the math:
tan (a) = dy / dx = f'(x0) dy = f(x0) dx = x0 - x1So, after substitutions:
f'(x) = f(x) / (x0 - x1)Or
x0 - x1 = f(x0) / f'(x0) x1 = x0 - f(x0) / f'(x0)In this case:
f(x) = x^2 - N f'(x) = 2.xResulting in:
x1 = x0 - (x0^2 - N) / 2.x0 = x0 - (x0 / 2) + N / 2.x0And that's it. You just pick a value for x0 and calculate x1. Now copy x1 to x0 and repeat. Within very few iterations you have reached the square root.
Cast it into code
Below is the Oberon source code for the Newton Raphson method.
PROCEDURE sqrt (number : INTEGER) : INTEGER; VAR x1, x0, t : INTEGER; BEGIN x0 := number; REPEAT t := x0; x1 := x0 - x0 DIV 2 + number DIV (x0 + x0); Out.Int (x1, 5); Out.Char (Tab); x0 := x1 UNTIL x1 = t; RETURN x1 END sqrt;Now let's see how it works. First make a program around it:
MODULE sq; IMPORT In, Out; VAR n : INTEGER; PROCEDURE sqrt (number : INTEGER) : INTEGER; VAR x1, x0, t : INTEGER; BEGIN x0 := number; REPEAT t := x0; x1 := x0 - x0 DIV 2 + number DIV (x0 + x0); Out.Int (x1, 5); Out.Char (09X); x0 := x1 UNTIL x1 >= t; RETURN x1 END sqrt; BEGIN LOOP Out.String ("Enter a number : "); In.Int (n); IF n < 0 THEN EXIT END; n := sqrt (n); Out.Ln; Out.String (" The integer root is "); Out.Int (n, 6); Out.Ln END; Out.Ln END sq.As you see I have changed the iteration algorithm. The new version has two divisions, no multiplications. I also changed the UNTIL condition of 'sqrt'. It now says 'x1 >= t'. At certain values the iterations get into an infinite flipflop between values n and n+1. Now the loop stops when x1 = x0 or when x1 rises above x0.
Remember, this is an integer version of REAL number math. Rounding errors are inevitable. Still, an astonishing result is achieved and within very few iterations a remarkably good value is reached.
The current algorithm is:
x1 := x0 - x0 DIV 2 + number DIV (x0 + x0);It is better to write this as
x1 := (x0 + number DIV x0) DIV 2;'number DIV x0' is a 'difficult division' and 'DIV 2' is merely a shift right one position. But now 5 operations are reduced to 3 operations. In integer math, less operations is more accurate since there are less rounding errors.
But is it better?
Left is the old algorithm, right is the improved algorithm. Not much of an improvement. There is less flipflopping at the stop condition. But in the new algorithm, 6 as a root pops up at n=35, in the old version it popped up at n=30. So the improvement is scrapped again....
See it run
This is what sq does. The red lines show flipflopping results that require the 'fuzzy UNTIL condition'.
oxygen:~/Oberon/prime$ ./sq Enter a number : 16 8 5 4 4 The integer root is 4 Enter a number : 36 18 10 6 6 The integer root is 6 Enter a number : 40 20 11 7 6 6 The integer root is 6 Enter a number : 41 21 11 7 6 6 The integer root is 6 Enter a number : 44 22 12 7 7 The integer root is 7 Enter a number : 300000
150000 75001 37502 18754 9384 4707 2385 1255 747 574 548 547 548The integer root is 548 Enter a number : 29 15 8 5 5 The integer root is 5 Enter a number : 30 15 9 6 5 6The integer root is 6 Enter a number : 10000000 5000000 2500001 1250002 625004 312509 156270 78166 39146 19700 10103 5546 3674 3197 3162 3162 The integer root is 3162 Enter a number : -9 oxygen:~/Oberon/prime$
Page created 27 Jan 2017,