LINUX.ORG.RU

Как равнить float/double с максимальной точностью и быстро? ::)

 , , , ,


0

4

UDP: не сравнивать float с double, а float c float и double c double. Мухи отдельно, котлеты отдельно.

Пишу маленькую, ущербную, но гордую штуку дабы тестировать «С» код, ну и подумалось что порой надо проверять выхлоп функций float/double c точностью до последнего бита.
Ну для тех кто не в курсе ::)

float a = 0.000001;
float b = 0.000001;
if(a == b)
{
   puts("Радость\n");
}else{
   puts("Мааам они чё тупые вообще чтоль маааам\n");
}
Я помню как отреагировал на подобное когда узнал про такое «сравнение», вернее не понимал что происходит ::)

Как железно сравнить float/double без math.h? будет ли такое, эмм нормальным?

if(0.0 == (a - b))
{
   puts("Гип-гип, все флоаты братья!\n");
}else{
   puts("Хммм, в наших рядах таиться враг!Эркюль Пуаро идёт на помощь!");
};

Или просто побитово сравнить? Но что-то это слишком топорно вроде

Ну и по традиции cast сишников i-rinat, beastie, ncrmnt, Iron_Bug. Здрасти снова :D Ну смысл поста ещё в том что может какие ньюансы есть компиляторов или ещё чего, в том плане вдруг я не в курсе, а есть палки какие?

Deleted

Ответ на: комментарий от telikan

Да я их почти всегда кастую, люди проверенные годами так сказать ) Хотя возможно это их раздражает ::)

Deleted ()

if(a == b)

Сравнивать float не зная допустимой погрешности - дурная идея. Зачем оно тебе вообще понадобилось?

peregrine ★★★★★ ()
Последнее исправление: peregrine (всего исправлений: 1)
Ответ на: комментарий от peregrine

Зачем оно тебе вообще понадобилось?

Спортивный интерес, как сказал выше пишу юнит тестирование и хочу такое

//я не знаю что выплюнет get_float(), но хочу точного сравнения
//если calculate_float() выдаст 0.63563546540001 то ошибка
compare_float(0.6356354654,calculate_float());
//или/и диапазоны точности, не более 0.0010 и не менее 0.00001
//тоесть любое число в этих пределах
range_float(0.0010,0.000001,calculate_float());

Сравнивать float не зная допустимой погрешности - дурная идея

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

Deleted ()

Что не так с первым куском кода? Он сделает побитово одинаковые значения и всё будет хорошо.

В целом для сравнения плавающих значений можно взять какой-нибудь epsilon и проверять что разница значений меньше чем он.

будет ли такое, эмм нормальным?

Нет, это ничем не отличается от первого варианта. Если медленный вариант, то if( fabs(a - b) < eps )

а есть палки какие?

Ну, например, на разных платформах разный float, и точность везде нужна разная.

alexanius ★★ ()

Никак. У float мантисса 23 бита и 8 бит на экспоненту, у double 52 и 11 бит на экспоненту.

sign = 0
exponent = 123
mantissa = 8181376
sign = 0
exponent = 1019
mantissa = 4392342702175935
sign = 0
exponent = 1019
mantissa = 4392342794534912
#include <stdio.h>
#include <stdint.h>

typedef union {
  float f;
  struct {
    uint32_t mantissa : 23;
    uint16_t exponent : 8;
    uint8_t sign : 1;
  } data;
} float_data;

typedef union {
    double d;
    struct {
        uint64_t mantissa : 52;
        uint16_t exponent : 11;
        uint8_t sign : 1;
    } data;
} double_data;

