Численная точность pow (a / b, x) и pow (b / a, -x) - PullRequest
Новая
       13

Численная точность pow (a / b, x) и pow (b / a, -x)

8 голосов
/ 09 апреля 2019

Есть ли разница в точности между pow(a/b,x) и pow(b/a,-x)? Если да, то дает ли повышение числа меньше 1 положительной степени или число больше 1 отрицательной степени, чтобы получить более точный результат?

Редактировать: Предположим, процессор x86_64 и компилятор gcc.

Редактировать: я пытался сравнивать, используя некоторые случайные числа. Например: 

printf("%.20f",pow(8.72138221/1.761329479,-1.51231)) // 0.08898783049228660424
printf("%.20f",pow(1.761329479/8.72138221, 1.51231)) // 0.08898783049228659037

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

Ответы [ 4 ]

2 голосов
/ 09 апреля 2019

В общем, форма с положительной силой немного лучше, хотя при малом количестве она, скорее всего, не будет иметь практического эффекта. Конкретные случаи можно выделить. Например, если a или b является степенью двойки, ее следует использовать в качестве знаменателя, поскольку в этом случае деление не имеет ошибки округления.

В этом ответе я предполагаю двоичную с плавающей точкой IEEE-754 с округлением до ближайшего числа, связанного с четностью, и что соответствующие значения находятся в нормальном диапазоне формата с плавающей запятой.

Дано a, b и x со значениями a , b и x , а также реализация pow который вычисляет представимое значение, ближайшее к идеальному математическому значению (фактические реализации обычно не так хороши), pow(a/b, x) вычисляет ( a / b • (1 + e * 1025) * 0 )) x • (1 + e 1 ), где e 0 - ошибка округления, которая возникает в делении, а e 1 - ошибка округления, возникающая в pow, и pow(b/a, -x) вычисляет ( * +1046 * б * 1 047 * / а * 1 049 * • (1 + е * ** 1 052 тысяча пятьдесят один * 2 * 1 053 *)) х • (1 + e 3 ), где e 2 и e 3 ошибки округления в этом делении и pow соответственно.

Каждая из ошибок e 0 … e 3 лежит в интервале [- u / 2, u / 2], где u - единица наименьшей точности (ULP), равная 1 в формате с плавающей запятой. (Обозначение [ p , q ] - это интервал, содержащий все значения от p до q , включая p и q .) В случае, если результат близок к краю бинарной области (где показатель степени с плавающей точкой изменяется, а значение около 1), нижняя граница может быть - u / 4. В настоящее время я не буду анализировать этот случай.

Перезапись, это ( a / b ) x • (1 + e 0 ) x • (1 + e 1 ) и ( a / б * * одна тысяча сто двадцать шесть) * ** +1128 одна тысяча сто двадцать семь * х * ** 1130 тысяча сто двадцать девять * • (1 + * +1131 * е * * 2 тысячи сто тридцать три ) х • (1 + е 3 ). Это показывает, что основная разница заключается в (1 + e 0 ) x против (1 + e 2 ) х . 1 + e 1 против 1 + e 3 также является разницей, но это только окончательное округление. [Я могу рассмотреть дальнейший анализ этого позже, но пока опущу его.]

Рассмотрим (1 + e 0 ) x и (1 + e 2 ) - x . Потенциальные значения первого диапазона выражений [(1− u / 2) x , (1 + u / 2) x ], в то время как второй охватывает [(1 + u / 2) ) - x , (1− u / 2) - x ]. Когда x > 0, второй интервал длиннее первого:

  • Длина первого: (1 + u / 2) x - (1 + u / 2) х .
  • Длина секунды равна (1 / (1− u / 2)) x - (1 / (1 + u / 2)) х .
  • Умножение последнего на (1− u 2 / 2 2 ) x производит (( 1- и * * 2 тысяча двести пятьдесят четыре / 2 2 ) / (1- * * и одна тысяча двести пятьдесят восемь / 2)) х - ((1− u 2 / 2 2 ) / (1 + u / 2)) x = (1 + u / 2) x - (1 + u / 2) x , то есть длина первого интервала.
  • 1− u 2 / 2 2 <1, сo (1− <em>u 2 / 2 2 ) x <1 для положительного <em>x .
  • Поскольку первая длина равна второй длине, умноженной на число меньше единицы, первый интервал короче.

Таким образом, форма, в которой показатель степени положителен,лучше в том смысле, что он имеет более короткий интервал потенциальных результатов.

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

2 голосов
/ 09 апреля 2019

... между pow(a/b,x) и pow(b/a,-x) ... действительно ли повышение числа меньше 1 до положительной степени или числа больше 1 до отрицательной степени дает более точный результат?

Какое бы деление ни было более дугообразным.

Рассмотрим z = x y = 2 y * log2 (x) .

Грубо говоря: ошибка в y * log2(x) увеличивается на значение z, чтобы сформировать ошибку в z. x y очень чувствителен к ошибке в x. Чем больше значение |log2(x)|, тем больше беспокойство.

В случае OP и pow(a/b,p), и pow(b/a,-p), как правило, имеют одинаковые y * log2(x) и одинаковые z и аналогичные ошибки в z. Речь идет о том, как x, y образуются:

