Блог

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

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

У меня была возможность ознакомиться с реализацией 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** вычисления и. И ни в том, ни в другом нет петель. Таким образом, несмотря на сложность примитивов с плавающей запятой, это выглядит как алгоритм с асимптотически постоянным временем.

Эксперты с плавающей запятой (к которым я не отношусь) могут прокомментировать. 🙂

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

Например, разложение в ряд Тейлора/Маклорена для 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^x на остаток, который довольно легко вычислить.

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

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

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

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

Как правило, он логарифмичен по размеру показателя. Алгоритм основан на инварианте «a^b*результат = 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
 

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

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

e^x = (1 дробь) * (2^показатель степени), 1 <= 1 дробь

x * log2(e) = log2(1 дробь) показатель степени, 0 <= log2(1 дробь)

показатель = этаж(x * log2(e))

1 дробь = 2^(x * log2(e) — показатель) = e^((x * log2(e) — показатель) * ln2) = e^(x — показатель * ln2), 0 <= x — показатель * ln2