Приближенное вычисление чисел с плавающей точкой, sin(x), cos(x) и sqrt(x*x + y*y)
06.12.2014 22:55
7
Потребовалось мне тут реализовать дискретное преобразование Фурье на микроконтроллере.
Пришлось вспоминать его формулу и математику. А там синусы, косинусы, да квадратные корни. А процессор умеет только складывать, вычитать, умножать и делить, да еще и только 16-битные целые числа.
Вспомнив, что там нам про тригонометрические функции рассказывали в институте, удалось сделать быстрое их вычисление с помощью совсем небольшой таблички (ибо памяти в микропроцессоре тоже кот наплакал).
А порывшись по закромам интернета, да по научным сайтам удалось наковырять и приближенное вычисление для sqrt(x2 + y2), причем такое, что аж завидовать стал математикам.
Плавающая точка
По скольку, что такое плавающая точка наш процессор никогда не слышал, эмулируем ее с помощью очень простого хака - умножаем все числа на фиксированную величину, например 1000, дробную часть откидываем - вот и все - дальше можно работать как с обычными числами!Единственное но - при выводе куда-либо нужно не забыть провернуть обратную операцию.
Подбирать множитель нужно так, чтобы и точность сохранить (после точки) и не вылезти за пределы 16 бит при больших числах.
Вычисление cos(x) и sin(x)
Аргумент - в радианах.Для градусов элементарно преобразовывается: αрад= α° * π / 180
Алгоритм прост до безобразия - вычисляем таблицу только для одной четверти круга (0 - π/2) и только для cos, а затем просто приводим аргумент к этому интервалу, подставляя нужные знаки и беря результат из нашей таблички.
Если нужно вычислить sin - просто сдвигаем угол на π/2 - это и будет sin (sin и cos можно выражать через друг - друга и отличаются они как раз на этот угол).
Вычисление очень быстрое, потребляет мало памяти и дает достаточно точный (да еще и настраиваемый) результат.
Пример данной функции на JavaScript:
//Big Math Multiplier - наш множитель //чем больше - тем выше точность, но и выше шанс на переполнение var bmm = 1000; //Константы var pi = Math.floor( Math.PI * bmm ); var dpi = 2 * pi; var hpi = Math.floor( pi / 2 ); var thpi = Math.floor( ( 3 * pi ) / 2 ); //сколько вычислений-точек на всю четверть круга будем хранить? //чем больше - тем выше точность, но больший расход памяти var points = 128; var cosSector = new Array( points ); var rpoints = points - 1; for ( var i = 0; i <= rpoints; i++ ) { cosSector[i] = Math.floor( Math.cos( ( ( i * Math.PI ) / 2 ) / rpoints ) * bmm ); } function cos( r, sin ) { //приводим аргумент к значению от 0 - 2*pi (весь тригонометрический круг) r = r % dpi; //если хотим вычислить sin (а табличка у нас для cos) - то просто //сдвигаем на pi/2 - это и будет sin if ( sin ) r = hpi - r; //если угол отрицателен, прибавляем 2*pi, //чтобы преобразовать его в положительный if ( r < 0 ) r += dpi; var m, a; //для угла от 0 <= a < pi if ( r >= 0 && r < hpi ) { m = false; a = r; } //для угла pi/2 <= a < pi else if ( r >= hpi && r < pi ) { m = true; a = hpi - ( r - hpi ); } //для угла pi <= a < 3pi/2 else if ( r >= pi && r < thpi ) { m = true; a = r - pi; } //для угла 3pi/2 <= a < 0 else if ( r >= thpi && r < dpi ) { m = false; a = hpi - ( r - thpi ); } //естественно, при реализации на микроконтроллере //вся математика целочисленная и никакой floor не нужен var v = Math.floor( ( a * rpoints ) / hpi ); v = cosSector[ v ]; return m ? -v : v; } function sin( r ) { return cos( r, true ); }
Использовать так: cos( значение в радианах * bmm ) и точно также sin.
Естественно, результат будет в единицах bmm:
cos( значение в радианах * bmm ) / bmm ≈ Math.cos( значение в радианах )
sin( значение в радианах * bmm ) / bmm ≈ Math.sin( значение в радианах )
sqrt(x2 + y2) он же sqrt(x*x + y*y)
Достаточно часто применяемая формула (например, для вычисления расстояний).Можно решить "в лоб", использовав приближенную формулу для квадратного корня, однако есть проблема: x*x + y*y - достаточно большое число и почти наверняка вызовет переполнение.
Оказывается, все уже давно решено за нас - формула поражает своей простотой и тем, насколько она экономит тактов процессора.
Если вдруг кто-то не знает >> 1 - сдвиг на один бит вправо - равносилен целочисленному делению на 2
function sqrtXXYY( x, y ) { //убираем знак у аргументов if ( x < 0 ) x = -x; if ( y < 0 ) y = -y; return x > y ? x + ( y >> 1 ) : y + ( x >> 1 ); }
Вот так вот все просто :)
06.12.2014, Protocoder
Конечно, можно и Math.abs() юзать и даже быстрее будет чем if, но потом все равно в проц укладывать, а там его нет :)
Я не пойму как с помощью твоего алгоритма вычислить cos(x) и sin(x).
У меня есть x=0,3959267453839191 и собственно нужен косинус и синус этого x.
Куда вставлять этот х в твой код или как им вообще пользоваться???
До этого я пробовал стандартные методы Math.cos() и Math.sin(), но точность они выдают никудышную, а как применить твой чудо код так и не понял, а главное как потом увидеть результат вычислений?
Скорее всего вопрос банальный, но я почти ноль, самоучка и спросить не у кого.
Если не в радианах - его требуется перевести в радианы.
Во-вторых, если Вас не устраивает качество cos/sin в JavaScript, то моя функция тем более Вам не поможет - она использует их для создания таблицы, по которой считает + огрубляет результат (т.к. идет интерполяция).
В-третьих, код использовать очень просто (добавил в статью) - например для cos это будет:
cos( значение в радианах * bmm )
Значение вернется также в bmm, если хотите перевести его в обычное с плавающей точкой - просто разделите его на bmm.
Результат можно увидеть с помощью команды alert примерно так (если x - радианы):
alert( cos( x * bmm ) )
Можно ознакомиться с тем, как их отлаживал автор на javascriptе (http://model.exponenta.ru/k2/Jigrein/JS/fwlink.htm#B2C7), а затем переписал программу для ARM-процессора (http://model.exponenta.ru/k2/Jigrein/md_106d.htm).
В частности, макросы вычисления функций IQsin и IQcos основаны на нескольких членах степенного ряда.
то точность вычисления увеличится почти в 2 раза.