Какие алгоритмы возведения в степень используют языки CPU/программирования?

Question

Я узнал о более быстрых алгоритмах возведения в степень (k-ary, раздвижные двери и т. Д.) И задавался вопросом, какие из них используются в процессорах/языках программирования? (Я неясен в отношении того, происходит ли это в CPU или через компилятор)

И только для пинков, что является самым быстрым?

Редактировать в отношении широты: это намеренно широкий, потому что я знаю, что для этого существует множество различных приемов. У проверенного ответа было то, что я искал.

Answer 1

Я предполагаю, что ваш интерес заключается в реализации функций возведения в степень, которые могут быть найдены в стандартных математических библиотеках для HLL, в частности, C/C++. Они включают в себя функции exp(), exp2(), exp10() и pow(), а также с одинарной точностью аналогов expf(), exp2f(), exp10f() и powf().

Методы возведения в степень, которые вы упоминаете (например, k-ary, скользящее окно), обычно используются в криптографических алгоритмах, таких как RSA, которые основаны на экспоненциальности. Они обычно не используются для функций возведения в степень, предоставляемых через math.h или cmath. Детали реализации для стандартных математических функций, как exp() различаются, но общая схема следует трехступенчатый процесс:

1. уменьшение аргумента функции к первичной аппроксимации интервала
2. приближения подходящей базовой функции на первичный интервал аппроксимации
3. отображение обратно результат для первичного интервала ко всему диапазону функции

вспомогательный шаг часто обработка специальных КАС эс. Они могут относиться к специальным математическим ситуациям, таким как log(0.0), или специальным операндам с плавающей запятой, таким как NaN (Not a Number).

