STM32. About sine

Instead of an introduction

Trigonometric functions, characterized by high consumption of processor time, can negatively affect the choice of budget microcontrollers (without FPU module) for tasks where counting speed is important, for example, attitude control.

Double-precision library trigonometric functions are twice as slow as single-precision trigonometric functions. This is an unfortunate fact.

In the illustration, the operating time of “sin()” and “sinf()”, measured in CPU cycles, on the range [-420° .. +420°] in 5° increments, plus low noise.

The performance of the library “sin()” and “sinf()” is high near the point “0” and decreases by a factor of two in proportion to the modulus of the argument.

There is a need to calculate trigonometry somewhat faster, without FPU, but allowing for worse accuracy:

argument – within the sensitivity of available electronic position sensors, for example, LSM6DSO;
the result is within the accuracy of available stepper motors and servos, for example, HS-5082MG.

Algorithm Requirements

Source Code Requirements:
a) Source code without third party intellectual ballast;
b) Self-contained compiled module based on standard libraries;
c) Cross-platform – no assembler inserts.
Requirements for the executable code:
a) Limited amount of binary code – the less, the better;
b) Limited amount of static and dynamic data – the less the better.
Interface requirements (usual):
a) float sint( float x );
b) float cost( float x ).
Argument requirements:
a) Entity: geometric flat corner;
b) Unit of measure: Angular degree Decimal: 1.0°;
c) Measurement accuracy: +/- 0.05°;
d) Operating range: +/- 5400°;
e) Data type: float, according to IEEE 754.
Result requirement:
a) Entity: trigonometric functions “sin” and “cos”;
b) Counting accuracy: +/- 0.005;
c) Data type: float, according to IEEE 754.

Applicable methodology

Comparative analysis of methods for approximating trigonometric functions is beyond the scope of this publication.

The applicable method considers the sine of a flat angle in the operating range from 0 to ¼π.

The sines of angles outside the working range are found through axial symmetry in relation to the result from the working range.

The calculation of the sine of the angle in the working range is replaced by the search for a point on a piece of a one-dimensional spline:

Under certain conditions, spline and sine exhibit similar behavior.

The approximating spline is represented by a Bezier shape. The spline point is calculated geometrically by the de Castelljou reduced method.

The algorithm is optimized to match the sine function.

The cosine is considered as a sine shifted in phase by ¾ π.

Achieved result

Characteristics of the binary code on the Cortex M0 platform:

binary code size – 168 bytes;
stack size – 40 bytes (local variables – 24 bytes);
static and dynamic memory – unnecessarily.

Result with three exact decimal places and a guaranteed error in the fourth decimal place — | 0.0003 |.

The speed of the newly created function on the Cortex M0 platform exceeds the library “sin” and “sinf” by 10 or more times.

The algorithm has options selected via conditional compilation:

slow counting – guaranteed error |0.0003| and less;
fast counting – gives an advantage of ~20 cycles on the Cortex M0, worsening the error to |0.0004|.

The advantage of 15-25 cycles ( ~1us ) is a significant performance boost for low-end controllers at low frequencies.

On the Cortex M4 platform, the fast counting method gives a dubious advantage of 3-5 cycles, commensurate with the measurement error and a deterioration in accuracy to |0.0004|.

In general, testing the algorithm on the Cortex M4 showed an acceptable result with a slight superiority in speed of the “sinf” library function, optimized for the existing FPU.

There is no “sin” function on the chart [ double sin ( double x ) ]showing the worst time in decimal order.

On the Cortex M4 platform, each “sin” call [ double sin ( double x ) ] takes about 2500 cycles.

The final algorithm accepts arguments in degrees of angle and in radians.
The choice of the argument unit is done through conditional compilation.

Instead of a conclusion

The algorithm meets the technical requirements, showing a satisfactory speed and the expected accuracy of the calculation.

There is some potential to improve the performance of the algorithm. From a practical point of view, it is advisable to postpone this potential for the future.

Application Library

File header

/******************************************************************************
 * 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_ */

Source

/******************************************************************************
 * 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

452