3

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;
kludg
  • 27,213
  • 5
  • 67
  • 118
  • just to note Extended on x64 is a 64bit floating point value http://docwiki.embarcadero.com/RADStudio/XE3/en/Delphi_Considerations_for_Cross-Platform_Applications – Sir Rufo Mar 09 '13 at 08:00
  • 2
    Can you explain why you think these assertions should hold? Also what are the if tests for? – David Heffernan Mar 09 '13 at 08:33
  • Why do you want to truncate the results of a Sqrt()? Just to see if the SQRT() is ever too large by an amount that renders it innacurate by more than +1.0 ? Why not measure the real errors and report accurate maximums? – Warren P Mar 09 '13 at 20:41

1 Answers1

2

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.
jachguate
  • 16,976
  • 3
  • 57
  • 98
  • 1
    Your tests are incomplete. Upper limits should be $FFFFFFFF for test1 and $FFFFFFFFFFFFFFFF for test2, with upper values included. Also compiler version and bitness (32- or 64-bit compiler) is needed. – kludg Mar 09 '13 at 03:45
  • 1
    A test with XE3 (changing the limits based on @Serge's comment) shows this fails (Win64, XE3 compiler version 17.0.4770.56661). It fails at Test1 with value 4294836225, which seems reasonable. +1 for providing a great basis from which to start (and correct within the range provided in the answer). – Ken White Mar 09 '13 at 05:12
  • The implementation in XE3 (indicated 64-bit compiler version, defined with `{$IFDEF CPUX64}`) is `SQRTSD XMM0, XMM0`, which addresses the `FPU function` portion of the question as well. (The 32-bit version uses `FSQRT`.) – Ken White Mar 09 '13 at 05:20
  • I've updated the test code. Test1 passes in Delphi XE, Test2 cannot be executed in reasonable time and commented. – kludg Mar 09 '13 at 05:54
  • @KenWhite the fail value you mentioned is $FFFE0001, very interesting. I consider it as a bug in `SQRT` function implementation. – kludg Mar 09 '13 at 06:01
  • @Serg: I posted both of the implementations (`FSQRT` for Win32, `SQRTSD` for Win64). Perhaps you should send the bug report to Intel? :-) As far as the question you asked ("Can I assume that the next assertions will never fail?"), my answer would be the same as always - one should never **assume** anything. ;-) – Ken White Mar 09 '13 at 06:04
  • @KenWhite does both 32- and 64-bit tests fail in XE3? As I said 80-bit FPU implementation used in Delphi XE passes the test. – kludg Mar 09 '13 at 06:11
  • @Serg: Don't know. Your question was "Can I assume it never fails?", and the answer is "no, you can't". Whether both fail or not is irrelevant; the answer to what you asked remains the same. :-) – Ken White Mar 09 '13 at 06:20
  • @KenWhite yes it answers the question and you are welcome to post a detailed answer. It also raises new question(s). – kludg Mar 09 '13 at 06:24
  • 1
    @Serg: "Raises new questions" means you should post a "new question". ;-) You should probably post it as such. As your last comment says, this post "answers the question". This seems to be turning into a chat, and it's not appropriate here. I still think this answer deserves my upvote, and I'll leave it at that for now. I hope you get an answer to what it really is you're hoping to learn here. Good luck. :-) – Ken White Mar 09 '13 at 06:29
  • @KenWhite, I cannot confirm the error you mentioned. The program given with the answer triggers no assertion on both platforms - at least on my system. – Uwe Raabe Mar 09 '13 at 15:22
  • @Uwe: I tested with the pre-edited code (before the 4th comment above by user24608) with the mentioned version of XE3, 64-bit, on Win64, with only the modified upper limits as I described. I haven't re-tested since. It's entirely possible that the *different code* may work now. :-) – Ken White Mar 09 '13 at 16:03
  • Great discussion here during my sleep, maybe I'm on the incorrect side of the world. :D. @user246408 I don't get what is the intent of your last edit: define the DEBUG symbol if it is not already defined? – jachguate Mar 09 '13 at 16:29
  • Great discussion it is. I've lost hours pondering over the problem that does not exist. Good thing is that now I will never take SO seriously. It sucks. – kludg Mar 09 '13 at 16:34
  • @user246408 I really don't get that, but just in case, I've added a new project with a different approach to test. It quickly shows failure of the assertion for 64 bits on my machine, just in case you're still interested. – jachguate Mar 09 '13 at 17:16