Ответ на »сообщение 2394« (Trurl)
___________________________

Возможно, люди забыли разницу между абсолютной и относительной погрешностью, ка это было в примере с корнем.

Именно этот факт я и назвал в »сообщение 1815« "падением культуры вычислений".
(Я думал, это очевидно.)
Можно ошибиться в деталях, забыть какие-то мелочи, но не различие между абсолютной и относительной погрешностью.
Проблема же состоит в том, что выбор "эпсилон" для расчетов плохо освещен в литературе (ИМХО).
Исправим ошибку:

float sqrt(float x)
{
    float y = 1, y0;
    do {
        y0 = y;
        y = (y0 + x / y0) * 0.5;
    } while (fabs(y-y0) >= EPS*y);
    return y;
}

Вопрос, как выбрать значение EPS?
(Следует учесть, что sqrt -- библиотечная функция.)

Ответ на »сообщение 2386« (AVC)
_____________________


Все вычисления над значениями типа REAL являются приближенными. Но некоторые вычисления приближеннее остальных. (c) :)

Точно!!! Это правило объясняет все возникшие проблемы :)

Почему-то сравнение с eps не помогает, а вот 100*eps проходит (???)

PROCEDURE Do2*;
VAR x, dt: REAL;
BEGIN
  dt :=0.1;
  x := 10/dt;
  ASSERT(dt=0.1);
  ASSERT(ABS(x-10/dt)<Math.Eps()) (* TRAP 0 *)
END Do2;

Ответ на »сообщение 2392« (AVC)
___________________________
Допустим, машинный эпсилон известен (Math.Eps() в ББ, DBL_EPSILON в Си etc).
Что дальше?
Я недавно приводил несколько странный пример (»сообщение 2319«) с выбором эпсилон для вычислений. Я так и не понял, почему там выбран множитель 100...
Возможно, люди забыли разницу между абсолютной и относительной погрешностью, ка это было в примере с корнем.

Ответ на »сообщение 2390« (Антон Григорьев)
___________________________

Цель же моего вопроса вот какая. Я сильно сомневаюсь, что таким очевидным вопросом до меня никто не задавался. И мне просто хотелось бы знать - есть ли какие-то работы самого Вирта или других разработчиков Оберона, где затрагивалась бы эта проблема? Наверняка они уже продвинулись в этом вопросе дальше, чем "Арифметические операции для чисел типа REAL выполняются приближенно", только вот я об этом пока не знаю. Короче, кинтье ссылку, кто знает.

Пока могу сказать совсем немного.
В "Compiler construction" Вирт говорит (совсем немного) о форматах IEEE для чисел с плавающей точкой.
В описании компилятора Oberon-SA, кажется, приводится реализация операций с плавающей точкой в виде библиотечного модуля на Обероне.
В недавней (2004) книге "Programming in Oberon" Вирт приводит несколько примеров вычислений с плавающей точкой, в частности -- суммы первых N чисел гармонического ряда. Там отмечается, что суммировать надо от меньших членов ряда к большим:
According to the rules of arithmetic, the two sums ought to be equal. However, if values are truncated (or even rounded), the sums will differ for sufficiently large n. The correct way is evidently to start with the small terms.

Но, AFAIK, Вирт нигде работу с плавающей точкой подробно не разбирает, полагая, вероятно, что этот вопрос относится не к Оберону и даже не к базовому курсу программирования.

Хочу заметить, что этот вопрос вообще мало освещается в учебниках программирования, даже "продвинутых" и "специализированных".
Например, в известной книге "Numerical recipes in C" говорится:
Pretty much any arithmetic operation among floating numbers should be thought of as introducing an additional fractional error of at least epsm This type of error is called roundoff error.
epsm -- это как раз машинный эпсилон.
Меня здесь не устраивает выражение at least.
Из at least кашу не сваришь.
Вот если бы at most, другое дело... :)

Ответ на »сообщение 2391« (Комбриг)
___________________________

Еще в семидесятых годах я видел программы на Фортране, которые принципиально были расчитаны на то, что могут выполняться на самых разных машинах и с самой разной разрядностью. И решалась эта проблема предельно просто - прежде всего, вычислялось "машинное эпсилон" и все округления и проверки делались с ним...

Ну, наконец-то! Есть все-таки знающие люди!
Допустим, машинный эпсилон известен (Math.Eps() в ББ, DBL_EPSILON в Си etc).
Что дальше?
Я недавно приводил несколько странный пример (»сообщение 2319«) с выбором эпсилон для вычислений. Я так и не понял, почему там выбран множитель 100...

