AiratGl Jun 5 at 12:50

Векторы в RISC-V на практике: вычисление softmax

Medium

13 min

1.4K

Tutorial

С 10 по 14 апреля 2025 года прошел первый онлайн RISC-V хакатон, организованный Ассоциацией RISC-V. Участникам на выбор давались 2 задачи. Одна задача от Codasip -доработать программу и кастомный процессор для вычисления LLM трансформера. Другая от Andes - улучшить вычисление функции softmax. Для демонстрации работы векторного расширения RISC-V задача с softmax мне показалась более подходящей.

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

Функция softmax

Функция softmax на вход принимает набор дискретных значений $\mathbf{X}$ и на выходе выдает функцию плотности вероятности. Выход функции softmax нормализован и в сумме дает 1, как и должна любая плотность вероятности. Используется функция в машинном обучении, в различных классификаторах, в частности в нейронных сетях, чтобы нормализовать выход, и чтобы он показывал вероятность принадлежности объекта к одному из классов.

Демонстрация применения функции softmax к выходу свёрточной нейронной сети для классификации. Источник.

Функция softmax определяется с помощью экспоненциальной функции следующим образом:

$\sigma(\mathbf{X})_i = \frac{e^{x_i}}{\sum\limits_{j=0}^{N}e^{x_j}}$

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

При вычислении на компьютере и использовании стандартных типов float или double для чисел с плавающей точкой необходимо помнить о стабильности вычислений. Возведение в степень какого-то большого числа может привести к “переполнению” и результат будет “бесконечностью”, специальным ненормальным float значением inf. Это испортит все дальнейшие вычисления, т.к. результат арифметических операций с inf будет либо inf, либо другим ненормальным значением NaN.

Чтобы сохранить стабильность вычисления softmax из входных данных вычитают максимальное значение. Так диапазон входных данных будет не положительным и экспоненты будут ограничены сверху 1-ей.

Это также поможет с проблемой деления на 0 и на близкие к 0 значения. Т.к. после вычитания максимума одно из значений входных данных будет 0 и экспонента от этого значения 1, то в знаменателе функции softmax будет значение не меньше 1. При этом вычитание максимума не влияет на результат, благодаря свойствам экспоненциальной функции.

$\sigma(\mathbf{X})_i = \frac{e^{x_i-\max(\mathbf{X})}}{\sum\limits_{j=0}^{N}e^{x_j-\max(\mathbf{X})}} = \frac{e^{x_i}e^{-\max(\mathbf{X})}}{e^{-\max(\mathbf{X})}\sum\limits_{j=0}^{N}e^{x_j}} = \frac{e^{x_i}}{\sum\limits_{j=0}^{N}e^{x_j}}$

В железе обычно есть только инструкции для простой арифметики: сложения, вычитания, умножения, деления, и инструкции битовых сдвигов. Для вычисления сложных функций вроде экспонент, логарифмов и тригонометрии используются различные аппроксимации. Подходов к приближенным методом множество, например такие как CORDIC метод, таблицы, аппроксимация полиномами и использование битового представление IEEE 754 float чисел [Muller, 2020].

В этой статье аппроксимация функции softmax будет основана на комбинации метода с битовыми манипуляциями с float числами и с дополнительными поправками с помощью метода аппроксимации полиномами.

Трюки с битами IEEE 754 float

Этот способ приближенного вычисления функций использует особенности представления чисел с плавающей точкой по стандарту IEEE 754. Из-за этих особенностей способ хорошо работает с функциями, связанными с экспонентой. Небезызвестный быстрый обратный корень (fast inverse square root) из движка Quake 3 как раз использует этот способ вычисления.

По стандарту IEEE 754 биты в записи 32-битного float числа определены подобно экспоненциальной форме записи чисел. Число состоит из бита знака, экспоненты и мантиссы. Мантисса - дробная часть числа, занимает биты 0..22. Поскольку мантисса принимает значения от 1 до 2, единицу в начале можно не хранить и значение непосредственно записанное в битах будет в диапазоне от 0 до 1. Мантисса закодирована как дробное двоичное число, т.е. бит 22 - половины $2^{-1}$ , бит 21 - четверти $2^{-2}$ и так далее. Экспонента лежит в битах 23..30 и значение в поле экспоненты минус 127 равно показателю степени двойки, на которую умножается мантисса. Значение экспоненты 255 специальное и предназначено для представления бесконечностей и NaN значений. А если экспонента 0, то представляются субнормальные числа, и мантисса начинается с 0, а не с 1. Значение 31-ого бита указывает на знак числа, если он в 1, то число отрицательное.

