procedure rsqrt(var x1,x2: hreal);
var
  n: integer;
  i,j: integer32;
  x,xt: hreal;
begin { rsqrt }
  trunc:=false;
  if x1.si<0 then writeln('rsqrt: error - negative argument')
  else if x1.si=0 then x2.si:=0
  else begin
    i:=10000*x1.ma[0]+x1.ma[1]+1;
    j:=1000;
    for n:=1 to 6 do j:=(j+(i div j)) div 2;
    while j*j<i do j:=j+1;
    with x2 do begin                      { crude upper bound on sqrt(x1) }
      si:=1;
      if odd(x1.ex) then begin
        ma[0]:=j;
        for n:=1 to hdim do ma[n]:=0;
        if x1.ex>0 then ex:=x1.ex div 2
        else ex:=(x1.ex div 2)-1;
        if j>9999 then begin
          ma[0]:=j div 10000;
          ma[1]:=j mod 10000;
          ex:=ex+1;
        end;
      end
      else begin
        ma[0]:=j div 100;
        ma[1]:=100*(j mod 100);
        for n:=2 to hdim do ma[n]:=0;
        ex:=x1.ex div 2;
      end;
    end;
    repeat                                                 { Newton steps }
      rquot(x1,x2,x);
      if trunc then rup(x);
      rdiff(x,x2,xt);
      if trunc then rup(xt);
      irdiv(2,xt);
      if trunc then rup(xt);
      x:=x2;
      rsum(x,xt,x2);
      if trunc then rup(x2);
    until irdiff(x2,x)>=0;
    repeat
      rprod(x2,x2,x);
      j:=irdiff(x,x1);
      if j>0 then rdown(x2);
    until j<=0;
    if j<0 then trunc:=true
    else if trunc then rdown(x2);
  end;
end { rsqrt };

