ISqrt Function
ISqrt (I) computes INT (SQRT (I)), that is, the integral part of the square root of integer I.

Unit
QESBPCSMath

Declaration
Function ISqrt(const I: LongWord): Longword;

Description
Code originally developed by Marcel Martin, used with permission. Rory Daulton introduced a faster routine (based on Marcel's) for most occassions and this is now used with Permission.

Parameters
Positive Integer Value to process.

Category
Arithmetic Routines for Integers

Implementation

function ISqrt (const I: LongWord): LongWord;
const
     Estimates: array [0..31] of Word = (
          // for index i, constant := Trunc(Sqrt((Int64(2) shl i) - 1.0)), which is
          //   the largest possible ISqrt(n) for  2**i <= n < 2**(i+1)
          1, 1, 2, 3, 5, 7, 11, 15,
          22, 31, 45, 63, 90, 127, 181, 255,
          362, 511, 724, 1023, 1448, 2047, 2896, 4095,
          5792, 8191, 11585, 16383, 23170, 32767, 46340, 65535);
     // eax  // ebx  // ecx  // edx
asm  // entry:   // eax = I
  // calc the result quickly for zero or one (sqrt equals the argument)
	   cmp     eax, 1
	   jbe     @@end
  // save registers and the argument
	   push    ebx
	   mov     ebx, eax                // ebx = I
  // use the logarithm base 2 to load an initial estimate, which is greater
  //   than or equal to the actual value
	   bsr     eax, ebx        // eax = ILog2(I) (note upper WORD is now zero)
	   mov     ax, [word ptr Estimates + eax * 2]
						  // eax = X
  // repeat ...
@@repeat:
  // --  save the last estimate [ X ]
	   mov     ecx, eax                        // ecx = X
  // --  calc the new estimate [ (I/X + X) / 2 ; the Newton-Raphson formula ]
	   xor     edx, edx                                // edx = 0
	   mov     eax, ebx        // eax = I (so edx:eax = I)
	   div     ecx             // eax = I/X            // edx = I mod X
	   add     eax, ecx        // eax = I/X + X
	   shr     eax, 1          // eax = XNew = (I/X+X)/2
  // until the new estimate >= the last estimate
  //   [which can never happen in exact floating-point arithmetic, and can
  //   happen due to truncation only if the last estimate <= Sqrt(I) ]
	   cmp     eax, ecx
	   jb      @@repeat
  // use the next-to-last estimate as the result
	   mov     eax, ecx        // eax = X
  // restore registers
	   pop     ebx
@@end:              //exit:     // eax = Result
End;


HTML generated by Time2HELP
http://www.time2help.com