Я пытаюсь протестировать блок simulink с его кодом c ++, блок simulink содержит некоторые алгебраические, тригонометрические функции и интеграторы. В моей тестовой процедуре генератор случайных чисел используется от входа блока simulink, и вход и выход записываются в файл mat (используя MatIO), который будет считываться кодом C ++, а результаты сравниваются с вычисляемыми значениями C ++. для сигналов, содержащих только алгебраические функции, результаты являются точными и разница равна нулю, для путей, содержащих тригонометрические функции, разница составляет около 10e-16.
Сообщество Matlab утверждает, что они верны, а glibc - нет.
Недавно я обнаружил, что выходные значения тригонометрических функций, реализованных в glibc, не совпадают со значениями, полученными в matlabs, согласно старым вопросам 1 2 3 и в моих экспериментах различия были связаны с точностью glibc. для большей части блока эта ошибка 10e-16 не имеет особого смысла, но на выходе интегратора 10e-16 накапливается все больше и больше, и окончательная ошибка интегратора будет около 1e-3, что немного выше и неприемлемо для такого типа блоков.
После долгих исследований этой проблемы я решил использовать другие подходы для вычисления функций sin / cos, чем те, которые представлены в glibc.
Я реализовал эти аппорации,
1 - ряд Тейлора с длинными двойными переменными и -O2 (что заставляет использовать FPU x87 с его арифметикой с плавающей запятой 80 бит)
2-ряд Тейлора с библиотекой квадратов GNU (точность 128 бит)
3 - библиотека MPFR (128 бит)
4- CRLibm (правильно округленный libm)
5- LibMCR от Sun (как и CRLibm)
6- X86 FSIN / FCOS с различными режимами округления
7- Java.lang.math через JNI (как я думаю, использует matlab)
8 - fdlibm (согласно одному из постов, которые я видел)
9 - openlibm
10 - вызов функции matlab через mex / matlab engine
Ни один из вышеперечисленных экспериментов, кроме последнего, не мог генерировать значения, равные matlab. я протестировал все эти библиотеки и подходы для широкого диапазона входных данных, некоторые из них, такие как libmcr и fdlibm, будут выдавать значения NAN для некоторых входных данных (похоже, они не имеют хороших проверок диапазона), а остальные - с ошибка 10е-16 и выше.
Только последнее выдает правильные значения по сравнению с matlab, как это было исключено, но вызов функции matlab не эффективен и намного медленнее, чем собственные реализации.
также я удивился, почему ряды MPFR и Тейлора с длинным двойным и четверным квадратами ошибаются.
Это ряд Тейлора с длинными двойными переменными (точность 80 бит)
и должен быть скомпилирован с -O2, который предотвращает сохранение значений из стека FPU в регистры (от 80 бит до 64 бит = потеря точности), также перед выполнением любых вычислений режим округления x87 будет установлен в ближайшее
typedef long double dt_double;
inline void setFPUModes(){
unsigned int mode = 0b0000111111111111;
asm(
"fldcw %0;"
: : "m"(mode));
}
inline dt_double factorial(int x) //calculates the factorial
{
dt_double fact = 1;
for (; x >= 1 ; x--)
fact = x * fact;
return fact;
}
inline dt_double power(dt_double x, dt_double n) //calculates the power of x
{
dt_double output = 1;
while (n > 0)
{
output = (x * output);
n--;
}
return output;
}
inline double sin(double x) noexcept //value of sine by Taylors series
{
setFPUModes();
dt_double result = x;
for (int y = 1 ; y != 44; y++)
{
int k = (2 * y) + 1;
dt_double a = (y%2) ? -1.0 : 1.0;
dt_double c = factorial(k);
dt_double b = power(x, k);
result = result + (a * b) / c;
}
return result;
}
Подход серии Тейлора, протестированный со всеми четырьмя режимами округления x87, лучший из которых имеет ошибку 10e-16
Это X87 fpu one
double sin(double x) noexcept
{
double d;
unsigned int mode = 0b0000111111111111;
asm(
"finit;"
"fldcw %2;"
"fldl %1;"
"fsin;"
"fstpl %0" :
"+m"(d) : "m"(x), "m"(mode)
);
return d;
}
также код x87 fpu не более точен, чем предыдущий
Вот код для MPFR
double sin(double x) noexcept{
mpfr_set_default_prec(128);
mpfr_set_default_rounding_mode(MPFR_RNDN);
mpfr_t t;
mpfr_init2(t, 128);
mpfr_set_d(t, x, MPFR_RNDN);
mpfr_t y;
mpfr_init2(y, 128);
mpfr_sin(y, t, MPFR_RNDN);
double d = mpfr_get_d(y, MPFR_RNDN);
mpfr_clear(t);
mpfr_clear(y);
return d;
}
Я не могу понять, почему версия MPFR не работает должным образом
также коды для всех других подходов, которые я тестировал, такие же, и все они имеют ошибки по сравнению с Matlab.
все коды проверены на широкий диапазон чисел, и я обнаружил простые случаи, которые они терпят неудачу. например:
в matlab следующий код выдает 0x3fe1b071cef86fbe, но в этих apporoches я получил 0x3fe1b071cef86fbf (разница в последнем бите)
format hex;
sin(0.5857069572718263)
ans = 0x3fe1b071cef86fbe
Чтобы прояснить вопрос,
Как я описал выше, эта однобитовая неточность важна, когда она подается в интегратор, и я ищу решение, чтобы получить значения, точно такие же, как у matlab. есть идеи?
Update1:
1 Ошибка Ulp никак не влияет на вывод алгоритма, но предотвращает проверку результатов с помощью matlab, особенно при выводе интеграторов.
Как сказал @Джон Боллинджер, ошибки не накапливаются в прямом пути с несколькими арифметическими блоками, но не при подаче в дискретный интегратор
Update2:
Я посчитал количество неравных результатов для всех вышеперечисленных подходов, ясно, что openlibm будет давать меньше неравных значений по сравнению с matlab, но это не ноль.