Introduction
Although some may think logarithms are rarely needed, there are occasions when their fast and (sufficiently) precise computation becomes essential.
The CMSIS DSP software library provides various processing functions for ARM Cortex processor based devices.
In this library, there are functions for calculating natural logarithms based on a Fast Binary Logarithm Algorithm. One such function is arm_vlog_q31(). In this post, we aim to explore how this function works under the hood.
Theoretical Basics
Logarithms follow specific rules that can be used to transform expressions, making it easier to calculate and manipulate them. We will denote some log rules that will be used. The base of two is used in the algorithm because it aligns with the binary system, making calculations simpler.
$$y = \log_2(x) \Rightarrow x = 2^y$$
As we know, ARM CMSIS-DSP Vector Log functions are used for computing natural logarithms (ln). So, what do base 2 logarithms have to do with this? This transformation is essential for the Fast Binary Logarithm Algorithm:
$$\log_a(x) = \frac{\log_2(x)}{\log_2(a)}$$
$$\log_e(x) = \frac{\log_2(x)}{\log_2(e)}$$
$$\ln(x) = \frac{\log_2(x)}{\log_2(e)}$$
$$\ln(x) = \log_2(x) \cdot \frac{1}{\log_2(e)}$$
Here are some values to keep in mind:
$$\log_2(1) = 0$$ $$\log_2(2) = 1$$
In the following text, the input will be referred to as x, and the output will be referred to as y.
Source code
Since the iterative part of our log algorithm assumes that $$1 \leq x < 2$$ we may need to “normalize” x. This normalization is performed by a simple succession of divides/multiplies by two. It is implemented by shifting the value bits. But how many places to shift? This function is used to calculate the number of leading zero bits.
int32_t c = __CLZ(src);
Next we need to normalize and convert to q30 format by shifting this value. This part of code literally serves to shift the input value to the left by as many fields so that the most significant bit is 0, and the first next bit is 1. That number is the value of the variable c.
/* Normalize and convert to q30 format */
x = src;
if ((c-1) < 0)
{
x = x >> (1-c);
}
else
{
x = x << (c-1);
}
normalization = c;
The number of divisions/multiplications determines the characteristic of the logarithm result. This characteristic will be added to the mantissa later. More details about the normalization variable and its purpose will be discussed further.
We now turn to the iterative part of the logarithm calculation, which is the core of the algorithm. It is necessary to explain the mathematical background that underlies the implementation of this section of the code.
Probably the first question that likely arises is why the value of x is squared in each iteration. Let’s try to explain. As mentioned earlier:
$$y = \log_2(x) \Rightarrow x = 2^y$$
The binary representation of the y value can be shown as follows:
$$y = y_1 \cdot 2^{-1} + y_2 \cdot 2^{-2} + y_3 \cdot 2^{-3} + y_4 \cdot 2^{-4} + …$$
It will be useful to express this in nested parentheses:
$$y = 2^{-1} ( y_1 + 2^{-1} ( y_2 + 2^{-1} ( y_3 + 2^{-1} ( y_4 + … ) ) ) )$$
The following notation can now be used:
$$x = 2^{2^{-1} ( y_1 + 2^{-1} ( y_2 + 2^{-1} ( y_3 + 2^{-1} ( y_4 + … ) ) ) )}$$
Finally, we reach the part where the value of x is squared: $$x^2 = 2^{2y}$$
$$x^2 = 2^{2[ 2^{-1} ( y_1 + 2^{-1} ( y_2 + 2^{-1} ( y_3 + 2^{-1} ( y_4 + … ) ) ) ) ]}$$
$$x^2 = 2^{y_1 + 2^{-1} ( y_2 + 2^{-1} ( y_3 + 2^{-1} ( y_4 + … ) ) )}$$
$$x^2 = 2^{y_1} \cdot 2^{2^{-1} ( y_2 + 2^{-1} ( y_3 + 2^{-1} ( y_4 + … ) ) ) ]}$$
Considering the following dependencies, we can easily determine whether the y_1 value is 1 or 0: $$y = \log_2(x) \Rightarrow x = 2^y$$
$$\text{if } x^2 \geq 2 \Rightarrow 2y \geq 1 \Rightarrow y_1 = 1$$ $$\text{if } x^2 < 2 \Rightarrow 2y < 1 \Rightarrow y_1 = 0$$
Extracting data from the remaining bits follows the same approach within a for loop.
for( i = 0; i < LOG_Q31_ACCURACY ; i++ )
{
x = ((int64_t)x*x) >> (LOG_Q31_ACCURACY - 1);
if (x >= LOQ_Q31_THRESHOLD)
{
y += inc ;
x = x >> 1;
}
inc = inc >> 1;
}
Next, we come across this line of code:
tmp = (int32_t)y - (normalization << (LOG_Q31_ACCURACY - LOG_Q31_INTEGER_PART));
The LOG_Q31_ACCURACY - LOG_Q31_INTEGER_PART shifting is needed to have output values in q5.26 format as stated in the documentation. But why the normalization value is subtracted from the previously calculated value of y?
Remember that at the beginning, the value of x was divided/shifted to the right for normalization. Now, it is necessary to subtract from the value of y because the following applies:
$$\frac{x}{2^n} = 2^{y - n}$$
At the end, it is necessary to multiply the value with a previously defined constant in order to calculate the natural logarithm:
$$\ln(x) = \log_2(x) \cdot \frac{1}{\log_2(e)}$$
i.e., in the code:
y = ((int64_t)tmp * LOG_Q31_INVLOG2EXP) >> 31;
Testing
The library is released in source form. We downloaded the files and compiled them, but without using -Ofast flag.
We used STM32F373 MCU and IAR Embedded Workbench for Arm, version 9.50.
We tested the execution time of the arm_vlog_q31() function for 32 values. It was about 500us.
DBG_PIN_LOW
__disable_interrupt();
DBG_PIN_HIGH
arm_vlog_q31( input_values_q31, output_values_q31, 32 );
DBG_PIN_LOW
__enable_interrupt();
DBG_PIN_HIGH
Below are the input and output values used during testing, along with the natural logarithms calculated for those input values, and the differences between them and the output values of the arm_vlog_q31() function.
The Q-format Converter & Calculator is a useful tool for performing arithmetic operations on fixed-point values.