Математический сопроцессор

Вещественные числа и их запись

Несколько слов о представлении вещественных чисел.

«Настоящее» вещественное число не может быть представлено в памяти классического компьютера исключительно как структура данных. Например, иррациональные вещественные числа в десятичном представлении апериодичны и содержат бесконечно много цифр. Рациональные числа можно представить в виде отношения двух целых или в виде конечной дроби с периодом, но такие представления весьма сложны аппаратно, и не решают основной задачи — смоделировать произвольное вещественное число с заданной точностью.

Таким образом, вещественное в памяти компьютера число неизбежно должно моделироваться рациональным и некоторым представлением о точности этой модели. В учебниках вещественное число записывается в виде десятичной дроби — либо в формате с «фиксированной точкой», например 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:

Подобный формат кладется и в основу представления вещественного числа в памяти компьютера. Естественно вместо представления десятичного числа с плавающей точкой используется соответствующее двоичное число с плавающей точкой. Здесь уже начинаются сложности, особенно в вопросе точности представления числа. Допустим, мы храним мантиссу в виде 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

Понятно, что число может быть либо неотрицательным, либо отрицательным, поэтому один бит на знак нам придется использовать, и это самый старший бит. 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 используется хитрость, позволяющая сэкономить ещё один бит мантиссы:

Итак, в нормализованном и в денормализованном виде нет смысла записывать старший бит мантиссы в отведенные 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 приведены максимальные ошибки) на всякий случай запасена тут

В RARS есть хороший инструмент для изучения числа в формате IEEE 754: «Tools → Floating Point Representation». В частности, легко проверить, что когда в мантиссу записаны одни только нули, это соответствует близкому к нулю числу в денормализованной форме.

FrBrGeorge/MyDict/speech_balloon_question.png Чему соответствует мантисса, состоящая из всех 1? NaN:

Нелишне заметить, что 0 00000000 0 0000000000000000000000 — это просто 0 (он также называется «+0», единственное вещественное число, совпадающее с его округлением побитно), а 1 00000000 0 0000000000000000000000 — это «-0», и они не равны!

Сопроцессоры

FPU RISC-V

  1. Расширение F поддерживает числа с плавающей точкой одинарной и двойной точности в IEEE-формате.

    • Есть также расширения D (double), Q (quadruple) и ZfH (half)

    • В RARS есть "F" и "D"
  2. Архитектура предусматривает наличие тридцати двух регистров для чисел с плавающей точкой "f0–f31". Это не те же самые регистры, которые мы использовали до сих пор.

    • Если в архитектуре поддерживается "D" или "Q", для представления в этих регистрах чисел меньшей разрядности используется т. н. «NaN boxing» — заполнение старших частей регистра специальной константой NaN — «Not A Number» (см. ниже)

  3. В соответствие с соглашением о вызове подпрограмм, у этих регистров разное назначение:

    Общее название

    Мнемоника ABI

    Назначение

    Переживают ли вызов подпрограммы?

    f0 - f7

    ft0 - ft7

    Временные

    Нет

    f8 - f9

    fs0 - fs1

    Сохраняемые при вызове

    Да

    f10 - f17

    fa0 - fa7

    Параметры и возвращаемые значения

    Нет

    f18 - f27

    fs2 - fs11

    Сохраняемые при вызове

    Да

    f28 - f31

    ft8 - ft11

    Временные

    Нет

TODO Далее требуется вставить больше мелких примеров из лекции и расписать конспективные утверждения

Инструкции FPU-сопроцессора

Условные операторы и fcsr

Помимо непосредственных сравнений, FPU отображает своё состояние в управляющем регистре fcsr.

Блок управляющих регистров (введение)

Подробнее про блок CSR

Регистр fcsr (0x003):

Условные операторы

В любом случае условный переход при сравнение вещественных неатомарен, потому что сравнение делает математический сопроцессор, а анализ результата и вычисление перехода — арифметико-логическое устройство центрального процессора:

FrBrGeorge/MyDict/speech_balloon_question.png Почему нет инструкции «прочесть флаг № такой-то»?

Пример. Неравенство треугольника (функция проверки теперь совсем простая — это в чистом виде само неравенство):

   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

Замечание: стоит вспомнить, что равенство вещественных — довольно сомнительное сравнение.

FrBrGeorge/MyDict/speech_balloon_question.png Приведите пример двух равных вещественных чисел, ни один десятичный знак которых не совпадает! (нажмите «Комментарии» в шапке страницы, чтобы прочитать спойлер)

Для эффективности можно хранить значения регистров 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

Д/З

  1. Почитать про представление вещественных чисел в формате IEEE754, потыкать в «RARS Floating Point Representation Tool»

  2. EJudge: CubicRoot 'Кубический корень'

    Написать программу, которая вводит два вещественных числа одинарной точности, 1⩽A⩽1000000 и 0.00001⩽e⩽0.01. Программа вычисляет и выводит кубический корень из A с точностью e (округлять до фиксированного знака не обязательно). Подсказка: что-что, а уж возвести число в куб просто!

    Input:

    1000
    0.0001
    Output:

    10.000086
  3. EJudge: DoubleTrunc 'Неточная дробь'

    Ввести вещественное число двойной точности A и целое неотрицательное N: -230 ⩽ A * 10N ⩽ 230 . Округлить A до N-го знака после запятой путём отбрасывания лишних цифр и незначащих нулей (при выводе может возникнуть ***.0)

    • Алгоритм округления: домножить на 10N, перевести в целое с округлением RTZ, затем обратно, а затем поделить на 10N .

    • Округление оформить в виде подпрограммы, которой в a0 подаётся количество знаков после запятой, а в fa0 — число, после чего в fa0 возвращается округление.

    Input:

    -12.77859
    4
    Output:

    -12.7785
    • Функция округления ещё пригодится!
  4. EJudge: TaylorSin 'Вычисление синуса'

    Написать программу, которая вводит вещественное число -π/2 ⩽ x ⩽ π/2, затем — неотрицательное целое n ⩽ 9 и вычисляет sin(x) с точностью n знаков после запятой, а остальные цифры отбрасывает (можно воспользоваться подпрограммой из ../Homework_DoubleTrunc). Синус удобно вычислять с помощью ряда Тейлора. Использовать двойную точность.

    Input:

    0.1
    9
    Output:

    0.099833416
  5. EJudge: FRFlags 'Ошибки вычислений'

    Ввести три числа одинарной точности A, B и C, и вычислить формулу A/(B*C), именно в таком порядке — сначала умножение, затем деление. Для каждого действия (умножения и деления) проверять состояние четырёх флагов fcsr — NX, UF, OF и DZ, в указанном порядке, и выводить соответствующую строку из двух букв, если флаг ненулевой. Не забыть очистить fcsr перед вторым действием! Затем вывести результат.

    Input:

    12
    5e-41
    0.000000000003
    Output:

    NX
    UF
    DZ
    Infinity

LecturesCMC/ArchitectureAssembler2024/04_FloatingPoint (последним исправлял пользователь FrBrGeorge 2024-11-09 23:10:44)