Представление 32-битного float числа в битах. Источник.

В виде формулы значение float числа будет равно:

$x = (-1)^s 2^{e - 127} (1 + m)$

где - значение бита знака, - значение поля экспоненты, - значение поля мантиссы.

Компиляторы С реализуют float тип как раз как 32-битный IEEE 754 формат. Но манипуляции с битами над объектами типа float в C запрещены. Битовые сдвиги и побитовые AND, OR, NOT, XOR можно производить только над переменными типа integer. В C (но не в C++) привести битовое представление float числа к целочисленному 32-битному числу и обратно - прочитать последовательность бит целого число как float число - можно с помощью union-а.

typedef union {
	int32_t i32;
	float f32;
} u_float_int_t;

Благодаря такому представлению float чисел вычисление экспоненциальной функции можно сделать довольно легко. Сперва представим float число как сумму целой и дробной частей , причем дробная часть лежит в диапазоне от 0 до 1.

В C выделение целой и дробной частей можно сделать так:

int x_i = (int) x;
float x_f = x - x_i + (x < 0) ? 1 : 0;

Округление float числа к int через преобразование типа производится отбрасыванием дробной части. Поэтому для положительных x выражение x - x_i будет лежать в диапазоне от 0 до 1, а для отрицательных чтобы получить дробную часть в нужном диапазоне надо добавить 1. Т.к. в поставленной задаче вычисления softmax входные данные будут не положительны после вычитания максимума для стабильности вычислений, далее в статье для взятия дробной части всегда будет добавляться 1.

Вычисление двойки в степени x тогда можно также разбить на вычисление двух экспонент с основанием 2: $2^x = 2^{x_i} 2^{x_f}$ .

Вычислить первую двойку в степени целого числа очень легко. Достаточно прибавить к 127, и сдвинуть побитово влево на 23 бита (что эквивалентно умножению на ), а затем прочитать получившуюся последовательность бит как float число.

$2^{x_i} \xleftarrow{(float)} (1 \ll 23)\cdot(x_i + 127)$

Биты мантиссы результата $2^{x_i}$ после этого вычисления будут заполнены нулями. Чтобы вычислить надо заполнить эту мантиссу каким-то значением, которое зависит от дробной части. Из формулы для float чисел выходит, что мантисса должна быть равна:

$1 + m = 2 ^{x_f} \\ m = 2^{x_f} - 1$

Введем функцию поправки , которую нужно прибавить к мантиссе $2^{x_i}$ , чтобы учесть дробную часть.

$2^x = 2^{x_i} 2^{x_f} = \\ = 2^{e-127}(1+m) \xleftarrow{(float)}(1 \ll 23) \cdot (x_i + 127) + g(x_f)$

Битовая последовательность значения такой поправки должна совпадать с битовой последовательностью в мантиссе результата вычисления .

$1 + m \xrightarrow{(int)} g(x_f)$

Т.е. поправка равна целому числу, которое получится, если взять биты мантиссы и интерпретировать их как целое число. Поскольку мантисса это дробное двоичное число и последний бит мантиссы соответствует степени двойки $2^{-23}$ , то это преобразование эквивалентно умножению float числа равного на $2^{23}$ и округления до целого.

$g(x_f) = \text{trunc}((1 \ll 23)\cdot(1+m)) = \text{trunc}((1 \ll 23)\cdot2^{x_f})$

Можно внести эту поправку перед умножением на $2^{23}$ . Тогда можно совместить её с дробной частью в виде функции $\delta(x_f)$ .

$(1 \ll 23)(x_i +127) + g(x_f) = (1 \ll 23)(x_i +127) + (1 \ll 23)2^{x_f} = \\ =(1 \ll 23)(x _i + 2^{x_f} + 127) = (1 \ll 23) (x - x_f - 1 + 2^{x_f} +127) = \\ =(1 \ll 23)(x - \delta(x_f) + 127) \\ \delta(x_f) = 1 + x_f - 2^{x_f}$

