Вероятно, в сети можно найти множество статей, посвященных вычислениям с фиксированной точкой. По роду своей деятельности, я разработчик электроники и программист микроконтроллеров. Математика, с фиксированной точкой в моих задачах практически не встречалась.
Однако.. Имеется фрезерный станок ЧПУ. Было принято решение, переделать его систему управления, используя ядро Cortex M3. В общем и целом, существуют программные пакеты, позволяющие сделать это. Например GRBL. Однако, у меня своя специфика. Для реализации ПО, требовались математические операции sin(x), cos(x), sqrt(x). И я задумался над их реализацией в формате фиксированной точки.
Для начала, приступил к реализации функций sin(x), cos(x). Да, конечно, можно использовать табличные методы, но я пошел другим путем и решил реализовать аппроксимацию, используя ряд Тейлора. Оценил участок, на котором буду раскладывать в ряд и приступил к реализации. Для разложения в ряд был выбран отрезок [0..PI/4], так как на этом отрезке, значение аргумента функции укладывается в промежуток [0..1], так и результат функции укладывается в отрезок от 0 до 1. Если быть более точным, то [0..(2^32)-1]/2^32. Таким образом, со всеми числами можно работать как с 32-битными положительными дробными числами на интервале [0..0.9999(9)]. Единственное, значение cos(0), принимает значение 1.0., но я ограничился значением 0.9999(9). Таким образом удалось расширить диапазон значений на два бита. Бит знака и бит единичного разряда. Т.е. при вычислений sin(x), cos(x), на этом отрезке, эффективно используются все 32-бита, а не 31 или 30. Откуда 31 или 30? Знаковое значение -1.0, скушало бы один разряд, но максимум был бы 2^31/2^32. А значение 1.0, забрало бы еще один разряд. При этом, младшие биты были бы нулями. Т.е. в сумме, минус 6 децибел динамического диапазона (если я не ошибся). Зачем их терять из за двух значений?
Поскольку sin(x) и cos(x), должны быть реализованы на интервале [-PI..PI] или [0..2*PI], задумался над их склеиванием в единую функцию. В ходе размышлений, появилась одна идея. Все равно, при большинстве вычислений, происходит масштабирование аргумента функций до 2*PI. почему бы не сделать это сразу, внутри функции. А в качестве аргумента использовать число в диапазоне [0..(2^31)-1]/32. Т.к. sin(x)/cos(x) функции периодические, то обрезанное значение 1.0 или 2^32, как раз попадет в 0, т.к. старший бит будет отброшен.
Подобное масштабирование приводит к одному интересному эффекту. Старшие биты числа можно использовать для выбора интервала на котором будет происходить вычисление функции. Т.е. старшие три бита просто определяют, какая из функций будет выполнена. Это позволяет избавиться от операции сравнения. Кроме того, аргумент функции может принимать как отрицательные так и положительные значения, и это практически не повлияет результат вычислений. Так устроен формат чисел.
...
Затем, приступил к реализации вычисления корня квадратного. При вычислении корня, тоже подошел к вопросу нестандартно, и стал аппроксимировать функцию на отрезках кратных степени двойки. Обнаружил ошибка аппроксимации, носит квадратичный характер. Однако, оказалось, что аппроксимирующая ошибку парабола сдвинута относительно центра. Я нашел коэффициенты, используя квадратичную аппроксимацию сдвига по оси X, и получил точность аппроксимации порядка 6.x, при извлечении квадратного корня из 32-битного числа.
После этого, решил повысить точность аппроксимации, аппроксимировав остаточную ошибку. И это привело к точности порядка 0.15, для 32-битного аргумента. Т.е. извлечение квадратного корня из целого 32-битного числа, до целого 16-битного результата происходит без использования итераций. Причем расчеты ошибки могут быть распараллелены (например в FPGA). Поскольку, каждая из ошибок имеет свой порядок, это позволяет обойтись определенным динамическим диапазоном в рамках определенного окна.
И наконец, добавил итеративную часть, которая использует бинарный поиск от 7 до 12 бит из 32-х бит результата (16-целая часть, 16-дробная).