int main(void) {
  float f1 = 0.123456;
  double d1 = 0.123456;
  float_data d = { .f = f1 };
  printf("sign = %u\n", d.data.sign);
  printf("exponent = %u\n", d.data.exponent);
  printf("mantissa = %u\n", d.data.mantissa);

  double_data x = { .d = d1 };
  printf("sign = %u\n", x.data.sign);
  printf("exponent = %u\n", x.data.exponent);
  printf("mantissa = %llu\n", x.data.mantissa);

  double_data y = { .d = f1 };
  printf("sign = %u\n", y.data.sign);
  printf("exponent = %u\n", y.data.exponent);
  printf("mantissa = %llu\n", y.data.mantissa);

}
xpahos ★★★★★ ()
Ответ на: комментарий от peregrine

Блин, да ты float с double сверять собрался

Неее, флоаты отдельно, дубли отдельно. Всё по отдельности в том числе и инты :)

Deleted ()
Ответ на: комментарий от xpahos

Извиняюсь,ввёл в заблуждение, я неточно выразился float и double отдельно, обновил шапку поста.

Deleted ()
Ответ на: комментарий от Deleted

Для слушателей ночного радио «Разработка на ЛОРе» даю специальный мастер-класс по сокращению примеров. Здесь суть проблемы такова. Переменная а имеет тип float, и туда записана константа соответствующего типа. А значение с которым она сравнивается интерпретируется как double. Из ассемблера мы можем видеть что инструкция cvtss2sd конвертирует значение из float в double, где и происходит потеря точности с последующим, казалось бы, неправильно работающим сравнением.

В данном случае надо было писать if( a == 0.000001f ) return 0; и тогда было бы автору счастье.

PS. Разумеется, реальный пример не соответствует тому что было написано ТС в первом посте.

alexanius ★★ ()

Как уже сказали выше, с твоим первым куском кода проблем нет.

Как равнить float/double с максимальной точностью и быстро? ::)

Как ты и написал в первом посте: «a == b». Это и есть сравнение с максимальной точностью.

Проблема такого подхода не в том, что сравнение чисел с плавающей точкой каким-то магическим образом отличается от сравнения целых чисел (а оно не отличается, кроме случая с NaN), а в том, что в реальном коде у тебя значения переменных «a» и «b» перед сравнением будут вычислены разным способом. И если чисто математически (== при вычислении с абсолютной точностью) они должны быть равны, то из-за ограниченной точности чисел с плавающей точкой у тебя в коде они неизбежно будут отличаться. Насколько отличаться - это другой вопрос, ответ на который зависит от конкретного кода.

Deleted ()

Определиться какая точность нужна (сколько знаков после запятой), домножить на 10^x, привести к целым числам и их сравнивать?

Nirdosh ()
Последнее исправление: Nirdosh (всего исправлений: 1)
Ответ на: комментарий от Nirdosh

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

peregrine ★★★★★ ()
Ответ на: комментарий от Nirdosh

Определиться какая точность нужна (сколько знаков после запятой), домножить на 10^x, привести к целым числам и их сравнивать?

Идеально, если бы я знал какую точность выплёвывает функция, у меня могут сравниваться выхлопы двух функций к примеру, и обработать их выхлопы нужно автоматом :)

Deleted ()

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

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

Могу посоветовать эту «Вычислительные методы для инженеров» группы авторов: https://libgen.pw/item/detail/id/8858?id=8858 Там прочитай главу 2. «Введение в элементарную теорию погрешностей.»

Особенно параграф 2.5 «Особенности машинной арифметики»

Пока же могу посоветовать сравнивать с машинным эпсилон.

if (fabs(a-b)<=epsilon) { }

Машинное эпсилон - это такое число:

double epsilon=1;
  while((epsilon+1)!=1)
     epsilon=epsilon/2;

Для каждого типа эпсилон своё. У меня это получилось 1.110223e-16

anonymous_incognito ★★★★★ ()
Последнее исправление: anonymous_incognito (всего исправлений: 1)
Ответ на: комментарий от anonymous_incognito

Это имело бы смысл, если сравниваемые величины нормированы к 1. Но обычно это не так. Для сравнения уже тысяч эта точность вроообще недостижима, для сравнения малых одного порядка (скажем, 1e-20) может быть сильно избыточна.