Для вычисления $2^{x_f}$ придется заменить какой-то аппроксимацией. К счастью, ограничен от 0 до 1, поэтому $2^{x_f}$ будет хорошо аппроксимироваться полиномом, а вычислять полиномы на процессоре с операциями сложения и умножения целых чисел и float чисел уже скорее всего возможно. Если целевой процессор архитектуры RISC-V, то для арифметики с целыми числами процессор должен реализовать расширение с базовыми инструкциями RV32I/RV64I. Для инструкций арифметики с float числами нужно расширение “F”.

Воспользуемся для примера полином 5-ой степени для аппроксимации функции $\delta(x_f)$ [Perini, 2018]:

$\delta(x_f) = 1 + x_f - 2^{x_f} = s_5 x^5 + s_4 x^4 + s_3 x^3 + s_2 x^2 + s_1 x + s_0$

Найти коэффициенты полинома можно приравняв значение функции и производных на концах диапазона . Для полинома 5 степени нужно найти 6 коэффициентов, значит нужно 6 уравнений. Для этого достаточно будет 2-x производных.

$\delta’(x_f) = 1 - \ln2 \cdot 2^{x_f} = 5 s_5 x^4 + 4 s_4 x^3 + 3 s_3 x^2 + 2 s_2 x + s_1 \\ \delta’’(x_f) = - (\ln2)^2 2^{x_f} = 20 s_5 x^3 + 12 s_4 x^2 + 6 s_3 x + 4 s_2$ $\delta(0) = 0 = s_0 \\ \delta’(0) = 1 - ln2 = s_1 \\ \delta’’(0) = - (\ln2)^2 = 2s_2 \\ \delta(1) = 0 = s_5 + s_4 + s_3 + s_2 + s_1 +s_0 \\ \delta’(1) = 1 - 2\ln2 = 5s_5 + 4s_4 +3s_3 +2s_2 + s_1 \\ \delta’’(1) = -2 (\ln2)^2 = 20 s_5 + 12 s_4 + 6 s_3 + 4 s_2 \\$

Решение этой системы уравнений дает коэффициенты полинома:

   const float s1 = 3.06852819440055e-1f;
   const float s2 = -2.40226506959101e-1f;
   const float s3 = -5.57129652016652e-2f;
   const float s4 = -9.01146535969578e-3f;
   const float s5 = -1.90188191959304e-3f;

Вычисление значения полинома 5-ой степени требует много умножений. Можно пожертвовать точностью и воспользоваться полиномом 3-ей степени или использовать наиболее простую для вычислений константную поправку. Константную поправку можно оптимизировать под разные критерии точности. Для сравнения на графике ниже даны сравнения точности экспоненциальной функции для полиномов 5-ой и 3-ей степеней, коэффициенты которых получены уравниванием производных и константная поправка, дающая наименьшую среднеквадратичную ошибку [Schraudolph, 1999].

Точность аппроксимации экспоненциальной функции с различными поправками.

Скалярная реализация в C

Используя приведенные выше математические выкладки можно реализовать аппроксимакцию softmax функции в простом C. Возведение двойки в степень реализованное на C выглядит так:

static float fast_pow2(const float in) {
   const float yf = in - (int32_t)(in) + 1;
   u_float_int_t value;
   value.i32 = (int32_t)((1 << 23) * (in - delta(yf) + 127.0f));
   return value.f32;
}

Где delta(x) вычисляет полином пятой степени с указанными выше коэффициентами.

Из экспоненты с показателем 2 легко получить экспоненту с показателем , надо просто поделить входное значение на $\ln2$

Вычисления softmax на C выглядит так:

int32_t scalar_softmax_f32(const float* in_vec, uint32_t size, float* out_vec) {
   float max = in_vec[0];
   float sum = 0;
   for (int i = 1; i < size; ++i) {
       if (in_vec[i] > max) {
           max = in_vec[i];
       }
   }
   for (int i = 0; i < size; ++i) {
       float value = (in_vec[i] - max) * k1OverLn2;
       out_vec[i] = fast_pow2(value);
       sum += out_vec[i];
   }
   for (int i = 0; i < size; ++i) {
       out_vec[i] = out_vec[i] / sum;
   }
   return 0;
}

