Реализовать вещественную арифметику (FPU) на базе целых C++

[ZYV]

Уважаемые коллеги,

Так случилось, что на экономическом (!!!) возникла задача написать программу на C++, реализующую класс вещественных чисел на основе целых. Образования по CS у меня, к сожалению, нет, поэтому обращаюсь к помощи зала

Конкретно вопрос заключается в алгоритме перевода дробной части десятичного вещественного числа в двоичный вид и наоборот, двоичного - в десятичный, причем только целочисленными вычислениями. Стандарт IEEE754 я прочитал и как они хранятся понял, всю оболочку сделал. Мне бы что-то, что из строки дробной части делало мне строку "ноликов и единичек", и, соответственно, наоборот. В битмаски потом рассую. Никаких ограничений на память и скорость работы нет.

Перевод "туда", хоть и завязанный на число позиций и с большими потерями точности я догадался как сделать "по-рабоче-крестьянски". Обратно - никак не могу понять. Умножать всё на $10^n$ - не хватит даже unsigned long long для эмуляции float. Помогите, пожалуйста, либо алгоритмом, либо исходником на C/C++/Pascal/PHP... К сожалению, время не бесконечно, учитывая ещё и собственную загруженность.

Спасибо!

[ZYV]

В продолжение темы про результатам работы за ночь: дана отсрочка до выходных. За ночь успел порыться в исходниках GCC на тему как же там эмулируется soft-FPU для архитектур, в которых не предусмотрено операций для работы с вещественными числами. Оказалось, меняется сама постановка задачи, о чем как-то не додумался. Ключевой момент тут в том, что мы должны обойтись без вещественных чисел, которые мы как раз реализуем, но любые операции над полем вещественных чисел, которые мы уже реализовали использовать имеем полное право.

Таким образом, soft-FPU, как я понял, там (в GCC / libc) реализуется следующим способом:

• Сначала реализуют элементарные операции с вещественными числами, а именно + и /.
• Затем реализуют операции сравнения <, > и = (>= и <= через <, >, =).
• Через (1) + (2) реализуют - и *. Примеры: 1 - 2 = 1 + (-2); 1 * 2 = 1 / (1 / 2 ).
• Дальше реализуют загрузку целых чисел в контейнер - это тривиально и может быть сделано только целочисленными операциями. Причем, чтобы сильно не думать, т.е. не делать самим перевод из десятичной в двоичную систему, берут готовую реализацию загрузки десятичных целых чисел и протыкивают 4-байтный контейнер int по битам, а затем копируют в контейнер IEEE754 с соответствующей экспонентой.
• Теперь можно получить любое наперед заданное вещественное число простым делением двух целых чисел и результат окажется автоматически в контейнере. Пример: 0.1 = 1 / 10.

Чтобы загрузить произвольное вещественное число из строки в контейнер, мы нормализуем его (приводим к виду +/-0.XXXXE+/-Y), а дальше кладем каждую цифру в контейнер как целое число, домножаем его на $10^n$ ($10 \times 10 \times 10 \dots \times 10$ или $1/10/10/10 \dots /10$) и складываем их все в результат.

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

Как реализовать сложение и вычитание однознаковых двоичных дробных чисел я себе более ли менее представляю - привести к одной экспоненте и сделать в столбик. С делением хуже - мне сейчас вообще не приходит в голову, как это реализовать. Помню, в школе, когда-то уголком делили дробные на целые...

Поэтому прошу всех участников, если у них есть такая возможность, поделиться простыми исходниками на любом языке или по крайней мере формализованными алгоритмами для этих операций. Скорость их выполнения и требуемые ресурсы совершенно не важны, поэтому можно делать самым "глупым" и простым способом.

Задача, по-моему, интересная, исходный код я без проблем в итоге могу опубликовать для потомков, я бы с удовольствием посвятил ей гораздо больше времени самостоятельно, если бы оно вообще у меня было.

[nestoklon]

