Search code examples
delphifpcsqrt

How accurate is SQRT function in Delphi and Free Pascal?


SQRT is implemented as a FPU function on 80-bit float value in Delphi XE; not sure how it is implemented in 64-bit compilers. Floating point functions are known to be approximate.

Can I assume that the next assertions will never fail?

procedure Test1(Value: Cardinal);
var
  Root: Cardinal;

begin
  Root:= Trunc(Sqrt(Value));
  Assert(Root * Root <= Value);
  if Root < $FFFF then
    Assert((Root + 1) * (Root + 1) > Value);
end;

procedure Test2(Value: UInt64);
var
  Root: UInt64;

begin
  Root:= Trunc(Sqrt(Value));
  Assert(Root * Root <= Value);
  if Root < $FFFFFFFF then
    Assert((Root + 1) * (Root + 1) > Value);
end;

Solution

  • More practice than theory:

    Perform a test on all numbers, like this:

    program Project1;
    
    {$APPTYPE CONSOLE}
    
    {$R *.res}
    
    uses
      System.SysUtils;
    
    {$IFNDEF DEBUG}
      {$DEFINE DEBUG}
    {$ENDIF}
    
    procedure Test1(Value: Cardinal);
    var
      Root: Cardinal;
    
    begin
      Root:= Trunc(Sqrt(Value));
      Assert(Root * Root <= Value);
      if Root < $FFFF then
        Assert((Root + 1) * (Root + 1) > Value);
    end;
    
    procedure Test2(Value: UInt64);
    var
      Root: UInt64;
    
    begin
      Root:= Trunc(Sqrt(Value));
      Assert(Root * Root <= Value);
      if Root < $FFFFFFFF then
        Assert((Root + 1) * (Root + 1) > Value);
    end;
    
    var
      VCar: Cardinal;
      VUInt: UInt64;
    const
      Limit1: Cardinal = $FFFFFFFF;
      Limit2: UInt64 = $FFFFFFFFFFFFFFFF;
    begin
      try
        for VCar := 0 to Limit1 do
        begin
          if (VCar mod 10000000) = 0 then
            Writeln('VCarTest ', VCar, ' ', (VCar / Limit1 * 100):0:2, '%');
          Test1(VCar);
        end;
        Writeln('VCarTest 0 .. $', IntToHex(Limit1, 8), ' passed');
    { commented because cannot be executed in a reasonable time
        VUInt := 0;
        while (VUInt <= Limit2) do
        begin
          if (VUInt mod 2500000) = 0 then
            Writeln('VUIntTest ', VUInt, ' ', (VUInt / Limit2 * 100):0:2, '%');
          Test2(VUInt);
          Inc(VUInt);
        end;
        Writeln('VUIntTest ', VUInt);
        Writeln('All passed');
    }
    
      except
        on E: Exception do
          Writeln(E.ClassName, ': ', E.Message);
      end;
      Readln;
    end.
    

    Since it really takes ages to test the whole range for UInt64, I changed the test a bit to test all perfect squares, the number before and the number after each one, just to make it faster and have a better idea. I personally ran the test for 32 bits for a while without failure (a 1% of the whole test), and on 64 bits it shows failure very fast. I'm still looking closer to this, but I posted the code just in case you're interested:

    program Project1;
    
    {$APPTYPE CONSOLE}
    
    {$R *.res}
    
    uses
      System.SysUtils;
    
    {$IFNDEF DEBUG}
      {$message error 'change your build configuration to Debug!'}
    {$ENDIF}
    
    procedure Test2(Value: UInt64);
    var
      Root: UInt64;
    begin
    //try/except block only for 64 bits, since in 32 bits it makes the process much slower
    {$ifdef CPUX64}
      try
    {$endif}
        Root:= Trunc(Sqrt(Value));
        Assert(Root * Root <= Value);
        if Root < $FFFFFFFF then
          Assert((Root + 1) * (Root + 1) > Value);
    {$ifdef CPUX64}
      except
        Writeln('Fails for value: ', Value, ' root: ', Root
          , ' test: ', (Root + 1) * (Root + 1));
        raise;
      end;
    {$endif}
    end;
    
    var
      RUInt, VUInt: UInt64;
    
    const
      Limit2: UInt64 = $FFFFFFFFFFF00000;
    begin
      try
        RUInt := 1;
        repeat
          Inc(RUInt);
          VUInt := RUInt * RUInt;
          if (RUInt mod 2500000) = 0 then
            Writeln('VUIntTest ', VUInt, ' ', (VUInt / Limit2 * 100):0:4, '%');
          Test2(VUInt - 1);
          Test2(VUInt);
          Test2(VUInt + 1);
        until (VUInt >= Limit2);
        Writeln('VUIntTest ', VUInt);
        Writeln('All passed');
      except
        on E:EAssertionFailed do
          Writeln('The assertion failed for value ', VUInt, ' root base ', RUInt);
        on E: Exception do
          Writeln(E.ClassName, ': ', E.Message);
      end;
      Readln;
    end.