Ваш пример не обязательно дает правильные результаты с оптимизацией
уровни -O1 или выше. Смотрите эту Годболт ссылку :
компилятором генерируется только одно дополнение vaddsd.
с оптимизацией
level -O0 сборка выглядит нормально, но это может привести к неэффективному коду.
Более того, звонить fegetround и fesetround относительно дорого,
по сравнению со стоимостью нескольких операций с плавающей запятой.
(Самоочевидный) код ниже, вероятно, является интересной альтернативой.
Он использует известные алгоритмы 2Sum и 2ProdFMA. В системах без аппаратного fma или эмуляции fma вы можете использовать алгоритм 2Prod вместо 2ProdFMA,
См., Например, Точный продукт с плавающей точкой и экспонента ,
Стеф Граиллат.
/*
gcc -m64 -Wall -O3 -march=haswell round_ex.c -lm
or with fma emulation on systems without hardware fma support, for example:
gcc -m64 -Wall -O3 -march=nehalem round_ex.c -lm
*/
#include<math.h>
#include<float.h>
#include<stdio.h>
int add_is_not_exact(double operand1, double operand2){
double a = operand1;
double b = operand2;
double s, t, a_1, b_1, d_a, d_b;
/* Algorithm 2Sum computes s and t such that a + b = s + t, exactly. */
/* Here t is the error of the floating-point addition s = a + b. */
/* See, for example, On the robustness of the 2Sum and Fast2Sum algorithms */
/* by Boldo, Graillat, and Muller */
s = a + b;
a_1 = s - b;
b_1 = s - a_1;
d_a = a - a_1;
d_b = b - b_1;
t = d_a + d_b;
return (t!=0.0);
}
int sub_is_not_exact(double operand1, double operand2){
return add_is_not_exact(operand1, -operand2);
}
int mul_is_not_exact(double operand1, double operand2){
double a = operand1;
double b = operand2;
double s, t;
/* Algorithm 2ProdFMA computes s and t such that a * b = s + t, exactly. */
/* Here t is the error of the floating-point multiplication s = a * b. */
/* See, for example, Accurate Floating Point Product and Exponentiation */
/* by Graillat */
s = a * b;
t = fma(a, b, -s);
if (s!=0) return (t!=0.0); /* No underflow of a*b */
else return (a!=0.0)&&(b!=0.0); /* Underflow: inexact if s=0, but (a!=0.0)&&(b!=0.0) */
}
int div_is_not_exact(double operand1, double operand2){
double a = operand1;
double b = operand2;
double s, t;
s = a / b;
t = fma(s, b, -a); /* fma(x,y,z) computes x*y+z with infinite intermediate precision */
return (t!=0.0);
}
int main(){
printf("add_is_not_exact(10.0, 1.0) = %i\n", add_is_not_exact(10.0, 1.0));
printf("sub_is_not_exact(10.0, 1.0) = %i\n", sub_is_not_exact(10.0, 1.0));
printf("mul_is_not_exact( 2.5, 2.5) = %i\n", mul_is_not_exact( 2.5, 2.5));
printf("div_is_not_exact( 10, 2.5) = %i\n", div_is_not_exact( 10, 2.5));
printf("add_is_not_exact(10.0, 0.1) = %i\n", add_is_not_exact(10.0, 0.1));
printf("sub_is_not_exact(10.0, 0.1) = %i\n", sub_is_not_exact(10.0, 0.1));
printf("mul_is_not_exact( 2.6, 2.6) = %i\n", mul_is_not_exact( 2.6, 2.6));
printf("div_is_not_exact( 10, 2.6) = %i\n", div_is_not_exact( 10, 2.6));
printf("\n0x1.0p-300 = %20e, 0x1.0p-600 = %20e \n", 0x1.0p-300 , 0x1.0p-600 );
printf("mul_is_not_exact( 0x1.0p-300, 0x1.0p-300) = %i\n", mul_is_not_exact( 0x1.0p-300, 0x1.0p-300));
printf("mul_is_not_exact( 0x1.0p-600, 0x1.0p-600) = %i\n", mul_is_not_exact( 0x1.0p-600, 0x1.0p-600));
}
Вывод:
$ ./a.out
add_is_not_exact(10.0, 1.0) = 0
sub_is_not_exact(10.0, 1.0) = 0
mul_is_not_exact( 2.5, 2.5) = 0
div_is_not_exact( 10, 2.5) = 0
add_is_not_exact(10.0, 0.1) = 1
sub_is_not_exact(10.0, 0.1) = 1
mul_is_not_exact( 2.6, 2.6) = 1
div_is_not_exact( 10, 2.6) = 1
0x1.0p-300 = 4.909093e-91, 0x1.0p-600 = 2.409920e-181
mul_is_not_exact( 0x1.0p-300, 0x1.0p-300) = 0
mul_is_not_exact( 0x1.0p-600, 0x1.0p-600) = 1
1020 *
*
<Ч />
Как отмечено в комментариях, также можно напрямую прочитать
регистр управления и состояния:
#include <fenv.h>
#pragma STDC FENV_ACCESS ON
int add_is_not_exact_v2(double a, double b)
{
fexcept_t excepts;
feclearexcept(FE_ALL_EXCEPT);
double c = a+b;
int tst = fetestexcept(FE_INEXACT);
return (tst!=0);
}
Обратите внимание, что это может не работать с уровнем оптимизации компилятора -O1 или выше.
В этом случае команда двойного добавления addsd иногда оптимизируется полностью,
приводя к неправильным результатам.
Например, с gcc 8.2 gcc -m64 -O1 -march=nehalem:
add_is_not_exact_v2:
sub rsp, 8
mov edi, 61
call feclearexcept
mov edi, 32
call fetestexcept
test eax, eax
setne al
movzx eax, al
add rsp, 8
ret
С уровнем оптимизации -O0, с 2 вызовами функций и с относительно
Обширные инструкции по изменению контроля и регистра состояния, это не обязательно самое эффективное решение.