Ответ 1
Я смог воспроизвести проблему с одной версией компилятора.
Мина - это MinGW g++ 4.6.2.
Когда я скомпилирую программу как g++ -g -O2 bugflt.cpp -o bugflt.exe, я получаю 720720.
Это разборка main():
_main:
pushl %ebp
movl %esp, %ebp
andl $-16, %esp
subl $16, %esp
call ___main
movl $720720, 4(%esp)
movl $__ZSt4cout, (%esp)
call __ZNSolsEi
movl %eax, (%esp)
call __ZSt4endlIcSt11char_traitsIcEERSt13basic_ostreamIT_T0_ES6_
xorl %eax, %eax
leave
ret
Как вы можете видеть, значение вычисляется во время компиляции.
Когда я скомпилирую его как g++ -g -O2 -fno-inline bugflt.cpp -o bugflt.exe, я получаю 720719.
Это разборка main():
_main:
pushl %ebp
movl %esp, %ebp
andl $-16, %esp
subl $32, %esp
call ___main
movl $1, 4(%esp)
movl $13, (%esp)
call __ZSt3powIiiEN9__gnu_cxx11__promote_2INS0_11__enable_ifIXaasrSt15__is_arithmeticIT_E7__valuesrS3_IT0_E7__valueES4_E6__typeES6_E6__typeES4_S6_
fmuls LC1
fnstcw 30(%esp)
movw 30(%esp), %ax
movb $12, %ah
movw %ax, 28(%esp)
fldcw 28(%esp)
fistpl 4(%esp)
fldcw 30(%esp)
movl $__ZSt4cout, (%esp)
call __ZNSolsEi
movl $__ZSt4endlIcSt11char_traitsIcEERSt13basic_ostreamIT_T0_ES6_, 4(%esp)
movl %eax, (%esp)
call __ZNSolsEPFRSoS_E
xorl %eax, %eax
leave
ret
...
LC1:
.long 1196986368 // 55440.0 exactly
Если я заменил вызов на exp() загрузкой 13.0 следующим образом:
_main:
pushl %ebp
movl %esp, %ebp
andl $-16, %esp
subl $32, %esp
call ___main
movl $1, 4(%esp)
movl $13, (%esp)
// call __ZSt3powIiiEN9__gnu_cxx11__promote_2INS0_11__enable_ifIXaasrSt15__is_arithmeticIT_E7__valuesrS3_IT0_E7__valueES4_E6__typeES6_E6__typeES4_S6_
fildl (%esp)
fmuls LC1
fnstcw 30(%esp)
movw 30(%esp), %ax
movb $12, %ah
movw %ax, 28(%esp)
fldcw 28(%esp)
fistpl 4(%esp)
fldcw 30(%esp)
movl $__ZSt4cout, (%esp)
call __ZNSolsEi
movl $__ZSt4endlIcSt11char_traitsIcEERSt13basic_ostreamIT_T0_ES6_, 4(%esp)
movl %eax, (%esp)
call __ZNSolsEPFRSoS_E
xorl %eax, %eax
leave
ret
Я получаю 720720.
Если я установил те же поля управления округлением и точностью управления управляющим словом x87 FPU на время exp(), как и для команды fistpl 4(%esp), как это:
_main:
pushl %ebp
movl %esp, %ebp
andl $-16, %esp
subl $32, %esp
call ___main
movl $1, 4(%esp)
movl $13, (%esp)
fnstcw 30(%esp)
movw 30(%esp), %ax
movb $12, %ah
movw %ax, 28(%esp)
fldcw 28(%esp)
call __ZSt3powIiiEN9__gnu_cxx11__promote_2INS0_11__enable_ifIXaasrSt15__is_arithmeticIT_E7__valuesrS3_IT0_E7__valueES4_E6__typeES6_E6__typeES4_S6_
fldcw 30(%esp)
fmuls LC1
fnstcw 30(%esp)
movw 30(%esp), %ax
movb $12, %ah
movw %ax, 28(%esp)
fldcw 28(%esp)
fistpl 4(%esp)
fldcw 30(%esp)
movl $__ZSt4cout, (%esp)
call __ZNSolsEi
movl $__ZSt4endlIcSt11char_traitsIcEERSt13basic_ostreamIT_T0_ES6_, 4(%esp)
movl %eax, (%esp)
call __ZNSolsEPFRSoS_E
xorl %eax, %eax
leave
ret
Я получаю 720720.
Отсюда я могу только заключить, что exp() не вычисляет 13 1 точно как 13.0.
Возможно, стоит посмотреть на исходный код этого __gnu_cxx::__promote_2<__gnu_cxx::__enable_if<(std::__is_arithmetic<int>::__value)&&(std::__is_arithmetic<int>::__value), int>::__type, int>::__type std::pow<int, int>(int, int), чтобы увидеть, как ему удается зачеркнуть экспоненцию целыми числами (см., в отличие от C exp(), вместо двух doubles требуется два ints)..
Но я бы не стал винить exp() за это. С++ 11 определяет float pow(float, float) и long double pow(long double, long double) в дополнение к C double pow(double, double). Но в стандарте нет double pow(int, int).
Тот факт, что компилятор предоставляет версию для целых аргументов, не дает дополнительной гарантии точности результата. Если exp() вычисляет a b как
a b= 2 b * log 2 (a)
или
a b= e b * ln (a)
для значений с плавающей запятой, определенно могут быть ошибки округления в процессе.
Если "целочисленная" версия exp() делает что-то подобное и несет аналогичную потерю точности из-за ошибок округления, она по-прежнему выполняет свою работу правильно. И он делает это, даже если потеря точности обусловлена некоторой глупой ошибкой, а не из-за обычных ошибок округления.
Как ни удивительно, это может показаться, это правильно. Или я верю, пока не доказано, что неправильно.