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

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

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

Хотелось бы напомнить, что мы имеем дело не с реальными вещественными числами, а с их (рациональным) представлением, в этом смысле любое вещественное число может быть представлено в привычном нам виде ( в виде числа с "фиксированной точкой"), например 1234.5 или 0.012345 или в научном (экспонециальном) виде (в виде числа с "плавающей точкой"). Например, 1234.5 = 1.2345*10³, 0.012345 = 1.2345*10^-2. Часто эти числа записываются в виде 1,2345·exp₁₀³ и 1,2345·exp₁₀^-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·2⁷ +0·2⁶+0·2⁵+1·2⁴+1·2³+0·2²+1·2¹+1·2⁰+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.625₁₀ = 10011011,101₂ - число в десятичной и в двоичной системе с плавающей точкой

Полученное число в экспоненциальном виде в двоичной системе будет иметь вид: 1,55625·exp₁₀⁺² = 1,0011011101·exp₂⁺¹¹¹

То есть его мантисса имеет вид 1.0011011101, а порядок — exp₂= +111

В действительности все несколько хитрее. Согласно стандарту IEEE 754, в памяти вещественные числа одинарной точности представляются в следующем виде (стандарт IEEE 754)

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 бита. Это обеспечивает довольно низкую точность: 2²³=8388608, все числа, по модулю большие, будут терять значимые знаки до запятой. В IEEE 754 используется хитрость, позволяющая сэкономить ещё один бит мантиссы:

• Если число достаточно большое, оно представляется в нормализованном виде: мантисса обязана быть в диапазоне от 1 до 2 (не включая 2). Двоичное представление позволяет для любого допустимого числа подобрать нормализованную мантиссу и соответствующий порядок. Так вот, в нормализованном виде старший бит мантиссы (а это всегда 1) не хранится вообще!

• Близкие нулю вещественные числа используются часто, например, для задания точности или шага. Если число настолько мало, что в его нормализованном представлении порядок состоит одних нулей (-127), оно хранится в денормализованном виде: мантисса при этом обязана быть не больше 0.5. Старший бит такой мантиссы (всегда 0) так что он тоже не хранится. Порядок денормализованного числа на 1 больше порядка нормализованного , следовательно, в денормализованном виде можно хранить в два раза меньшее число, чем в нормализованном.

Итак, в нормализованном и в денормализованном виде нет смысла записывать старший бит мантиссы в отведенные 23 бита, и по этой причине в отведенные 23 бита записывают «остаток от мантиссы». В случае нашего примера остаток от мантиссы будет 00110111010000 0000 0000.

Второй пример — это близкое к нулю число, мантисса которого хранится в денормализованом виде.

Вещественные числа двойной точности занимают 64 бита. Поле S в них такое же однобитное, на поле E (порядок) отводится 11 битов, а остальные 52 — на поле M (мантисса). Вещественные числа четверной точности — соответственно, 128 всего, 15 — E и 112 — M.

Дело в том, что 24 бита на мантиссу — это очень мало. Для чисел из обыденного диапазона сохраняется примерно 7-8 десятичных знаков в худшем случае. Например десять миллиардов представлены с точностью плюс-минус полтысячи.

Более подробно смотрите основную статью Стандарт IEEE 754 (в таблице 1 приведены максимальные ошибки).

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

А чему соответствует мантисса, состоящая из всех 1?

• Стандарт IEEE 754 (большая, но весьма толковая статья) на всякий случай запасена тут

  • Подробно о недостатках формата см.ту же статью, §9

• Манифест сообщества ненавистников IEEE 754 (разумеется, я к этому отношение не имею, просто спас из вебархива любопытный текст по ссылке из той же статьи☺)

NaN:

• 0 11111111 X **********************

  • x=0 — «NaN по-тихому», для вычислений, возвращающих NaN
  • x=1 — «сигнальный NaN», для исключений

  • RARS: ± 11111111 0 0000000000000000000000 — ±∞

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

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

• Различные, почти не пересекающиеся задачи ⇒ сопроцессоры:
  • Сопроцессор - FPU
    • Другая математика
    • Свои регистры
    • Почти не смешивается с целочисленной арифметикой
    • «Тяжёлые» вычисления
  • Сопроцессор 0 - управления (см. следующую лекцию)

  • RARS -> Tools -> Floating Point Representation

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-сопроцессора

