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

Но что-то это слишком топорно вроде

Пойдет. Бывает и хуже.

Чего так мало сишников нашли?

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

if(a == b)

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

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

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

//я не знаю что выплюнет 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 не зная допустимой погрешности - дурная идея

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

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

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

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

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

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

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

Никак. У 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);

}

По подводным камням, можешь погуглить про catastrophic cancellation.

Блин, да ты float с double сверять собрался (или нет?), это ещё хуже.

Что не так с первым куском кода?

У меня 0.000001 != 0.000001 И я расстроен ::)

if( fabs(a - b) < eps )

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

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

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

А можно показать вот здесь полный участок примера? Потому что такое поведение... несколько нетипично.

тыцлистай в самый низ

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

Там ошибки, но не суть.

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

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

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

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

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

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

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

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

А значение с которым она сравнивается интерпретируется как double

Пасиба большое, приведение типов... лицо, рука, пойду покурю с горя ;

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

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

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

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

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

Могу посоветовать эту «Вычислительные методы для инженеров» группы авторов: 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

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

Спасибо, почитаю

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

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

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

https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers...

Quick'n'dirty:

int equal(float a, float b) {
    return (a-b) < 1e-6 && (b-a) < 1e-6
}

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.
  };

Да, то есть 1e-7 равно 9е-7

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

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

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 :)

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

*((uint*)(fSmth0)) == *((uint*)(fSmth1))

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

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

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

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

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

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

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

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

И до сих пор никто не посоветовал результат разницу двух чисел сравнивать с 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;
}

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

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

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

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

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

Правильное я тоже не помню, нужно искать.

Ну нате, 15 теорем и прочие следствия: https://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html

Самое смешное, что для зачади ТСа «сравнить до последнего бита» ничего кроме (a==b) не надо делать.

Нужна погрешность, say 1e-8.

http://wiki.seas.harvard.edu/geos-chem/index.php/Floating_point_math_issues#E...

https://randomascii.wordpress.com/category/floating-point/

https://randomascii.wordpress.com/2017/06/19/sometimes-floating-point-math-is...

Можно больше деталей про сравнение с нулём?

Машинный 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 передавать в этом месте оценку погрешности.