a/b и b/a, в общем, оба имеют одинаковую ошибку +/- 0,5 * на последнем месте , и поэтому оба подхода имеют аналогичная ошибка .

Тем не менее, при выбранных значениях a/b против b/a один фактор будет более точным, и именно этот подход с меньшей ошибкой pow().

pow(7777777/4,-p) может быть более точным, чем pow(4/7777777,p).

При отсутствии уверенности в ошибке в делении применяется общее правило: без существенной разницы.

2 голосов
/ 09 апреля 2019

Вот один из способов ответить на такие вопросы, чтобы увидеть, как ведет себя плавающая точка. Это не на 100% правильный способ анализа такого вопроса, но он дает общее представление.

Давайте сгенерируем случайные числа. Вычислите v0=pow(a/b, n) и v1=pow(b/a, -n) в точности с плавающей точкой. И вычислите ref=pow(a/b, n) с двойной точностью, и округлите его, чтобы плавать. Мы используем ref в качестве справочного значения (мы предполагаем, что double имеет гораздо большую точность, чем float, поэтому мы можем верить, что ref можно считать наилучшим из возможных значений. Это верно для IEEE-754 в большинстве случаев) , Затем сложите разницу между v0-ref и v1-ref. Разница должна рассчитываться как «число чисел с плавающей запятой между v и ref».

Обратите внимание, что результаты могут зависеть от диапазона a, b и n (и от качества генератора случайных чисел. Если он действительно плохой, он может дать необъективный результат). Здесь я использовал a=[0..1], b=[0..1] и n=[-2..2]. Кроме того, этот ответ предполагает, что алгоритм float / double Division / pow одного типа, имеет одинаковые характеристики.

Для моего компьютера суммированные различия: 2604828 2603684, это означает, что между ними нет существенной разницы в точности.

Вот код (обратите внимание, этот код предполагает арифметику IEEE-754): 

#include <cmath>
#include <stdio.h>
#include <string.h>

long long int diff(float a, float b) {
    unsigned int ai, bi;
    memcpy(&ai, &a, 4);
    memcpy(&bi, &b, 4);
    long long int diff = (long long int)ai - bi;
    if (diff<0) diff = -diff;
    return diff;
}

int main() {
    long long int e0 = 0;
    long long int e1 = 0;
    for (int i=0; i<10000000; i++) {
        float a = 1.0f*rand()/RAND_MAX;
        float b = 1.0f*rand()/RAND_MAX;
        float n = 4.0f*rand()/RAND_MAX - 2.0f;

        if (a==0||b==0) continue;

        float v0 = std::pow(a/b, n);
        float v1 = std::pow(b/a, -n);
        float ref = std::pow((double)a/b, n);

        e0 += diff(ref, v0);
        e1 += diff(ref, v1);
    }

    printf("%lld %lld\n", e0, e1);
}
0 голосов
/ 09 апреля 2019

Для оценки ошибок округления, как в вашем случае, может быть полезно использовать некоторую библиотеку с множественной точностью, такую ​​как Boost.Multiprecision. Затем вы можете сравнить результаты для различной точности, например, с помощью следующей программы: 

#include <iomanip>
#include <iostream>
#include <boost/multiprecision/cpp_bin_float.hpp>
#include <boost/multiprecision/cpp_dec_float.hpp>

namespace mp = boost::multiprecision;

template <typename FLOAT>
void comp() {
  FLOAT a = 8.72138221;
  FLOAT b = 1.761329479;
  FLOAT c = 1.51231;

  FLOAT e = mp::pow(a / b, -c);
  FLOAT f = mp::pow(b / a, c);

  std::cout << std::fixed << std::setw(40) << std::setprecision(40) << e << std::endl;
  std::cout << std::fixed << std::setw(40) << std::setprecision(40) << f << std::endl;
}

int main() {
  std::cout << "Double: " << std::endl;
  comp<mp::cpp_bin_float_double>();
  td::cout << std::endl;

  std::cout << "Double extended: " << std::endl;
  comp<mp::cpp_bin_float_double_extended>();
  std::cout << std::endl;

  std::cout << "Quad: " << std::endl;
  comp<mp::cpp_bin_float_quad>();
  std::cout << std::endl;

  std::cout << "Dec-100: " << std::endl;
  comp<mp::cpp_dec_float_100>();
  std::cout << std::endl;
}

На моей платформе вывод: 

Double: 
0.0889878304922865903670015086390776559711
0.0889878304922866181225771242679911665618

Double extended: 
0.0889878304922865999079806265115166752366
0.0889878304922865999012043629334822725241

Quad: 
0.0889878304922865999004910375213273866639
0.0889878304922865999004910375213273505527

Dec-100: 
0.0889878304922865999004910375213273881004
0.0889878304922865999004910375213273881004

Демонстрация в реальном времени: https://wandbox.org/permlink/tAm4sBIoIuUy2lO6

Для double первый расчет был более точным, однако, я думаю, здесь нельзя сделать какие-либо общие выводы.

Также обратите внимание, что ваши входные числа не могут быть точно представлены с помощью типа с плавающей точкой двойной точности IEEE 754 (ни один из них). Вопрос в том, заботишься ли ты о точности расчетов с этими точными числами их ближайших представлений.

...