Intereting Posts

matlab и c отличаются функцией cos

У меня есть программа, реализованная в 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 8 double So, b differ a bit (very slightly, but enough to make my debuging task difficult) b = -0.96928123535654842068964853751822374761104583740234 c b = -0.96928123535654830966734607500256970524787902832031 matlab 

Я использую тот же компьютер, Intel Core 2 Duo и 8-байтовый двойной тип данных.

Почему результат отличается?

не использует ли Matlab аппаратную функцию cos, встроенную в Intel?

Есть ли простой способ использовать одну и ту же функцию cos в matlab и c (с точными результатами), даже если немного медленнее, чтобы я мог смело сравнивать результаты моей программы matlab и c?


Обновить:

Большое спасибо за ваши ответы!

Итак, как вы указали, функция cos для matlab и c отличается. Это восхитительно! Я думал, что они используют функцию cos, встроенную в микропроцессор Intel.

Косо версия MATLAB равна (по крайней мере, для этого теста) одной из матриц. вы также можете попробовать из Matlab: b = java.lang.Math.cos (a)

Затем я сделал небольшую функцию MEX, чтобы использовать версию cos c из matlab, и она отлично работает; Это позволяет мне отлаживать мою программу (ту же, что реализована в Matlab и c), и посмотреть, в какой момент они отличаются, что было целью этой публикации.

Единственное, что вызов MEX c cos версии от matlab слишком медленный.

Теперь я пытаюсь вызвать функцию Java cos из c (так же, как и у matlab), посмотрите, идет ли это быстрее.

Использование сценария по адресу http://www.mathworks.com/matlabcentral/fileexchange/1777-from-double-to-string

разница между двумя числами находится только в последнем бите:

 octave:1> bc = -0.96928123535654842068964853751822374761104583740234; octave:2> bm = -0.96928123535654830966734607500256970524787902832031; octave:3> num2bin(bc) ans = -.11111000001000101101000010100110011110111001110001011*2^+0 octave:4> num2bin(bm) ans = -.11111000001000101101000010100110011110111001110001010*2^+0 

Один из них должен быть ближе к «правильному» ответу, если предположить, что значение, заданное для a является точным.

 >> be = vpa('cos(2.89308776595231886830)',50) be = -.96928123535654836529707365425580405084360377470583 >> bc = -0.96928123535654842068964853751822374761104583740234; >> bm = -0.96928123535654830966734607500256970524787902832031; >> abs(bc-be) ans = .5539257488326242e-16 >> abs(bm-be) ans = .5562972757925323e-16 

Таким образом, результат библиотеки C является более точным.

Однако в целях вашего вопроса вы не должны ожидать получить тот же ответ в Matlab и в какой бы библиотеке C вы не связались.

Числа с плавающей запятой хранятся в двоичном, а не в десятичном значении. Палитра с double точностью имеет 52 бита точности, что соответствует примерно 15 значащим десятичным разрядам. Другими словами, первых 15 отличных от нуля десятичных цифр double печати в десятичной форме достаточно, чтобы однозначно определить, какая double была напечатана.

Являясь диадическим рациональным , double имеет точное представление в десятичном значении , которое занимает гораздо больше десятичных знаков, чем 15, чтобы представлять (в вашем случае 52 или 53 места, я считаю). Тем не менее, стандарты для printf и подобных функций не требуют, чтобы цифры, предшествующие 15-м, были правильными; они могут быть полной ерундой. Я подозреваю, что одна из двух сред печатает точное значение, а другая печатает плохое приближение и что в действительности оба соответствуют точному двоящему двоичному значению.

Результат будет таким же, что и до 15 знаков после запятой, я подозреваю, что этого достаточно для почти всех приложений, и если вам требуется больше, вы, вероятно, должны реализовать свою собственную версию косинуса так, чтобы вы контролировали специфику, а ваш код переносится для разных компиляторов C.

Они будут отличаться, потому что они, несомненно, используют разные методы для вычисления аппроксимации результата или повторения нескольких раз. Поскольку косинус определяется как бесконечный ряд терминов, для его реализации программного обеспечения необходимо использовать аппроксимацию. Алгоритм CORDIC является одной общей реализацией.

