BKM algorithm

The BKM algorithm is a shift-and-add algorithm for computing elementary functions, first published in 1994 by Jean-Claude Bajard, Sylvanus Kla, and Jean-Michel Muller. BKM is based on computing complex logarithms (L-mode) and exponentials (E-mode) using a method similar to the algorithm Henry Briggs used to compute logarithms. By using a precomputed table of logarithms of negative powers of two, the BKM algorithm computes elementary functions using only integer add, shift, and compare operations.

BKM is similar to CORDIC, but uses a table of logarithms rather than a table of arctangents. On each iteration, a choice of coefficient is made from a set of nine complex numbers, 1, 0, −1, i, −i, 1+i, 1−i, −1+i, −1−i, rather than only −1 or +1 as used by CORDIC. BKM provides a simpler method of computing some elementary functions, and unlike CORDIC, BKM needs no result scaling factor. The convergence rate of BKM is approximately one bit per iteration, like CORDIC, but BKM requires more precomputed table elements for the same precision because the table stores logarithms of complex operands.

As with other algorithms in the shift-and-add class, BKM is particularly well-suited to hardware implementation. The relative performance of software BKM implementation in comparison to other methods such as polynomial or rational approximations will depend on the availability of fast multi-bit shifts (i.e. a barrel shifter) or hardware floating point arithmetic.

Overview

In order to solve the equation

ln ⁡ ( x ) = y {\displaystyle \ln(x)=y}

the BKM algorithm takes advantage of a basic property of logarithms

ln ⁡ ( a b ) = ln ⁡ ( a ) + ln ⁡ ( b ) {\displaystyle \ln(ab)=\ln(a)+\ln(b)}

Using Pi notation, this identity generalizes to

ln ⁡ ( ∏ k = 0 n a k ) = ∑ k = 0 n ln ⁡ ( a k ) {\displaystyle \ln \left(\prod _{k=0}^{n}a_{k}\right)=\sum _{k=0}^{n}\ln(a_{k})}

Because any number can be represented by a product, this allows us to choose any set of values a k {\displaystyle a_{k}} which multiply to give the value we started with. In computer systems, it's much faster to multiply and divide by multiples of 2, but because not every number is a multiple of 2, using a k = 1 + 2 m {\displaystyle a_{k}=1+2^{m}} is a better option than a more simple choice of a k = 2 m {\displaystyle a_{k}=2^{m}}. Since we want to start with large changes and get more accurate as k {\displaystyle k} increases, we can more specifically use a k = 1 + 2 − k {\displaystyle a_{k}=1+2^{-k}}, allowing the product to approach any value between 1 and ~4.768, depending on which subset of a k {\displaystyle a_{k}} we use in the final product. At this point, the above equation looks like this:

ln ⁡ ( ∏ k ∈ K ⊂ Z 0 + 1 + 2 − k ) = ∑ k ∈ K ⊂ Z 0 + ln ⁡ ( 1 + 2 − k ) {\displaystyle \ln \left(\prod _{k\in K\subset \mathbb {Z} _{0}^{+}}1+2^{-k}\right)=\sum _{k\in K\subset \mathbb {Z} _{0}^{+}}\ln(1+2^{-k})}

This choice of a k {\displaystyle a_{k}} reduces the computational complexity of the product from repeated multiplication to simple addition and bit-shifting depending on the implementation. Finally, by storing the values ln ⁡ ( 1 + 2 − k ) {\displaystyle \ln(1+2^{-k})} in a table, calculating the solution is also a simple matter of addition. Iteratively, this gives us two separate sequences. One sequence approaches the input value x {\displaystyle x} while the other approaches the output value ln ⁡ ( x ) = y {\displaystyle \ln(x)=y}:

