Как рассчитываются показатели?

Question

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

Я специально ищу алгоритм pow (), используемый для чисел с плавающей запятой двойной точности.

Answer 1

У меня была возможность взглянуть на реализацию fdlibm. Комментарии описывают используемый алгоритм:

 *                    n
 * Method:  Let x =  2   * (1+f)
 *      1. Compute and return log2(x) in two pieces:
 *              log2(x) = w1 + w2,
 *         where w1 has 53-24 = 29 bit trailing zeros.
 *      2. Perform y*log2(x) = n+y' by simulating muti-precision
 *         arithmetic, where |y'|<=0.5.
 *      3. Return x**y = 2**n*exp(y'*log2)

с последующим перечислением всех обработанных особых случаев (0, 1, inf, nan).

Наиболее интенсивные разделы кода, после всей обработки специального случая, включают вычисления log2 и 2**. И нет петель ни в одном из них. Таким образом, несмотря на сложность примитивов с плавающей точкой, он выглядит как асимптотически постоянный алгоритм времени.

Эксперты с плавающей точкой (из которых я не один) приглашаются комментировать. : -)

Answer 2

Если они не нашли лучшего способа сделать это, я считаю, что приблизительные значения для тригонометрических, логарифмических и экспоненциальных функций (например, для экспоненциального роста и затухания) обычно рассчитываются с использованием арифметических правил и ряд Тейлора расширения для получения приблизительного результата с точностью до запрошенной точности. (Подробную информацию о степенных рядах, рядах Тейлора и рядах Маклаурина можно найти в любой книге по исчислению.) Обратите внимание, что с тех пор, как я это сделал, прошло довольно много времени, поэтому я не мог вам сказать, например, как именно рассчитать количество терминов в серии, которое необходимо включить, гарантирует ошибку, достаточно малую, чтобы ее можно было пренебречь в вычислениях с двойной точностью.

Например, разложение в ряд Тейлора / Маклорина для е ^ х таково:

      +inf [ x^k ]           x^2    x^3      x^4        x^5
e^x = SUM  [ --- ] = 1 + x + --- + ----- + ------- + --------- + ....
      k=0  [  k! ]           2*1   3*2*1   4*3*2*1   5*4*3*2*1

Если вы возьмете все члены (k от 0 до бесконечности), это расширение будет точным и полным (без ошибок).

Однако, если вы не принимаете все термины, идущие в бесконечность, но останавливаетесь после, скажем, 5 терминов или 50 терминов или чего-то еще, вы получите приблизительный результат, который отличается от фактического e ^ x значение функции по остатку, который довольно легко вычислить.

Хорошей новостью для экспонент является то, что они хорошо сходятся, и условия его полиномиального расширения довольно легко кодировать итеративно, поэтому вы могли бы (повторяю, MIGHT - помните, это Некоторое время назад) даже не нужно было предварительно рассчитывать, сколько терминов нужно, чтобы гарантировать, что ваша ошибка меньше точности, потому что вы можете проверить размер вклада на каждой итерации и остановиться, когда он станет достаточно близким к нулю. На практике я не знаю, является ли эта стратегия жизнеспособной или нет - я должен был бы попробовать это. Есть важные детали, о которых я давно забыл. Такие вещи, как: точность машины, ошибка машины и ошибка округления и т. Д.

Также обратите внимание, что если вы не используете e ^ x, но вы делаете рост / затухание с другой базой, такой как 2 ^ x или 10 ^ x, аппроксимирующая функция многочлена изменяется.

Answer 3

Обычный подход, чтобы поднять a до b для целочисленной степени, выглядит примерно так:

result = 1
while b > 0
  if b is odd
    result *= a
    b -= 1
  b /= 2
  a = a * a

Обычно логарифмичен по величине показателя степени. Алгоритм основан на инварианте "a ^ b * result = a0 ^ b0", где a0 и b0 - начальные значения a и b.

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

Edit: так как, кажется, есть некоторый интерес, вот версия без дополнительного умножения.

result = 1
while b > 0
  while b is even
    a = a * a
    b = b / 2
  result = result * a
  b = b - 1

Answer 4

Вы можете использовать exp (n * ln (x)) для вычисления x ⁿ . И x, и n могут быть числами с плавающей запятой двойной точности. Натуральный логарифм и экспоненциальная функция могут быть вычислены с использованием рядов Тейлора. Здесь вы можете найти формулы: http://en.wikipedia.org/wiki/Taylor_series

Answer 5

Если бы я писал функцию pow, предназначенную для Intel, я бы возвратил exp2 (log2 (x) * y). Микрокод Intel для log2, безусловно, быстрее, чем все, что я смогу закодировать, даже если бы я мог вспомнить численный анализ исчисления в первом курсе и аспирантуры.