bormant ★★★★★ ()
Ответ на: комментарий от bormant

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

anonymous_incognito ★★★★★ ()
Последнее исправление: anonymous_incognito (всего исправлений: 1)
Ответ на: комментарий от anonymous_incognito

Ну какой же адекватный, если он недостижимый? Выше считали его для 1, теперь нужно сравнивать числа порядка 1e6, пересчитайте epsilon для этого случая, оно будет примерно в 1e6 раз больше. А в другую сторону он наоборот слишком свободен.

Отчего б не взять за eps мантиссу 1<<грубость с порядком, равным минимальному порядку сравниваемых величин?

bormant ★★★★★ ()
Последнее исправление: bormant (всего исправлений: 4)

if( fabs(a - b) < eps )

Нуу, если ничего другого не придумаю то да с максимальным esp

В GoogleTest функция AlmostEquals используется для ASSERT_DOUBLE_EQ и ASSERT_FLOAT_EQ

// This template class represents an IEEE floating-point number
// (either single-precision or double-precision, depending on the
// template parameters).
//
// The purpose of this class is to do more sophisticated number
// comparison.  (Due to round-off error, etc, it's very unlikely that
// two floating-points will be equal exactly.  Hence a naive
// comparison by the == operation often doesn't work.)
//
// Format of IEEE floating-point:
//
//   The most-significant bit being the leftmost, an IEEE
//   floating-point looks like
//
//     sign_bit exponent_bits fraction_bits
//
//   Here, sign_bit is a single bit that designates the sign of the
//   number.
//
//   For float, there are 8 exponent bits and 23 fraction bits.
//
//   For double, there are 11 exponent bits and 52 fraction bits.
//
//   More details can be found at
//   http://en.wikipedia.org/wiki/IEEE_floating-point_standard.
//
// Template parameter:
//
//   RawType: the raw floating-point type (either float or double)
template <typename RawType>
class FloatingPoint {
 public:
  // Defines the unsigned integer type that has the same size as the
  // floating point number.
  typedef typename TypeWithSize<sizeof(RawType)>::UInt Bits;

  // Constants.

  // # of bits in a number.
  static const size_t kBitCount = 8*sizeof(RawType);

  // # of fraction bits in a number.
  static const size_t kFractionBitCount =
    std::numeric_limits<RawType>::digits - 1;

  // # of exponent bits in a number.
  static const size_t kExponentBitCount = kBitCount - 1 - kFractionBitCount;

  // The mask for the sign bit.
  static const Bits kSignBitMask = static_cast<Bits>(1) << (kBitCount - 1);

  // The mask for the fraction bits.
  static const Bits kFractionBitMask =
    ~static_cast<Bits>(0) >> (kExponentBitCount + 1);

  // The mask for the exponent bits.
  static const Bits kExponentBitMask = ~(kSignBitMask | kFractionBitMask);

  // How many ULP's (Units in the Last Place) we want to tolerate when
  // comparing two numbers.  The larger the value, the more error we
  // allow.  A 0 value means that two numbers must be exactly the same
  // to be considered equal.
  //
  // The maximum error of a single floating-point operation is 0.5
  // units in the last place.  On Intel CPU's, all floating-point
  // calculations are done with 80-bit precision, while double has 64
  // bits.  Therefore, 4 should be enough for ordinary use.
  //
  // See the following article for more details on ULP:
  // http://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
  static const size_t kMaxUlps = 4;

  // Constructs a FloatingPoint from a raw floating-point number.
  //
  // On an Intel CPU, passing a non-normalized NAN (Not a Number)
  // around may change its bits, although the new value is guaranteed
  // to be also a NAN.  Therefore, don't expect this constructor to
  // preserve the bits in x when x is a NAN.
  explicit FloatingPoint(const RawType& x) { u_.value_ = x; }

  // Static methods

  // Reinterprets a bit pattern as a floating-point number.
  //
  // This function is needed to test the AlmostEquals() method.
  static RawType ReinterpretBits(const Bits bits) {
    FloatingPoint fp(0);
    fp.u_.bits_ = bits;
    return fp.u_.value_;
  }