x k = { 1 if k = 0 x k − 1 ⋅ ( 1 + 2 − k ) if x k would be ≤ x x k − 1 otherwise {\displaystyle x_{k}={\begin{cases}1&{\text{if }}k=0\\x_{k-1}\cdot (1+2^{-k})&{\text{if }}x_{k}{\text{ would be}}\leq x\\x_{k-1}&{\text{otherwise}}\end{cases}}}

Given this recursive definition and because x k {\displaystyle x_{k}} is strictly increasing, it can be shown by induction and convergence that

lim k → ∞ x k = x {\displaystyle \lim _{k\to \infty }x_{k}=x}

for any 1 ≤ x ≲ 4.768 {\displaystyle 1\leq x\lesssim 4.768}. For calculating the output, we first create the reference table

A k = ln ⁡ ( 1 + 2 − k ) {\displaystyle A_{k}=\ln(1+2^{-k})}

Then the output is computed iteratively by the definition y k = { 0 if k = 0 y k − 1 + A k if x k would be ≤ x y k − 1 otherwise {\displaystyle y_{k}={\begin{cases}0&{\text{if }}k=0\\y_{k-1}+A_{k}&{\text{if }}x_{k}{\text{ would be}}\leq x\\y_{k-1}&{\text{otherwise}}\end{cases}}} The conditions in this iteration are the same as the conditions for the input. Similar to the input, this sequence is also strictly increasing, so it can be shown that

lim k → ∞ y k = y {\displaystyle \lim _{k\to \infty }y_{k}=y}

for any 0 ≤ y ≲ 1.562 {\displaystyle 0\leq y\lesssim 1.562}.

Because the algorithm above calculates both the input and output simultaneously, it's possible to modify it slightly so that y {\displaystyle y} is the known value and x {\displaystyle x} is the value we want to calculate, thereby calculating the exponential instead of the logarithm. Since x becomes an unknown in this case, the conditional changes from

… if x k would be ≤ x {\displaystyle \dots {\text{if }}x_{k}{\text{ would be}}\leq x}

to

… if y k would be ≤ y {\displaystyle \dots {\text{if }}y_{k}{\text{ would be}}\leq y}

Logarithm function

To calculate the logarithm function (L-mode), the algorithm in each iteration tests if x n ⋅ ( 1 + 2 − n ) ≤ x {\displaystyle x_{n}\cdot (1+2^{-n})\leq x}. If so, it calculates x n + 1 {\displaystyle x_{n+1}} and y n + 1 {\displaystyle y_{n+1}}. After N {\displaystyle N} iterations the value of the function is known with an error of Δ ln ⁡ ( x ) ≤ 2 − N {\displaystyle \Delta \ln(x)\leq 2^{-N}}.

Example program for natural logarithm in C++ (see A_e for table):

Logarithms for bases other than e can be calculated with similar effort.

Example program for binary logarithm in C++ (see A_2 for table):

The allowed argument range is the same for both examples (1 ≤ Argument ≤ 4.768462058...). In the case of the base-2 logarithm the exponent can be split off in advance (to get the integer part) so that the algorithm can be applied to the remainder (between 1 and 2). Since the argument is smaller than 2.384231..., the iteration of k can start with 1. Working in either base, the multiplication by s can be replaced with direct modification of the floating point exponent, subtracting 1 from it during each iteration. This results in the algorithm using only addition and no multiplication.

Exponential function

To calculate the exponential function (E-mode), the algorithm in each iteration tests if y n + ln ⁡ ( 1 + 2 − n ) ≤ y {\displaystyle y_{n}+\ln(1+2^{-n})\leq y}. If so, it calculates x n + 1 {\displaystyle x_{n+1}} and y n + 1 {\displaystyle y_{n+1}}. After N {\displaystyle N} iterations the value of the function is known with an error of Δ exp ⁡ ( x ) ≤ 2 − N {\displaystyle \Delta \exp(x)\leq 2^{-N}}.

Example program in C++ (see A_e for table):

Tables for examples

Notes

• Bajard, Jean-Claude; Kla, Sylvanus; Muller, Jean-Michel (August 1994). (PDF). IEEE Transactions on Computers. 43 (8): 955–963. doi:. ISSN . Zbl . (PDF) from the original on 2018-11-03.
• Skaf, Ali; Muller, Jean-Michel; Guyot, Alain (20–22 September 1994). . ESSCIRC '94: Twientieth European Solid-State Circuits Conference. Ulm, Germany. CiteSeerX . ISBN 2-86332-160-9. (5 pages)
• Bajard, Jean-Claude; Imbert, Laurent (1999-11-02). (PDF). In Luk, Franklin T. (ed.). Advanced Signal Processing Algorithms, Architectures, and Implementations IX. SPIE Proceedings. Vol. 3807. Society of Photo-Optical Instrumentation Engineers (SPIE). pp. 2–9. Bibcode:. doi:. S2CID . (PDF) from the original on 2020-06-09. (8 pages)
• Imbert, Laurent; Muller, Jean-Michel; Rico, Fabien (2006-05-24) [2000-06-01, September 1999]. . Journal of VLSI Signal Processing (Research report). 25 (2). Kluwer Academic Publishers / Institut national de recherche en informatique et en automatique (INRIA): 179–186. doi:. ISSN . S2CID . RR-3754. INRIA-00072908. Theme 2 ISSN . from the original on 2018-07-11. (1+15 pages)
• Didier, Laurent-Stéphane; Rico, Fabien (2002-01-21). (PDF). S2CID . lip6.2002.009. hal-02545612. (PDF) from the original on 2021-08-23. (1+11 pages)
• Didier, Laurent-Stéphane; Rico, Fabien (2004-12-01). . Numerical Algorithms. SCAN'2002 International Conference. 37 (1–4 [4]). Springer Science+Business Media, LLC: 113–125. Bibcode:. doi:. eISSN . ISSN . S2CID .
• Muller, Jean-Michel (2006). Elementary Functions: Algorithms and Implementation (2 ed.). Boston, MA, USA: Birkhäuser. ISBN 978-0-8176-4372-0. LCCN .
• Muller, Jean-Michel (2016-12-12). Elementary Functions: Algorithms and Implementation (3 ed.). Boston, MA, USA: Birkhäuser. ISBN 978-1-4899-7981-0. ISBN 1-4899-7981-6.

Further reading

• Jorke, Günter; Lampe, Bernhard; Wengel, Norbert (1989). (in German) (1 ed.). Berlin, Germany: VEB Verlag Technik. pp. 280–282. ISBN 3-34100515-3. ISBN 978-3-34100515-6. EAN . MPN 5539165. License 201.370/4/89.
• Meggitt, John E. (1961-08-29). . IBM Journal of Research and Development. 6 (2). Riverton, New Jersey, USA: IBM Corporation (published April 1962): 210–226, 287. doi:.
• Chi Chen, Tien (July 1972). . IBM Journal of Research and Development. 16 (4). San Jose, California, USA; Riverton, New Jersey, USA: IBM San Jose Research Laboratory; IBM Corporation: 380–388. doi:.
• Revol, Nathalie; Yakoubsohn, Jean-Claude (2000-05-01). (PDF). Reliable Computing. 6 (2). Boston, USA: Laboratoire d'Analyse Numérique et d'Optimisation (ANO) de l'Université des Sciences et Technologies de Lille; Kluwer Academic Publishers: 193–205. doi:. eISSN . ISSN . OCLC . S2CID . (PDF) from the original on 2021-08-23. (14 pages)