Een integer square root algoritme

Delphi kent de operators div en mod voor integer divide en rest berekening.
Maar worteltrekken betreft alleen floating point getallen.
Dit project voorziet in worteltrekken uit gehele getallen, waarbij ook de rest beschikbaar komt.

Voor geheel getal N schrijven we:
    N = wortel2 + rest.
De plaatjes hierna tonen het programma aan het werk:



Het algoritme

Het integer wortel programma werkt met positieve getallen van 32 bit (cardinal, dword) voor N en
16 bit positieve integers (word) voor de wortel en de rest.

Voor de duidelijkheid volgt hier een voorbeeld waarbij N een 8 bit getal is.
We berekenen de wortel uit 121 = $79 = [0111 1001]2.

Vooraf verdelen we N in groepjes van 2 bits, vanaf de rechterzijde gerekend.
Variabelen N, root en b zijn bytes.

b is een 1 bit waarde (rood getekend) die staat in het begin rechts in het meest linkse groepje van 2 bits van N.
De variable root is aanvankelijk nul (zwart getekend).
In het plaatje hierna hebben lege cellen de waarde nul.

N wordt verminderd met root + b als N >= root + b.



Na een poging tot vermindering van N:
  • Als aftrekking plaatsvond dan b bit kopieren naar root
  • root 1 bit positie naar rechts schuiven
  • het b bit 2 posities naar rechts schuiven
Bovenstaande stappen herhalen totdat b buiten het byte is geschoven.
N bevat dan de rest.

Waarom dit werkt

Beschouw een tweetallig getalletje van 2 bits a,b : [ab]2 = 2a+b.
Het kwadraat is (2a+b)2 = 4a2 + 4ab + b2
Dat is een vier bits getalletje.

Nu is a2 = a, b2 = b omdat a en b enkele bits zijn.
Als a = 1:
4a2=[100]2 en 4ab+b2=[101]2 als b=1.
Als N >= [100] dan wordt [100] afgetrokken van N en a=1.
Nu nog de waarde van b vinden.
Als a=1 en ook b=1 dan kunnen we 4ab + b= [101] aftrekken van N.
Dit betreft maar 4 bits maar bedenk:
deze bits kunnen de meest linkse bits zijn van een groter getal.
Na elke poging tot vermindering van N vinden we een waarde van bit b en kunnen we a vervangen door [ab]2.
N is al verminderd met a2 dus een volgend bit b vinden we door een test of N verminderd
kan worden met [a01]2.

Het programma

N is een 32 bit positief getal waarvan we de wortel berekenen.
Aan het eind staat de rest in N.
Variabele root bevat de waarde van a.
Variabele b is b in de uitleg hiervoor.
Vooraf zetten we a op nul, b is de eerste waarde die we trachten af te trekken van N om a te vinden.
Hieronder staat het programma: 
function IntRoot(N : dword; var rem : word) : word;
//return integer square root of N and remainder rem
var b,root,sub : dword;
begin
 b := $40000000;
 root := 0;
 repeat
  sub := root or b;
  if sub <= N then
   begin
    N := N - sub;
    root := (root shr 1) or b;
   end else root := root shr 1;
  b := b shr 2;
 until b = 0;
 result := root;
 rem := N;
end;
procedure TForm1.GoBtnClick(Sender: TObject);
var N,root : dword;
    rem : word;
begin
 if length(numberEdit.Text) = 0 then numberEdit.Text := '0';
 N := strtoInt(numberEdit.Text);
 root := IntRoot(N,rem);
 rootLabel.Caption := inttostr(root);
 remainderlabel.Caption := inttostr(rem);
end;

procedure TForm1.FormKeyPress(Sender: TObject; var Key: Char);
var OK,mt : boolean;
begin
 OK := false;
 with numberEdit do
  begin
   mt := length(text) = 0;
   case key of
    #8       : OK := true;
    '1'..'9' : OK := true;
    '0'      : OK := not(mt);
    '$'      : OK := mt;
    'a'..'f' : OK := text[1] = '$';
   end;
   if OK = false then key := #0;
  end;
end;


 download root programma

 download Delphi(7) project