Насколько точна функция SQRT в Delphi и Free Pascal?

SQRT реализован как функция FPU для 80-битного значения с плавающей запятой в Delphi XE; не уверен, как это реализовано в 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;

person kludg    schedule 09.03.2013    source источник
comment
просто отметим, что Extended на x64 — это 64-битное значение с плавающей запятой RADStudio/XE3/ru/   -  person Sir Rufo    schedule 09.03.2013
comment
Можете ли вы объяснить, почему, по вашему мнению, эти утверждения должны быть верными? И для чего нужны тесты if?   -  person David Heffernan    schedule 09.03.2013
comment
Почему вы хотите обрезать результаты Sqrt()? Просто чтобы убедиться, что SQRT() когда-либо слишком велик на величину, которая делает его неточным более чем на +1,0? Почему бы не измерить реальные ошибки и не указать точные максимумы?   -  person Warren P    schedule 10.03.2013


Ответы (1)


Больше практики, чем теории:

Выполните тест на всех числах, например:

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.
person jachguate    schedule 09.03.2013
comment
Ваши тесты неполные. Верхние пределы должны быть $FFFFFFFF для test1 и $FFFFFFFFFFFFFFFF для test2, включая верхние значения. Также необходима версия и разрядность компилятора (32- или 64-битный компилятор). - person kludg; 09.03.2013
comment
Тест с XE3 (изменение ограничений на основе комментария @Serge) показывает, что это не удается (Win64, версия компилятора XE3 17.0.4770.56661). Он терпит неудачу в Test1 со значением 4294836225, что кажется разумным. +1 за отличную основу для начала (и исправить в пределах диапазона, указанного в ответе). - person Ken White; 09.03.2013
comment
Реализация в XE3 (указана 64-битная версия компилятора, определенная с помощью {$IFDEF CPUX64}) — это SQRTSD XMM0, XMM0, что также касается FPU function части вопроса. (В 32-разрядной версии используется FSQRT.) - person Ken White; 09.03.2013
comment
Я обновил тестовый код. Test1 проходит в Delphi XE, Test2 не может быть выполнен за разумное время и прокомментирован. - person kludg; 09.03.2013
comment
@KenWhite значение ошибки, которое вы упомянули, равно $FFFE0001, очень интересно. Я считаю это ошибкой в ​​реализации функции SQRT. - person kludg; 09.03.2013
comment
@Serg: я опубликовал обе реализации (FSQRT для Win32, SQRTSD для Win64). Возможно, вам следует отправить отчет об ошибке в Intel? :-) Что касается вопроса, который вы задали (Могу ли я предполагать, что следующие утверждения никогда не подведут?), мой ответ будет таким же, как всегда - никогда не следует предполагать что-либо. ;-) - person Ken White; 09.03.2013
comment
@KenWhite не проходят ли 32- и 64-битные тесты в XE3? Как я уже сказал, 80-битная реализация FPU, используемая в Delphi XE, проходит тест. - person kludg; 09.03.2013
comment
@Серг: Не знаю. Ваш вопрос был «Могу ли я предположить, что он никогда не выйдет из строя?», и ответ — нет, вы не можете. Потерпят неудачу оба или нет, не имеет значения; ответ на то, что вы спросили, остается прежним. :-) - person Ken White; 09.03.2013
comment
@KenWhite да, это отвечает на вопрос, и вы можете опубликовать подробный ответ. Это также поднимает новые вопросы. - person kludg; 09.03.2013
comment
@Serg: поднимает новые вопросы означает, что вы должны опубликовать новый вопрос. ;-) Вы, вероятно, должны опубликовать это как таковое. Как говорится в вашем последнем комментарии, этот пост отвечает на вопрос. Кажется, это превращается в чат, и здесь это неуместно. Я все еще думаю, что этот ответ заслуживает моего одобрения, и я пока оставлю его на этом. Я надеюсь, вы получите ответ на то, что на самом деле вы надеетесь узнать здесь. Удачи. :-) - person Ken White; 09.03.2013
comment
@KenWhite, я не могу подтвердить упомянутую вами ошибку. Программа, указанная в ответе, не вызывает никаких утверждений на обеих платформах — по крайней мере, в моей системе. - person Uwe Raabe; 09.03.2013
comment
@Uwe: я тестировал предварительно отредактированный код (до 4-го комментария пользователя 24608 выше) с упомянутой версией XE3, 64-разрядной, на Win64, только с измененными верхними пределами, как я описал. С тех пор я не проводил повторное тестирование. Вполне возможно, что теперь другой код может работать. :-) - person Ken White; 09.03.2013
comment
Отличная дискуссия здесь во сне, может быть, я не на той стороне света. : Д. @ user246408 Я не понимаю, какова цель вашего последнего редактирования: определить символ DEBUG, если он еще не определен? - person jachguate; 09.03.2013
comment
Отличная дискуссия. Я потерял часы, размышляя над проблемой, которой не существует. Хорошо, что теперь я никогда не буду воспринимать ТАК серьезно. Это отстой. - person kludg; 09.03.2013
comment
@user246408 user246408 Я действительно этого не понимаю, но на всякий случай я добавил новый проект с другим подходом к тестированию. Он быстро показывает несостоятельность утверждения для 64-битной версии на моей машине, на случай, если вы все еще заинтересованы. - person jachguate; 09.03.2013