Неприводимость сравнения с плавающей точкой

Question

Я и мой доктор философииСтудент столкнулся с проблемой в контексте анализа физических данных, о которой я мог бы рассказать.У нас есть код, который анализирует данные одного из экспериментов на LHC и дает невоспроизводимые результаты.В частности, результаты вычислений, полученных из того же двоичного файла , запущенного на той же машине , могут различаться при последовательном выполнении.Мы знаем о многих различных источниках невоспроизводимости, но исключили обычных подозреваемых.

Мы отследили проблему до невоспроизводимости операций сравнения с плавающей запятой (двойной точности) при сравнении двух чисел, которые номинально имеют одинаковое значение.Это может происходить иногда в результате предыдущих шагов в анализе.Пример Мы только что нашли пример, который проверяет, является ли число меньше 0,3 (обратите внимание, что мы НИКОГДА не проверяем на равенство между плавающими значениями).Оказывается, что из-за геометрии детектора, расчет мог иногда давать результат, который был бы точно 0,3 (или его ближайшее представление с двойной точностью).

Мы хорошо осведомлены о подводных камнях при сравнении чисел с плавающей запятой, а также о возможности чрезмерной точности в FPU, чтобы повлиять на результаты сравнения.На вопрос, на который я хотел бы ответить, «почему результаты невоспроизводимы?»Это потому, что загрузка регистра FPU или другие инструкции FPU не очищают лишние биты и, таким образом, «оставшиеся» биты из предыдущих вычислений влияют на результаты?(это кажется маловероятным) На другом форуме я видел предложение о том, что переключение контекста между процессами или потоками может также вызвать изменение результатов сравнения с плавающей запятой из-за того, что содержимое FPU хранится в стеке и, следовательно, усекается.Будем благодарны за любые комментарии к этим = или другим возможным объяснениям.

Answer 1

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

Однако, когда происходит переключение контекста, состояние FPU должно быть сохранено и восстановлено - и есть, по крайней мере, довольно большой шанс, что эти дополнительные биты не сохраняются и восстанавливаются впереключение контекста.Когда это произойдет, это, вероятно, не приведет к серьезным изменениям, но если (например) вы позже вычтете фиксированную сумму из каждого и умножите то, что осталось, разница также будет умножена.

Чтобы бытьЯсно: я сомневаюсь, что «оставшиеся» биты были бы виновником.Скорее, это будет потеря лишних битов, вызывающая округление в немного разных точках вычисления.

Answer 2

Какая платформа?

Большинство FPU могут внутренне хранить больше точности, чем двойное представление ieee - чтобы избежать ошибки округления в промежуточных результатах.Часто компилятор переключается на скорость / точность торговли - см. http://msdn.microsoft.com/en-us/library/e7s85ffb(VS.80).aspx

Answer 3

Является ли программа многопоточной?

Если да, я бы заподозрил состояние гонки.

Если нет, выполнение программы детерминировано. Наиболее вероятная причина получения разных результатов при одинаковых входных данных - неопределенное поведение, то есть ошибка в вашей программе. Чтение неинициализированной переменной, устаревшего указателя, перезапись младших битов некоторого числа FP в стеке и т. Д. Возможности безграничны. Если вы работаете в Linux, попробуйте запустить его под valgrind и посмотрите, не обнаружит ли он некоторые ошибки.

Кстати, как вы сузили проблему до сравнения FP?

(Долгосрочный вариант: аппаратное сбой? Например, сбой чипа ОЗУ может привести к тому, что данные будут считываться по-разному в разное время. Хотя это, вероятно, приведет к быстрому падению ОС.)

Любое другое объяснение неправдоподобно - ошибки в ОС или в HW не остались бы незамеченными надолго.

Answer 4

Я сделал это:

#include <stdio.h>
#include <stdlib.h>

typedef long double ldbl;

ldbl x[1<<20];

void hexdump( void* p, int N ) {
  for( int i=0; i<N; i++ ) printf( "%02X", ((unsigned char*)p)[i] );
}