Вычисления производятся в 3 этапа. Первый этап - поиск максимального значения в данных и вычитание его из данных. Далее вычисление экспоненты от входных данных и суммирование. И последнее действие - деление экспоненты на сумму экспонент от других элементов входного вектора.

Эта имплементация скалярная и не использует в явном виде параллелизацию через SIMD. Компилятор может попробовать автовекторизовать, но вручную получится лучше, что и будет сделано далее для архитектуры RISC-V с его реализацией SIMD в векторном расширении “V”.

RVV имплементация с интринзиками

Основной особенностью SIMD в RISC-V это использование идеологии векторного процессора. В векторном процессоре количество элементов обрабатываемых за одну инструкцию определяется динамически. Это количество обозначается как vl и для выставления значения vl используется инструкция vsetvl. Максимальное возможное значение vl - vlmax - определяется размером векторных регистров VLEN, размером обрабатываемых элементов SEW и множителем количества векторных регистров на один вектор LMUL. Значение VLEN определяется имплементацией процессора и повлиять на него мы не можем. SEW равен 32 битам, потому что в поставленной задаче softmax вычисляется для 32-битных float данных. А LMUL можно задать программно, инструкцией vsetvl. Поскольку векторный процессор эффективнее работает с векторами большей длины, LMUL будет равен наибольшему доступному по спецификации RISC-V значению - 8.

В C компиляторах gcc и clang есть заголовочный файл <riscv_vector.h> с интринзиками для векторных инструкций RISC-V. Это дает удобный инструмент для работы с векторным расширением без необходимости писать ассемблер вручную.

Функция-интринзик vsetvl принимает желаемое количество элементов в векторе и возвращает реально установленный vl. Принцип обработки данных с помощью векторного процессора заключается в разбиении в цикле данных на куски длиной vlmax. Причем оставшийся хвост данных не кратный vlmax обработается естественным образом в последней итерации цикла, без необходимости каких-то дополнительных манипуляций. Потому что vsetvl на значение меньшее vlmax установит vl на это запрашиваемое значение. Такое удобство это ещё одно преимущество векторной архитектуры.

С учетом векторного подхода поиск максимума с помощью векторных инструкций будет выглядеть так:

static float rvv_find_max(const float *in_vec, uint32_t size) {
   float *cur_pos = (float *)in_vec;
   size_t rem_size = size;
   size_t vl = __riscv_vsetvl_e32m8(rem_size);
   vfloat32m8_t vInput = __riscv_vle32_v_f32m8(cur_pos, vl);


   cur_pos += vl;
   rem_size -= vl;
   while (rem_size > 0) {
       vl = __riscv_vsetvl_e32m8(rem_size);
       const vfloat32m8_t vNext = __riscv_vle32_v_f32m8(cur_pos, vl);
       cur_pos += vl;
       rem_size -= vl;
       vInput = __riscv_vfmax_vv_f32m8(vInput, vNext, vl);
   }


   vl = __riscv_vsetvl_e32m1(size);
   vfloat32m1_t vMax = __riscv_vfmv_v_f_f32m1(0.0f, vl);


   vl = __riscv_vsetvl_e32m8(size);
   vMax = __riscv_vfredmax_vs_f32m8_f32m1(vInput, vMax, vl);
   float max = __riscv_vfmv_f_s_f32m1_f32(vMax);
   return max;
}

Векторы в C коде представляются через специальные типы vfloat32m8_t.

Загрузка данных из памяти в векторные регистры осуществляется инструкцией vle. В C коде соответствующая ей функция-интринзик принимает в качестве аргументов указатель на начальную позицию последовательно расположенных элементов данных в памяти и текущее значение vl. Возвращает объект типа vfloat32m8_t который представляет вектор с этими данными в элементах.

В интризиках вшиты тип элементов и LMUL, а vl всегда передается как параметр, поскольку он меняется в процессе работы.

Для поиска максимума используется интринзик __riscv_vfmax_vv_f32m8 соответствующий ассемблерной инструкции. Эта инструкция принимает на вход два вектора и кладет в первый вектор поэлементно максимальное из элементов этих векторов.

После отработки цикла, вектор vInput хранит максимальные поэлементные значения кусков по vlmax изначального вектора и теперь остается найти максимальное значение внутри vInput чтобы найти максимум в данных. Для этого используется инструкция vfredmax свертки вектора применением функции max.