Хм. А почему вы пытаетесь сделать класс чисел с плавающей запятой?
Вещественных чисел нет как таковых. Всё равно оные в компе представляются как 1) рациональные 2) с фиксированной запятой 3) с плавающей запятой.
Третьи делаются сложнее всего. Сделайте рациональные или с фиксированной точкой.
Рациональныое число -- это дробь, то есть пара. Тонкости две -- 1) надо научится сокращать сократимые дроби. 2) для того, чтобы их складывать, надо научится сначала приводить их к общему знаменателю. После этого все арифметические операции делаются тривиально.
С фиксированной запятой -- совсем элементарно. При вводе/выводе делим на 10,100,1000 или сколько душа пожелает. При умножении на то же число. Это совсем тривиально. Можно к тому же сделать этот финт через сдвиги -- будет к тому же чуть ли не быстрее арифметики с плавающей точкой. Во времена моего детства так игрушки писали.

[ZYV]

nestoklon
Хм..., т.к. таково задание. Разумеется, реализация должна быть самой примитивной - ни о каких-там распространениях NaN, округлениях и т.п. речи не идет, т.к. важен сам принцип.

[worm2]

Концептуально деление "столбиком" в двоичном коде просто. Правда, детали мешаются. Если ограничитья только случаем беззнакового деления нормализованных [примечание: под нормализацией я понимаю приведение чисел к такому виду, что их старшие биты равны 1] ненулевых беззнаковых чисел, с выдачей только нормализованной мантиссы результата, то алгоритм может выглядеть примерно так:

unsigned Мантисса_Частного(unsigned Мантисса_Делимого, unsigned Мантисса_Делителя)
начало
  unsigned A = Мантисса_Делимого
  unsigned B = Мантисса_Делителя
  unsigned Result = 0
  for i = 1 to K // K - требуемая разрядность результата
    Result = Result << 1 // поразрядный сдвиг влево
    если A >= B
      Result = Result + 1
      A = A - B;
    конец если
    B = B >> 1 // поразрядный сдвиг вправо
  конец for
  // Нормализация результата
  если <Старший бит Result> = 0
    Result = Result << 1 // предыдущий по старшинству обязан быть=1
    // Показатель в этом случае нужно уменьшить на 1!
  конец если
  return Result
конец

В этом псевдокоде я предполагаю, что K --- требуемая разрядность мантиссы --- совпадает с разрядностью используемых целых чисел (Мантисса_Делимого, Мантисса_Делителя, A, B, Result). Если разрядность целых чисел меньше, тогда караул Если больше, то под <Старший бит Result> нужно понимать (K-1)-й бит, и ещё надо побитово сдвинуть Result влево на недостающее число бит.

К этому ещё нужно присобачить правильную обработку показателей (которая в вышенаписанном коде только обозначена), знаков и нулевых значений.

Добавлено спустя 13 минут 54 секунды:

Я тут упустил ещё один нюанс: при поразрядном сдвиге B вправо теряются значимые цифры делителя, что вносит искажения в результат. Вероятно, искажения будут не слишком сильными, но всё-же, по-хорошему, чтобы от них избавиться, разрядности A и B надо бы увеличить в 2 раза (со сдвигом значащих цифр на соотв. число бит влево) по сравнению с разрядностью делимого и делителя.

Добавлено спустя 3 минуты 8 секунд:

Имхо, для любого факультета, кроме Computer Science (ну, ещё математиков-прикладников), это задание неоправданно сложное. Поневоле в голову лезут нехорошие мысли

[ZYV]

worm2
Спасибо большое за подробный комментарий по самому важному вопросу. Нули и знаки мне теперь реализовать не составит труда.

В IEEE754 single мантисса хранится как раз в нормализованном виде, за исключением подразумеваемого бита, т.е. обеспечивает (с его учетом) точность 24 бита, поэтому я всё уже нормализовал unsigned int обеспечивает мне 32 бита для использования по моему усмотрению, поэтому, я так понимаю, проблемы тут нет. Для того, чтобы убрать искажения можно использовать unsigned long int - это даст мне 64 бита, что больше 48 бит. Правильно ли я понимаю, что при подготовке A и B мне надо будет сдвинуть их до упора влево?

Ещё мне совершенно непонятно как в вашем методе учитываются экспоненты. В нормализованном виде у чисел же не обязательно будет двоичная экспонента, равная нулю.

Кстати, подумал тут немного над сложением и вычитанием. На самом деле, видимо, не всё так просто. Надо реализовать отдельно сложение положительных чисел и вычитание положительных чисел, когда одно число больше другого, а потом уже на основе этого добавить обертки, которые разбирались бы с нулями и с тем, что второе больше первого и т.п.

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

Добавлено спустя 48 секунд:

P.S. Мысли мне тоже по этой ситуации приходят разные, но, видимо, пока их лучше держать при себе.

[nestoklon]

ZYV в сообщении #204819 писал(а):

Хм..., т.к. таково задание.

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

[ZYV]

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

[PAV]

Я не очень следил за темой, но почему бы просто не представить вещественное число в виде пары (целое+степень 10). Например, число $10.34$ представить в виде $1034$ и $2$. Какие здесь возникают сложности?

[worm2]

ZYV в сообщении #204837 писал(а):

Правильно ли я понимаю, что при подготовке A и B мне надо будет сдвинуть их до упора влево?

Считая, что исходные данные нормализованы, достаточно сдвига на 24 бита (для точности float).

ZYV писал(а):

Ещё мне совершенно непонятно как в вашем методе учитываются экспоненты. В нормализованном виде у чисел же не обязательно будет двоичная экспонента, равная нулю.

А в моём методе и не учитываются экспоненты (в моём предыдущем сообщении они называются "показатели"), я писал только для мантисс.
При делении из экспоненты делимого вычитается экспонента делителя, а затем результат уменьшается на 1, если нормализованная мантисса делимого меньше нормализованной мантиссы делителя (в предыдущем сообщении я в коде указал это место комментарием "// Показатель в этом случае нужно уменьшить на 1!").

ZYV писал(а):

Но теперь уже и со сложением ситуация перестала быть понятной - чтобы сложить в столбик не думая, надо сдвинуть нормализованные мантиссы вправо до предела и оставить 1 бит на переполнение (или не один?!). Но где взять столько памяти? Надо опять как-то учитывать экспоненту... Sad

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

• представлять каждое число в виде числа с фикс. точкой очень большой разрядности (256+24=280 бит, множитель $2^{-128}$)
• отдельно рассматривать случаи сложения чисел с одинаковым знаком, с разными знаками (нужно запомнить знак, из большего по модулю вычитать меньшее, затем восстанавливать знак) и с нулём
• для каждого слагаемого сдвигать мантиссу вправо на (128-Exp) битов (т.о. число x будет представлено в виде числа $A = x\cdot 2^{128}$, у которого экспонента = 128)
• осуществить длинное сложение или вычитание
• не забыть отдельно проверить, не равен ли результат точно нулю
• сдвигать результат влево, пока старший бит не станет равным 1
• Экспонента результата будет равна 128-<число сдвигов>

Конечно, можно оптимизировать эту концептуальную модель. Легко видеть, что если экспоненты слагаемых отличаются более, чем на 24, меньшим по модулю можно просто пренебречь, а в остальных случаях хватит 64-битных операций.
Но всё же не исключаю, что реализовать эту модель, как есть (с длинным сложением/вычитанием) будет проще.

[ZYV]

Уважаемый worm2!

Спасибо большое за квалифицированную помощь. Вчера закончил длинное сложение и вычитание (сделал через обратный код) по вашему алгоритму, ноль учел при денормализации, знаки - два варианта - когда одинаковые и когда разные, у разных - 4 подварианта. Всё вроде бы работает.

Могу отметить, что зря я беспокоился по поводу памяти, ничто не помешало взять и сделать bitset размера REAL_SFD + REAL_BIAS + REAL_BIAS - 1 и с ним оперировать, а bitset, между прочим, он на то и bitset, что для хранения каждого бита действительно использует 1 bit. Только насчет 280 бит у вас немного завышенная оценка, т.к. на самом деле экспоненты могут меняться только от -126 до 127, ещё два значения - все нули и все единицы зарезервированы для специальных значений - нули для нуля, а единицы, как я понял, для NaN и переполнений.

Сейчас буду пытаться осознать и реализовать ваш алгоритм деления. Надо бы закончить сегодня-завтра, а то в воскресение я улетаю. Скажите, пожалуйста, а тут нельзя провернуть такой же фокус как со сложением и вычитанием, чтобы не учитывать экспоненты в алгоритме? Памяти и времени не жалко. Если можно, то как?

Правильно ли я понимаю, что при реализации умножения через деление никаких подводных камней нет, кроме деления на ноль?

Спасибо!

[worm2]

ZYV в сообщении #205556 писал(а):

а тут нельзя провернуть такой же фокус как со сложением и вычитанием, чтобы не учитывать экспоненты в алгоритме? Памяти и времени не жалко. Если можно, то как?

