Delphi和Free Pascal中的SQRT函数有多精确?

3

SQRT在Delphi XE中作为80位浮点值的FPU函数实现;不确定64位编译器中是如何实现的。浮点数函数被认为是近似的。

我能否假设下面的断言永远不会失败?

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;

只是需要注意的是,在x64上扩展是一个64位浮点值。有关跨平台应用程序的Delphi注意事项,请参见http://docwiki.embarcadero.com/RADStudio/XE3/en/Delphi_Considerations_for_Cross-Platform_Applications。 - Sir Rufo
2
你能解释一下为什么你认为这些断言应该成立吗?另外,这些if测试是用来做什么的? - David Heffernan
为什么要截断Sqrt()的结果?只是为了查看SQRT()是否过大,导致其精度误差超过+1.0吗?为什么不测量实际误差并报告准确的最大值呢? - Warren P
1个回答

2

实践比理论更重要:

对所有数字进行测试,例如:

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.

由于测试UInt64的整个范围确实需要很长时间,因此我稍微修改了测试方法,改为只测试所有完全平方数及其前后的数字,以加快速度并获得更好的结果。我个人对32位进行了一段时间的测试,没有出现故障(整个测试的1%),而在64位上则很快就会出现故障。我仍在仔细研究这个问题,但我已经发布了代码,以防您有兴趣:

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.

1
你的测试不完整。test1 的上限应为 $FFFFFFFF,test2 的上限应为 $FFFFFFFFFFFFFFFF,包括上限值。此外,需要编译器版本和位数(32位或64位编译器)。 - kludg
1
使用XE3进行测试(根据@Serge的评论更改限制)显示此测试失败(Win64,XE3编译器版本17.0.4770.56661)。它在值为4294836225的Test1处失败,这似乎是合理的。提供了一个很好的基础,可以在答案提供的范围内开始(并进行更正),加一分。 - Ken White
@KenWhite 你提到的失败值是$FFFE0001,非常有趣。我认为这是SQRT函数实现中的一个错误。 - kludg
@Serg:我发布了两个实现(Win32的FSQRT,Win64的SQRTSD)。也许你应该将错误报告发送给英特尔? :-) 至于你提出的问题(“我可以假设下一个断言永远不会失败吗?”),我的回答始终如一 - 一个人永远不应该假设任何事情。;-) - Ken White
1
@Serg: "Raises new questions" 意味着你应该发布一个“新问题”。;-) 你应该将其作为新问题发布。正如你最后的评论所说,这篇文章“回答了问题”。这似乎正在变成一次聊天,在这里不合适。我仍然认为这个答案值得我的点赞,现在就这样吧。希望你能在这里学到真正想学的东西。祝好运。 :-) - Ken White
显示剩余10条评论

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接