Porównaj podwójnie do zera za pomocą epsilon

Dzisiaj przeglądałem kod C ++ (napisany przez kogoś innego) i znalazłem tę sekcję:

double someValue = ...
if (someValue <  std::numeric_limits<double>::epsilon() && 
    someValue > -std::numeric_limits<double>::epsilon()) {
  someValue = 0.0;
}

Próbuję dowiedzieć się, czy to w ogóle ma sens.

Dokumentacja epsilon()mówi:

Funkcja zwraca różnicę między 1 a najmniejszą wartością większą niż 1, która jest reprezentowalna [podwójnie].

Czy dotyczy to również 0, tj epsilon(). Czy najmniejsza wartość jest większa od 0? Czy też są liczby pomiędzy 0i 0 + epsilonktóre mogą być reprezentowane przez double?

Jeśli nie, to czy porównanie nie jest równoważne someValue == 0.0?

Zakładając, że 64-bitowy podwójny IEEE ma 52-bitową mantysę i 11-bitowy wykładnik. Podzielmy to na kawałki:

1.0000 00000000 00000000 00000000 00000000 00000000 00000000 × 2^0 = 1

Najmniejsza reprezentowalna liczba większa niż 1:

1.0000 00000000 00000000 00000000 00000000 00000000 00000001 × 2^0 = 1 + 2^-52

W związku z tym:

epsilon = (1 + 2^-52) - 1 = 2^-52

Czy są jakieś liczby od 0 do epsilon? Dużo ... Np. Minimalna liczba reprezentatywna (normalna) to:

1.0000 00000000 00000000 00000000 00000000 00000000 00000000 × 2^-1022 = 2^-1022

W rzeczywistości istnieją (1022 - 52 + 1)×2^52 = 4372995238176751616liczby od 0 do epsilon, co stanowi 47% wszystkich dodatnich liczb reprezentowalnych ...

Test na pewno nie jest taki sam jak someValue == 0. Cała idea liczb zmiennoprzecinkowych polega na tym, że przechowują one wykładnik i znaczenie. Reprezentują zatem wartość z pewną liczbą binarnych znaczących liczb precyzji (53 w przypadku podwójnego IEEE). Reprezentatywne wartości są znacznie gęstiej upakowane w pobliżu 0 niż w pobliżu 1.

Aby użyć bardziej znanego systemu dziesiętnego, załóżmy, że przechowujesz wartość dziesiętną „do 4 cyfr znaczących” z wykładnikiem wykładniczym. Następnie następna reprezentowalna wartość większa niż 1jest 1.001 * 10^0i epsilonjest 1.000 * 10^-3. Ale 1.000 * 10^-4jest również reprezentowalny, zakładając, że wykładnik może przechowywać -4. Możesz wierzyć mi na słowo, że podwójne IEEE może przechowywać wykładniki mniejsze niż wykładnik epsilon.

Nie można stwierdzić na podstawie samego kodu, czy ma sens, czy nie należy go używać epsilonjako granicy, należy spojrzeć na kontekst. Może to być epsilonuzasadnione oszacowanie błędu w obliczeniach, które się pojawiło someValue, i może być tak, że nie jest.

Istnieją liczby, które istnieją między 0 a epsilon, ponieważ epsilon jest różnicą między 1 a następną najwyższą liczbą, którą można przedstawić powyżej 1, a nie różnicą między 0 a następną najwyższą liczbą, którą można przedstawić powyżej 0 (jeśli tak, to że kod zrobiłby bardzo mało): -

#include <limits>

int main ()
{
  struct Doubles
  {
      double one;
      double epsilon;
      double half_epsilon;
  } values;

  values.one = 1.0;
  values.epsilon = std::numeric_limits<double>::epsilon();
  values.half_epsilon = values.epsilon / 2.0;
}

Za pomocą debugera zatrzymaj program na końcu main i spójrz na wyniki, a zobaczysz, że epsilon / 2 różni się od epsilon, zero i jeden.

Ta funkcja przyjmuje więc wartości pomiędzy +/- epsilon i czyni je zerami.

Przybliżenie epsilon (najmniejsza możliwa różnica) wokół liczby (1,0, 0,0, ...) można wydrukować za pomocą następującego programu. Wyświetla następujące dane wyjściowe:
epsilon for 0.0 is 4.940656e-324
epsilon for 1.0 is 2.220446e-16
Trochę myślenia wyjaśnia, że ​​epsilon staje się mniejszy, im mniejszą liczbę używamy do sprawdzania jej wartości epsilon, ponieważ wykładnik może dostosować się do wielkości tej liczby.

#include <stdio.h>
#include <assert.h>
double getEps (double m) {
  double approx=1.0;
  double lastApprox=0.0;
  while (m+approx!=m) {
    lastApprox=approx;
    approx/=2.0;
  }
  assert (lastApprox!=0);
  return lastApprox;
}
int main () {
  printf ("epsilon for 0.0 is %e\n", getEps (0.0));
  printf ("epsilon for 1.0 is %e\n", getEps (1.0));
  return 0;
}

Załóżmy, że pracujemy z zabawkowymi liczbami zmiennoprzecinkowymi, które pasują do 16-bitowego rejestru. Jest bit znaku, 5-bitowy wykładnik i 10-bitowa mantysa.