int main( int argc, char** argv ) {

  printf( "sizeof(long double)=%i\n", sizeof(ldbl) );

  if( argc<2 ) return 1;

  int i;
  ldbl a = ldbl(1)/atoi(argv[1]);

  for( i=0; i<sizeof(x)/sizeof(x[0]); i++ ) x[i]=a;

  while(1) {
    for( i=0; i<sizeof(x)/sizeof(x[0]); i++ ) if( x[i]!=a ) {
      hexdump( &a, sizeof(a) );
      printf( " " );
      hexdump( &x[i], sizeof(x[i]) );
      printf( "\n" );
    }
  }

}

, скомпилированный с IntelC с использованием / Qlong_double, так что он произвел это:

;;;     for( i=0; i<sizeof(x)/sizeof(x[0]); i++ ) if( x[i]!=a ) {

        xor       ebx, ebx                                      ;25.10
                                ; LOE ebx f1
.B1.9:                          ; Preds .B1.19 .B1.8
        mov       esi, ebx                                      ;25.47
        shl       esi, 4                                        ;25.47
        fld       TBYTE PTR [?x@@3PA_TA+esi]                    ;25.51
        fucomp                                                  ;25.57
        fnstsw    ax                                            ;25.57
        sahf                                                    ;25.57
        jp        .B1.10        ; Prob 0%                       ;25.57
        je        .B1.19        ; Prob 79%                      ;25.57
[...]
.B1.19:                         ; Preds .B1.18 .B1.9
        inc       ebx                                           ;25.41
        cmp       ebx, 1048576                                  ;25.17
        jb        .B1.9         ; Prob 82%                      ;25.17

и запустил 10 экземпляров с разными "семенами".Как вы можете видеть, он сравнивает 10-байтовые длинные дубликаты из памяти с одним в стеке FPU, поэтому в случае, когда ОС не сохраняет полную точность, мы наверняка увидим ошибку.И они по-прежнему работают, ничего не обнаруживая ... что неудивительно, поскольку в x86 есть команды для одновременного сохранения / восстановления всего состояния FPU, и в любом случае ОС, которая не сохранит полную точность, будет полностью сломана..

Таким образом, либо какая-то уникальная ОС / процессор / компилятор, либо отличающиеся результаты сравнения создаются после изменения чего-либо в программе и перекомпиляции, или из-за ошибки в программе, например.переполнение буфера.

Answer 5

Внутренний FPU ЦП может хранить числа с плавающей запятой с большей точностью, чем double или float.Эти значения должны быть преобразованы всякий раз, когда значения в регистре должны храниться где-то еще, в том числе, когда память выгружается в кеш (это я знаю точно), и переключение контекста или прерывание ОС на этом ядре звучит как другой простой источник,Конечно, синхронизация прерываний ОС или переключений контекста или замена не горячей памяти совершенно непредсказуема и неуправляема приложением.

Конечно, это зависит от платформы, но ваше описание звучит так, как будто вы работаете насовременный рабочий стол или сервер (т. е. x86).

Answer 6

Вы, безусловно, нажимаете Ошибка GCC n ° 323 , что, как и другие признаки, связано с чрезмерной точностью FPU.

Решения:

• Использование SSE (или AVX, это 2016 ...) для выполнения вычислений
• Использование переключателя компиляции -ffloat-store.Из GCC документы.

Не хранить переменные с плавающей запятой в регистрах и запрещать другие параметры, которые могут изменить то, берется ли значение с плавающей запятой из регистра или памяти.
Эта опция предотвращает нежелательныеизбыточная точность на машинах, таких как 68000, где плавающие регистры (из 68881) сохраняют большую точность, чем предполагалось в двойном.Аналогично для архитектуры x86.Для большинства программ избыточная точность приносит только пользу, но некоторые программы полагаются на точное определение плавающей запятой IEEE.Используйте -ffloat-store для таких программ, изменив их, чтобы сохранить все соответствующие промежуточные вычисления в переменных.

Answer 7

Я просто объединю некоторые комментарии Дэвида Родригеса и Бо Перссона и сделаю дикое предположение.

Может ли это быть переключением задач при использовании инструкций SSE3? Основываясь на этой статье Intel об использовании инструкций SSE3 , команды для сохранения статуса регистров FSAVE и FRESTOR были заменены FXSAVE и FXRESTOR, которые должны обрабатывать полная длина аккумулятора.

Я полагаю, что на компьютере x64 «неправильная» инструкция может содержаться в некоторой внешней скомпилированной библиотеке.