MIPS: алгоритм деления (деление значений в формате IEEE-754) дает неверный ответ для последних 4-5 бит (LSB)

Question

Для случая, когда делитель> дивиденд, вычисленная мантисса дает неверный результат для последних 4-5 младших разрядов.

Я пытаюсь разделить значения и мантиссу двух чисел с плавающей точкой IEEE-754. Я использовал этот алгоритм деления

Когда делитель> дивиденд, мантисса нормализуется, но все еще неверна для последних 4-5 бит.

DivisionAlgorithm: # Делит Мантиссу Дивиденда (в $ a0) на Мантиссу Дивизора (в $ a1) за 25 Итераций # # Возвращает частное значение в $ v0 #

    addi $sp, $sp, -12      #Decrement in $sp (Stack Pointer)
    sw $s0, 0($sp)      #Pushing $s0 into Stack
    sw $ra, 4($sp)      #Pushing $ra into Stack (Function Call Exists)
    sw $s1,8($sp)       #Pushing $s1 into stack 

    move $t0, $a0       #$t0 = Mantissa of Dividend/Remainder
    move $t1, $a1       #$t1 = Mantissa of Divisor

    add $s0, $0, $0     #$s0 = Initialization
    add $v1, $0, $0     #$v1 = 0 (Displacement of Decimal Point Initialized)
    addi $s1,$0,1       #$s1 = 1 (initialize loop variable to 1)
    add $t3,$0,33

loop:   
    beq $t1, $0, check      #If Divisor = 0, Branch to check
    sub $t0, $t0, $t1       #Dividend = Dividend - Divisor
    sll $s0, $s0, 1     #Quotient Register Shifted Left by 1-bit
    slt $t2, $t0, $0
    bne $t2, $0, else       #If Dividend < 0, Branch to else
    addi $s0, $s0, 1        #Setting Quotient LSb to 1
    j out
else:   add $t0, $t0, $t1       #Restoring Dividend Original Value

out:    srl $t1, $t1, 1     #Divisor Register Shifted Right by 1-bit
    j loop

check:  slt $t2, $a0, $a1       #If Dividend < Divisor, Call Function 'Normalization'
    beq $t2, $0, exit       #If Dividend > Divisor, Branch to exit
    move $a0, $s0       #$a0 = Quotient

    jal Normalization       #Function Call 
    j return

exit:   move $v0, $s0       #$v0 = Calculated Mantissa

return: lw $ra, 4($sp)      #Restoring $ra
    lw $s0, 0($sp)      #Restoring $s0
    lw $s1, 8($sp)      #restoring $s1
    addi $sp, $sp, 8        #Increment in $sp (Stack Pointer)
    jr $ra          #Return

Нормализация: # Нормализует Мантиссу (в $ a0) и считает десятичные разряды, сдвинутые на десятичную точку # #Returns: # i) $ v0 = нормализованная мантисса # # ii) $ v1 = Количество десятичных знаков #

    lui $t0, 0x40       #$t0 = 0x40 (1 at 23rd-bit)
    addi $t2, $0, 1     #$t2 = 1 (Initialization)

loop2:  and $t1, $a0, $t0       #Extracting 23rd-bit of Mantissa 
    bne $t1, $0, else2      #If 23rd-bit = 1; Branch to else2
    addi $t2, $t2, 1        #Increment in Count of Decimal Places Moved
    sll $a0, $a0, 1     #Mantissa Shifted Left (To Extract Next Bit)
    j loop2         

else2:  sll $a0, $a0, 1     #Setting 24th-bit = 1 (Implied)
    move $v0, $a0       #$v0 = Normalized Mantissa
    move $v1, $t2       #$v1 = Displacement of Decimal Point    
    jr $ra          #Return

Например, я ожидаю, что выход 2.75 / 6.355 будет 00111110110111011000111011001110, но фактический вывод будет 00111110110111011000111011010110.

Answer 1

Ваш алгоритм неверен.

Правильный алгоритм для восстанавливающего деления был бы

qu=0
rem=dividend
repeat N times
  rem = rem - divisor
  if rem < 0
    rem = rem + divisor
    qu = (qu<<1)
  else
    qu = (qu<<1) + 1
  end
  rem = rem << 1
end

пока ваш

qu=0
rem=dividend
repeat N times
  rem = rem - divisor
  if rem < 0
    rem = rem + divisor
    qu = (qu<<1)
  else
    qu = (qu<<1) + 1
  end
  divisor = divisor >> 1  // instead of left shift remainder
end

Поскольку алгоритм основан только на сравнении divisor и rem, он кажется эквивалентным сдвигу вправо divisor или сдвигу влево rem. Но это не так.
При смещении делителя вправо вы теряете младшие значащие биты делителя.
Это может повлиять на сравнение и, следовательно, частное.
Если вы напечатаете остаток, вы увидите, что он оказывает существенное влияние и может иметь двукратную разницу в его значении по сравнению с правильным результатом.

