Подтвердить что ты не робот

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

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

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;

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

4b9b3361

Ответ 1

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

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

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

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

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


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

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

Возьмем пару x = 1000 и y = 0,001. Пусть возьмем "стартовую" относительную ошибку у при 10 -6.

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

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

1/y составляет около 1000. Его относительная ошибка равна 10 -6. ( "1" не имеет ошибок)

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

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

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

  • Основываясь на х, у имеет большее значение абс и делит 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 работает даже для поплавков.

Ответ 2

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

 1/5: 0.0011001100110011001100110011001100110011...

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

   x: 0.0011001100110011001100110011001

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

 1/y: 0.0011001100110011001100110011001100110011001100110011
   x: 0.0011001100110011001100110011001000000000000000000000

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

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

   x: 0.0011001100110011001100110011001
   z: 0.0011001100110011001100110011001

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

Ответ 3

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

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

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

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

Ответ 4

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

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

Одно дополнительное предположение: - он получает фактические поплавки, которые * будет умножаться (т.е. умножение поплавков не будет использовать более высокие значения точности, которые математическое аппаратное обеспечение/библиотека может использовать внутри)

#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

Ответ 5

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

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

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