我有一个在matlab中实现的程序和c中的相同程序,结果不同.
我有点困惑的是cos函数没有返回完全相同的结果.
在这两种情况下,我使用相同的计算机,Intel Core 2 Duo和8字节双数据类型.
为什么结果不同？
这是测试:
c:
double a = 2.89308776595231886830;
double b = cos(a);
printf("a = %.50f\n", a);
printf("b = %.50f\n", b);
printf("sizeof(a): %ld\n", sizeof(a));
printf("sizeof(b): %ld\n", sizeof(b));
a = 2.89308776595231886830106304842047393321990966796875
b = -0.96928123535654842068964853751822374761104583740234
sizeof(a): 8
sizeof(b): 8
matlab:
a = 2.89308776595231886830
b = cos(a);
fprintf('a = %.50f\n', a);
fprintf('b = %.50f\n', b);
whos('a')
whos('b')
a = 2.89308776595231886830106304842047393321990966796875
b = -0.96928123535654830966734607500256970524787902832031
Name Size Bytes Class Attributes
a 1x1 8 double
Name Size Bytes Class Attributes
b 1x1 …我试图用它的c ++代码测试一个simulink块,simulink块包含一些algebratic,三角函数和积分器.在我的测试过程中,从simulink块输入使用随机数生成器,输入和输出都被记录到mat文件中(使用MatIO),将由C++代码读取,输出与C++计算得到.对于仅包含代数函数的信号,结果是精确的,差值为零,对于包含三角函数的路径,差值约为10e-16.matlab社区声称他们是正确的而glibc不是.
最近我发现在glibc中实现的三角函数的输出值不等于matlabs中产生的值,根据旧问题1 2 3和我的实验,这些差异与1 glp> glibc的准确性有关.对于大多数块而言,这个10e-16误差感觉不大,但是在积分器的输出中,10e-16积累的越来越多,积分器的最终误差将是大约1e-3,这有点高,这种阻止是不可接受的.
经过对该问题的大量研究后,我决定使用其他方法来计算sin/cos函数,而不是glibc中提供的函数.
我实施了这些apporaches,
1-泰勒系列具有长双变量和-O2(强制使用x87 FPU及其80位浮点运算)
2-taylor系列与GNU quadmath库(128位精度)
3个MPFR库(128位)
4- CRLibm(正确舍入的libm)
5- Sun的LibMCR(就像CRLibm一样)
6- X86 FSIN/FCOS具有不同的舍入模式
7- Java.lang.math到JNI(我认为matlab使用)
8-fdlibm(根据我见过的一篇博文)
9- openlibm
10-通过mex/matlab引擎调用matlab函数
除了最后一个以外的实验不能生成等于matlab的值.我测试了所有这些库和方法的广泛输入,其中一些像libmcr和fdlibm将为一些输入产生NAN值(看起来他们没有良好的范围检查),其余的产生值错误10e-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 …