Подтвердить что ты не робот

(0.3) ^ 3 == (0.3) * (0.3) * (0.3) возвращает false в matlab?

Я пытаюсь понять ошибку округления для основных арифметических операций в MATLAB, и я наткнулся на следующий любопытный пример.

(0.3)^3 == (0.3)*(0.3)*(0.3)

ans = 0

Я хотел бы точно знать, как вычисляется левая часть. Документация MATLAB предполагает, что для целых степеней используется алгоритм "возведения в степень возведения в квадрат".

"Матричная мощность X ^ p равна X мощности p, если p - скаляр. Если p - целое число, мощность вычисляется путем повторного квадратирования."

Итак, я предположил, что (0.3)^3 и (0.3)*(0.3)^2 вернут одно и то же значение. Но это не так. Как объяснить разницу в ошибке округления?

4b9b3361

ОТВЕТЫ

Ответ 1

Благодаря @Dougal я нашел это:

#include <stdio.h>
int main() {
  double x = 0.3;
  printf("%.40f\n", (x*x*x));

  long double y = 0.3;
  printf("%.40f\n", (double)(y*y*y));
}

который дает:

0.0269999999999999996946886682280819513835
0.0269999999999999962252417162744677625597

Случай странный, потому что вычисление с большим количеством цифр дает худший результат. Это связано с тем, что в любом случае начальное число 0,3 приближается к нескольким цифрам и, следовательно, мы начинаем с относительно большой ошибки. В этом конкретном случае происходит то, что вычисление с несколькими цифрами дает еще одну "большую" ошибку, но с противоположным знаком... поэтому компенсирует начальную. Вместо этого вычисление с большим количеством цифр дает вторую меньшую ошибку, но первая остается.

Ответ 2

Я ничего не знаю о MATLAB, но я попробовал его в Ruby:

irb> 0.3 ** 3
  => 0.026999999999999996
irb> 0.3 * 0.3 * 0.3
  => 0.027

В соответствии с исходным кодом Ruby оператор экспоненции передает правый операнд в float, если левый операнд является поплавком, а затем вызывает стандартная функция C pow(). Вариант float функции pow() должен реализовывать более сложный алгоритм обработки нецелых показателей, который будет использовать операции, которые приводят к ошибке округления. Возможно, MATLAB работает аналогично.

Ответ 3

Интересно, что скаляр ^, по-видимому, реализован с использованием pow, а матрица ^ реализована с использованием квадрата и умножения. К остроумию:

octave:13> format hex
octave:14> 0.3^3
ans = 3f9ba5e353f7ced8
octave:15> 0.3*0.3*0.3
ans = 3f9ba5e353f7ced9
octave:20> [0.3 0;0 0.3]^3
ans =

  3f9ba5e353f7ced9  0000000000000000
  0000000000000000  3f9ba5e353f7ced9

octave:21> [0.3 0;0 0.3] * [0.3 0;0 0.3] * [0.3 0;0 0.3]
ans =

  3f9ba5e353f7ced9  0000000000000000
  0000000000000000  3f9ba5e353f7ced9

Это подтверждается запуском октавы под gdb и установкой точки останова в pow.

То же самое, скорее всего, верно в matlab, но я не могу проверить.

Ответ 4

Здесь небольшая тестовая программа, которая следует за тем, что система pow() из Source/Intel/xmm_power.c, в Apple Libm-2026, в этом случае

#include <stdio.h>
int main() {
    // basically lines 1130-1157 of xmm_power.c, modified a bit to remove
    // irrelevant things

    double x = .3;
    int i = 3;

    //calculate ix = f**i
    long double ix = 1.0, lx = (long double) x;

    //calculate x**i by doing lots of multiplication
    int mask = 1;

    //for each of the bits set in i, multiply ix by x**(2**bit_position)
    while(i != 0)
    {
        if( i & mask )
        {
            ix *= lx;
            i -= mask;
        }
        mask += mask;
        lx *= lx; // In double this might overflow spuriously, but not in long double
    }

    printf("%.40f\n", (double) ix);
}

Это выводит 0.0269999999999999962252417162744677625597, что согласуется с результатами, которые я получаю для .3 ^ 3 в Matlab и .3 ** 3 в Python (и мы знаем последний просто вызывает этот код). Напротив, .3 * .3 * .3 для меня получает 0.0269999999999999996946886682280819513835, это то же самое, что вы получаете, если просто попросите распечатать 0.027 для того, что много десятичных знаков и, следовательно, является ближайшим двойным.

Итак, есть алгоритм. Мы могли точно определить, какое значение задано на каждом шаге, но это не слишком удивительно, что он будет округлять до очень немного меньшего числа, учитывая другой алгоритм для этого.