  // Returns the floating-point number that represent positive infinity.
  static RawType Infinity() {
    return ReinterpretBits(kExponentBitMask);
  }

  // Returns the maximum representable finite floating-point number.
  static RawType Max();

  // Non-static methods

  // Returns the bits that represents this number.
  const Bits &bits() const { return u_.bits_; }

  // Returns the exponent bits of this number.
  Bits exponent_bits() const { return kExponentBitMask & u_.bits_; }

  // Returns the fraction bits of this number.
  Bits fraction_bits() const { return kFractionBitMask & u_.bits_; }

  // Returns the sign bit of this number.
  Bits sign_bit() const { return kSignBitMask & u_.bits_; }

  // Returns true iff this is NAN (not a number).
  bool is_nan() const {
    // It's a NAN if the exponent bits are all ones and the fraction
    // bits are not entirely zeros.
    return (exponent_bits() == kExponentBitMask) && (fraction_bits() != 0);
  }

  // Returns true iff this number is at most kMaxUlps ULP's away from
  // rhs.  In particular, this function:
  //
  //   - returns false if either number is (or both are) NAN.
  //   - treats really large numbers as almost equal to infinity.
  //   - thinks +0.0 and -0.0 are 0 DLP's apart.
  bool AlmostEquals(const FloatingPoint& rhs) const {
    // The IEEE standard says that any comparison operation involving
    // a NAN must return false.
    if (is_nan() || rhs.is_nan()) return false;

    return DistanceBetweenSignAndMagnitudeNumbers(u_.bits_, rhs.u_.bits_)
        <= kMaxUlps;
  }

 private:
  // The data type used to store the actual floating-point number.
  union FloatingPointUnion {
    RawType value_;  // The raw floating-point number.
    Bits bits_;      // The bits that represent the number.
  };

fsb4000 ★★★★★ ()
Ответ на: комментарий от fsb4000
 // Given two numbers in the sign-and-magnitude representation,
  // returns the distance between them as an unsigned number.
  static Bits DistanceBetweenSignAndMagnitudeNumbers(const Bits &sam1,
                                                     const Bits &sam2) {
    const Bits biased1 = SignAndMagnitudeToBiased(sam1);
    const Bits biased2 = SignAndMagnitudeToBiased(sam2);
    return (biased1 >= biased2) ? (biased1 - biased2) : (biased2 - biased1);
  }

  // Converts an integer from the sign-and-magnitude representation to
  // the biased representation.  More precisely, let N be 2 to the
  // power of (kBitCount - 1), an integer x is represented by the
  // unsigned number x + N.
  //
  // For instance,
  //
  //   -N + 1 (the most negative number representable using
  //          sign-and-magnitude) is represented by 1;
  //   0      is represented by N; and
  //   N - 1  (the biggest number representable using
  //          sign-and-magnitude) is represented by 2N - 1.
  //
  // Read http://en.wikipedia.org/wiki/Signed_number_representations
  // for more details on signed number representations.
  static Bits SignAndMagnitudeToBiased(const Bits &sam) {
    if (kSignBitMask & sam) {
      // sam represents a negative number.
      return ~sam + 1;
    } else {
      // sam represents a positive number.
      return kSignBitMask | sam;
    }
  }

fsb4000 ★★★★★ ()
Ответ на: комментарий от peregrine

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

Бывали времена, когда работать с большими целыми числами было быстрее, если они были таки double. Главное за их целый диаппазон не велезти. Та же PDP с двойными словами работала медленно и печально, а уж делить и умножать вообще можно было только через функции, а double работали через быстрый сопроцессор (если он был) и ещё и параллельно.

vodz ★★★★★ ()

1. Выкинуть floating-point, перейти на fixed-point

2. Выкинуть floating-point, перейти на fixed-point

3. Если все-таки руки шаловливые, почитать что-то типа https://www.python.org/dev/peps/pep-0485/ для прозрения, почему floating-point плохо. (Ответ: потому что сложно, и в 99% случаях бизнеслогика отлично ложится на fixed point, что избавляет от тонн геморроя). Ну и куча других ссылок про fp, типа https://rsdn.org/forum/cpp/2640596 или сам найди.

4. Выкинуть floating-point, перейти на fixed-point

PS. вопрос и тред читал по-диагонали. Сагрился на слово float :)