Визуализация поиска максимума с помощью векторного расширения RISC-V.

Аппроксимация экспоненциальной функции будет работать следующим образом:

static vfloat32m8_t rvv_fast_pow2(const vfloat32m8_t vIn, size_t vl) {
   const vfloat32m8_t vFrac = rvv_get_frac(vIn, vl);
   const vfloat32m8_t vDelta = rvv_delta(vFrac, vl);
   const int32_t pow2_23 = 1 << 23;


   vfloat32m8_t vTmp;
   vTmp = __riscv_vfsub_vv_f32m8(vIn, vDelta, vl);
   vTmp = __riscv_vfadd_vf_f32m8(vTmp, 127.0f, vl);
   vTmp = __riscv_vfmul_vf_f32m8(vTmp, pow2_23, vl);
   vint32m8_t vTmp_i32 = __riscv_vmv_v_x_i32m8(0, vl);
   vTmp_i32 = __riscv_vfcvt_rtz_x_tu(vTmp_i32, vTmp, vl);
   const vfloat32m8_t vResult = __riscv_vreinterpret_v_i32m8_f32m8(vTmp_i32);
   return vResult;
}

для взятия дробной части float чисел в векторе можно использовать интринзик __riscv_vfcvt_f_x_v_f32m8 - ассемблерная инструкция конвертации float в int.

vfloat32m8_t rvv_get_frac(const vfloat32m8_t vIn, size_t vl) {
   vint32m8_t vInt = __riscv_vmv_v_x_i32m8(0, vl);
   vInt = __riscv_vfcvt_rtz_x_tu(vInt, vIn, vl);
   const vfloat32m8_t vIntF = __riscv_vfcvt_f_x_v_f32m8(vInt, vl);
   vfloat32m8_t vResult = __riscv_vfsub_vv_f32m8(vIn, vIntF, vl);
   vResult = __riscv_vfadd_vf_f32m8(vResult, 1.0f, vl);
   return vResult;
}

Для интерпретации последовательности бит целого числа как float число необходимо использовать интринзик __riscv_vreinterpret_v_i32m8_f32m8, который не отображается напрямую в ассемблерную инструкцию, но позволяет интерпретировать последовательность бит в элементах вектора с 32-х битными целыми числами как элементы вектора float чисел.

Визуализация аппроксимации экспоненциальной функции с помощью векторного расширения RISC-V.

Имея функцию для вычисления экспоненты можно дописать векторный softmax. В нем присутствуют те же 3 части: поиск максимума, вычитание максимума, масштабирование и вычисление экспоненты, деление результата на сумму экспонент. Поскольку для softmax нужно делить на сумму всех входных значений приходится сохранять промежуточные результаты вычисления экспоненты в память инструкцией vse. Затем промежуточные данные вчитываются, чтобы поделить на сумму. В векторе сумм vSum находятся промежуточные результаты суммирования элементов векторов по vlmax . Вектор затем редуцируется суммой, чтобы получить значение суммы всех элементов.

int32_t rvv_pure_softmax_f32(const float *in_vec, uint32_t size,
                            float *out_vec) {
    float max = rvv_find_max(in_vec, size);


   float *cur_pos = (float *)in_vec;
   float *cur_out_pos = (float *)out_vec;
   size_t rem_size = size;


   size_t vl = __riscv_vsetvl_e32m8(rem_size);


   vfloat32m8_t vSum = __riscv_vfmv_v_f_f32m8(0.0f, vl);


   while (rem_size > 0) {
       vl = __riscv_vsetvl_e32m8(rem_size);
       vfloat32m8_t vData = __riscv_vle32_v_f32m8(cur_pos, vl);


       cur_pos += vl;
       rem_size -= vl;
       // substract max and scale
       vData = __riscv_vfsub_vf_f32m8(vData, max, vl);
       vData = __riscv_vfmul_vf_f32m8(vData, k1OverLn2, vl);


       //--- pow2_f32 begin ---
       vData = rvv_fast_pow2(vData, vl);
       //--- pow2_f32 end ---


       vSum = __riscv_vfadd_vv_f32m8(vSum, vData, vl);  // accumulate


       // save value to memory
       __riscv_vse32_v_f32m8(cur_out_pos, vData, vl);


       cur_out_pos += vl;
   }


   vl = __riscv_vsetvl_e32m1(1);
   vfloat32m1_t vResult = __riscv_vfmv_v_f_f32m1(0.0f, vl);
   vl = __riscv_vsetvl_e32m8(size);
   vResult = __riscv_vfredosum_vs_f32m8_f32m1(vSum, vResult, vl);
   float sum = __riscv_vfmv_f_s_f32m1_f32(vResult);


   // divide the output by sum
   rem_size = size;


   while (rem_size > 0) {
       vl = __riscv_vsetvl_e32m8(rem_size);
       vfloat32m8_t vData = __riscv_vle32_v_f32m8(out_vec, vl);
       vData = __riscv_vfdiv_vf_f32m8(vData, sum, vl);
       __riscv_vse32_v_f32m8(out_vec, vData, vl);
       out_vec += vl;
       rem_size -= vl;
   }

   return 0;
}

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