Код C99 для expf(float) ниже показывает примерным образом, что эти шаги выглядят для конкретного примера. Аргумент a сначала разделить таким образом, что exp(a) = е ^г * 2 ^я, где i представляет собой целое число, и r находится в [LOG (SQRT (0.5), журнал (SQRT (2.0)], первичный интервал аппроксимации. На втором этапе мы теперь приближаем e ^r с полиномом. Такие аппроксимации могут быть сконструированы в соответствии с различными критериями проектирования, такими как минимизация абсолютной или относительной погрешности. Полином можно оценивать по-разному, включая схему Хорнера и схему Эстрина.

Приведенный ниже код использует очень общий подход, используя минимаксную аппроксимацию, которая минимизирует максимальную ошибку во всем интервале аппроксимации. rd алгоритм вычисления таких аппроксимаций является алгоритмом Ремеза. Оценка осуществляется по схеме Хорнера; численная точность этой оценки повышается за счет использования fmaf().

Эта стандартная математическая функция реализует то, что известно как сплавленное многократное добавление или FMA. Это вычисляет a*b+c, используя полный продукт a*b во время добавления и нанесения единственного округления в конце. На большинстве современных аппаратных средств, таких как GPU, IBM Power CPU, последние процессоры x86 (например,Haswell), последние процессоры ARM (в качестве дополнительного расширения), это прямое указание на аппаратную инструкцию. На платформах, которым не хватает такой инструкции, fmaf() будет отображать довольно медленный код эмуляции, и в этом случае мы не хотим использовать его, если нас интересует производительность.

Окончательное вычисление представляет собой умножение на 2 ⁱ, для которого C и C++ предоставляют функцию ldexp(). В коде библиотеки «промышленная прочность» обычно используется машинная специфика, в которой используется двоичная арифметика IEEE-754 для float. Наконец, код очищает случаи переполнения и переполнения.

x86 FPU внутри процессоров x86 имеет инструкцию F2XM1, которая вычисляет 2 ^x -1 на [-1,1]. Это можно использовать для второго этапа вычисления exp() и exp2(). Существует инструкция FSCALE, которая используется для умножения на 2 ⁱ на третьем шаге. Обычный способ реализации F2XM1 сам по себе является микрокодом, который использует рациональное или полиномиальное приближение. Обратите внимание, что в настоящее время x90 FPU поддерживается в основном для поддержки старых версий. На современных платформах платформы x86 обычно используются чистые программные реализации на основе SSE и алгоритмы, аналогичные тем, которые показаны ниже. Некоторые объединяют небольшие таблицы с полиномиальными приближениями.

pow(x,y) может быть концептуально реализован как exp(y*log(x)), но страдает от значительной потери точности при x близка к единице и y в большом по величине, а также неправильное обращение многочисленных особых случаев, указанная в C/C++ стандартов C. Один из способов обойти проблему точности - вычислить log(x) и продукт y*log(x)) в некоторой форме расширенной точности. Детали заполнили бы полный, длинный отдельный ответ, и у меня нет кода, который бы мог его продемонстрировать. В различных математических библиотеках C/C++ pow(double,int) и powf(float, int) вычисляются с помощью отдельного пути кода, который применяет метод «квадрат-и-умножить» с битовым сканированием двоичного представления целочисленной экспоненты.

#include <math.h> /* import fmaf(), ldexpf() */ 

/* Like rintf(), but -0.0f -> +0.0f, and |a| must be < 2**22 */ 
float quick_and_dirty_rintf (float a) 
{ 
    float cvt_magic = 0x1.800000p+23f; 
    return (a + cvt_magic) - cvt_magic; 
} 

/* Approximate exp(a) on the interval [log(sqrt(0.5)), log(sqrt(2.0))]. 
    Maximum ulp error = 0.70951 
*/ 
float expf_poly (float a) 
{ 
    float r; 

    r =    0x1.6ab95ep-10f; 
    r = fmaf (r, a, 0x1.126d28p-07f); 
    r = fmaf (r, a, 0x1.5558f8p-05f); 
    r = fmaf (r, a, 0x1.55543ap-03f); 
    r = fmaf (r, a, 0x1.fffffap-02f); 
    r = fmaf (r, a, 0x1.000000p+00f); 
    r = fmaf (r, a, 0x1.000000p+00f); 
    return r; 
} 

/* Approximate exp2() on interval [-0.5,+0.5] 
    Maximum ulp error = 0.79961 
*/ 
float exp2f_poly (float a) 
{ 
    float r; 

    r =    0x1.4308f2p-13f; 
    r = fmaf (r, a, 0x1.5f0722p-10f); 
    r = fmaf (r, a, 0x1.3b2bd4p-07f); 
    r = fmaf (r, a, 0x1.c6af80p-05f); 
    r = fmaf (r, a, 0x1.ebfbdep-03f); 
    r = fmaf (r, a, 0x1.62e430p-01f); 
    r = fmaf (r, a, 0x1.000000p+00f); 
    return r; 
} 

/* Approximate exp10(a) on [log(sqrt(0.5))/log(10), log(sqrt(2.0))/log(10)]. 
    Maximum ulp error = 0.77879 
*/ 
float exp10f_poly (float a) 
{ 
    float r; 

    r =    0x1.a65694p-3f; 
    r = fmaf (r, a, 0x1.158a00p-1f); 
    r = fmaf (r, a, 0x1.2bda9ap+0f); 
    r = fmaf (r, a, 0x1.046f72p+1f); 
    r = fmaf (r, a, 0x1.535248p+1f); 
    r = fmaf (r, a, 0x1.26bb1cp+1f); 
    r = fmaf (r, a, 0x1.000000p+0f); 
    return r; 
} 

/* Compute exponential base e. Maximum ulp error = 0.88790 */ 
float my_expf (float a) 
{ 
    float t, r; 
    int i; 

    t = a * 0x1.715476p+0f;   // 1/log(2) 
    t = quick_and_dirty_rintf (t); 
    i = (int)t; 
    r = fmaf (t, -0x1.62e430p-1f, a); // log_2_hi 
    r = fmaf (t, 0x1.05c610p-29f, r); // log_2_lo 
    t = expf_poly (r); 
    r = ldexpf (t, i); 
    if (a < -105.0f) r = 0.0f; 
    if (a > 105.0f) r = 1.0f/0.0f; // +INF 
    return r; 
} 

/* Compute exponential base 2. Maximum ulp error = 0.87922 */ 
float my_exp2f (float a) 
{ 
    float t, r; 
    int i; 

    t = quick_and_dirty_rintf (a); 
    i = (int)t; 
    r = a - t; 
    t = exp2f_poly (r); 
    r = ldexpf (t, i); 
    if (a < -152.0f) r = 0.0f; 
    if (a > 152.0f) r = 1.0f/0.0f; // +INF 
    return r; 
} 

/* Compute exponential base 10. Maximum ulp error = 0.95588 */ 
float my_exp10f (float a) 
{ 
    float r, t; 
    int i; 

    t = a * 0x1.a934f0p+1f; 
    t = quick_and_dirty_rintf (t); 
    i = (int)t; 
    r = fmaf (t, -0x1.344136p-2f, a); // log10(2)_hi 
    r = fmaf (t, 0x1.ec10c0p-27f, r); // log10(2)_lo 
    t = exp10f_poly (r); 
    r = ldexpf (t, i); 
    if (a < -46.0f) r = 0.0f; 
    if (a > 46.0f) r = 1.0f/0.0f; // +INF 
    return r; 
}