🔮 Как считать экспоненту с помощью умножения и сложения 🔮

const int L = 1 << 23; // == 2^23
const float a = L / std::log(2);
const float b = L * 127;

float fastExp(float x) {
    x = a * x + b;
    return std::bit_cast<float>((int)x);
}

Выше упрощенная версия кода, который считает приблизительное значение экспоненты, который я недавно встретил на GitHub.

Я был супер удивлен, что тут умножение + сложение + преобразование типов, и всё — экспонента подсчитана.
В моей голове экспонента всегда была нетривиальной операцией, на которую требуется много тактов. Как же так, давайте разбираться.

Целые и дробные числа хранятся в битах по-разному

Числа в формате int хранятся ожидаемым образом, например, для числа 123:
123: int = 0000000000000000000000000 1111011

Числа в формате float хранятся в виде маски битов:
s eeeeeeee mmmmmmmmmmmmmmmmmmmmmmm

s — 1 бит знака. Определяет, положительное число или отрицательное.
e — 8 битов экспоненты. Показывают степень двойки.
m — 23 бита мантиссы. Хранят дробную часть числа.

Само число = (-1)^{s} * 2^{e} * 1.{m}

Например для нашего числа 123:

123: float = 0 00000110 11101100000000000000000

e = 00000110b = 6
m = 0.1110110b = 0.921875

Проверяем: 2^{e} x 1.{m} = 2^6 x 1.921875 = 123. Сходится!

Переводить представления из int в float и обратно несложно, и процессор умеет делать это очень эффективно.

Чтобы понять суть алгоритма, вместо e^{x} будем считать 2^{x}

Считаем быстро 2^{x} если x - int

Заметим, что если мы знаем двоичное представление x, то если мы его засунем в биты экспоненты float, то задача уже решена (e = x, m = 0). Это ключевое наблюдение для понимания алгоритма!

В int представлении все биты числа x у нас хранятся в правом конце числа. Чтобы сдвинуть их в область экспоненты, мы умножим наш int x на 2^23. Тем самым просто допишем 23 нуля справа.

В этих двух идеях и есть вся суть. Как перейти от 2^{x} к e^{х} в других деталях можно разобраться самому или заглянуть в оригинальную статью.

Cчитать экспоненту в десятичной системе мы (человеки) тоже можем с помощью одного умножения (wow!)

Например, чтобы оценить exp(36) мы можем переписать

exp(36) = 10^{36 * log10(e)} = 10^{15.634} ~ 10^16

Эта неожиданная похожесть алгоритмов, получается из-за того что научная запись чисел вроде 5.34e7 очень похожа на то, как хранятся float в компьютере.

P. S. На самом деле степень двойки в экспоненте float сдвинута на 127, чтобы мы могли записывать как большие числа, так и очень маленькие. Мы учитываем этот сдвиг, добавляя b = L * 127. Я убрал его, чтобы он не отвлекал от основной идеи.

P. P. S. Для тех кто всё понял: почему умножение на 2^32 в коде происходит в float, а не в int и что происходит с мантиссой в этом случае?

@lovesyuk