ошибка точности gcc?

Question

Я могу только предположить, что это ошибка. Первое утверждение проходит, а второе не проходит:

double sum_1 =  4.0 + 6.3;
assert(sum_1 == 4.0 + 6.3);

double t1 = 4.0, t2 = 6.3;

double sum_2 =  t1 + t2;
assert(sum_2 == t1 + t2);

Если не ошибка, то почему?

Answer 1

Вы сравниваете числа с плавающей запятой. Не делайте этого, числа с плавающей запятой имеют присущие ошибки точности в некоторых обстоятельствах. Вместо этого возьмите абсолютное значение разности двух значений и подтвердите, что это значение меньше некоторого небольшого числа (эпсилон).

void CompareFloats( double d1, double d2, double epsilon )
{
    assert( abs( d1 - d2 ) < epsilon );
} 

Это не имеет ничего общего с компилятором, и все, что связано со способом реализации чисел с плавающей запятой. вот спецификация IEEE:

http://www.eecs.berkeley.edu/~wkahan/ieee754status/IEEE754.PDF

Answer 2

Это то, что меня укусило тоже.

Да, числа с плавающей точкой никогда не должны сравниваться на равенство из-за ошибки округления, и вы, вероятно, знали это.

Но в этом случае вы вычисляете t1+t2, а затем вычисляете его снова. Неужели что должно дать идентичный результат?

Вот что, вероятно, происходит. Держу пари, что вы используете это на процессоре x86, правильно? FPU x86 использует 80 бит для своих внутренних регистров, но значения в памяти хранятся как 64-битные двойные числа.

Итак, t1+t2 сначала вычисляется с точностью 80 битов, затем - я полагаю - сохраняется в памяти в sum_2 с точностью 64 бита - и происходит некоторое округление. Для утверждения он загружается обратно в регистр с плавающей запятой, и t1+t2 вычисляется снова, снова с 80 битами точности. Итак, теперь вы сравниваете sum_2, который ранее был округлен до 64-битного значения с плавающей запятой, с t1+t2, который был вычислен с более высокой точностью (80 бит) - и поэтому значения не совсем идентичны .

Редактировать Так почему первый тест проходит? В этом случае компилятор, вероятно, оценивает 4.0+6.3 во время компиляции и сохраняет его как 64-битную величину - как для назначения, так и для утверждения. Таким образом, сравниваются идентичные значения, и утверждение проходит.

Второе редактирование Вот код сборки, сгенерированный для второй части кода (gcc, x86), с комментариями - в значительной степени соответствует сценарию, описанному выше:

// t1 = 4.0
fldl    LC3
fstpl   -16(%ebp)

// t2 = 6.3
fldl    LC4
fstpl   -24(%ebp)

// sum_2 =  t1+t2
fldl    -16(%ebp)
faddl   -24(%ebp)
fstpl   -32(%ebp)

// Compute t1+t2 again
fldl    -16(%ebp)
faddl   -24(%ebp)

// Load sum_2 from memory and compare
fldl    -32(%ebp)
fxch    %st(1)
fucompp

Интересное примечание: это было скомпилировано без оптимизации. Когда он скомпилирован с -O3, компилятор оптимизирует весь кода.

Answer 3

При сравнении чисел с плавающей запятой для близости вы обычно хотите измерить их относительную разницу, которая определяется как

if (abs(x) != 0 || abs(y) != 0) 
   rel_diff (x, y) = abs((x - y) / max(abs(x),abs(y))
else
   rel_diff(x,y) = max(abs(x),abs(y))

Например,

rel_diff(1.12345, 1.12367)   = 0.000195787019
rel_diff(112345.0, 112367.0) = 0.000195787019
rel_diff(112345E100, 112367E100) = 0.000195787019

Идея состоит в том, чтобы измерить количество старших значащих цифр, которые имеют общие числа; если вы возьмете -log10 0,000195787019, вы получите 3,70821611, что соответствует числу старших десятичных цифр во всех примерах.

Если вам нужно определить, равны ли два числа с плавающей запятой, вы должны сделать что-то вроде

if (rel_diff(x,y) < error_factor * machine_epsilon()) then
  print "equal\n";

где машинный эпсилон - это наименьшее число, которое может содержаться в мантиссе используемого аппаратного обеспечения с плавающей запятой. Большинство компьютерных языков имеют вызов функции, чтобы получить это значение. error_factor должен основываться на количестве значащих цифр, которые, по вашему мнению, будут использованы при округлении ошибок (и других) при вычислениях чисел x и y. Например, если бы я знал, что x и y были результатом около 1000 суммирований, и не знал никаких границ суммируемых чисел, я бы установил error_factor равным примерно 100.

Пытался добавить их как ссылки, но не смог, так как это мой первый пост:

• en.wikipedia.org / вики / Relative_difference
• en.wikipedia.org / вики / Machine_epsilon
• ru.wikipedia.org / wiki / Significand (мантисса)
• en.wikipedia.org / вики / Rounding_error

Answer 4

Я продублировал вашу проблему на моем Intel Core 2 Duo и посмотрел код сборки. Вот что происходит: когда ваш компилятор оценивает t1 + t2, он делает

load t1 into an 80-bit register
load t2 into an 80-bit register
compute the 80-bit sum

Когда он сохраняется в sum_2, он делает

round the 80-bit sum to a 64-bit number and store it

Затем сравнение == сравнивает 80-битную сумму с 64-битной суммой, и они отличаются, в основном потому, что дробная часть 0.3 не может быть точно представлена ​​с помощью двоичного числа с плавающей запятой, поэтому вы сравниваете «повторяющийся десятичный» (фактически повторяющийся двоичный), который был усечен до двух разных длин.

Что действительно раздражает, так это то, что если вы компилируете с gcc -O1 или gcc -O2, gcc сделает неправильную арифметику во время компиляции, и проблема исчезнет. Может быть, это нормально в соответствии со стандартом, но это еще одна причина, по которой gcc не мой любимый компилятор.

---

P.S. Когда я говорю, что == сравнивает 80-битную сумму с 64-битной, конечно, я действительно имею в виду, что она сравнивает расширенную версию 64-битной суммы. Вы могли бы хорошо подумать

sum_2 == t1 + t2

разрешается до

extend(sum_2) == extend(t1) + extend(t2)

и

sum_2 = t1 + t2

разрешается до

sum_2 = round(extend(t1) + extend(t2))

Добро пожаловать в удивительный мир с плавающей точкой!

Answer 5

Может случиться так, что в одном из случаев вы закончите сравнение 64-битного двойного с 80-битным внутренним регистром Может быть поучительно взглянуть на инструкции по сборке, которые GCC испускает для двух случаев ...

Answer 6

Сравнения чисел двойной точности по своей сути неточны. Например, вы часто можете найти 0.0 == 0.0, возвращающий false . Это связано с тем, как FPU хранит и отслеживает номера.

Википедия говорит :

Проверка на равенство проблематична. Две вычислительные последовательности, которые математически равны, могут давать разные значения с плавающей запятой.

Вам потребуется использовать дельту, чтобы дать допуск для ваших сравнений, а не точное значение.

Answer 7

Эту «проблему» можно «исправить» с помощью следующих опций:

-msse2 -mfpmath = sse

как объяснено на этой странице:

http://www.network -theory.co.uk / документы / gccintro / gccintro_70.html

Как только я использовал эти опции, оба утверждения прошли.