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; |