Deleted ()
Ответ на: комментарий от anonymous

Это float, ничего не поделаешь. Где-нибудь приходится искать компромисс.

beastie ★★★★★ ()
Ответ на: комментарий от anonymous

развели тут философию на пустом месте

Ну да, ну да. И сравнили только 32 бита на 64 битной тачке.

vodz ★★★★★ ()
Ответ на: комментарий от bormant

Это имело бы смысл, если сравниваемые величины нормированы к 1. Но обычно это не так. Для сравнения уже тысяч эта точность вроообще недостижима, для сравнения малых одного порядка (скажем, 1e-20) может быть сильно избыточна.

int feql(double x, double y) {
  return x == 0 ? y == 0 : fabs(y / x - 1) < epsilon;
}

Ноль --- специальный случай, и сравнивать с ним надо более аккуратно.

Jini ★★ ()
Ответ на: комментарий от Jini

Вопрос не в этом, вопрос в величине epsilon — для 1e20 и для 1e-20 оно не одинаковое. Попытка избавиться делением может легко привести к переполнению.

Тот же самый дефект на сравнении x или y как результата вычислений с нулем.

bormant ★★★★★ ()
Последнее исправление: bormant (всего исправлений: 2)

И до сих пор никто не посоветовал результат разницу двух чисел сравнивать с FLT_EPSILON или DBL_EPSILON:

int cmpd(double a, double b){
  double d = a - b;
  if(d > DBL_EPSILON) return 1;
  if(d < -DBL_EPSILON) return -1;
  return 0;
}

anonymous ()
Ответ на: комментарий от anonymous

И правильно сделал, что не посоветовал, потому, что это не работает для чисел порядка 1e16 и больше.

d просто будет равнятся нулю для double a = 2e16 и double b = 2e16 - 1, а (a == b), вернёт единичку.

Самое странное, что вопрос давно много где рассмотрен, но в интернете полно неверных решений :( Правильное я тоже не помню, нужно искать.

grem ★★★★★ ()
Последнее исправление: grem (всего исправлений: 1)

Если fpu/vfp не богаты, и именно точное равенство (без погрешности) - скастовать/запихнуть в union c uint64_t и сравнить их.

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

ncrmnt ★★★★★ ()
Последнее исправление: ncrmnt (всего исправлений: 2)
Ответ на: комментарий от ZERG

Машинный epsilon задаёт не абсолютную, а относительную погрешность. Абсолютную погрешность для произвольного float можно оценить как epsilon * abs(float). Ноль здесь является исключением, так как погрешность у него есть, но оценить её в рамках функции feql невозможно --- для этого надо знать как этот ноль получается в процессе расчёта. К примеру, если x = 0, а y --- результат вычитания двух чисел, близких к 1, то погрешность y будет порядка epsilon, но epsilon * abs(x) = 0, а epsilon * abs(y) много меньше epsilon.

Следует ещё иметь в виду, что редко когда нужно сравнивать сферические числа в вакууме; чаще это результат какого-то расчёта, и эти числа являются физическими величинами. В этом случае у них есть не только машинная погрешность, но и физическая, и они обе накапливаются в процессе вычислений. Бессмысленно класть epsilon равным 1e-20 если вы используете приближённую формулу с погрешностью на уровне 1%. Поэтому лучше epsilon сделать параметром и при вызове feql передавать в этом месте оценку погрешности.

Jini ★★ ()
Вы не можете добавлять комментарии в эту тему. Тема перемещена в архив.