Как проверить зависимости поплавков

Question

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

float x=5,y=0.2;
if(x==(1/y)) cout<<"They are the multiplicative inverse of eachother"<<endl;
else cout<<"They are NOT the multiplicative inverse of eachother"<<endl;

выведет: «они не ...», что неверно, и этот код:

float x=5,y=0.2,z;
z=1/y;
if(x==z) cout<<"They are the multiplicative inverse of eachother"<<endl;
else cout<<"They are NOT the multiplicative inverse of eachother"<<endl;

выведет: "они ...", что правильно.
Почему это происходит?

Answer 1

Проблема точности поплавка

У вас есть две проблемы здесь, но оба происходят из одного корня

Вы не можете сравнивать поплавки точно. Вы не можете вычесть или разделить их точно. Вы не можете рассчитывать ничего для них точно. Любая операция с ними может (и почти всегда) приводит к некоторой ошибке в результате. Даже a=0.2f не является точной операцией. Более глубокие причины этого очень хорошо объяснены авторами других ответов здесь. _{(Моя благодарность и голос за них.)}

Вот ваша первая и более простая ошибка. Никогда, никогда , никогда , никогда , НИКОГДА не используйте на них == или его эквивалент в любой язык.

Вместо a==b используйте Abs(a-b)<HighestPossibleError.

---

Но это не единственная проблема в вашей задаче.

Abs(1/y-x)<HighestPossibleError тоже не сработает. По крайней мере, это не будет работать достаточно часто. Зачем?

Давайте возьмем пару х = 1000 и у = 0,001. Давайте возьмем «начальную» относительную ошибку y для 10 ^-6 .

_{(Относительная ошибка = ошибка / значение).}

Относительные погрешности значений складываются при умножении и делении.

1 / y составляет около 1000. Его относительная ошибка такая же, как 10 ^-6 . («1» не содержит ошибок)

Это делает абсолютную ошибку = 1000 * 10 ^-6 = 0,001. Когда вы вычтете x позже, этой ошибкой будет все, что останется. (Абсолютные ошибки добавляются к при сложении и вычитании, а ошибка x пренебрежимо мала.) Конечно, вы не рассчитываете на такие большие ошибки, HighestPossibleError, несомненно, будет установлен ниже, и ваша программа выдаст хорошую пару x у

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

Существует два простых способа избежать этой проблемы.

• Найдя, что из x, у имеет большее значение abs и делит 1 на большее и только позже вычитает меньшее.

• Если вы хотите сравнить 1/y against x, пока вы еще работаете с буквами, а не со значениями, и ваши операции не дают ошибок , умножьте обе стороны сравнения на y и у вас есть 1 against x*y. _{(Обычно вы должны проверять знаки в этой операции, но здесь мы используем значения abs, поэтому он чистый.)} Сравнение результатов вообще не имеет деления.

Короче:

1/y V x   <=>   y*(1/y) V x*y   <=>   1 V x*y 

Мы уже знаем, что такое сравнение как 1 against x*y должно быть сделано так:

