Обнаружение ошибок округления

Question

У меня есть два целых числа n и d. Они могут быть точно представлены двойным dn (n) и двойным dd (d). Есть ли надежный способ в C ++ проверить, если

double result = dn/dd

содержит ошибку округления? Если бы это была просто целочисленная проверка, будет ли работать (n/d) * d==n, но выполнение с арифметикой двойной точности может скрыть ошибки округления.

Редактировать: Вскоре после публикации это меня поразило, что изменение режима округления на round_down заставило бы тест (n/d)*d==n работать на удвоение. Но если есть более простое решение, я все еще хотел бы услышать это.

Answer 1

Если частное q=dn/dd является точным, оно разделит dn ровно на dd раз.
Поскольку у вас dd целое число, вы можете проверить точность с целочисленным делением.
Вместо того, чтобы проверять частное, умноженное на dd с (dn/dd)*dd==dn, где ошибки округления могут компенсировать, вы должны скорее проверить остаток.
Действительно, std:remainder всегда точно:

if(std:remainder(dn,dn/dd)!=0)
    std::cout << "Quotient was not exact." << std::endl;

Answer 2

Если доступен аппаратный FMA, то в большинстве случаев (в случаях, когда n ожидается не малым, как указано ниже), самый быстрый тест может быть:

#include <cmath>
…
double q = dn/dd;
if (std::fma(-q, dd, dn))
    std::cout << "Quotient was not exact.\n";

Это может не сработать, если nd - q • dd настолько мало, что округляется до нуля, что происходит при округлении до ближайших связей. в режим четности, если его величина меньше половины наименьшего представимого положительного значения (обычно 2 ^-1074 ). Это может произойти, только если dn само по себе мало. Я ожидаю, что мог бы рассчитать некоторую оценку для dn для этого, если это необходимо, и, учитывая, что dn = n и n - целое число, это не должно происходить.

Игнорирование границ экспонент, способ проверки значений и делимости:

#include <cfloat>
#include <cmath>
…
int sink; // Needed for frexp argument but will be ignored.
double fn = std::ldexp(std::frexp(n, &sink), DBL_MANT_DIG);
double fd = std::frexp(d, &sink);
if (std::fmod(fn, fd))
    std::cout << "Quotient will not be exact.\n";

Учитывая, что n и d являются целыми числами, которые точно представимы в типе с плавающей запятой, я думаю, мы могли бы показать, что их показатели не могут быть такими, чтобы вышеприведенный тест не удался. Есть случаи, когда n - это маленькое целое число, а d - большое (значение от 2 ¹⁰²³ до 2 ¹⁰²⁴ −2 ⁹⁷² включительно), о котором мне нужно подумать.

Answer 3

Есть ли в C ++ надежный способ проверить, содержит ли double result = dn/dd ошибку округления?

Если ваша система разрешает доступ к различным флагам FP, проверьте на FE_INEXACT после разделения.

Если код FP стоит дорого, то по крайней мере этот код можно использовать для проверки решений только для целых чисел.

---

Далее следует решение C (у меня нет доступа к совместимому компилятору C ++ для тестирования прямо сейчас)

#include <fenv.h>
// Return 0: no rounding error
// Return 1: rounding error
// Return -1: uncertain
#pragma STDC FENV_ACCESS ON
int Rounding_error_detection(int n, int d) {
  double dn = n;
  double dd = d;

  if (feclearexcept(FE_INEXACT)) return -1;
  volatile double result = dn/dd;
  (void) result;
  int set_excepts = fetestexcept(FE_INEXACT);
  return set_excepts != 0;
}

Тестовый код

void Rounding_error_detection_Test(int n, int d) {
  printf("Rounding_error_detection(%d, %d) --> %d\n", 
    n, d, Rounding_error_detection(n,d));
}

int main(void) {
  Rounding_error_detection_Test(3, 6);
  Rounding_error_detection_Test(3, 7);
}

выход

Rounding_error_detection(3, 6) --> 0
Rounding_error_detection(3, 7) --> 1

Answer 4

Если вы игнорируете переполнение и недостаток (что вы должны делать, если целочисленные типы, представляющие d и n очень широкие), то (двоичное) деление с плавающей точкой dn/dd точно тогда и только тогда d - это делитель, умноженный на n на двойную степень.

Алгоритм проверки этого может выглядеть следующим образом:

assert(d != 0);
while (d & 1 == 0) d >>= 1; // extract largest odd divisor of d
int exact = n % d == 0;

Это дешевле, чем изменение режима округления FPU, если вы хотите, чтобы режим округления был «ближайшим» в остальное время, и, вероятно, существуют хитрые трюки, которые могут ускорить извлечение наибольшего нечетного делителя d.