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

IEEE 754

Базовая лекция

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

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

Хотелось бы напомнить, что мы имеем дело не с реальными вещественными числами, а с их (рациональным) представлением, в этом смысле любое вещественное число может быть представлено в привычном нам виде ( в виде числа с "фиксированной точкой"), например 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.2345-E2 соответственно. Таким образом вещественное число состоит из двух частей: мантиссы (в примере 1.2345) и показателя степени — порядка (в примере 3 и -2 соответственно). Одно и то же число можно записать несколькими способами, в зависимости от того, чему равна мантисса. Например, все эти записи соответствуют одному и тому же числу 155.625:

• 155.625E0
• 1.55625E2
• 0.001555625E5
• 155625E-3

  Это представление кладется в основу представления вещественного числа в памяти компьютера. Естественно вместо представления десятичного числа с плавающей точкой используется соответствующее двоичное число с плавающей точкой. Здесь уже начинаются сложности, особенно в вопросе точности представления числа. Допустим, мы храним мантиссу в виде 16-битного числа. То есть точность представления — до 16-го двоичного знака. Но в десятичном представлении 16-битное число — это 0…65535. Стало быть, за пятый десятичный знак мы ручаться не можем, только за четвёртый.

Хранение вещественных чисел в памяти

Представим число 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 на картинке)

В IEEE754 порядок (который может быть и положительным, и отрицательным, и нулём) хранится в виде всегда неотрицательного числа, для чего к нему добавляют «смещение» размером в половину допустимого диапазона. Например, в стандарте на 32-разрядное целое порядку отведено 8 битов (область E), в которых можно хранить числа 0…255. Допустимый порядок числа, таким образом, колеблется от -127 до 127, а при записи в ячейку к нему добавляется 127. В примере порядок = +7 (+111 в двоичной), соответственно, смещённый порядок = 7+127=134 ( 10000110 ). А если бы порядок был -7 , то смещенный порядок = 127-7 =120. Чтобы узнать истинный порядок, надо вычесть 127 из содержимого области E.

Если вспомнить, что «вещественные» числа IEE754 — на самом деле целочисленные двоичные дроби, каждое число можно представить себе в виде мантиссы, сдвинутой влево на соответствующее порядку число разрядов. Сдвиг на 0 разрядов будет соответствовать наименьшему по модулю числу, а сдвиг на 254 — наибольшему. Таким образом, для сравнения и проведения арифметических операций это представление наиболее удобно.

В 32 битном представлении на мантиссу (область M) приходится 23 бита. Это обеспечивает довольно низкую точность: 2²³=8388608, все числа, по модулю большие, будут терять значимые знаки до запятой. В IEE754 используется хитрость, позволяющая сэкономить ещё один бит мантиссы:

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

• Если число очень малое (близкие нулю вещественные числа используются часто, например, для задания точности или шага), оно хранится в денормализованном виде: мантисса обязана быть в диапазоне от 0.5 до 1 (не включая 1). Это также возможно, причем старший бит такой мантиссы (всегда 0) также не хранится. Порядок денормализованного числа на 1 больше порядка нормализованного, следовательно, в денормализованном виде можно хранить в два раза меньшее число, чем в нормализованном. Итак, в нормализованном и в денормализованном виде нет смысла записывать старший бит мантиссы в отведенные 23 бита, и по этой причине в отведенные 23 бита записывают «остаток от мантиссы». В случае нашего примера остаток от мантиссы будет 00110111010000 0000 0000.

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

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

Дело в том, что 24 бита на мантиссу — это очень мало. Например TODO

Более подробно смотрите основную статью Стандарт IEEE 754

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

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

FPU, он же C1

Базовая лекция

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

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

• Стандарт IEEE 754 (большая, но весьма толковая статья)

  • Mars -> Tools -> Floating Point Representation

  • О недостатках формата (см. ту же статью)

  • Манифест сообщества ненавистников IEEE 754

• FPU MIPS
  1. MIPS поддерживает числа с плавающей точкой одинарной и двойной точности в IEEE-формате.
  2. Архитектура MIPS предусматривает наличие тридцати двух 32-битных регистров для чисел с плавающей точкой "$f0–$f31". Это не те же самые регистры, которые мы использовали до сих пор.
  3. Числа двойной точности (64-битные) размещаются в парах 32-битных регистров, так что для операций с такими числами доступны только 16 четных регистров "$f0, $f2, $f4, ..., $f30".
  4. В соответствии с соглашением, принятым в архитектуре MIPS, у этих регистров разное назначение:

     Название	Номер	Назначение
     $fv0–$fv1	0, 2	Значение, возвращаемое функцией
     $ft0–$ft3	4, 6, 8, 10	Временные переменные
     $fa0–$fa1	12, 14	Параметры функции
     $ft4–$ft5	16, 18	Временные переменные
     $fs0–$fs5	20, 22, 24, 26, 28, 30	Сохраняемые переменные

  Примечание: это соглашение используется нечасто, например в Mars мнемонические имена регистров вещественного сопроцессора отсутствуют

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

