Связанные полномочия с плавающей запятой

Question

Учитывая переменную времени выполнения double значения x и y, если мне нужно вычислить:

double c1 = pow(x, y)
double c2 = pow(x, y + 1.0)

Безопасно ли с точки зрения числовой точности реализовать это как:

double c1 = pow(x, y)
double c2 = c1 * x;

Вы можете предположить, что y> = -1,0.

Answer 1

Каждый шаг математики FP, подобный этому, влечет потенциальную ошибку округления.

Рассмотрим, как y = pow(a,b) приближает y = (a + error_a) ^{b + error_b} + error_y. Ошибка error_y очень чувствительна к error_b, особенно когда y большое.

Добавление может вызывать проблемы, когда положительные y и 1.0 имеют существенно разные величины. Когда они сильно различаются, double c2 = pow(x, y) * x; лучше.

y + 1.0 отменяет биты, когда y находится рядом с -1.0, double c2 = pow(x, y + 1.0) лучше.

Когда бы y + 1.0 точно, используйте pow(x, y + 1.0). Когда y < 0.0, я бы go за pow(x, y + 1.0). Когда y > 1.0, используйте pow(x, y) * x.

В общем, я бы использовал pow(x, y) * x, так как в худшем случае ошибка, вероятно, всего несколько ULP s. (Несколько штук за хороший pow() и еще ½ за *). Численно рискованно вводить ошибку в b из pow(a,b).

Примечания:

Учтите, что худший случай - это когда y велико, а не бесконечно. Это действительно небольшое подмножество всех возможных a,b в pow(a,b), поскольку результат часто равен 0,0, ∞ или NaN.

Answer 2

В общем, вы можете получить разные результаты, используя оба выражения, поскольку pow() использует логарифмы для вычисления окончательного результата, а во втором случае используется смесь логарифмов (в использовании pow()), смешанная с продуктами.

Математически оба выражения

pow(x, y+1.0) ~= exp((y+1.0)*log(x))    /* ~= here means approximately equal, not the C operator */

и

pow(x, y+1.0) ~= pow(x, y)*x ~= exp(y*log(x))*x

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

Обе функции (exp() и log()) имеют альтернативы для случаев, когда их аргументы близки к 0.0 и 1.0 соответственно. Итак, если это так (x слишком близко к 1.0 и / или y*log(x) слишком близко к 0.0), тогда будет лучше использовать варианты log1p() и / или expm1() )

В любом случае для своего эксперимента вы можете попробовать обе формулы и сравнить результаты. Даже вы можете сделать это с помощью pow() (и, вероятно, вы получите разные результаты --- более точные, поскольку, вероятно, pow() будет решать при каждом вызове, какую функцию лучше использовать)

Другой Подсказка состоит в том, что ошибки округления относятся к числам, которые умножаются ... поэтому, если у вас есть относительная ошибка 1.0E-10% в двух числах перед их умножением, относительные ошибки сохранятся (вы получите 2.0E-10%) Но когда вы вычитаете две приблизительно равные величины, относительная ошибка возникает полностью (1.0000000000003 - 1.0 == 3.0E-13, но абсолютная ошибка сохраняется, как и до добавления, на 1.0E-13), что означает, что теперь ошибка составляет 1.0E-13 для величины, которая составляет 3.0E-13, относительная ошибка возросла до 33.0%. Если позже вы будете использовать это число в продукте, относительная ошибка будет go для продукта, что сделает весь ваш расчет бесполезным.

Исходя из этого принципа, кажется, что второй подход лучше (поскольку он имеет без добавлений)

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

/* fea stands for (f)ast
 * (e)xponential (a)lgorithm (d)ouble */
double fead(double x, unsigned n)
{
    double result = 1.0;
    while (n) {
        if (n & 1)
            result *= x;
        x *= x; /* square x */
        n >>= 1; /* divide n by 2 */
    }
    return result;
}

Ниже приведен полный реализация некоторых тестовых примеров:

#include <stdio.h>

/* fea stands for (f)ast
 * (e)xponential (a)lgorithm (d)ouble */
double fead(double x, unsigned n)
{
    double result = 1.0;
    while (n) {
        if (n & 1)
            result *= x;
        x *= x; /* square x */
        n >>= 1; /* divide n by 2 */
    }
    return result;
}

/* test cases to probe */
struct test_case {
   double base; /* base */
   unsigned exp; /* exponent */
   double expctd; /* expected result */
} test_cases_start [] = {
    { 10.0,     45, 1.0E45 },
    { 2.0,      63, 9.223372036854775808E18 },
    { 3.141593, 15, 2.865819336962266259E7 },
    { 0.1,      15, 1.0E-15 },
};

int main()
{
    struct test_case *test_cases_end
        = (struct test_case *)(&test_cases_start + 1);

    struct test_case *p;

    for (p = test_cases_start;
         p < test_cases_end;
         p++)
    {
        double calc = fead(p->base, p->exp);
        double diff = calc - p->expctd;
        printf("fead(%20.17g, %u)\n"
            " calculated = %20.17g\n"
            "   expected = %20.17g\n"
            " difference = %20.17g\n"
            " rel. error = %20.17g\n",
            p->base, p->exp,
            calc,
            p->expctd,
            diff,
            diff/p->expctd);
    }
}

и пробный запуск:

$ fead
fead(                  10, 45)
 calculated = 1.0000000000000001e+45
   expected = 9.9999999999999993e+44
 difference = 1.5845632502852868e+29
 rel. error = 1.584563250285287e-16
fead(                   2, 63)
 calculated = 9.2233720368547758e+18
   expected = 9.2233720368547758e+18
 difference =                    0
 rel. error =                    0
fead(  3.1415929999999999, 15)
 calculated =   28658193.369622648
   expected =   28658193.369622663
 difference = -1.4901161193847656e-08
 rel. error = -5.199616389511387e-16
fead( 0.10000000000000001, 15)
 calculated = 1.0000000000000017e-15
   expected = 1.0000000000000001e-15
 difference = 1.5777218104420236e-30
 rel. error = 1.5777218104420234e-15
$ _