Математический сопроцессор
Вещественные числа и их запись
Несколько слов о представлении вещественных чисел.
«Настоящее» вещественное число не может быть представлено в памяти классического компьютера исключительно как структура данных. Например, иррациональные вещественные числа в десятичном представлении апериодичны и содержат бесконечно много цифр. Рациональные числа можно представить в виде отношения двух целых или в виде конечной дроби с периодом, но такие представления весьма сложны аппаратно, и не решают основной задачи — смоделировать произвольное вещественное число с заданной точностью.
Таким образом, вещественное в памяти компьютера число неизбежно должно моделироваться рациональным и некоторым представлением о точности этой модели. В учебниках вещественное число записывается в виде десятичной дроби — либо в формате с «фиксированной точкой», например 1234.5 или 0.012345, либо в научном (экспоненциальном) формате — с "плавающей точкой", то есть в виде десятичной дроби, содержащей все цифры числа, степени 10 на которую надо эту дробь домножить. В простых случаях эти записи эквивалентны, например, 1234.5 = 1.2345*103, 0.012345 = 1.2345*10-2. Часто эти числа записываются в виде 1,2345·exp103 и 1,2345·exp10-2 , а в тексте программы (и при вводе-выводе) — 1.2345E3 и 1.2345E-2 соответственно. В этой записи состоит из двух частей: мантиссы (в примере 1.2345) и показателя степени — порядка (в примере 3 и -2 соответственно).
Одно и то же число можно записать несколькими способами, в зависимости от того, чему равна мантисса. Например, все эти записи соответствуют одному и тому же числу 155.625:
- 155.625E0
- 1.55625E2
- 0.001555625E5
- 155625E-3
Подобный формат кладется и в основу представления вещественного числа в памяти компьютера. Естественно вместо представления десятичного числа с плавающей точкой используется соответствующее двоичное число с плавающей точкой. Здесь уже начинаются сложности, особенно в вопросе точности представления числа. Допустим, мы храним мантиссу в виде 16-битного числа. То есть точность представления — до 16-го двоичного знака. Но в десятичном представлении 16-битное число — это 0…65535. Стало быть, за пятый десятичный знак мы ручаться не можем, только за четвёртый.
Хранение вещественных чисел в памяти по IEEE 754
Представим число 155.625 в виде двоичного числа с плавающей точкой, для этого разложим заданное число по двоичным разрядам:
155.625 = 1·27 +0·26+0·25+1·24+1·23+0·22+1·21+1·20+1·2-1+0·2-2+1·2-3 155.625 =128 + 0 + 0 + 16 + 8 + 0 + 2 + 1 + 0.5 + 0 + 0.125 155.62510 = 10011011.1012 - число в десятичной и в двоичной системе с плавающей точкой
Полученное число в экспоненциальном виде в двоичной системе будет иметь вид: 1.55625·exp10+2 = 1.0011011101·exp2+111
То есть его мантисса имеет вид 1.0011011101, а порядок — exp2= +111
В действительности все несколько хитрее. Согласно стандарту IEEE_754-2008, в памяти вещественные числа одинарной точности представляются в следующем виде:
S[ E ][ M ] 01110101111100010110101110101111
- S - бит знака
E - смещенный порядок двоичного числа; для 32-битного представления — 8 битов
M - остаток мантиссы двоичного числа с плавающей точкой, для 32-битного представления — 23 бита
Понятно, что число может быть либо неотрицательным, либо отрицательным, поэтому один бит на знак нам придется использовать, и это самый старший бит. 0 - число неотрицательное, 1 - число отрицательное (область S на картинке); никакого дополнительного кода, как при хранении целых чисел, не используется — от него никакого толку.
В IEEE 754 порядок (который может быть и положительным, и отрицательным, и нулём) хранится в виде всегда неотрицательного числа, для чего к нему добавляют «смещение» размером в половину допустимого диапазона. Например, в стандарте на 32-разрядное целое порядку отведено 8 битов (область E), в которых можно хранить числа 0…255. Допустимый порядок числа, таким образом, колеблется от -127 до 127, а при записи в ячейку к нему добавляется 127. В примере порядок = +7 (+111 в двоичной), соответственно, смещённый порядок = 7+127=134 ( 10000110 ). А если бы порядок был -7 , то смещенный порядок = 127-7 =120. Чтобы узнать истинный порядок, надо вычесть 127 из содержимого области E.
Если вспомнить, что «вещественные» числа IEEE 754 — на самом деле целочисленные двоичные дроби, каждое число можно представить себе в виде мантиссы, сдвинутой влево на соответствующее порядку число разрядов. Сдвиг на 0 разрядов будет соответствовать наименьшему по модулю числу, а сдвиг на 254 — наибольшему. Таким образом, для сравнения и проведения арифметических операций это представление наиболее удобно.
В 32-битном представлении на мантиссу (область M) приходится 23 бита. Это обеспечивает довольно низкую точность: 223=8388608, все числа, по модулю большие, будут терять значимые знаки до запятой. В IEEE 754 используется хитрость, позволяющая сэкономить ещё один бит мантиссы:
Если число достаточно большое, оно представляется в нормализованном виде: мантисса обязана быть в диапазоне от 1 до 2 (не включая 2). Двоичное представление позволяет для любого допустимого числа подобрать нормализованную мантиссу и соответствующий порядок. Так вот, в нормализованном виде старший бит мантиссы (а это всегда 1) не хранится вообще!
Близкие нулю вещественные числа используются часто, например, для задания точности или шага. Если число настолько мало, что в его нормализованном представлении порядок состоит одних нулей (-127), оно хранится в денормализованном виде: мантисса при этом обязана быть не больше 0.5. Старший бит такой мантиссы (всегда 0) так что он тоже не хранится. Порядок денормализованного числа на 1 больше порядка нормализованного , следовательно, в денормализованном виде можно хранить в два раза меньшее число, чем в нормализованном.
- В таком представлении оказывается, что ячейка, содержащая все нулевые биты, оказывается и целочисленным нулём, и нулём IEEE 754!
Итак, в нормализованном и в денормализованном виде нет смысла записывать старший бит мантиссы в отведенные 23 бита, и по этой причине в отведенные 23 бита записывают «остаток от мантиссы». В случае 155.625 остаток от мантиссы будет 00110111010000 0000 0000.
Число |
1 бит |
8 бит |
23 бит |
Шестнадцатеричное |
|
знак числа |
смещённый порядок |
мантисса |
|
155.625 |
0 |
10000110 |
00110111010000000000000 |
431BA000 |
-5.23E-39 |
1 |
00000000 |
01110001111001100011101 |
8038f31d |
- Второй пример — это близкое к нулю число, мантисса которого хранится в денормализованом виде.
Вещественные числа двойной точности занимают 64 бита. Поле S в них такое же однобитное, на поле E (порядок) отводится 11 битов, а остальные 52 — на поле M (мантисса). Вещественные числа четверной точности — соответственно, 128 всего, 15 — E и 112 — M.
Дело в том, что 24 бита на мантиссу — это очень мало. Для чисел из обыденного диапазона сохраняется примерно 7-8 десятичных знаков в худшем случае. Например десять миллиардов представлены с точностью плюс-минус полтысячи.
Более подробно смотрите основную статью Стандарт IEEE 754 (в таблице 1 приведены максимальные ошибки) на всякий случай запасена тут
- Подробно о недостатках формата см. ту же статью, §9
Манифест сообщества ненавистников IEEE 754 (разумеется, я к этому отношение не имею, просто спас из вебархива любопытный текст по ссылке из той же статьи☺)
В RARS есть хороший инструмент для изучения числа в формате IEEE 754: «Tools → Floating Point Representation». В частности, легко проверить, что когда в мантиссу записаны одни только нули, это соответствует близкому к нулю числу в денормализованной форме.
Чему соответствует мантисса, состоящая из всех 1? ⇒ NaN:
± 11111111 0 0000000000000000000000 — ±∞
* 11111111 X ********************** в остальных случаях — NaN
- x=0 — «NaN по-тихому», для вычислений, возвращающих NaN
- x=1 — «сигнальный NaN», для исключений
Нелишне заметить, что 0 00000000 0 0000000000000000000000 — это просто 0 (он также называется «+0», единственное вещественное число, совпадающее с его округлением побитно), а 1 00000000 0 0000000000000000000000 — это «-0», и они не равны!
Сопроцессоры
- Различные, почти не пересекающиеся задачи ⇒ сопроцессоры:
- Сопроцессор - FPU
- Другая математика
- Свои регистры
- Почти не смешивается с целочисленной арифметикой
- «Тяжёлые» вычисления
- Сопроцессор 0 - управления (см. следующую лекцию)
RARS -> Tools -> Floating Point Representation
- Сопроцессор - FPU
FPU RISC-V
Расширение F поддерживает числа с плавающей точкой одинарной и двойной точности в IEEE-формате.
Архитектура предусматривает наличие тридцати двух регистров для чисел с плавающей точкой "f0–f31". Это не те же самые регистры, которые мы использовали до сих пор.
Если в архитектуре поддерживается "D" или "Q", для представления в этих регистрах чисел меньшей разрядности используется т. н. «NaN boxing» — заполнение старших частей регистра специальной константой NaN — «Not A Number» (см. ниже)
В соответствие с соглашением о вызове подпрограмм, у этих регистров разное назначение:
Общее название
Мнемоника ABI
Назначение
Переживают ли вызов подпрограммы?
f0 - f7
ft0 - ft7
Временные
Нет
f8 - f9
fs0 - fs1
Сохраняемые при вызове
Да
f10 - f17
fa0 - fa7
Параметры и возвращаемые значения
Нет
f18 - f27
fs2 - fs11
Сохраняемые при вызове
Да
f28 - f31
ft8 - ft11
Временные
Нет
TODO Далее требуется вставить больше мелких примеров из лекции и расписать конспективные утверждения
Инструкции FPU-сопроцессора
Регистры f* не имеют отношения к регистрам x*, но количество их такое же, и значит, они могут встречаться в командах типа S (запись в память, fs*), I (чтение из памяти, fl*) и R (вычисления).
Процессор ничего не знает о формате IEEE754, так что в командах работы с памятью используется терминология «word» / «double word»: flw, fsw, fld, fsd, (а также *q и *h, если они реализованы)
Точность арифметических команд сопроцессора указывается суффиксом инструкции (s, d, q, h)
Арифметические команды: fCMD.P, CMD — мнемоника инструкции, P — точность
- CMD: add, sub, mul, div, sqrt, min, max
- P: s и d (а также q и h)
- Вычисление полусуммы:
Что это за t0 (регистр общего назначения) в псевдоинструкции flw? Посмотрим на раскрытие псевдоинструкции:
0x00400000 0x0fc10297 auipc x5,0x0000fc10 6 flw ft0 a t0 0x00400004 0x0002a007 flw f0,0(x5)
Полный формат flw — flw f-регистр, смещение(x-регистр), а в этом x-регистре (в примере — t0) как раз и формируется адрес слова для загрузки (с помощью уже известного нам auipc)
Заодно освоили ecall вывода вещественного числа
Четырёхместные команды «умносложения» f[n]madd/f[n]msub по формуле a*b+c: новый тип команд R4:
Используются, например, для подсчёта многочленов в форме a0 + a1(x + a2(x + …))
- См. пример ниже
- Можно обмениваться не с памятью, а с регистрами общего назначения
В мнемонике используются два суффикса
Перемещать машинное слово из одного регистра в другой умеет центральный процессор: fmv.s.x и fmv.x.s (для двойной точности такой инструкции нет)
- Не преобразованное целое, лежащее в регистре сопроцессора и вещественное, лежащее в регистре общего назначения, нельзя осмысленно обработать, но можно зачем-то там хранить
Но преобразовывать из вещественного формата в целый и обратно (а также из двойного в одинарный и обратно) умеет только FPU: fcvt.d.s fcvt.s.d, fcvt.P.w[u] и fcvt.w[u].P
(x - 1)² подсчитанный с раскрытием скобок:
В раскрытии псевдоинструкции fcvt присутствует таинственный параметр dyn — это константа, означающая, что используется стратегия округления по умолчанию
Перемещение между f-регистрами fmv f1 f2 — в действительности псевдоинструкция на базе инструкции расширения знака fsgnj.P
Во многих архитектурах (в частности, MIPS) недооценили важность и частоту операции «смены знака по образцу», а вычислительно, как ни странно, эта операция непростая, особенно для вещественных чисел. Поэтому в RISC-V есть специальная операция «взять знак из fA, а мантиссу и порядок из fB, и всё это положить в fC».
Понятно, что тогда fmv раcrрывается так:
0x00400028 0x222100d3 fsgnj.d f1,f2,f2 13 fmv.d ft1 ft2
f{eq|lt|le}.P x1 f1 f2 — записывает 0 или 1 в x1
- остальные сравнения — псевдоинструкции
fclass.P x1 f1 — рассматривает f1 и записывает в x1 результат исследования:
бит
свойства f1
0
f1 — это −∞
1
f1 — это отрицательное нормализованное число
2
f1 — это отрицательное денормализованное число
3
f1 — это −0
4
f1 — это +0
5
f1 — это положительное денормализованное число
6
f1 — это положительное нормализованное число
7
f1 — это +∞.
8
f1 — это сигнализирующий NaN.
9
f1 — это тихий NaN.
Условные операторы и fcsr
Помимо непосредственных сравнений, FPU отображает своё состояние в управляющем регистре fcsr.
Блок управляющих регистров (введение)
- 4K (12 битов) пространство управляющих счётчиков-флагов-масок-прочего
Атомарные инструкции типа I (фактически — обмен значениями между регистром общего назначения и CSR)
csrrw[i], причём часть i тут нестандартная: 12-битное поле IMM занято номером CSR, так что число можно задать совсем маленькое, 5 битов вместо номера регистра-источника
а также csrrs[i] / csrrc[i] включение/выключение битов
- 5 битов флагов
- NX потеря точности
- UF сверхмалое число
- OF переполнение
- DZ деление на 0
- NV недопустимая операция
- 3 бита: тип округления
- RNE ближайшее, лучше чётное
- RTZ ближайшее к нулю
- RDN ближайшее к -∞
- RUP ближайшее к +∞
- RMM ближайшее, лучше с бо́льшим модулем
- …
- …
- DYN умолчания не менять (используется в инструкциях явного округления)
псевдоинструкции f[sr]csr
доступ только к отдельным битам: frm и fflags (это другие адреса CSR!)
отдельные Инструкции frflags и fsrm для атомарного чтения/записи битов только из этого поля
- RARS:
1 .data 2 numb: .float 7.65 3 4 .text 5 flw ft0 numb t0 6 fcvt.w.s a0 ft0 # По умолчанию RNE (ближайшее, чётное) — 8 7 jal outn 8 fcvt.w.s a0 ft0 rtz # Ближайшее к 0 — 7 9 jal outn 10 li t0 2 # RDN — ближайшее к -∞ 11 fsrm t0 12 fcvt.w.s a0 ft0 # теперь по умолчанию RDN — 7 13 jal outn 14 li a7 10 15 ecall 16 17 outn: li a7 1 18 ecall 19 li a0 '\n' 20 li a7 11 21 ecall 22 ret
В этом примере используется округление по умолчанию; в сторону 0 (явно); в сторону минус бесконечности (изменяется умолчание с помощью fsrm);
Условные операторы
В любом случае условный переход при сравнение вещественных неатомарен, потому что сравнение делает математический сопроцессор, а анализ результата и вычисление перехода — арифметико-логическое устройство центрального процессора:
fCMP.P / fclas.P:
- Сравнение
- Условный переход по содержимому x-регистра (0 / не 0)
fcsr*/frflags:
- Чтение всех флагов в x-регистр
Выделение конкретного флага с помощью инструкции andi, например
- Условный переход по содержимому x-регистра (0 / не 0)
Почему нет инструкции «прочесть флаг № такой-то»?
- (предположительный!) ответ: придётся делать AND — а это целая операция
- ⇒ сделать её явно одной инструкцией лучше, чем «подпольно» (противоречит RISC)
- ⇒ аппаратная и программная оптимизации тоже будут за это благодарны
Пример. Неравенство треугольника (функция проверки теперь совсем простая — это в чистом виде само неравенство):
1 .data
2 yes: .asciz "Is a triangle\n"
3 no: .asciz "Not a triangle\n"
4 .text
5 jal input
6 fmv.d fs2 fa0
7 jal input
8 fmv.d fs1 fa0
9 jal input
10 fmv.d fa1 fs1
11 fmv.d fa2 fs2
12 jal check
13 bnez a0 true
14 la a0 no
15 b output
16 true: la a0 yes
17 output: li a7 4
18 ecall
19 li a7 10
20 ecall
21
22 .data
23 prompt: .ascii "Enter triangle side: "
24 .text
25 input: la a0 prompt
26 li a7 4
27 ecall
28 li a7 7
29 ecall
30 ret
31
32 check: fadd.d ft0 fa1 fa2
33 flt.d t0 fa0 ft0
34 fadd.d ft1 fa2 fa0
35 flt.d t1 fa1 ft1
36 fadd.d ft2 fa1 fa0
37 flt.d t2 fa2 ft2
38 and a0 t0 t1
39 and a0 a0 t2
40 ret
Пример. (a + b) / c с предупреждением о точности вычислений:
1 .data
2 Exact: .asciz " , однозначно!\n"
3 Inex: .asciz " , но это не точно…\n"
4
5 .text
6 li a7 6
7 ecall
8 fmv.s fs0 fa0 # A
9 li a7 6
10 ecall
11 fmv.s fs1 fa0 # B
12 li a7 6
13 ecall
14 fmv.s fs2 fa0 # C
15
16 fadd.s ft0 fs0 fs1
17 frflags t0 # Флаги FPU
18 andi t0 t0 1 # Потеря точности?
19 fdiv.s fa0 ft0 fs2
20 frflags t1 # Флаги FPU
21 andi t1 t1 1 # Потеря точности?
22 or s1 t0 t1 # …хотя бы раз
23 li a7 2
24 ecall
25 la a0 Inex
26 bnez s1 out
27 la a0 Exact
28
29 out: li a7 4
30 ecall
31 li a7 10
32 ecall
При C = 0 приезжает флаг DZ (деление на 0), но мы его игнорируем!
Замечание: стоит вспомнить, что равенство вещественных — довольно сомнительное сравнение.
Приведите пример двух равных вещественных чисел, ни один десятичный знак которых не совпадает! (нажмите «Комментарии» в шапке страницы, чтобы прочитать спойлер)
Для эффективности можно хранить значения регистров FPU в регистрах общего назначения, но в силу разности форматов только число 0 не нуждается в преобразовании. Попытаемся вычислить квадратный корень из целого, заодно проверим его на <0 уже в регистре FPU:
1 li a7 5
2 ecall # Ввод целого
3 fcvt.s.w ft1 a0 # Преобразование в вещественное
4 fmv.s.x ft0 zero # 0 не надо преобразовывать!
5 flt.s t0 ft1 ft0 # t0 = ft1 < 0 ?
6 fmv.s fa0 ft0 # результат пока 0…
7 bnez t0 negat # …и останется 0, если число отрицательное
8 fsqrt.s fa0 ft1 # вычислим корень
9 negat: li a7 2 # выведем результат
10 ecall
11 li a7 10
12 ecall
Пример. Вычисление числа e как бесконечной суммы 1/(n!) при n→∞, т. е. Σ∞n=01/(n!)
1 .data
2 one: .double 1
3 ten: .double 10
4
5 .text
6 fld fs1 one t0 # Это 1
7 fld fs10 ten t0 # Это 10
8 fcvt.d.w fs2 zero # Это текущее n, пока 0
9 fmv.d fs0 fs1 # Это текущий n!, пока 1
10 fmv.d fs4 fs1 # Здесь будет e, пока 1
11 li a7 5 # количество знаков после запятой K
12 ecall
13
14 fmv.d fs3 fs1 # Здесь будет P = 10**(K + 1)
15 enext: fmul.d fs3 fs3 fs10 # P *= 10
16 addi a0 a0 -1
17 bgtz a0 enext
18 fdiv.d fs3 fs1 fs3 # А теперь тут ε = 1 / P
19
20 loop: fadd.d fs2 fs2 fs1 # n = n + 1
21 fmul.d fs0 fs0 fs2 # n! = (n - 1)! * n
22 fdiv.d ft0 fs1 fs0 # 1/n!
23 fadd.d fs4 fs4 ft0 # e += 1/n!
24 flt.d t0 ft0 fs3 # 1/n! < ε?
25 beqz t0 loop
26
27 fmv.d fa0 fs4 # Выведем e
28 li a7 3
29 ecall
30 li a7 10
31 ecall
- Для того, чтобы гарантировать попадание приближения в ε-окрестность результата, надо убедиться, что всякое очередное приближение будет отстоять не больше, чем на ε от текущего. В нашем случае это просто:
- все слагаемые ряда положительные
- все слагаемые строго убывающие
- значение слагаемого убывает более, чем экспоненциально,
Следовательно, достаточно, чтобы очередное слагаемое оказалось < ε