• Код операции у всех команд с плавающей точкой равен 17 (10001₂ ). Для определения типа команды в них должны присутствовать поля funct и cop (от англ. coprocessor). Ниже показан формат команд типа F, используемый для чисел с плавающей точкой.

  • op	cop	ft	fs	fd	funct
    6bits	5bits	5bits	5bits	5bits	6bits

    Поле cop равно 16 (10000₂ ) для команд одинарной точности и 17 (10001₂ ) для команд двойной точности. Аналогично командам типа R, у команд типа F есть два операнда-источника (fs и ft) и один операнд-назначение (fd). Формат команд типа F

• Команды одинарной и двойной точности различаются суффиксом в мнемонике (.s и .d соответственно).

  • Общий вид: CMD.T $f2 $f4 $f6

    • где CMD = add, sub, div, mul; T — s или d

  • Общий вид: CMD.T $f2 $f4

    • где cmd = neg, abs, mov, sqrt, movf, movt; T — s или d

    • movf и movt — условное присваивание (см. далее)
  • Например:
    • mov.d $f4 $f6
    • mul.s $f1 $f2 $f8
    • sqrt.d $f0 $f4
• Работа с памятью и регистрами общего назначения (во второй колонке даны варианты написания)

  Загрузить слово из памяти по адресу offs($rn) в регистр $fs сопроцессора1
  lwc1 $fs offs($rn)	l.s $fs offs($rn)	$fs ← offs($fn)
  Записать слово из регистра $fs сопроцессора1 в память по адресу offs($rn)
  swc1 $fs offs($rn)	s.s $fs offs($rn)	$fs ← offs($fn)
  Загрузить двойное слово из памяти по адресу offs($rn) в чётный регистр $fd сопроцессора1
  ldc1 $fd offs($rn)	l.d $fd offs($rn)	$fd ← offs($fn)
  Записать двойное слово из чётного регистра $fd сопроцессора1 в память по адресу offs($rn)
  sdc1 $fd offs($rn)	s.d $fd offs($rn)	$fd ← offs($fn)
  Перенести слово из регистра $rn в регистр сопроцессора1 $fw
  mtc1 $rn $fw		$fw ← $rn
  Перенести слово из регистра сопроцессора1 $fw в регистр $rn
  mfc1 $rn $fw		$rn ← $fw
  Перенести двойное слово из регистров $rn и $rn+1 в чётный регистр сопроцессора1 $fd
  mtc1.d $rn $fd		$fd ← $rn
  Перенести двойное слово из четного регистра сопроцессора1 $fd в регистры $rn и $rn+1
  mfc1.d $rn $fd		$rn ← $fd

  Никакого преобразования эти инструкции не производят. Они просто перекладывают слово из регистра FPU в регистр или память и обратно. Инструкции ldc1/sdc1 работают с двойным словом и чётными регистрами FPU

• Преобразование форматов возможно между тремя представлениями: s (одинарной точности), d (двойной точности, чётные регистры) и w (целое);

  cvt.s.d $fs $fd	$fs ← float($fd)	convert d to s
  cvt.s.w $fs $fw	$fs ← float($fw)	convert w to s
  cvt.d.s $fd $fs	$fd ← double($fs)	convert s to d
  cvt.d.w $fd $fw	$fd ← double($fw)	convert w to d
  floor.w.s $fw $fs	$fw ← int($fs)	floor s to w
  floor.w.d $fw $fd	$fw ← int($fd)	floor d to w
  trunc.w.s $fw $fs	$fw ← int($fs)	truncate s to w
  trunc.w.d $fw $fd	$fw ← int($fd)	truncate d to w
  round.w.s $fw $fs	$fw ← round($fs)	round s to w
  round.w.d $fw $fd	$fw ← round($fd)	round d to w

  При преобразовании в/из целого это целое лежит в регистре FPU и ни к чему не пригодно, кроме преобразования в s или d или выгрузки в целый регистр MIPS с помощью mfc1.

Вопрос: чем плохо было бы объединить инструкции преобразования плавающего в целое и обратно с инструкциями обмена между регистрами общего назначения и регистрами FPU?

Особенности работы с сопроцессором

Сопроцессор FPU в MIPS использует биты состояния, которые хранит в специальном регистре. Этот регистр доступен побитно основному процессору.

Так же, как, например, в УМ-2, операции условного перехода становятся неатомарными:

• Сначала выполняется инструкция сравнения регистров вида "c.CMD.F $fn, $f", где F — .s или .d, а CMD — одна из команд le, lt или eq. Эта инструкция устанавливает флаг 0 в регистре состояния FPU в 1, если сравнение истинно (в формате "c.CMD.F b, $fn, $f" можно задать другой флаг с помощью b, предыдущая форма — по сути часто используемая псевдоинструкция).

• Затем выполняется условный переход: инструкция "bc1t метка" (общий вид "bc1t бит метка") или "bc1f метка" (общий вид "bc1f бит метка"), что расшифровывается как «branch if coprocessor 1 true» (и «false» соответственно) проверяет соответствующий бит регистра состояния FPU (по умолчанию нулевой).

