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

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

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

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

      +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 ^ х значение функции от остатка, который довольно легко вычислить.

Хорошая новость для экспонент состоит в том, что они хорошо сходятся, и условия его полиномиального расширения довольно легко кодировать итеративно, поэтому вы могли бы (повторить, 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.

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

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

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

Answer 4

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

Answer 5

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