К сожалению, я не знаю специфики реализации в любом случае, действительно, C будет зависеть от используемой вами реализации стандартной библиотеки C.

Как объяснили другие, когда вы вводите этот номер непосредственно в исходном коде, будут использоваться не все цифры цифр, так как вы получаете только 15/16 десятичных знаков для точности. Фактически, они преобразуются в ближайшее двойное значение в двоичном формате (все, что выходит за пределы фиксированного предела цифр, отбрасывается).

Чтобы ухудшить ситуацию, и согласно @R , IEEE 754 допускает ошибку в последнем бите при использовании функции косинуса. Я на самом деле столкнулся с этим при использовании разных компиляторов.

Чтобы проиллюстрировать, я протестировал со следующим файлом MEX, после компиляции с компилятором LCC по умолчанию, а затем с помощью VS2010 (я на 32-разрядной версии WinXP).

В одной функции мы непосредственно вызываем функции C ( mexPrintf является просто макросом #define как printf ). В другом случае мы называем mexEvalString для того, чтобы выкачать вещи в MATLAB-движке (что эквивалентно использованию командной строки в MATLAB).

prec.c

 #include  #include  #include  #include "mex.h" void c_test() { double a = 2.89308776595231886830L; double b = cos(a); mexPrintf("[C] a = %.25Lf (%16Lx)\n", a, a); mexPrintf("[C] b = %.25Lf (%16Lx)\n", b, b); } void matlab_test() { mexEvalString("a = 2.89308776595231886830;"); mexEvalString("b = cos(a);"); mexEvalString("fprintf('[M] a = %.25f (%bx)\\n', a, a)"); mexEvalString("fprintf('[M] b = %.25f (%bx)\\n', b, b)"); } void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { matlab_test(); c_test(); } 

скопирован с LCC

 >> prec [M] a = 2.8930877659523189000000000 (4007250b32d9c886) [M] b = -0.9692812353565483100000000 (bfef045a14cf738a) [C] a = 2.8930877659523189000000000 ( 32d9c886) [C] b = -0.9692812353565484200000000 ( 14cf738b) <--- 

скомпилирован с VS2010

 >> prec [M] a = 2.8930877659523189000000000 (4007250b32d9c886) [M] b = -0.9692812353565483100000000 (bfef045a14cf738a) [C] a = 2.8930877659523189000000000 ( 32d9c886) [C] b = -0.9692812353565483100000000 ( 14cf738a) <--- 

Я скомпилирую вышеуказанное, используя: mex -v -largeArrayDims prec.c , и переключаюсь между mex -setup компиляторами, используя: mex -setup

Заметьте, что я также попытался напечатать шестнадцатеричное представление чисел. Мне удалось показать только нижнюю половину двоичных двоичных чисел в C (возможно, вы можете получить другую половину, используя какие-то бит-манипуляции, но я не уверен, как!)

Наконец, если вам нужна более высокая точность вычислений, рассмотрите возможность использования библиотеки для арифметики с переменной точностью. В MATLAB, если у вас есть доступ к Symbolic Math Toolbox , попробуйте:

 >> a = sym('2.89308776595231886830'); >> b = cos(a); >> vpa(b,25) ans = -0.9692812353565483652970737 

Таким образом, вы можете видеть, что фактическое значение находится где-то между двумя разными приближениями, которые я получил выше, и на самом деле все они равны до 15-го знака после запятой:

 -0.96928123535654831.. # 0xbfef045a14cf738a -0.96928123535654836.. # <--- actual value (cannot be represented in 64-bit) -0.96928123535654842.. # 0xbfef045a14cf738b ^ 15th digit --/ 

ОБНОВИТЬ:

Если вы хотите правильно отобразить шестнадцатеричное представление чисел с плавающей запятой в C, используйте вместо этого эту вспомогательную функцию (аналогично функции NUM2HEX в MATLAB):

 /* you need to adjust for double/float datatypes, big/little endianness */ void num2hex(double x) { unsigned char *p = (unsigned char *) &x; int i; for(i=sizeof(double)-1; i>=0; i--) { printf("%02x", p[i]); } }