Приближенное вычисление чисел с плавающей точкой, sin(x), cos(x) и sqrt(x*x + y*y)


Потребовалось мне тут реализовать дискретное преобразование Фурье на микроконтроллере.

Пришлось вспоминать его формулу и математику. А там синусы, косинусы, да квадратные корни. А процессор умеет только складывать, вычитать, умножать и делить, да еще и только 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 ) / bmmMath.cos( значение в радианах )
sin( значение в радианах * bmm ) / bmmMath.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
Дрозд Бункерный07.12.2014 11:21:40#ответить
Ты такой хардкорный программист, что даже Math.abs() не используешь :)
Protocoder07.12.2014 14:25:43#ответить
С ассемблера привычка осталась :D

Конечно, можно и Math.abs() юзать и даже быстрее будет чем if, но потом все равно в проц укладывать, а там его нет :)
JJ15.04.2015 00:09:08#ответить
Слушай, Protocoder, ну помоги начинающему, "программисту".

Я не пойму как с помощью твоего алгоритма вычислить cos(x) и sin(x).
У меня есть x=0,3959267453839191 и собственно нужен косинус и синус этого x.

Куда вставлять этот х в твой код или как им вообще пользоваться???

До этого я пробовал стандартные методы Math.cos() и Math.sin(), но точность они выдают никудышную, а как применить твой чудо код так и не понял, а главное как потом увидеть результат вычислений?

Скорее всего вопрос банальный, но я почти ноль, самоучка и спросить не у кого.
Protocoder15.04.2015 02:41:52#ответить
Ну, во-первых, в чем у вас x? в градусах? радианах?
Если не в радианах - его требуется перевести в радианы.

Во-вторых, если Вас не устраивает качество cos/sin в JavaScript, то моя функция тем более Вам не поможет - она использует их для создания таблицы, по которой считает + огрубляет результат (т.к. идет интерполяция).

В-третьих, код использовать очень просто (добавил в статью) - например для cos это будет:
cos( значение в радианах * bmm )

Значение вернется также в bmm, если хотите перевести его в обычное с плавающей точкой - просто разделите его на bmm.

Результат можно увидеть с помощью команды alert примерно так (если x - радианы):
alert( cos( x * bmm ) )
Борис25.04.2017 10:23:41#ответить
Функции: IQsin, IQcos, IQatan2, IQsqrt, IQdiv необходимы, в том числе, и при написании цифровых систем управления электропривода.

Можно ознакомиться с тем, как их отлаживал автор на javascriptе (http://model.exponenta.ru/k2/Jigrein/JS/fwlink.htm#B2C7), а затем переписал программу для ARM-процессора (http://model.exponenta.ru/k2/Jigrein/md_106d.htm).

В частности, макросы вычисления функций IQsin и IQcos основаны на нескольких членах степенного ряда.
Protocoder25.04.2017 20:23:38#ответить
Спасибо - действительно интересно.
Написать комментарий