Кажется опасным умножать на 2 остаток, так как может быть переполнение. Но если мы посмотрим на алгоритм, мы увидим, что этого не произойдет.
Первоначально дивиденд и делитель являются мантиссой некоторого FP, и, следовательно, 1 & le; делитель, дивиденд <2 и то же самое относится к rem, который изначально является копией дивиденда. Обратите внимание, что это означает, что rem <2 * div <br> Теперь мы делаем первый шаг вычисления
⚫ если rem

⚫ если rem & ge; div, затем мы вычисляем rem & minus; div и умножаем его на 2.
Как изначально rem <2 * div, rem (= 2 * (rem & minus; div)) <2 * (2 * div & minus; div), а свойство rem <2 * div все еще имеет значение true. </p>

Таким образом, на каждом шаге у нас всегда есть свойство rem <2 * div, и, если 2 * div можно закодировать, это гарантирует, что rem никогда не будет переполнен. </p>

С точки зрения реализации, вы можете закодировать эти числа в 24 младших бита целочисленного регистра. Это в значительной степени достаточно, поскольку точность остается неизменной.
В вашей реализации вы выполняете цикл 32 раза. Если вы хотите записать результат в мантиссу IEEE, он бесполезен, и количество итераций можно уменьшить. Достаточно петли
24 раза (размер мантиссы)
+ 1 (в случае деления <делитель, первый бит будет равен 0, но второй бит гарантированно будет равен 1) </p>

Если вы хотите выполнить некоторое округление результата, вы должны сделать два дополнительных шага, чтобы получить биты округления и защиты. После этих двух шагов вычисления, если остаток & ne; 0, установите липкий бит на 1, в противном случае установите его на 0. Затем выполните округление с обычными правилами.

Answer 2

я. В описании самого алгоритма деления в той форме, с которой вы начали (и, насколько я вижу, код ассемблера) отсутствует главная часть. При подготовке к запуску цикла деления, когда вы сдвигаете делитель вправо на 1 бит на каждой итерации, вы должны были сначала сдвинуть его влево, насколько это возможно и максимально полезно.

Позвольте мне объяснить это на десятичных числах: представьте, что вы делите 10001 на 73 на 6-значном аппарате (так что это 010001 на 000073). Вам следует сдвинуть 73 влево до тех пор, пока он не перестанет соответствовать (таким образом, максимальное смещение равно 4, а смещенный делитель равен 730000), или его наиболее значимая цифра (MSD) не станет выше MSD дивидендов (поэтому мы можем остановиться на 73000). Затем цикл запускается с

• смещенный_дивизор = 73000, смещение = 3
• shift__disisor = 7300, Shift = 2
• shift__disisor = 730, Shift = 1
• shift__disisor = 73, Shift = 0

(С десятичной дробью для каждого сдвига должен выполняться вложенный цикл вычитания; с двоичным - нет необходимости.)

Если вы сдвинете на 73 вправо больше, вы потеряете значащие цифры в делителе, и поэтому конечный коэффициент будет больше точного.

Вы можете безоговорочно сместить левый делитель на максимальную ширину (730000, смещение = 4), но это приведет к потере циклов ЦП. Если процессор имеет CountLeadingZeros операцию, вы можете сэкономить время с ним.

II. Второй главный момент - вы делите мантиссы, полученные из чисел с плавающей точкой. Это означает, что в наиболее типичном случае двух нормализованных чисел MSB деления и делителя будет одинаковым, и вы получите только 2 варианта отношения - 0 и 1. (Для двоичного кода IEEE32 это означает, что установлен бит 23 и дивиденд, и делитель находятся в диапазоне 8388608..16777215.)

Чтобы получить больше реальных цифр, вы должны выполнить длинное деление. (См. Ниже для более экономичного подхода; теперь, просто сравните с разделением учебника.) Для двоичного IEEE32, это означает, что перед окончательным округлением необходимо получить как минимум 27 чисел (24 из полученной мантиссы + охранник + округление + наклейка). Сдвиньте делимое влево на 27 бит и продолжите деление 51-битного дивиденда на 24-битный делитель, как описано выше, на 28 циклов.

Версия, описанная @AlainMerigot, по сути такая же, как я описал - вы все равно вычисляете (dividend<<quotient_width)/divisor, но менее очевидным способом. Если сдвинуть влево остаток на 1, это, вероятно, совпадает с делителем правого сдвига на 1, если данные не потеряны. Итак, первая итерация (число 0) сравнивает равные позиции и выдает MSB результата, но в данный момент она находится на бите 0; затем вы «активируете» еще один бит остатка на каждой итерации. Это позволяет избежать операций с двумя словами. В этом алгоритме, если вам нужен последний остаток, он должен был быть восстановлен путем сдвига вправо на счетчик итераций; но самому плавающему делению это не нужно.

Но: для его правильной работы вы должны правильно выровнять дивиденд и делитель с самого начала. В десятичном примере это означает, что вы должны начинать не с divident = 10001, divisor = 73, а с дивиденда = 10001, divisor = 73000. Если ваши исходные мантиссы были нормализованы до получения целочисленных значений, это автоматически выполняется. В противном случае требуются дополнительные усилия, и они практически такие же, как я описал для сдвига делителей при разделении учебников. (Похоже, именно поэтому процессоры Intel десятилетиями страдали от неограниченного времени ненормальной обработки;))