Вместо вступления
В статье алгоритмическая оптимизация функции sin() для бюджетных микроконтроллеров stm32, повышающая производительность в 10 и более раз.
Тригонометрические функции, характеризующиеся высоким потреблением процессорного времени, могут негативно влиять на выбор бюджетных микроконтроллеров ( без модуля FPU ) для задач, где важна скорость счёта, например, контроль пространственного положения.
Библиотечные тригонометрические функции с двойной точностью в два раза медленнее, чем с одинарной. Это досадный факт.
На иллюстрации время работы «sin()» и «sinf()», измеренное в тактах CPU, на диапазоне [-420° .. +420°] с шагом 5°, плюс слабый шум.
Производительность библиотечных «sin()» и «sinf()» высока вблизи точки «0» и снижается в два раза пропорционально модулю аргумента.
Возникает потребность считать тригонометрию несколько быстрее, без FPU, но допуская худшую точность:
у аргумента — в пределах чувствительности доступных электронных датчиков положения, например, LSM6DSO;
у результата — в пределах точности доступных шаговых двигателей и сервомашин, например, HS-5082MG.
Требования к алгоритму
Требования к исходному коду:
a) Исходный код без интеллектуального балласта третьих лиц;
b) Самодостаточный компилируемый модуль на основе стандартных библиотек;
c) Кросс-платформенность — отсутствие вставок на ассемблере.Требования к исполняемому коду:
a) Ограниченный объём бинарного кода — чем меньше, тем лучше;
b) Ограниченный объём статических и динамических данных — чем меньше, тем лучше.Требования к интерфейсу (привычный):
a) float sint( float x );
b) float cost( float x ).Требования к аргументу:
a) Объектная сущность: геометрический плоский угол;
b) Единица измерения: угловой градус в десятичном виде: 1,0 °;
c) Точность измерения: +/- 0,05 °;
d) Рабочий диапазон: +/- 5400 °;
e) Тип данных: float, согласно IEEE 754.Требование к результату:
a) Объектная сущность: тригонометрические функции «sin» и «cos»;
b) Точность счёта: +/- 0,005;
c) Тип данных: float, согласно IEEE 754.
Применяемая методика
Сравнительный анализ способов приближения тригонометрических функций за границами настоящей публикации.
Применяемый метод считает синус плоского угла в рабочем диапазоне от 0 до ½π.
Синусы углов за пределами рабочего диапазона находятся через осевую симметрию применительно к результату из рабочего диапазона.
Вычисление синуса угла в рабочем диапазоне заменяется поиском точки на куске одномерного сплайна:
В некоторых условиях сплайн и синус демонстрируют схожее поведение.
Аппроксимирующий сплайн представляется формой Безье. Точка сплайна вычисляется геометрически, редуцированным методом де Кастельжо.
Алгоритм оптимизируется на соответствие функции синус.
Косинус считается как синус, сдвинутый по фазе на 1½ π.
Достигнутый результат
Характеристики бинарного кода на платформе Cortex M0:
размер бинарного кода — 168 байт;
размер стека — 40 байт (локальные переменные — 24 байта);
статическая и динамическая память — без надобности.
Результат с тремя точными знаками после запятой и гарантированной ошибкой в четвёртом знаке — | 0,0003 |.
Скорость работы вновь созданной функции на платформе Cortex M0 превосходит библиотечные «sin» и «sinf» в 10 и более раз.
У алгоритма есть опции, выбираемые через условную компиляцию:
медленный счёт — гарантирует ошибку |0,0003| и менее;
быстрый счёт — даёт преимущество в ~20 тактов на Cortex M0, ухудшая ошибку до |0,0004|.
Преимущество в 15-25 тактов ( ~1us ) — это существенный прирост производительности для бюджетных контроллеров на малых частотах.
На платформе Cortex M4 быстрый метод счёта даёт сомнительное преимущество в 3-5 тактов, соразмерное погрешности измерения, и ухудшение точности до |0,0004|.
В целом, тестирование алгоритма на Cortex M4 показало приемлемый результат с малым превосходством по скорости библиотечной функции «sinf», оптимизированной под имеющийся FPU.
На графике отсутствует функция «sin» [ double sin ( double x ) ], показывающая худшее на десятичный порядок время.
На платформе Cortex M4 каждый вызов «sin» [ double sin ( double x ) ] забирает порядка 2500 тактов.
Итоговый алгоритм допускает аргумент в угловых градусах и в радианах.
Выбор единицы измерения аргумента осуществляется через условную компиляцию.
Вместо заключения
Алгоритм соответствует техническим требованиям, показывая удовлетворительную скорость и ожидаемую точность счёта.
Есть некоторый потенциал для улучшения характеристик алгоритма. С практической точки зрения этот потенциал целесообразно отложить в будущее.
Прикладная библиотека
Файл заголовок
/******************************************************************************
* File: cosint.h Created on 31 мар. 2022 г.
*
* Copyright (c) 2022, "Nikolay E. Garbuz" <nik_garbuz@list.ru>
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License version 3 as
* published by the Free Software Foundation.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
* Authored by "Nikolay E. Garbuz" <nik_garbuz@list.ru>
* Modified by
*
* TAB Size .EQ 4
********************************************************************************/
#ifndef _COSINT_H_
#define _COSINT_H_
#ifdef __cplusplus
extern "C"
{
#endif
//**********************************************************
//---------------- module interface ------------------------
/**
* @brief Fast sine function for approximate calculations
* @param Geometric angle between -5400.00 and 5400.00 degrees
* with two correct numbers after the decimal separator
* @retval Sine with three correct numbers after the decimal separator
*/
float sint( float x );
/**
* @brief Fast cosine function for approximate calculations
* @param Geometric angle between -5400.00 and 5400.00 degrees
* with two correct numbers after the decimal separator
* @retval Cosine with three correct numbers after the decimal separator
*/
float cost( float x );
//**********************************************************
//---------------- module options --__----------------------
/*
* The [COSINT_RAD] directive sets the angle units.
* By default, module functions take a parameter measured in degrees.
*
* Remove or comment out the [#undef COSINT_RAD] line
* to switch the interface [cosint] from degrees to radiant.
*/
#define COSINT_RAD
#undef COSINT_RAD
/*
* The [COSINT_FAST] directive defines the linear proportion method.
* The default is the slow method with a maximum error modulus is |0.0003|.
*
* Remove or comment out the [#undef COSINT_FAST] line to switch module [cosint] to fast mode.
* About ~1.6% speedup, maximum error modulus is 0.0004
*/
#define COSINT_FAST
#undef COSINT_FAST
#ifdef __cplusplus
}
#endif
#endif /* _COSINT_H_ */Исходный код
/******************************************************************************
* File: cosint.c Created on 31 мар. 2022 г.
*
* Copyright (c) 2022, "Nikolay E. Garbuz" <nik_garbuz@list.ru>
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License version 3 as
* published by the Free Software Foundation.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
* Authored by "Nikolay E. Garbuz" <nik_garbuz@list.ru>
* Modified by
*
* TAB Size .EQ 4
********************************************************************************/
#pragma GCC push_options
#pragma GCC optimize ("O3")
#include "cosint.h"
#include <stdint.h>
//**********************************************************
//---------------- module implementation -------------------
typedef enum
{
coname = -1, siname
} EFName;
inline float cosint( float x, EFName fn );
//---------------- exported module functions ---------------
float sint( float x )
{
return cosint( x, siname );
}
float cost( float x )
{
return cosint( x, coname );
}
//**********************************************************
//------------- helper data and function -------------------
#ifndef M_PI
# define M_PI 3.14159265358979323846 /* pi */
#endif
#ifdef COSINT_RAD
# define R_UNIT M_PI
#else
# define R_UNIT 180.
#endif
#define R_TU ( R_UNIT / 2. ) // The real measure unit in use
#define M_PRC 0x0A // The precision of calculations
#define M_MAX 0x1F // The bit range of calculations
#define M_DAT ( M_MAX >> 1 ) // The bit range of a data
#define M_SGN 0x02 // The sign mask for a result of calculations
#define M_QTR 0x01 // The mask of a circle quarter
#define M_ASQ ( M_SGN | M_QTR ) // The mask of common characteristics for the angular data
#define M_ADV ( 1 << M_DAT ) // The allowed data value
#define M_ADM ( M_ADV - 1 ) // The mask of allowed data value
#define M_MTR (int32_t)( ( M_ADV << M_PRC ) / R_TU )
// Slow and precise, maximum error modulus is |0.0003|
#define DO_SLOW( t, pw, A, B ) ( ( ( t * ( B - A ) + ( ( 1 << pw ) - 1 ) ) >> pw ) + A )
// Fast and rough, about ~1.6% speedup, maximum error modulus is |0.0004|
#define DO_FAST( t, pw, A, B ) ( ( ( t * ( B - A ) ) >> pw ) + A )
#ifdef COSINT_FAST
# define LN_ DO_FAST
#else
# define LN_ DO_SLOW
#endif
// Control data for the spline
#define Ax 0x0000
#define Bx 0x320B
#define Cx 0x655D
#define Dx 0x8000
// helper function
float cosint( float x, EFName fn )
{
int32_t sq, tm;
int32_t A, B, C, D;
// Converting from real to integer
tm = M_MTR;
tm *= x;
tm = tm >> M_PRC;
// Choosing a function via a phase shift by it name
tm += fn & ( M_ADM );
// Extracting a quarter of the angle and a sign of the result
sq = ( tm >> M_DAT ) & M_ASQ;
// Removing periods and reverse from the angle value
// Transforming the angle value to the first quarter
if ( sq & M_QTR )
tm = ( tm & M_ADM ) ^ M_ADM;
else
tm &= M_ADM;
// Calculating a result thru a line proportion
A = Dx;
B = LN_( tm, M_DAT, Cx, Dx );
C = LN_( tm, M_DAT, Bx, Cx );
D = LN_( tm, M_DAT, Ax, Bx );
A = LN_( tm, M_DAT, B, A );
B = LN_( tm, M_DAT, C, B );
C = LN_( tm, M_DAT, D, C );
A = LN_( tm, M_DAT, B, A );
B = LN_( tm, M_DAT, C, B );
A = LN_( tm, M_DAT, B, A );
// Correcting a numeric sign of the result
x = ( sq & M_SGN ) ? -A : A;
// Converting from integer to real and return
x /= M_ADV;
return x;
}
#pragma GCC pop_options