§ Разбиение на части

В этой статье мне бы хотелось рассказать о том, как сделать умножение двух Half Precision Float Point (назвал FP16) так, чтобы все работало корректно и без ошибок. Для начала, необходимо определиться с модулем, его входами и выходами.
1module fpalu_mul
2(
3    input    [15:0]  a, // Операнд A
4    input    [15:0]  b, // Операнд B
5    output   [15:0]  c  // Результат C
6);
7endmodule
Здесь только 2 входа — операнды A и B, и результат C. Конечно же, на вход и выход подаются 16-битные float point. Про формат этих чисел я рассказывал в этой и этой статьях.
Snimok_ekrana_ot_2023-07-31_08-49-15.png
В выше приведенной картинке показана структура битов у числа Half Precision (FP16). Для более простого доступа я создам именованные wire для каждого битового ранга.
1wire        sa = a[   15], sb = b[   15]; // Знак +/-
2wire [4:0]  ea = a[14:10], eb = b[14:10]; // Порядок 5 бит
3wire [9:0]  ma = a[  9:0], mb = b[  9:0]; // Значение мантиссы 10 бит
Теперь у нас есть именованные провода для дальнейших операции. Теперь же надо проверить числа на денормализованность. Это делается достаточно просто, необходимо, чтобы в порядке (ea или eb) были нули.
1wire da = (ea == 1'b0),
2     db = (eb == 1'b0);
Таким образом, если в da или в db находится 1, то число является денормализованным.

§ Процедура умножения

Самая основная процедура. Формируются два числа, у которых в скрытой части может быть либо 1, либо 0.
1wire [21:0] m0 = {!da, ma} * {!db, mb};
Как можно догадаться, умножаются два числа вида:
1.mmmmmm или 1.mmmmmm или 0.mmmmmm или 0.mmmmmm
1.mmmmmm     0.mmmmmm     1.mmmmmm     0.mmmmmm
Получается 22-х битное число, где старшие 2 бита m[21:20] — целая часть (от 0 до 3), а биты m[19:0] — мантисса.
В целой части может получаться число от 0 до 3, и это тоже необходимо обработать.
  • Если в целой части получилось число 00, это значит, что результат умножения меньше 1, что значит, что один из операндов был денормализованным числом
  • Если получилось 01, то оба операнда были нормализованными числами и результат лежит в диапазоне от 1 до 1.(1)
  • Если 10 или 11 — результат больше или равен 2, необходимо скорректировать итоговую степень.
Известно правило умножения степеней с одинаковыми основаниями 2^x 2^y = 2^{x+y} . Известно, что если степень равна 15 (bias), то на самом деле, это степень 0, так что для нашего случая запишем умножение степеней так:
2^{x-15} 2^{y-15} = 2^{x-15+y-15}
При этом, необходимо не забыть добавить +15 обратно к результату. То есть, чтобы вычислить новую степень, необходимо сделать следующее:
1wire [5:0]  e1 = ea + eb - 15;
То есть, e1 = (ea - 15) + (eb - 15) + 15 = ea + eb - 15. В случае, если старший бит e1[5] получается не равным 0, то это значит, что есть переполнение, либо в положительную сторону, либо в отрицательную.
  • Превышение в положительную часть. Максимальный ea=30, eb=30, тогда 30 + 30 - 15 = 45. Здесь 45 <= 47 (максимальное возможное)
  • Минимальные ea=0, eb=0, 0 + 0 - 15 = 48, так что это определяется как переполнение в отрицательную сторону (т.к. 48 > 47)
Поскольку, помимо нормализованных чисел, существуют так же денормализованные, то степень -15 не существует, а только степень -14, так что, при сложении итоговых степеней, обязательно требуется скорректировать ее, при получении денормализованного числа.
1wire [5:0]  e1 = (ea + da) + (eb + db) - 15;
Поскольку, если ea = 0, то на самом деле, он должен быть равен 1, так что добавляется +da к нему, также как и eb + db.
Перейдем к следующей ситуации. Допустим, при умножении получилось в целой части число равное 2 или 3. Чтобы нормализовать число, потребуется сместить его на 1 бит вправо, и добавить в степень +1, чтобы восстановить порядок.
1wire [20:0] m1 = m0[21] ? m0[21:1] : m0[20:0];
2wire [ 5:0] e1 = (ea + da) + (eb + db) + m0[21] - 15;
Теперь в m1 находится число, не превышающее 1.(1), а также в e1 добавлена коррекция +1, если в старшем бите m0[21] появился бит 1. Если бы мы работали только с нормализованными числами, то процедура была бы почти закончена. Но не все так просто, нужно учесть работу с денормализованными числами.

§ Нормализация числа

Итак, по итогу, есть результат в m1 и e1. Рассмотрим такой случай, когда в итоговом результате вышло не 1.mmmm, а 0.mmmm, например m1=0.00011010101011111010. Здесь я поставил точку только для того, чтобы отделить целые от дробных. Как видно, полученное число можно вполне себе нормализовать:
0.00011010101011111010
1.1010101011111010____
Чтобы это сделать, потребуется 4 сдвига влево. Соответственно, из итоговой степени необходимо будет вычесть 4. Но, для того чтобы найти сколько надо вычитать, необходимо найти количество лидирующих нулей. Здесь количество таких нулей равно 4. Если же в целой части числа будет находится 1, то количество лидирующих нулей будет равно 0, то есть, двигать никуда не надо. Сначала я переделаю получение m1 с 21-х битного на 32-х битное.
1wire [31:0] m1 = {m0[21] ? m0[21:1] : m0[20:0], 11'b0000000000};
В старших битах находится результат. В младших — нули, для случая, если получится 0 в мантиссе. Я взял 32х битное число потому, что буду использовать бинарный поиск лидирующих нулей.
Алгоритм поиска довольно просто. Разбиваем 32х битное число на 2 части. Если в одной части 16 нулей (старшей), то ставится 1 в результат и двигается число вправо на 16 бит. Если нет — то ничего не трогаем и передаем на следующий этап. Там проверяем уже 8 битов, потом 4, потом 2, и в конце — 1 бит. Вот и весь алгоритм. Результатом будет 5-битное поле, исходное же число будет нормализовано. Если везде будут 0, то итоговое число будет равно 0, а количество сдвигов — 31, максимальное.
1wire [4:0]  np = {
2    ~|m1[31:16], // 4
3    ~|m2[31:24], // 3
4    ~|m3[31:28], // 2
5    ~|m4[31:30], // 1
6     ~m5[31]     // 0
7};
8
9wire [31:0] m2 = np[4] ? {m1[15:0], 16'b0} : m1; // 16 бит
10wire [31:0] m3 = np[3] ? {m2[23:0],  8'b0} : m2; // 8 бит
11wire [31:0] m4 = np[2] ? {m3[27:0],  4'b0} : m3; // 4 бит
12wire [31:0] m5 = np[1] ? {m4[29:0],  2'b0} : m4; // 2 бит
13wire [31:0] m6 = np[0] ? {m5[30:0],  1'b0} : m5; // 1 бит
К завершению всех стадии m6 будет содержать нормализованное число (если это возможно), а np — количество лидирующих нулей.
С учетом подобной нормализации, необходимо будет также поправить итоговый порядок:
1wire [ 6:0] e1 = (ea + da) + (eb + db) + m0[21] - 15 - np;
Тут можно отметить следующее. Возможно появление нулевого или отрицательного e1, когда при нормализации числа появилась степень, уходящая в денормализованное число. Если это произошло, то итоговая степень становится равной 0, а лишние биты скрываются путем сдвига их направо.

§ Финализация

Сначала проведем провод dn, где будет 1, если e1 меньше или равен 0.
1wire dn = e1[6] || e1 == 0;
Поэтому, если получилась степень 0, то тогда надо сдвинуть на 1 вправо, если степень -1, то тогда на 2 вправо двигаем и так далее. Здесь за количество сдвигов при получении денормализованного отвечает провод e2.
1wire [ 5:0] e2 = 1 - e1;
2wire [ 9:0] rm = dn ? m6[31:21] >> e2 : m6[30:21];
Так как скрытая единица не участвует, то она не будет показана при dn=0, результат берется без нее в m6[30:21]. В случае денормализованного числа, скрытая единица "переплывает" в мантиссу. Я намеренно сделал поле m6[31:21] именно 11-битным, а не 10 битным, чтобы не потерять эту единицу с 31-го бита при сдвиге. Она там все равно появится, если была.
Единственный случай, где скрытой единицы не будет — это результат 0.
1wire        vp = e1[6:5] == 2'b01 || e1 == 31 || &ea || &eb;
2assign c = vp ? {sn, 5'h1F, 10'b0} : {sn, dn ? 5'b00000 : e1[4:0], rm};
Вот собственно, и весь результат. При получении dn=1, в экспоненте записывается 0, иначе степень от -14 до +14. Результирующая степень не больше чем 64, но и не меньше 31 (должно быть меньше), то это означает переполнение сверху, потому результату устанавливается +/- INF. Также, если одно из чисел, входящих в операнды, равно inf, то результат будет тоже inf.

§ Функция перевода FP16 <=> FP32

Однажды, я подумал, что неплохо бы написать пару функции для тестирования перевода чисел. Вот функция перевода числа FLOAT32 в 16-битное.
1uint16_t getfp16(float a) {
2
3    // Определить знак
4    int s = a < 0 ? 1 : 0;
5    a = a < 0 ? -a : a;
6
7    // Вычисление экспоненты и мантиссы
8    int e = 15;
9    int m = 0;
10
11    // Нормализация мантиссы
12    if      (a == 0) { e = 0; }
13    else if (a >= 2) { while (a >= 2) { a /= 2; e++; } }
14    else if (a <  1) { while (a <  1) { a *= 2; e--; } }
15
16    // Итоговая мантисса
17    m = ((int)(a * 1024)) & 0x3FF;
18
19    // Денормализация числа
20    if (e <= 0 && a != 0) { m = (0x400 + m) >> (1 - e); e = 0; }
21
22    // Число бесконечности
23    if (e > 30) { m = 0; e = 31; }
24
25    return m + 1024*e + 32768*s;
26}
И наоборот тоже можно:
1float getfp32(uint16_t a) {
2
3    int s = (a & 0x8000);
4    int p = (a & 0x7C00) >> 10;
5    int m = (a & 0x03FF);
6
7    float fm =     ((float)m / 1024) + (p ? 1.0 : 0.0);
8    float em = p ? ((float)p - 15.0) : -14.0;
9
10    return (s ? -1.0 : 1.0) * pow(2.0, em) * fm;
11}
Код АЛУ для верилоге приложен тут https://cdn.neurofox.ru/ov/kh4slsym8c-fpualu.v.zip