• Регистры 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
• Вычисление полусуммы:

     1 .data
     2 a:      .float  123.456
     3 b:      .float  654.321
     4 _2:     .float  2
     5 .text
     6         flw     ft0 a t0
     7         flw     ft1 b t0
     8         flw     ft2 _2 t0
     9         fadd.s  ft3 ft2 ft1
    10         fdiv.s  fa0 ft3 ft2
    11         li      a7 2
    12         ecall
  

• Что это за 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)² подсчитанный с раскрытием скобок:

       1 .data
       2 x:      .double 12.34
       3 .text
       4         li              t2 2
       5         fcvt.d.w        ft2 t2                  # ft2 = 2.0
       6         li              t1 1
       7         fcvt.d.w        ft1 t1                  # ft1 = 1.0
       8         fld             ft0 x t0                # ft0 = x
       9         fnmsub.d        ft3 ft2 ft0 ft1         # ft3 = -(2 * x) + 1
      10         fmadd.d         fa0 ft0 ft0 ft3         # fa0 = x * x + ft3
      11         li              a7 3                    # Вывод числа двойной точности
      12         ecall
    

    • В раскрытии псевдоинструкции 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

    • остальные сравнения — псевдоинструкции

• Инструкция классфикации

  • fclas.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.

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

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

• 4K (12 битов) пространство управляющих счётчиков-флагов-масок-прочего

• Атомарные R/W инструкции типа I

  • csrrw[i], причём часть i тут нестандартная: 12-битное поле IMM занято номером CSR, так что число можно задать совсем маленькое, 5 битов вместо номера регистра-источника

  • а также csrrs[i] / csrrc[i] включение/выключение битов

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

• 5 битов флагов
  • NX потеря точности
  • UF сверхмалое число
  • OF переполнение
  • DZ деление на 0
  • NV недопустимая операция
• 3 бита: тип округления
  0. RNE ближайшее, лучше чётное
  1. RTZ ближайшее к нулю
  2. RDN ближайшее к -∞
  3. RUP ближайшее к +∞
  4. RMM ближайшее, лучше с бо́льшим модулем
  5. …
  6. …
  7. DYN умолчания не менять (используется в инструкциях явного округления)

• псевдоинструкции f[sr]csr

• RARS:
  • зачем-то отдельные регистры frm и fflags, хотя это просто поля

  • отдельные Инструкции frflags и fsrm для атомарного чтения/записи битов только из этого поля

       1 .data
       2 numb:   .float  7.65
       3 
       4 .text   
       5         flw     ft0 numb t0
       6         fcvt.w.s a0 ft0i        # По умолчанию 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:

    1. Сравнение
    2. Условный переход по содержимому x-регистра (0 / не 0)

  • fcsr*/frflags:

    1. Чтение всех флагов в x-регистр

    2. Выделение конкретного флага с помощью инструкции andi, например

    3. Условный переход по содержимому 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

• Для того, чтобы гарантировать попадание приближения в ε-окрестность результата, надо убедиться, что всякое очередное приближение будет отстоять не больше, чем на ε от текущего. В нашем случае это просто:
  • все слагаемые ряда положительные
  • все слагаемые строго убывающие
  • значение слагаемого убывает более, чем экспоненциально,

• Следовательно, достаточно, чтобы очередное слагаемое оказалось < ε

Д/З

• Почитать про представление вещественных чисел в формате IEEE754

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

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

  Input:

  1000
  0.0001

  Output:

  10.000086

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

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

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

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

  Input:

  -12.77859
  4

  Output:

  -12.7785

  • Функция округления ещё пригодится!

• EJudge: TaylorSin 'Вычисление синуса'

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

  Input:

  0.1
  9

  Output:

  0.099833416

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

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

  • Внимание! В задании и тестах используется ошибочная аббревиатура UN (вместо канонической UF), исправлять в этот раз не буду

  Input:

  12
  5e-41
  0.000000000003

  Output:

  NX
  UN
  DZ
  Infinity