Расчет максимальной ошибки с плавающей запятой

Question

pow(x,y) вычисляется как e^(y*log(x)) Обычно математические библиотеки вычисляют log(x) с четверной точностью (что отнимает много времени), чтобы избежать потери точности при вычислении y*log(x), поскольку ошибки точности будут увеличиваться при вычислении e^(y*log(x)).Теперь, на случай, если я захочу вычислить pow(x,y) в следующих шагах.

double pow(double x,double y) {
   return exp(y*log(x)); // Without any quad precision multiplication
}

Какова будет максимальная ошибка ULP этой функции.Я знаю, что стандарт IEEE-754 гласит, что любая операция с плавающей запятой должна иметь ошибку менее 0,5 ULP, т.е. 0.5*2^(-52).Так что, если моя операция y*log(x) страдает от ошибки 0,5 ULP, как мне вычислить максимально возможную ошибку ULP для e^(y*log(x))

---

Согласна, что вычисление pow(x,y) довольно сложно.Алгоритмы обычно вычисляют log(x) с более высокой точностью, и умножение между y и log(x) не является простым.Поскольку ошибка ULP зависит от y*log(x), максимальная ошибка будет для наибольшего значения Y*log(x), для которого e^(y*log(x)) не является бесконечностью.Правильно?Как вычислить количество ULP для такого случая?Каково максимальное количество битов мантиссы в формате двойной точности, которое будет меняться от фактического значения в случае наибольшего значения y*log(x)?

---

Обновлен вопрос.Спасибо за всю помощь!

Так что разница в 10 бит приведет к тому, сколько будет ошибка ULP?Я вычислил это как

    ULP = (actual - computed)/ 2^(e-(p-1)) 

, где e - показатель фактического числа, p = 53 для двойной точности.Я читал, что я ULP = 2 ^ (e- (p-1)) Давайте предположим,

    Actual = 1.79282279439444787915898270592 *10^308
    Computed = 1.79282279439451553814547593293 * 10^308 
    error= actual - computed = 6.7659e+294 

Сейчас

1 ULP = 2^(e- (p-1)) 
e = 1023 (exponent of the actual number - bias)
1 ULP = 2^(1023 - (53-1)) = 2^971
ULP error = error/ 1 ULP = 339

Это правильно?

Answer 1

максимальная ошибка с плавающей запятой (для exp(y*log(x)))

x == 0

Ошибка бесконечна;

x <0 </strong>

OP pow() завершается неудачей, однако, когда y не имеет дробной части, она должна завершиться.

x> 0

Пусть z = y*log(x).Ошибка в этом расчете может быть вполне разумной - всего несколько ULP или меньше.Давайте предположим, что 1 ULP или, возможно, z_error = z* 2 ^-53 с учетом обычного double.

Все же рассмотрим его влияние: exp(z + error_z) = exp(z)*exp(error_z).

С z около 710, exp(709.78) составляет около DBL_MAX, поэтому давайте рассмотрим, что гораздо большие значения приводят к бесконечности с OP, pow().

exp(some_small_value)о 1 + some_small_value.exp(error_z) - это тогда 1 + 710*pow(2,-53).Таким образом, окончательная pow() может потерять около log2(710) или 9, 10 плюс несколько битов точности.

Answer 2

У меня нет времени для более полного ответа, но вот что нужно начать:

Предположим, что ошибка в каждой отдельной операции - exp, * и log - равнаБольшая ½ ULP.Это означает, что вычисленное значение log(x) является точно логарифмом ( x ) • (1 + e ₀ ) для некоторых e ₀ удовлетворяющий -? / 2 ≤ e ₀ ≤ ? / 2, где ? - ULP 1 (2 ⁻⁵² IEEE-754основной 64-битный двоичный формат).И вычисленное значение y*log(x) равно y • log ( x ) • (1 + e ₀ ) • (1+ e ₁ ) для некоторых -? / 2 ≤ e ₀ ≤ ? / 2 и для некоторых -? / 2 ≤ e ₁ ≤ ? / 2.И вычисленное значение exp(y*log(x)) точно равно e ^{y • log ( x ) • (1 + e ₀ ) • (1 + e ₁ )} • (1 + e ₂ ) для некоторых -? / 2 ≤ e ₀ , e ₁ , e ₂ ≤ ? / 2.И ошибка, конечно, e ^{y • log ( x ) • (1 + e ₀ ) • (1+ e ₁ )} • (1 + e ₂ ) - e ^{y • log ( x )} .

Затем вы можете начать анализировать это выражение на предмет максимальных и минимальных значений по возможным значениям e ₀ , e ₁ и e ₂ .Если 0 <<em> y и 1 <<em> x , наибольшая ошибка возникает, когда e ₀ = e ₁ = e ₂ = ? / 2.

Обратите внимание, что общие реализации exp и log обычно не правильно округляются.Ошибки нескольких ULP могут быть распространены.Таким образом, вы не должны принимать значение ½ ULP, если не задокументировано, что подпрограмма правильно округлена.