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.

Newton Raphson

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 - x1
So, 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.x
Resulting in: 
x1 = x0 - (x0^2 - N) / 2.x0 = x0 - (x0 / 2) + N / 2.x0
And 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.
A good initial value for x0 is 'N' itself.

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.

Small improvement

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     548
 The 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       6
 The 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$