• Пример: сравнить $f1 с $f2, установить флаг состояния 0 в 1, если $f1<$f2; если флаг — истина, перейти на метку lesser

     1 c.lt.s    $f1 $f2
     2 bc1t      lesser
  

• Сравнения на gt и ge легко получить перестановкой операндов в le и lt соответственно, либо в последующей проверке использовать bc1f вместо bc1t

• Для условного вычисления выражения, можно воспользоваться условными пересылками:

  поместить значение $rs в $rt если флаг 0 вещественного сопроцессора — истина
  movt $rd $rs	$rd ← $rs
  поместить значение $rs в $rt если флаг n вещественного сопроцессора — истина
  movt $rd $rs n	$rd ← $rs
  поместить значение $rs в $rt если флаг 0 вещественного сопроцессора — ложь
  movf $rd $rs	$rd ← $rs
  поместить значение $rs в $rt если флаг n вещественного сопроцессора — ложь
  movf $rd $rs n	$rd ← $rs

  Условные пересылки могут иметь суффикс .s или .d, тогда пересылаются не регистры общего назначения, а регистры вещественного сопроцессора. Например, команда "movf.s $f0 $f1" означает пересылку вещественного регистра одинарной точности $f1 в регистр $f0 при условии, что флаг 0 сопроцессора ложен

• Замечание: стоит вспомнить, что равенство вещественных — довольно сомнительное сравнение Для эффективности можно хранить значения регистров FPU в регистрах общего назначения, но в силу разности форматов только число 0 не нуждается в преобразовании. Попытаемся вычислить квадратный корень из целого, заодно проверим его на >0 уже в регистре FPU

     1 .data
     2 src:    .word   100
     3 dst:    .float  0
     4 .text
     5         lw      $t0 src         # исходное число целое
     6         mtc1    $t0 $f2         # положим в FPU
     7         cvt.s.w $f2 $f2         # превратим в вещественное
     8         mtc1    $zero $f0       # ноль не надо преобразовывать
     9         c.lt.s  $f2 $f0         # если <0 …
    10         bc1t    nosqrt          # обойдём корень
    11         sqrt.s  $f2 $f2
    12 nosqrt: s.s     $f2 dst         # результат вещественный
  

  Замечание: не путать в текстах lt и 1t :(

Пример вычислений с плавающей точкой

Вычисление числа e как бесконечной суммы 1/(n!) при n→∞, т. е. Σ^∞_n=01/(n!)

   1 .data
   2 one:    .double 1
   3 ten:    .double 10
   4 .text
   5         l.d     $f2 one         # 1
   6         sub.d   $f4 $f4 $f4     # n
   7         mov.d   $f6 $f2         # n!
   8         mov.d   $f8 $f2         # здесь будет e
   9         l.d     $f10 ten        # здесь будет ε
  10         mov.d   $f0 $f2         # количество знаков после запятой K
  11         li      $v0 5
  12         syscall
  13 
  14 enext:  blez    $v0 edone       # посчитаем 10**(K+1)
  15         mul.d   $f0 $f0 $f10
  16         subi    $v0 $v0 1
  17         b       enext
  18 edone:  div.d   $f10 $f2 $f0    # посчитаем ε
  19 
  20 loop:   add.d   $f4 $f4 $f2     # n=n+1
  21         mul.d   $f6 $f6 $f4     # n!=(n-1)!*n
  22         div.d   $f0 $f2 $f6     # очередное слагаемое
  23         add.d   $f8 $f8 $f0
  24         c.lt.d  $f0 $f10        # очередное слагаемое < ε
  25         bc1f    loop
  26 
  27         mov.d   $f12 $f8        # выведем результат
  28         li      $v0 3
  29         syscall

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

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

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

Д/З

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

• TODO

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

  Ввести два натуральных числа — A, B и целое неотрицательное n. Вывести вещественное число двойной точности F, в дробной части которого ровно n знаков числа A/B, с учётом возможных нулей, которые выводить не обязательно. Округление двойного вещественного числа до n знаков оформить в виде функции (она ещё понадобится ). Подсказка: предполагается, что 10ⁿ*A/B помещается в одно машинное слово.

  Input:

  123
  456
  7

  Output:

  0.2697368

• EJudge: LeibPi 'Вычисление Пи'

  Воспользоваться рядом Лейбница для вычисления числа Пи с точностью до определённого знака. Ввести N — количество точных знаков после запятой, вывести вычисленное с помощью ряда число Пи, пользуясь решением задачи ../Homework_FractionTruncate. Поскольку ряд Лейбница знакочередующийся, но сходящийся и монотонно убывающий по модулю, достаточно, чтобы сумма модулей двух соседних элементов ряда была меньше 10**(-N) (на самом деле точность должна быть выше, потому что ряд вычисляет Пи/4, а не Пи; я ставил 10**(-N)/8) . Ряд сходится медленно, дальше пятого знака добраться трудно.

  Input:

  4

  Output:

  3.1415

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

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

  Input:

  1000
  0.0001

  Output:

  9.99995