Wartością tej liczby zmiennoprzecinkowej jest mantysa, interpretowana jako binarna wartość dziesiętna, razy dwa do potęgi wykładnika.

Około 1 wykładnik równa się zero. Najmniejsza cyfra mantysy to jedna część na 1024.

Blisko 1/2 wykładnika wynosi minus jeden, więc najmniejsza część mantysy jest o połowę mniejsza. Z pięciobitowym wykładnikiem może osiągnąć wartość ujemną 16, w którym to momencie najmniejsza część mantysy jest warta jedną część na 32m. Przy ujemnym wykładniku 16 wartość wynosi około jednej części na 32k, znacznie bliżej zera niż epsilon wokół tego, który obliczyliśmy powyżej!

Teraz jest to zabawkowy model zmiennoprzecinkowy, który nie odzwierciedla wszystkich dziwactw prawdziwego systemu zmiennoprzecinkowego, ale zdolność do odzwierciedlania wartości mniejszych niż epsilon jest dość podobna do rzeczywistych wartości zmiennoprzecinkowych.

Różnica między Xkolejną wartością a kolejną wartością Xróżni się w zależności od X.
epsilon()jest tylko różnicą między 1i następną wartością 1.
Różnica pomiędzy 0i następną wartością 0nie jest epsilon().

Zamiast tego możesz użyć std::nextafterdo porównania podwójnej wartości z 0następującymi:

bool same(double a, double b)
{
  return std::nextafter(a, std::numeric_limits<double>::lowest()) <= b
    && std::nextafter(a, std::numeric_limits<double>::max()) >= b;
}

double someValue = ...
if (same (someValue, 0.0)) {
  someValue = 0.0;
}

Myślę, że to zależy od precyzji twojego komputera. Spójrz na tę tabelę : możesz zobaczyć, że jeśli twój epsilon jest reprezentowany przez dwukrotność, ale twoja precyzja jest wyższa, porównanie nie jest równoważne z

someValue == 0.0

W każdym razie dobre pytanie!

Nie możesz zastosować tego do 0, z powodu części mantysy i wykładnika. Z powodu wykładnika możesz przechowywać bardzo małe liczby, które są mniejsze niż epsilon, ale gdy spróbujesz zrobić coś takiego (1.0 - „bardzo mała liczba”), otrzymasz 1.0. Epsilon jest wskaźnikiem nie wartości, ale precyzji wartości, która znajduje się w mantysie. Pokazuje, ile poprawnych kolejnych cyfr dziesiętnych liczby możemy zapisać.

Z zmiennoprzecinkową wartością IEEE, pomiędzy najmniejszą niezerową wartością dodatnią i najmniejszą niezerową wartością ujemną, istnieją dwie wartości: zero dodatnie i zero ujemne. Testowanie, czy wartość mieści się między najmniejszymi niezerowymi wartościami, jest równoważne testowaniu na równość z zerem; przypisanie może jednak mieć wpływ, ponieważ zmieniłoby ujemne zero na dodatnie zero.

Można sobie wyobrazić, że format zmiennoprzecinkowy może mieć trzy wartości między najmniejszymi skończonymi wartościami dodatnimi i ujemnymi: dodatni nieskończenie mały, zero bez znaku i ujemny nieskończenie mały. Nie znam żadnych formatów zmiennoprzecinkowych, które w rzeczywistości działałyby w ten sposób, ale takie zachowanie byłoby całkowicie rozsądne i prawdopodobnie lepsze niż IEEE (być może nie na tyle lepsze, aby warto było dodać dodatkowy sprzęt do obsługi, ale matematycznie 1 / (1 / INF), 1 / (- 1 / INF) i 1 / (1-1) powinny reprezentować trzy różne przypadki ilustrujące trzy różne zera). Nie wiem, czy jakikolwiek standard C nakazałby, aby podpisane nieskończenie małe, jeśli istnieją, musiałyby się równać zero. Jeśli nie, kod taki jak powyższy może pożytecznie zapewnić, że np

Powiedzmy, że system nie może odróżnić 1.000000000000000000000 i 1.000000000000000000001. to jest 1,0 i 1,0 + 1e-20. Czy uważasz, że nadal istnieją pewne wartości, które można przedstawić między -1e-20 a + 1e-20?

Ponadto dobrym powodem takiej funkcji jest usunięcie „denormałów” (te bardzo małe liczby, które nie mogą już używać domyślnego wiodącego „1” i mają specjalną reprezentację FP). Dlaczego chcesz to zrobić? Ponieważ niektóre maszyny (w szczególności niektóre starsze Pentium 4) działają naprawdę bardzo wolno podczas przetwarzania denormałów. Inne stają się nieco wolniejsze. Jeśli twoja aplikacja tak naprawdę nie potrzebuje tych bardzo małych liczb, dobrym rozwiązaniem jest spłukanie ich do zera. Dobrym miejscem do rozważenia tego są ostatnie kroki filtrów IIR lub funkcji rozpadu.

Zobacz także: Dlaczego zmiana 0,1f na 0 spowalnia działanie 10-krotnie?

i http://en.wikipedia.org/wiki/Denormal_number