Еще в семидесятых годах я видел программы на Фортране, которые принципиально были расчитаны на то, что могут выполняться на самых разных машинах и с самой разной разрядностью. И решалась эта проблема предельно просто - прежде всего, вычислялось "машинное эпсилон" и все округления и проверки делались с ним...

Ответ на »сообщение 2389« (AVC)
___________________________

Здесь причина - непредставимость 0.1 в виде конечной двоичной дроби.

Причина-то как раз понятна - в своё время я немного поразвлекался с такими примерами - см. http://www.delphikingdom.com/asp/viewitem.asp?catalogid=374

Мои сомнения вы поняли правильно: я не вижу способа объяснить это замкнутой математической моделью, не прибегая к объяснению того, как эта арифметика реализована.

Цель же моего вопроса вот какая. Я сильно сомневаюсь, что таким очевидным вопросом до меня никто не задавался. И мне просто хотелось бы знать - есть ли какие-то работы самого Вирта или других разработчиков Оберона, где затрагивалась бы эта проблема? Наверняка они уже продвинулись в этом вопросе дальше, чем "Арифметические операции для чисел типа REAL выполняются приближенно", только вот я об этом пока не знаю. Короче, кинтье ссылку, кто знает.

Ответ на »сообщение 2387« (Антон Григорьев)
___________________________


Вот простейший пример:

var
  I:Integer;
  S:Single;
begin
  S:=1.0;
  for I:=1 to 10 do
    S:=S-0.1;
  ASSERT(S=0);  // Здесь возникнет исключение
end;


Здесь причина - непредставимость 0.1 в виде конечной двоичной дроби.
Интересно, что некоторые программисты не сразу верят в этот факт. :)
(У меня был опыт объяснения возникшей из-за этого ошибки.)

Более известная причина - конечность числа разрядов мантиссы.
Классический пример - конечность суммы единиц.

PROCEDURE Sum* ;
VAR s, s0: SHORTREAL;
BEGIN
s := 0.0;
REPEAT
s0 := s;
s := s + 1;
UNTIL s = s0;
Log.Real(s); Log.Ln
END Sum;

Сумма бесконечного ряда равна в данном случае равна 16777216.0. :)

Ответ на »сообщение 2387« (Антон Григорьев)
___________________________

Итак, имеем тип Single, для которого гарантируется точность не менее 7 десятичных цифр. Имеем два числа: 1.0 и 0.1, каждое из которых влезает в отведённую точность. Как, исходя из этих условий, объяснить, почему S не будет равно нулю? Пока мы не объясним, что такое число с плавающей двоичной точкой - никак.

Поэтому я и не могу согласиться ни с точкой зрения уважаемого info21, ни с точкой зрения уважаемого Трурля (»сообщение 1822«).
ИМХО, знания численных методов здесь недостаточно, надо представлять особенности используемой реализации плавающей арифметики.
В конце концов, FPU представляет собой компромисс между точностью и быстродействием, причем иногда сильно за счет точности. :(

Ответ на »сообщение 2385« (info21)
___________________________

А если вот так:
"Арифметические операции для чисел типа REAL выполняются приближенно."
:))


А если вот так: "Тип integer имеет верхнюю и нижнюю границу" без указания, какой должны быть эти границы? Может, математически и можно построить на этом непротиворечивую модель, но будет ли удобно пользоваться таким языком? Про переносимость я вообще молчу. Но с границами integer'ов всё просто - мы их можем просто полстулировать, не заостряя внимания на том, что именно такие значения определяются аппаратной реализацией. А с вещественной арифметикой вы так просто не разберётесь. И назвать границу точности, не переходя к аппаратной реализации, здесь невозможно.

Вот простейший пример:

var
  I:Integer;
  S:Single;
begin
  S:=1.0;
  for I:=1 to 10 do
    S:=S-0.1;
  ASSERT(S=0);  // Здесь возникнет исключение
end;

Итак, имеем тип Single, для которого гарантируется точность не менее 7 десятичных цифр. Имеем два числа: 1.0 и 0.1, каждое из которых влезает в отведённую точность. Как, исходя из этих условий, объяснить, почему S не будет равно нулю? Пока мы не объясним, что такое число с плавающей двоичной точкой - никак.