Сравнение производительности

Не имея на руках процессора с RISC-V и векторным расширением, можно проверить программу кросс-компилированную под RISC-V с помощью эмулятора Spike. Это всего лишь эмулятор, поэтому по времени исполнения надежных выводов не сделать. Но работает по инструкциям, поэтому можно оценить и сравнить количество инструкций, которое понадобится для исполнения скалярной и векторной реализаций softmax.

В моем проекте представлены инструкции как кросс-компилировать код softmax-а под RISC-V и как запустить её в эмуляторе. Spike позволяет использовать псевдо-операционную систему proxy-kernel, что дает возможность программе из эмулятора читать файлы из системы-хоста. Скрипт Jupyther Notebook сгенерирует файл с набором входных данных и файл с эталонными значениями softmax от этих данных.

Для подсчета количества исполненных процессором инструкций при сборке под ОС Linux можно воспользоваться системными вызовами для измерения производительности perf.

В C коде можно определить, для какой архитектуры собирается программа и под какую ОС с помощью предопределенных макросов, полный список которых можно вывести вызвав командуecho | gcc -dM -E - .

Наличие полной стандартной библиотеки под Линукс, включая perf, можно проверить по макросу linux.В прокси ядре Spike-а perf-а нет, но в RISC-V можно посчитать количество выполненных инструкций с помощью специальной инструкции rdinstret, которая возвращает значение внутреннего счетчика выполненных инструкций. Интринзика под эту инструкцию нет, поэтому воспользуемся ассемблерной вставкой.

asm volatile("rdinstret %0" : "=r"(counter));

При сборке скалярной версии компилятору передается опция -march=rv64gc т.е. сборка под 64-битный RISC-V процессор с набором инструкций для общих задач (g - general) и сжатые инструкции (c - compressed). При сборке векторного решения используется архитектура -march=rv64gcv, добавляется векторное расширение V.

В моей имплементации использование векторных инструкций уменьшило количество исполняемых инструкций в 8 раз.

Результат скалярного решения:

spike --isa=rv64gcv pk ./bin/riscv-test -f ./bin/data -g ./bin/golden
Softmax size=2048
rdinstret 532637
rdinstret 616609
Used 83972 instructions
------------------
maxAbsDiff = 0.000000
sumSquareErr = 0.000000
squareGolden = 0.000528
snr_result = 114.316017
f_scores = 106.918671

--> Test PASSED! Your scores = 100

Результат векторного решения:

spike --isa=rv64gcv pk ./bin/rvv-test -f ./bin/data -g ./bin/golden
Softmax size=2048
rdinstret 522514
rdinstret 532074
Used 9560 instructions
------------------
maxAbsDiff = 0.000000
sumSquareErr = 0.000000
squareGolden = 0.000528
snr_result = 115.325989
f_scores = 107.389938

--> Test PASSED! Your scores = 100

P.S.

В прошлой статье рассматривался векторный сопроцессор Ara. В ara есть собственный пример softmax. Аппроксимация экспоненты в нем сделана как в библиотеке https://www.netlib.org/cephes/ и использует аппроксимацию Паде. Будет интересно посмотреть на то, что происходит внутри лейна при исполнении векторных инструкций, но для этого придется создать тестбенч и это задача для будущих статей.

Кроме упомянутых статей, материалов хакатона и спецификаций на RISC-V и RVV мне очень помогла статья от fprox https://fprox.substack.com/p/implementing-softmax-using-risc-v

Hubs: