Синяя Гусеница: А ну нас-то не собьешь. Мы себе сидим, знаем: ждут нас превращения. А что? А ничего! Сидим, курим, ждем…
Алиса- кукла: Чего?
Синяя Гусеница: Чего, чего! Превращений. Дом — в дым, дым- в даму, а дама — в маму. Вот так-то. Не мешайте, не заскакивайте вперед, а то сами еще превратитесь преждевременно в бабочку какую-нибудь.
Просматривая код на одном из форумов, посвященных Ардуино, обнаружил забавный способ работы с числом с плавающей точкой (ПТ). Второе общепринятое название для чисел в таком формате — с плавающей запятой, но возникающее при этом сокращение (ПЗ) лично у меня вызывает совсем другие ассоциации, так что будем пользоваться именно этим вариантом. Первое впечатление (от увиденного кода) — что за фигня тут написана (надо сказать, что и второе такое же, хотя есть нюансы, но об этом позже), но возникает вопрос — а как надо на самом деле — ответ на который и дается в дальнейшем тексте.
Часть первая — постановка вопроса
Формулируем задачу — нам нужно вывести на консоль (превратить в символьное представление) число с плавающей точкой, не используя опции печати, именно для этой цели предназначенные. Почему мы хотим это сделать самостоятельно —
- использование формата %f влечет за собой подключение библиотеки для работы с плавающей точкой и расширенного варианта функции prntf (вернее, делает невозможным использование ее усеченного варианта), что приводит к существенному увеличению размера исполняемого модуля,
- стандартное решение требует значительных временных затрат (оно всегда работает с числом двойной точности), что может быть неприемлемым в данной конкретной ситуации,
- ну и (last, but not least), это просто интересно.
Для начала рассмотрим тот вариант, который был предложен в упомянутом материале, что то вроде:
for (float Power10=10000.0; Power10>0.1; Power10/=10.0; )
{char c=(int)(Fdata/Power10); Fdata -=Power10*c; };
и согласимся, что задачу он вполне себе решает. Более того, это не самый плохой вариант, поскольку его быстродействие может быть вполне приемлемым. Вот этот момент давайте рассмотрим подробнее — мы видим деление чисел ПТ, но, если глубже вникнуть в суть вопроса, то выяснится, что оно по скорости работы почти не уступает делению целых чисел соответствующей разрядности. Вообще то, прежде чем оценивать быстродействие алгоритма, следует оценить быстродействие различных элементарных операций, чем мы и займемся.
Часть вторая — оценка быстродействия элементарных операций
Первая интересная операция — сложение (вычитание, в смысле затрат времени они эквивалентны) целах чисел и мы можем считать, что она занимает единицу времени (такт) со следующей оговоркой — это справедливо только для «нативных» данных. К примеру, для МК серии AVR это 8-битовое слово, для MSP430 это 16-битное слово (ну и, конечно, меньшее по размеру), для Cortex-M это 32-битное слово и так далее. Тогда операция сложения чисел длиной в Н раз больше нативной может быть оценена как Н тактов. Есть исключения, например, AddW в AVR контроллерах, но оно не отменяет правило.
Следующая операция — умножение целых чисел (но не деление, оно отличается в смысле быстродействия) и для него уже не все столь однозначно. Прежде всего, умножение может быть реализовано аппаратно и, например, в AVR MEGA она требует 2 тактов, а в улучшенном 51 целых 6 (для перемножения нативных чисел).
Но рассмотрим случай, когда аппаратной реализации нет и нам придется реализовывать умножение в виде подпрограммы. Поскольку при умножении Н разрядных чисел получается 2Н разрядное произведение, то оценку классического варианта со сдвигами можно найти следующим образом: нам потребуется Н сдвигов множителя с 1 тактом на сдвиг, Н сдвигов второго множителя длиной 2Н с 2 тактами на сдвиг, далее Н принятий решения и, в среднем, Н/2 сложений чисел длиной 2Н, в завершение организация цикла 2 такта. Итого Н+2Н+Н+2Н/2+2Н=7Н тактов, причем собственно выполнение арифметических операций из них занимает всего Н тактов (ничего себе КПД, хотя паровоз нам обойти удалось).
То есть для умножения двух 8р чисел на 8р МК потребуется 56 тактов, а для умножения 16р чисел уже 112 (чуть меньше, но пренебрегаем точным значением) тактов, что несколько больше, чем хотелось. К счастью, направление сдвигов можно модифицировать и существует, причем единственный, способ осуществления умножения, который потребует только Н сдвигов числа длиной 2Н разрядов и Н/2 сложений нативных чисел, что улучшает время работы алгоритма умножения до 0+2+1+1/2+2=5.5Н — конечно, с аппаратной реализацией не сравнить, но хоть какой-то выигрыш без потери функциональности. Есть улучшения данного алгоритма, например, анализ 2 битов за цикл, но кардинально оценку ситуации они не меняют — время умножения на порядки превосходит время сложения.
А вот с делением ситуация похуже — даже аппаратно реализованное деление проигрывает умножению почти в два раза, к тому же существуют МК с аппаратным умножением, но без аппаратного деления. При определенных условиях можно заменить деление умножением на обратную величину, но эти условия специфичны и дают аналогичный результат — требуется две итерации умножения с последующей суммой, так что проигрыш в 2 раза. Если же мы реализуем деление в виде подпрограммы, то потребуется Н сдвигов делителя длиной 2Н, Н вычитаний из делимого длиной 2Н, Н сдвигов результата, 2Н организации цикла, но всему этому предшествует выравнивание, которое займет еще 5Н тактов, так что общая цифра составит 2+2+1+2+5=12Н, что примерно в 2 раза хуже умножения.
Ну а теперь рассмотрим операции ПТ, и тут ситуация несколько парадоксальная — операция умножения требует почти столько же времени, как и для целых чисел (соответствующей разрядности, как правило, 24 бита), поскольку мы должны перемножить мантиссы и всего лишь сложить порядки, нормализация не требуется. С делением тоже хорошо, делим мантиссы и вычитаем порядки, нормализация опять не нужна. Поэтому для этих двух операций проигрыш по сравнению с целыми числами не слишком значителен, хотя и имеет место быть.
А вот операция сложения и вычитания требует, прежде всего, выравнивания порядков (а это сдвиги и их может быть много, хотя есть нюансы), потом собственно операции, и (при вычитании) нормализации (при сложении тоже, но там требуется не более 1 сдвига), что расточительно по времени, поэтому операции данного класса для ПТ существенно медленнее, чем для целых, особенно в относительном выражении.
Вернемся к нашим баранам и согласимся, что, исходя из ранее данных оценок, предложенный способ может оказаться не слишком длительным, тем более, что он сразу же дает результат, но у него есть существенное ограничение — он применим к весьма ограниченному диапазону входных значений ПТ. Поэтому будет искать универсальное (более-менее) решение.
Сразу оговоримся, что наше решение не должно использовать операции с плавающей точкой вообще (от слова совсем), чтобы подчеркнуть достоинства нашего варианта. А на недоуменный вопрос — откуда тогда вообще появится число такого типа, если операции недоступны, ответим — вполне может появиться, например, при считывании информации с датчика освещенности (как и было в исходном примере), который выдает данные именно в формате ПТ.
Как именно устроено число ПТ, Вы легко найдете на многочисленных сайтах, на Хабре недавно была статья, с этим проблем возникнуть не должно. Тем не менее, вызывает интерес ряд вопросов к формату ПТ в стиле «если бы директором был я» — почему именно так, а не иначе. Я дам свои ответы на некоторые из них, если кто знает более правильные — прошу в комменты.
Первый вопрос — почему мантисса хранится в прямом коде, а не в дополнительном? Мой вариант ответа — потому что так проще работать с нормализованной мантиссой со скрытым (дополнительным) битом.
Второй вопрос — почему порядок хранится со смещением, а не как-либо иначе? Мой ответ — потому что в таком случае легко провести сравнение модулей двух ПТ, как целых чисел, при других способах это сложнее.
Третий вопрос — почему отрицательный знак кодируется единицей, а не нулем, ведь тогда можно было бы просто сравнивать два ПТ, как целые? Мой ответ — не знаю, просто «здесь так принято».
Часть третья — необходимые пояснения
В предыдущем параграфе я мог дать непонятные термины, поэтому немного о представлении чисел. Разумеется, они разные, иначе бы не было необходимости их обсуждать. Сразу же отметим, что в памяти МК (насчет ЭВМ верно то же самое, хотя насчет самых современных архитектур я не столь категоричен — они настолько сложны, что всего можно ожидать) нет никаких чисел, там лежат только элементарные единицы хранения — биты, сгрупированные в байты и далее в слова. Когда мы говорим о представлении числа, то имеется в виду, что набор битов конкретной длины мы так или иначе интерпретируем, то есть задаем закон, по которому можем найти некоторое число, соответствующее данному набору бит, и больше ничего.
Таких законов можно придумать бесчисленное множество, но некоторые из них будут обладать рядом полезных свойств в плане проведения разнообразных операций, поэтому они будут чаще применяться в практике. Одним из таких свойств, которое подразумевается неявно, например, является детерменированность, а другим — независимость от окружения — свойства, на первый взгляд, очевидные, хотя есть нюансы. Другие же свойства типа взаимно однозначного соответствия уже являются предметом дискуссии и не всегда имеют место в конкретном представлении. Сама по себе тема представления чисел необычайно увлекательна, у Кнута (в томе втором) она раскрыта в полной мере, так что за глубинами туда, а мы пробежимся по поверхности.
В предположении, что набор битов имеет длину н (пронумеруем их подряд от 0 до н-1) и взвешен равномерно с шагом 2 и младший бит (с номером 0) имеет вес 1 (что, вообще то говоря, совершенно не обязательно, просто мы к таким вещам привыкли, и они нам кажутся очевидными) мы получаем двоичное представление числа, при котором формула приведения выглядит следующим образом: число, отображаемое набором битов (Ч2) = В(0)*2^0 + В(1)*2^1 + ... + В(н-1)*2^(н-1) или в каскадной форме Ч2(н) = В(0)+2*(В(1)+2*(...+2*(В(н-1))..))), здесь и далее В(к) обозначает бит с номером к. Обратим внимание, что подобное представление не накладывает никаких ограничений на расположение байтов числа в памяти, но будет логичнее младший байт располагать в младших адресах (вот как легко и непринужденно я разрешил «вечный спор славян между собой» относительно того, с какого конца удобнее разбивать яйцо).
При такой интерпретации набора битов длиной н(=8) мы получаем представление для чисел от 0 до (2^н)-1(=255) (здесь и далее в скобках будет конкретное значение для набора из 8 бит), которое обладает рядом замечательных и полезных свойств, почему оно и получило широкое распространение. К сожалению, оно обладает и рядом недостатков, одним из которых является то, что отрицательные числа мы представить в такой записи не можем в принципе.
Можно предложить разнообразные решения данной задачи (представление отрицательных чисел), среди которых есть и имеющие практическое значение, они перечислены ниже.
Представление со смещением описывается формулой Ч = Ч2(н) — смещение(С), где Ч2 — число, получаемой при двоичной записи с н битами, а С — некоторое заранее выбранное значение. Тогда мы представляем числа от 0-С до 2^(н)-1-С и, если мы выберем С=2^(н-1)-1(=127) (это совершенно не обязательно, но весьма удобно), то получим диапазон от 0-(2^(н-1)-1)(=-127) до 2^(н-1)(=128). Главное достоинство данного представления — монотонность (более того, возрастание) на всем интервале, есть и недостатки, среди которых выделим несимметричность (есть и другие, связанные со сложностью выполнения операций над числом в этом представлении), но разработчики стандарта IEEE 457 (это стандарт для ПТ) превратили этот недостаток в достоинство (использовав лишнее значение для кодировки ситуации nan), что еще раз подчеркивает верность классного высказывания: «Если ты выше противника, то это твое преимущество. Если противник выше тебя, то это тоже твое преимущество».
Обратим внимание, что поскольку общее количество возможных комбинаций любого количества бит четное (если у Вас нет запрещенных по религиозным мотивам комбинаций), то симметрия между положительными и отрицательными представимыми числами принципиально недостижима (вернее, достижима, но при определенных дополнительных условия, о чем далее).
Представление в виде прямого кода, когда один из битов (старший) представляет закодированный знак числа Ч = (-1)^В(н-1) * Ч2(н-1) имеет диапазон от 0-(2^(н-1)-1) (=-127) до 2^(н-1)-1 (=127). Интересно отметить, что я только что заявил о принципиальной невозможности симметрии, а тут она явно есть — максимальное представимое положительное число равно модулю минимального представимого отрицательного числа. Достигается такой результат путем наличия двух представлений для нуля (00...00 и 10...00), что обычно считают главным недостатком данного способа. Это действительно недостаток, но не такой страшный, как принято считать, поскольку есть более существенные, которые и ограничили его применение.
Представление в виде обратного кода, когда мы в прямом представлении инвертируем все биты значения для отрицательных чисел Ч=(1-В(н-1))*Ч2(н-1) + В(н-1)*(2^(н-1)-Ч2(н-1)) — это из определения, можно сделать гораздо более понятную формулу Ч=Ч2(н-1)-В(н-1)*(2^(н-1)-1), что позволяет представить числа от 0-2^(н-1)+1(=-127) до 2^(н-1)-1(=127). Видно, что это представление со смещением, но смещение меняется скачкообразно, что делает данное представление не монотонным. Опять имеем два нуля, что не очень страшно, гораздо хуже возникновение кругового переноса при сложении, что создает определенные проблемы в реализации АЛУ.
Устранить последний недостаток предыдущего представления оказывается необычайно просто, достаточно изменить смещение на единицу, тогда получаем Ч=Ч2(н-1)-В(н-1)*2^(н-1) и можем представить числа от 0-2^(н-1)(=-128) до 2^(н-1)-1(=127). Нетрудно видеть, что представление ассимеричное, зато ноль единственный. Значительно более интересно следующее свойство, «совершенно очевидно, что» кольцевого переноса для операции типа сложение не возникает, что и послужило причиной (вместе с другими приятными особенностями) всеобщего распространения именно этого способа кодирования отрицательных чисел.
Составим таблицу интересных значений для различных способов кодирования чисел, обозначив через Н значение 2^(н-1) (128)
| Биты | 00..00 | 01..11 | 10..00 | 11..11 |
|---|---|---|---|---|
| Ч(н) | 0 | Н-1(127) | Н(128) | 2*Н-1(255) |
| Ч(н-1) | 0 | Н-1(127) | 0 | Н-1(127) |
| Смещ. Н | -Н+1 (-127) | 0 | 1 | Н(128) |
| Прямой | 0 | Н-1(127) | 0 | -Н+1(-127) |
| Обратный | 0 | Н-1(127) | -Н+1(-127) | 0 |
| Дополн. | 0 | Н-1(127) | -Н(-128) | -1 |
Ну и в завершение темы приведем графики для перечисленных представлений, из которых сразу видны их достоинства и недостатки (конечно же, далеко не все, что заставляет вспомнить интересное изречение «Преимущество графического представления информации — в наглядности, и других преимуществ у него нет»).
Часть четвертая — собственно решение исходной задачи (лучше поздно, чем никогда).
Небольшое отступление
Для начала я хотел отпечатать ПТ в шестнадцатиричном формате (и в конечном счете я это сделал), но вполне ожидаемо/совершенно неожиданно (нужное подставить) наткнулся на следующий результат. Что, по Вашему, будет отпечатано в результате выполнения операторов:
printf("%f %x", 1.0,1.0); printf("%f %x",2.0,2.0);
printf("%x %d",1.0,1.0); printf("%x %d",2.0,2.0);, обратите также внимание на следующую конструкцию и ее результат:
printf("%x %x %f",1.0,1.0);Пояснений данному явлению давать не буду, «умному достаточно».
Тем не менее, как же нам правильно вывести на печать шестнадцатиричное представление ПТ? Первое решение очевидно — union, а вот второе предназначено для любителей однострочников printf("%x",*( (int *) (&f))); (прошу прощения, если кого обидели лишние скобки, но я никогда не мог, да и не собирался, запомнить приоритеты операций, особенно, если учесть, что скобки кода не порождают, поэтому буду продолжать так же и дальше). И вот оно, решение поставленной задачи — мы видим строку символов, 0х45678, которые однозначно определяют нам искомое число, но в такой форме, что мы (не знаю, как Вы, я то точно) ничего внятного об этом числе сказать не можем. Я думаю, что академик Карналь, который мог указать на ошибку в перфоленте с исходным кодом, и справился бы с данной задачей, но не все столь продвинуты, поэтому продолжим.
Попытаемся получить информацию в более понятной форме.
Для этого вернемся к формату ПТ (далее по тексту я рассматриваю только float), представляющем из себя набор битов, из которых можно извлечь (по определенным правилам) три набора битов для представления трех чисел — знака(з), мантиссы(м) и порядка(п), и закодированное этими числами искомое число будет определяться следующей формулой Чз * Чм * Чп. Здесь символами Ч обозначены числа, представленные соответствующим набором бит, поэтому, для того, чтобы найти искомое число, нам надо знать законы, по которым из исходного набора битов мы извлекаем эти три набора, а также вид кодирования для каждого из них.
В решении данной задачи обращаемся к стандарту IEEE и выясняем, что знак представляет собой один (старший) бит исходного набора и формула для кодирования Чз=(-1)^В(0). Порядок занимает следующие за старшим 8 бит, записан в коде со смещением 127, и представляет собой степень двух, тогда Чп=2^(Ч2(8)-127). Мантисса занимает следующие за порядком 23 разряда и представляет собой число Чм=1+Ч2(23)/2^23.
Теперь мы имеем все необходимые данные и вполне можем решить поставленную задачу — создать строку с символами, которая при определенном прочтении представит число, равное закодированному в ПТ. Для этого нам следует путем несложных операций извлечь вышеуказанные числа, а затем вывести их на печать, снабдив необходимыми атрибутами. Будем считать, что преобразовывать в строку символов целое число с не более, чем 32 битами мы умеем, тогда это совсем несложно.
К сожалению, мы находимся всего лишь в начале пути, поскольку мало кто из читателей этого поста в записи "+1.625*2^3" узнает несчастливое число, которое кодируется более привычным десятичным «13», а уж угадать в записи «1.953125*2^9» несложное «1Е3» или «1*10^3» или же совсем привычное «1000» способны вообще единицы людей, я к ним точно не отношусь. Странно как то получилось, ведь мы выполнили исходное задание, что лишний раз свидетельствует, как внимательно следует относиться к формулировкам. И дело не в том, что десятичная запись лучше либо хуже двоичной (в данном случае имеется в виду двойка в основании степени), а в том, что к десятичной мы привыкли с детства и переделать людей намного сложнее, чем программу, поэтому будем приводить нашу запись к более привычной.
С точки зрения математики нами предстоит несложная операция — имеется запись ПТ=(-1)^з * м * 2^п, а нам требуется преобразовать ее к виду ПТ = (-1)з' * м' * 10^п'. Приравниваем, преобразуем и получаем (один из возможных вариантов) решения з'=з', м'=м, п'=п*lg(2). Если мы оставим за скобками необходимость умножать на явно иррациональное число (это можно сделать, если число рационализировать, но об этом поговорим позже), то задача выглядит решенной до тех пор, пока мы не увидим ответ, потому что, если запись вида "+1.953125*2^9" представляется нам малопонятной, то запись "+1.953125*10^2.70927" еще менее приемлема, хотя и казалось, что хуже некуда.
Продолжаем улучшать решение и находим следующее решение — уравнения приведения к степени по основанию 10 м'=м * 10^{п * lg(2)}, п'= [п * lg(2)], где фигурными и квадратными скобками обозначены дробная и целая часть некоего числа соответственно. Тогда для рассматриваемого примера получим (1.953125*10^0.7 0927)*10^2=«10*10^2», что уже значительно более приемлемо, хоть и не идеально, но уже можно реализовывать.
Дело за малым, нам надо научиться:
- умножать целое число (п) на заранее известное иррациональное (lg(2)) (это несложно при определенных ограничениях по точности результата);
- брать целую и дробную часть числа с фиксированной точкой (это несложно);
- возводить известное целое (10) в иррациональную степень (мда...);
- умножать целое на произвольное иррациональное («мы упростим вычисления, говорили они ...»).
Тем не менее, попробуем двигаться в этом направлении и рассмотрим то, что сделать несложно, а именно пункт 1. Сразу отметим, что задача эта принципиально нерешаемая и мы не можем вычислить н * lg(2), от слова «совсем» не можем, за исключением тривиального случая н=0 (ну и очевидного случая н=к/lg(10)). Интересное заявление, особенно после утверждения «это несложно», но видимое противоречие снимается фразой «с определенной точностью». То есть мы все таки можем вычислить произведение произвольного целого на известное иррациональное и это несложно для результата с определенной точностью. Например, если нас интересует результат с точностью до одного процента, то представив искомый результат п' = п * lg(2) в виде п * [lg(2) *256 + 1/2] / 256 мы получим значение с нужной нам точностью, поскольку возможная относительная погрешность не может превзойти 1/2/77 = 1/144, что явно лучше, чем требуемая 1/100. Следует принять во внимание одно важное соображение — малая величина относительного отклонения ровно ничего не говорит о поведении функции при применении к ней нелинейного преобразования, а операция взятия целой части очевидно нелинейна. Приведем простой пример [4.501]=5, а [4.499]=4 и, несмотря на то, что относительное отклонение в исходных данных составит 0.002/4.5=0.04%, отклонение результата составит целых 1/4=25%. К сожалению, в общем виде задача не решается вообще, при применении любого алгоритма округления. Можно решить только частный случай, когда входные данные ограничены и, более того, принимают фиксированный набор значений, тогда подбором начального смещения и угла наклона можно получить абсолютно точную, в смысле округления, аппроксимацию.
Для нашего случая такой идеальной аппроксимацией будет функция п'=п*77/256
Прежде чем продолжить проектирование алгоритма, мы должны оценить необходимую нам точность. Поскольку мантисса 24 разрядная, то представимое ей число имеет относительную погрешность 2^-24=2^-4*2^-20=16^-1*(2^10)^-2~(10)^-1*(10^3)^-2=10^-7, что означает 7 точных десятичных цифр. Умножение двух 24 битовых чисел будет достаточно для удержания точности в этом диапазоне (ну почти достаточно). Сразу заметим, что переход к 32 битовым числам (оба сомножителя) уменьшает относительную ошибку более, чем в 100 (256) раз, этот факт нам дальше пригодится.
Самая неприятная в смысле точности формула вычисляет новую мантиссу и выглядит так
м' = м * 10^{п * lg(2)}
Почему она самая неприятная — 1) она содержит цепочку вычислений относительно п и погрешность будет накапливаться, 2) в ней есть крайне плохая с точки зрения точности операция, и это не взятие дробной части, поскольку, если делать ее одновременно со взятием целой части, все не так уж и плохо, а экспонента. Если остальные операции — это умножения и относительные погрешности при этом просто складываются, являются прогнозируемыми и зависят исключительно от длины разрядной сетки представления операндов, то для экспоненты все крайне плохо и совершенно очевидно, что относительная погрешность будет весьма велика при больших значениях аргумента.
«Ну да, совершенно очевидно, что»
q(10^x) = Δ(10^x)/10^x = (10^(x +Δx) — 10^x)/10^x = 10^Δx -1 = 10^(x*qx)-1,
10^(x*qx) >~ 10^(x*0) + (10^(x*0))'*qx = 1 + x*ln(10)*10^(0)*qx = 1+x*ln(10)*qx,
отсюда получаем
Что означает приведенное выражение — то, что на краю диапазона значений, при п=127, относительная погрешность вырастет в 292 раза и для того, чтобы удержать точность результата в требуемом пределе, нам надо точность аргумента существенно поднять.
Вспомним, что переход от 24 к 32 бит дает нам требуемое увеличение точности (не совсем дает, но весьма близко), отсюда понимаем, что первое умножение (п*lg(2)) следует проводить с 32 разрядными операндми, то есть с такой точностью следует выразить логарифм двойки, тогда он будет равен 1'292'914'005/2^32. Отметим, что в коде числитель этой константы следует записываться, как (int)((lg(2)*float(2^32))+0.5), но ни в коем случае не как загадочное 0х4d104d42, пусть даже и с комментарием относительно его вычисления, поскольку хорошо написанный код самодокументирован.
Далее нам потребуется целая часть результата, это несложно, поскольку мы точно знаем положение десятичной точки в обоих сомножителях и, соответственно, в результате.
А вот далее нам надлежит вычислить 10 в степени от 0 до 1 и здесь для получения требуемой точности прибегнем к небольшому трюку. Поскольку, согласно формуле для погрешности, точность на правом краю диапазона падает более, чем двое, то если мы представим значение аргумента как сумму логарифма десятичного двух и некоторого остатка, п''=lg(2)*i+(п''-lg(2)*i), то первый член суммы приведет к умножению на 2 в соответствующей степени, что легко осуществить с нулевой погрешностью (пока не наступит переполнение), а оставшийся будет ограничен значением lg(2) и мы не потеряем в точности при вычислении 10^п'' (при условии правильного вычитания).
Тем не менее, экспоненциальную функцию для ограниченного значением lg(2) аргумента все равно придется вычислять и единственный путь, который я вижу, это разложение в ряд Тэйлора. К сожалению, он не слишком быстро сходится на нашем диапазоне, и, например, для достижения точности в 10Е-7 нам потребуется 9 членов суммы, что приводит к необходимости осуществить 1+9*2=19 умножений 32-разрядных целых чисел, что несколько превышает желательные показатели. Также все еще остаются смутные сомнения относительно нашей способности провести вычисление п'=п*lg(2) с такой точностью, чтобы ее хватило при максимальном значении п.
Тем не менее агоритм получился вполне работоспособный и нам требуется для получения результата только 32-разрядное умножение в количестве 1+19+1=21 операций что и определяет вычислительную сложность алгорима
Можно ли уменьшить вычислительную сложность нашего преобразования — вроде бы нет, мы все аккуратно посчитали — но внезапно оказывается, что все таки можно. Несколько неожиданное заявление, ключ к пониманию такой возможности лежит в природе порядка ПТ — он принимает фиксированный (и относительно небольшой) набор значений, а мы с Вами при выводе формул преобразования это не учитывали, а неявно работали с непрерывной величиной.
Простейшее решение — меняем память на время — заранее вычислить для всех возможных (2^8=256) показателей степени п соответствующие значения [п'] (наиболее подходящего показателя степени 10) и {п'} (коректирующего множителя для мантиссы), занести их в таблицу и далее просто использовать в процессе вычислений. Формула получается совсем несложной — ПТ=м*2^п=м*10^п'*(2^п/10^п')=(м*(2^п/10^п'))*10^п'.
В простейшем случае нам потребуется 256*3 (корректирующий множитель из 24 разрядов, больше не нужно) + 256*1 (порядок по основанию 10 гарантированно меньше порядка по основанию 2) = 1кбайт констант. При этом нам остается сделать всего лишь одно умножение 24*24 разряда (скорее всего это будет 32*32), что существенно ускоряет работу по сравнению с вариантом настоящего вычисления.
Посмотрим, что можно сделать с точки зрения экономии памяти (при этом нам опять придется платить временем, так что ищем разумный компромисс). Прежде всего, если мы будем учитывать знак порядка отдельно, то можно обойтись только половиной требуемой памяти (из 256 байт для порядка 10) и проинвертировать результат в случае необходимости. К сожалению, с корректирующим множителем так просто не получится, поскольку
2^-п/10^-п' = 1/(2^п/10^п') != 2^п/10^п',
а жаль. Мы либо должны оставить длинную таблицу, либо для отрицательных показателей проводить деление на константу для положительных показателей. Конечно, деление — это не 18 умножений, но все равно по скорости оно двум умножениям точно эквивалентно, так что время точно возрастет в два раза, чтобы сэкономить память в два раза, до 512 байт. Стоит ли оно того — вопрос непростой, но, к счастью, у нас есть значительно более красивый способ, который позволяет избавиться от страданий выбора.
Способ этот в общем случае называется кусочно-линейная аппроксимация и заключается в задании констант не для каждой точки исходного значения, а лишь для некоторых и вычисления отсутствующих значений (с требуемой точностью) с использованием заданых значений по несложной формуле. Применительно к нашей задаче получается (не учитывая знак)
ПТ=м*2^п=м*2^(п0+п1)=м*10^п'*(2^(п0+п1)/10^п')=м*(2^п0/10^п')*2^п1*10^п',
где п0-некоторое опорное значение, а п1=п-п0. Тогда новая мантисса вычисляется путем умножение двух чисел с фиксированной точкой с последующими сдвигом результата, который точности не ухудшает.
Тогда возникает законный вопрос — а зачем нам вообще таблица, ведь можно в качестве п0 взять минимальный показатель и обойтись всего одним значением корректирующего множителя? К сожалению, такой подход контрпродуктивен в силу двух взаимодополняющих обстоятельств — необходимости получить наиболее подходящий показатель степени 10 и появлением очень длинных сдвигов при подобном подходе. Последнее обстоятельство наводит на мысль о границах применимости подобного метода — если мы проводим умножение 32*32, а исходная мантисса имеет 24 разряда, то сдвиг на 8 разрядов не приведет к переполнению и нам потребуется одна опорная точка на 8 значений двоичного порядка. Тогда общий объем требуемой памяти составит 256/8*4=32*4=128 байт — хорошая экономия памяти ценой времени выполнения из за необходимости сдвига целого результата произведения на максимум 8 разрядов.
Можно еще немного сократить объем констант за счет симметричности показателя степени относительно 0, о чем я говорил ранее, но экономия составит 32/2=16 байт, не уверен, что это оправдает усложнение (и увеличение длины кода) собственно программы.
Кстати, недавно смотрел код широко известной в узких кругах библиотеки adafruit и был слегка удивлен следующим фрагментом кода
const UINT8 Bits[] = {0x01, 0x02, 0x04, 0x08, 0x10, 0x20, 0x40, 0x80};
... data = data | Bits[n];с комментарием, что операция 1 << n на AVR исполняется долго. В своем посте я уже показывал, какие чудеса творит компилятор при константной параметре, но это не тот случай.
Мне показалось сомнительным, чтобы взятие битовой маски из массива было быстрее, нежели выполнение непосредственно операций сдвига, и последующий анализ кода (с помощью сайта godbolt, хотя и крайне маловероятно, чтобы его создатель читал Хабр, тем не менее в очередной раз приношу ему свою искреннюю благодарность) показал, что это действительно так.
Сгенерированный компилятором код для обоих вариантов (вот правильный вариант со сдвигами с учетом особенности задачи, ведь нам нужен только 1 бит)
ldi r18,lo8(4)
sbrs r25,1
ldi r18,lo8(1)
sbrc r25,0
lsl r18
sbrs r25,2
swap r18занял ровно одинаковое место в памяти, а если сделать все аккуратно на ассемблере, то вариант с индексом вырывается вперед 8:7 за счет лишних 8 байт программы (разумеется, если мы не воспринимаем всерьез воистину восхитительное решение с отдельным хранением инвертированной маски, которое обойдется в 16 байт — и ЭТО используют повсеместно — «я знал, что будет плохо, но не знал, что так скоро»). Ну да упомянутый пакет — это вообще отдельная песня, как нельзя лучше описываемый следующей цитатой из одной замечательной книги: «Этот замок просился на обложку труда по фортификации с подписью „Как не надо строить замки или найди 12 ошибок“ (»Последний кольценосец", если кто не читал, рекомендую).
Вернемся к нашим баранам с плавающей точкой и построим результирующую формулу
ПТ = м * 2^п=(м * пк[п/8]) * 2^(п%8) * 10^пп[п/8],
где квадратные скобки означают взятие элемента массивов пк-коррекция показателя и пп-порядок показателя. Сразу видна вычислительная сложность алгоритма, которая определяется умножением 32*32(24*24) и последующими сдвигами. Дальше можно учесть возможность объединения в одном 32 разрядном слове показателя степени 10 и корректирующего множителя, это оставим на долю пытливого (и терпеливого, ведь он дочитал до конца) читателя данного поста.
Единственное замечание напоследок — когда будем создавать таблицу констант, ни в коем случае нельзя это делать в следующем стиле
const uint32_t Data[32] PROGMEM = { 0xF82345,… }
и дело, конечно, не в атрибутах описания массива, а в самих данных в виде магических чисел. Как справедливо отмечали авторы, точно не глупее меня, хорошо написанный код самодокументирован и, если мы напишем вышеприведенную константу (и остальные) в виде
#define POWROUD(pow) ((uint8_t)((pow & 0x07)*log(2)+0.5))
#define MULT(pow) (2^pow / 10^POWROUND(pow))
#define MULTRAW(pow) (uint32_t((MULT(pow) << 24) +0.5))
#define BYTEMASK 0xFF
#define POWDATA(pow) ((POWROUND(pow) & BYTEMASK)|
(MULTRAW(pow) & (~BYTEMASK)))
const uint32_t Data[(BYTEMASK/8)+1] = { POWDATA(0x00),POWDATA(0x08), ..POWDATA(0xF8)}то нам никто не пришлет недоуменных вопросов, а если кто и пришлет, мы можем точно на них не отвечать, все равно бесполезно.
Можно предложить модификацию данного метода, в которой подходящая степень десяти будет вычисляться не для правого края сегмента, а для левого и тогда результат будет сдвигаться не вправо, чтобы учесть степень двойки, а влево. С точки зрения математики способы абсолютно равнозначные. Посмотрим на полученный результат:
1.953125*2^9=1.953125*2^(8+1)=1.953125*42949673/256/256/256(2.56)*2*10^2=10*10^2
тут уже очень легко узнать 1000. Конечно, надо бы еще преобразовать полученные мантиссу и порядок в строки, аккуратно выполнить округление, подогнать результат к требуемому формату, добавить знак, учесть особые случаи и так далее, но это все уже не так интересно, главную часть преобразования мы выполнили.
Автор: GarryC