const float HighestPossibleError=1e-10;
if(Abs(x*y-1.0)<HighestPossibleError){...

Вот и все.

---

P.S. Если вам действительно нужно все в одной строке, используйте:

if(Abs(x*y-1.0)<1e-10){...

Но это плохой стиль. Я бы не советовал.

P.P.S. Во втором примере компилятор оптимизирует код так, чтобы он установил z на 5, прежде чем запускать какой-либо код. Таким образом, проверка 5 против 5 работает даже для поплавков.

Answer 2

Проблема в том, что 0.2 нельзя представить точно в двоичном виде, поскольку его двоичное расширение имеет бесконечное число цифр:

 1/5: 0.0011001100110011001100110011001100110011...

Это похоже на то, как 1/3 не может быть представлено точно вдесятичный.Поскольку x хранится в float, который имеет конечное число битов, в какой-то момент эти цифры будут обрезаны, например:

   x: 0.0011001100110011001100110011001

Проблема возникает из-за того, что процессоры часто используют более высокиевнутренняя точность, поэтому, когда вы только что вычислили 1/y, результат будет иметь больше цифр, а когда вы загрузите x для их сравнения, x будет расширен, чтобы соответствовать внутренней точности ЦП.

 1/y: 0.0011001100110011001100110011001100110011001100110011
   x: 0.0011001100110011001100110011001000000000000000000000

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

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

   x: 0.0011001100110011001100110011001
   z: 0.0011001100110011001100110011001

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

Answer 3

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

0.2 не имеет точного двоичного представления.Если вы храните числа, которые не имеют точного представления с ограниченной точностью, вы не получите точных ответов.

То же самое происходит в десятичной форме.Например, 1/3 не имеет точного десятичного представления.Вы можете сохранить его как .333333.Но тогда у вас есть проблема.Являются ли 3 и .333333 мультипликативными инверсиями?Если вы умножите их, вы получите .999999.Если вы хотите, чтобы ответ был «да», вам нужно создать тест для мультипликативных инверсий, который не так прост, как умножение и проверка на равенство 1.

То же самое происходит с двоичным.

Answer 4

Обсуждения в других ответах великолепны, и поэтому я не буду повторять их, но кода нет. Вот небольшой код, чтобы фактически проверить, дает ли пара чисел с плавающей точкой ровно 1,0 при умножении.

Код делает несколько предположений / утверждений (которые обычно выполняются на платформе x86):
- float 32-битные двоичные (AKA single precision) IEEE-754
- либо int, либо long 32-битные (я решил не полагаться на наличие uint32_t)
- memcpy() копирует числа с плавающей запятой в целые / длинные значения, так что 8873283.0f становится 0x4B076543 (т.е. ожидается определенная "последовательность")

Еще одно допущение:
- он получает действительные числа с плавающей точкой, которые * умножит (то есть, при умножении чисел с плавающей запятой не будут использоваться значения более высокой точности, которые математическое оборудование / библиотека могут использовать для внутреннего использования)

#include <stdio.h>
#include <string.h>
#include <limits.h>
#include <assert.h>

#define C_ASSERT(expr) extern char CAssertExtern[(expr)?1:-1]

#if UINT_MAX >= 0xFFFFFFFF
typedef unsigned int uint32;
#else
typedef unsigned long uint32;
#endif
typedef unsigned long long uint64;

C_ASSERT(CHAR_BIT == 8);
C_ASSERT(sizeof(uint32) == 4);
C_ASSERT(sizeof(float) == 4);

int ProductIsOne(float f1, float f2)
{
  uint32 m1, m2;
  int e1, e2, s1, s2;
  int e;
  uint64 m;

  // Make sure floats are 32-bit IEE754 and
  // reinterpreted as integers as we expect
  {
    static const float testf = 8873283.0f;
    uint32 testi;
    memcpy(&testi, &testf, sizeof(testf));
    assert(testi == 0x4B076543);
  }

  memcpy(&m1, &f1, sizeof(f1));
  s1 = m1 >= 0x80000000;
  m1 &= 0x7FFFFFFF;
  e1 = m1 >> 23;
  m1 &= 0x7FFFFF;
  if (e1 > 0) m1 |= 0x800000;

  memcpy(&m2, &f2, sizeof(f2));
  s2 = m2 >= 0x80000000;
  m2 &= 0x7FFFFFFF;
  e2 = m2 >> 23;
  m2 &= 0x7FFFFF;
  if (e2 > 0) m2 |= 0x800000;

  if (e1 == 0xFF || e2 == 0xFF || s1 != s2) // Inf, NaN, different signs
    return 0;

  m = (uint64)m1 * m2;

  if (!m || (m & (m - 1))) // not a power of 2
    return 0;

  e = e1 + !e1 - 0x7F - 23 + e2 + !e2 - 0x7F - 23;
  while (m > 1) m >>= 1, e++;

  return e == 0;
}

const float testData[][2] =
{
  { .1f, 10.0f },
  { 0.5f, 2.0f },
  { 0.25f, 2.0f },
  { 4.0f, 0.25f },
  { 0.33333333f, 3.0f },
  { 0.00000762939453125f, 131072.0f }, // 2^-17 * 2^17
  { 1.26765060022822940E30f, 7.88860905221011805E-31f }, // 2^100 * 2^-100
  { 5.87747175411143754E-39f, 1.70141183460469232E38f }, // 2^-127 (denormalized) * 2^127
};

int main(void)
{
  int i;
  for (i = 0; i < sizeof(testData) / sizeof(testData[0]); i++)
    printf("%g * %g %c= 1\n",
           testData[i][0], testData[i][1],
           "!="[ProductIsOne(testData[i][0], testData[i][1])]);
  return 0;
}

Вывод (см. ideone.com ):

0.1 * 10 != 1
0.5 * 2 == 1
0.25 * 2 != 1
4 * 0.25 == 1
0.333333 * 3 != 1
7.62939e-06 * 131072 == 1
1.26765e+30 * 7.88861e-31 == 1
5.87747e-39 * 1.70141e+38 == 1

Answer 5

Что поразительно, так это то, что независимо от правила округления, вы ожидаете, что результаты двух версий будут одинаковыми (дважды неправильно или дважды правильно)!

Скорее всего, в первом случае повышение оценки в регистрах FPU происходит при оценке x == 1 / y, тогда как z = 1 / y действительно сохраняет результат с одинарной точностью.

Другие авторы объясняют, почему 5 == 1 / 0.2 может потерпеть неудачу, мне не нужно повторять это.