А зачем не учитывать экспоненты в алгоритме? Я же написал, их при делении, и тем более умножении, учитывать очень просто, гораздо проще, чем при сложении. Экспоненты просто вычитаются при делении и складываются при умножении. Никакого сдвига мантисс здесь не нужно, за исключением однократной корректировки при "переполнении" старшего разряда (умножение) или при мантиссе делимого < мантиссы делителя (деление).

ZYV в сообщении #205556 писал(а):

Правильно ли я понимаю, что при реализации умножения через деление никаких подводных камней нет, кроме деления на ноль?

Да, вроде бы, других нет. Ну, ещё может получиться так, что какое-то число представимо в Вашем типе, а обратное к нему --- нет. Например, в IEEE, насколько я знаю, есть денормализованное представление очень маленьких по модулю чисел, а попытка вычисления обратных к ним вызывает переполнение. Но такие маргинальные случаи можно со спокойной совестью проигнорировать (независимо от того, IEEE там у Вас или нет). Чай, не Mission Critical система

[ZYV]

Уважаемый worm2, спасибо большое за помощь!

Всё действительно получилось, умножение через деление работает, с экспонентами разобрался, наконец Остался только один (вернее два), но очень существенных вопроса, а именно, как реализовать загрузку чисел из строк в этот класс и выгрузку из класса в строку?

Сначала я подошел к решению просто - вводимое число разделяется на целую и дробные части, затем целая загружается через long int, а дробная преобразовывается разложением по обратным степеням двоек. Похоже и с выгрузкой в виде строки.

Однако оказалось, что вещественные числа за счет экспоненты могут хранить в себе очень большие числа (больше unsigned long int), пусть и с точностью всего 24 двоичных знака. Поэтому в моих алгоритмах возникает переполнение long'а, и, хотя числа-то там где-то внутри в битовой маске есть, вывести я их не могу, равно как и загрузить.

Меня вполне устроило бы отображение их как X.YYYYE+/-ZZ, ну и приводить любой ввод пользователя к такому виду я, наверное, тоже смогу. Не могли бы вы меня снабдить алгоритмами для загрузки и выгрузки, чтобы я мог получать из/в такой нотации бинарную экспоненту и сигнификанд?

Может быть есть какой-то другой способ всё это реализовать? Теперь у меня всё есть - и операции сравнения, и сложение, и умножение, и деление.

[worm2]

Поскольку производительность не критична, можно пойти по такому пути:

I. Преобразование из строки (из десятичной системы).
1) Вычисление $10^{-38}$, $10^{-37}$, $\dots$, $10^{38}$ (для float) и сохранение этих значений в таблице T.
2) Ввод только десятичной мантиссы M (например, 3141E+30 вводится как 3.141) с домножением j-го знака после запятой на $10^{-j}$ и с запоминанием десятичной экспоненты D (в данном случае будет 33).
3) Умножение полученной мантиссы M на $10^D$ из таблицы.

II. Преобразование в строку (в десятичную систему).
1) Нахождение (по таблице T) такого j, что $10^j\leqslant |x| < 10^{j+1}$.
2) Делением на $10^j$ получаем первую цифру десятичной мантиссы $d_0$.
3) Вычитаем из |x| число $d_0\times 10^j$, уменьшаем j на 1 и продолжаем процесс.
4) Заканчиваем процесс, когда вывели достаточное число десятичных знаков (7 для float).

Немного коряво, но в большинстве случаев работать будет. Разумеется, знаки нужно обрабатывать отдельно.

[ZYV]

worm2, спасибо огромное! Аналогичный алгоритм, как я выяснил, используется и в dietlibc, но не мог разобраться в нем сам до конца. Теперь, с вашими комментариями, всё получилось. Больше наводить красоту не буду, итак потратил почти 3 полных дня, т.к. пришлось для начала выучить C++ (никогда на нем не писал), а потом ещё, к тому же, сами видите, какие у меня проблемы с алгоритмами.

Если интересно, могу выложить финальную версию класса как приеду через две недели под лицензией MIT. Сейчас его не хотелось бы светить до сдачи, на случай, если другие студенты умеют пользоваться Google и, в то же время, не хотят платить $300

Можно задать вам нескромный вопрос? А откуда вы всё знаете?