BKM algorithm

In order to solve the equation

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

the BKM algorithm takes advantage of a basic property of logarithms

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

Using Pi notation, this identity generalizes to

${\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 ${\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 ${\displaystyle a_{k}=1+2^{m}}$ is a better option than a more simple choice of ${\displaystyle a_{k}=2^{m}}$. Since we want to start with large changes and get more accurate as ${\displaystyle k}$ increases, we can more specifically use ${\displaystyle a_{k}=1+2^{-k}}$, allowing the product to approach any value between 1 and ~4.768, depending on which subset of ${\displaystyle a_{k}}$ we use in the final product. At this point, the above equation looks like this:

${\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 ${\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 ${\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 ${\displaystyle x}$ while the other approaches the output value ${\displaystyle \ln(x)=y}$:

$${\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 ${\displaystyle x_{k}}$ is strictly increasing, it can be shown by induction and convergence that

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

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

${\displaystyle A_{k}=\ln(1+2^{-k})}$

Then the output is computed iteratively by the definition

$${\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

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

for any ${\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 ${\displaystyle y}$ is the known value and ${\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

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

to

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

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

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

 double log_e (const double Argument, const int Bits = 53) // 1 <= Argument <= 4.768462058
 {
    double  x = 1.0, y = 0.0, s = 1.0;

    for (int k = 0; k < Bits; k++) {
      double const  z = x + x*s;
      if (z <= Argument) {
          x  = z;
          y += A_e[k];
      }
      s *= 0.5;
    }
    return y;
 }

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

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

 double log_2 (const double Argument, const int Bits = 53) // 1 <= Argument <= 4.768462058
 {
    double  x = 1.0, y = 0.0, s = 1.0;

    for (int k = 0; k < Bits; k++) {
      double const  z = x + x*s;
      if (z <= Argument) {
          x  = z;
          y += A_2[k];
      }
      s *= 0.5;
    }
    return y;
 }

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.

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

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

 double exp (const double Argument, const int Bits = 54)	// 0 <= Argument <= 1.5620238332
 {
   double  x = 1.0, y = 0.0, s = 1.0;

   for (int k = 0; k < Bits; k++) {
      double const  z = y + A_e[k];
      if (z <= Argument) {
         y = z;
         x = x + x*s;
      }
      s *= 0.5;
   }
   return x;
 }

• Bajard, Jean-Claude; Kla, Sylvanus; Muller, Jean-Michel (August 1994). "BKM: A new hardware algorithm for complex elementary functions" (PDF). IEEE Transactions on Computers. 43 (8): 955–963. doi:10.1109/12.295857. ISSN 0018-9340. Zbl 1073.68501. Archived (PDF) from the original on 2018-11-03. Retrieved 2017-12-21.
• Skaf, Ali; Muller, Jean-Michel; Guyot, Alain (20–22 September 1994). On-Line Hardware Implementation for Complex Exponential and Logarithm. ESSCIRC '94: Twientieth European Solid-State Circuits Conference. Ulm, Germany. CiteSeerX 10.1.1.47.7521. ISBN 2-86332-160-9. Retrieved 2021-08-23. (5 pages)
• Bajard, Jean-Claude; Imbert, Laurent (1999-11-02). Luk, Franklin T. (ed.). "Evaluation of Complex Elementary Functions: A New Version of BKM" (PDF). SPIE Proceedings, Advanced Signal Processing Algorithms, Architectures, and Implementations IX. Advanced Signal Processing Algorithms, Architectures, and Implementations IX. Society of Photo-Optical Instrumentation Engineers (SPIE). 3807: 2–9. Bibcode:1999SPIE.3807....2B. doi:10.1117/12.367631. S2CID 121001818. Archived (PDF) from the original on 2020-06-09. Retrieved 2020-06-09. (8 pages)
• Imbert, Laurent; Muller, Jean-Michel; Rico, Fabien (2006-05-24) [2000-06-01, September 1999]. "Radix-10 BKM Algorithm for Computing Transcendentals on Pocket Computers". Journal of VLSI Signal Processing (Research report). Kluwer Academic Publishers / Institut national de recherche en informatique et en automatique (INRIA). 25 (2): 179–186. doi:10.1023/A:1008127208220. ISSN 0922-5773. S2CID 392036. RR-3754. INRIA-00072908. Theme 2 ISSN 0249-6399. Archived from the original on 2018-07-11. Retrieved 2018-07-11. (1+15 pages)
• Didier, Laurent-Stéphane; Rico, Fabien (2002-01-21). "High radix BKM algorithm with Selection by Rounding" (PDF). S2CID 17750192. lip6.2002.009. hal-02545612. Archived (PDF) from the original on 2021-08-23. Retrieved 2021-08-23. (1+11 pages)
• Didier, Laurent-Stéphane; Rico, Fabien (2004-12-01). "High Radix BKM Algorithm". Numerical Algorithms. SCAN'2002 International Conference. Springer Science+Business Media, LLC. 37 (1–4 [4]): 113–125. Bibcode:2004NuAlg..37..113D. doi:10.1023/B:NUMA.0000049459.69390.ff. eISSN 1572-9265. ISSN 1017-1398. S2CID 2761452.
• Muller, Jean-Michel (2006). Elementary Functions: Algorithms and Implementation (2 ed.). Boston, MA, USA: Birkhäuser. ISBN 978-0-8176-4372-0. LCCN 2005048094.
• 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.

• Jorke, Günter; Lampe, Bernhard; Wengel, Norbert (1989). Arithmetische Algorithmen der Mikrorechentechnik (in German) (1 ed.). Berlin, Germany: VEB Verlag Technik. pp. 280–282. ISBN 3-34100515-3. ISBN 978-3-34100515-6. EAN 9783341005156. MPN 5539165. License 201.370/4/89. Retrieved 2015-12-01.
• Meggitt, John E. (1961-08-29). "Pseudo Division and Pseudo Multiplication Processes". IBM Journal of Research and Development. Riverton, New Jersey, USA: IBM Corporation (published April 1962). 6 (2): 210–226, 287. doi:10.1147/rd.62.0210. Retrieved 2015-12-01.
• Chi Chen, Tien (July 1972). "Automatic computation of exponentials, logarithms, ratios and square roots". IBM Journal of Research and Development. San Jose, California, USA; Riverton, New Jersey, USA: IBM San Jose Research Laboratory; IBM Corporation. 16 (4): 380–388. doi:10.1147/rd.164.0380. Retrieved 2015-12-01.
• Revol, Nathalie; Yakoubsohn, Jean-Claude (2000-05-01). "Accelerated Shift-and-Add algorithms" (PDF). Reliable Computing. Boston, USA: Laboratoire d'Analyse Numérique et d'Optimisation (ANO) de l'Université des Sciences et Technologies de Lille; Kluwer Academic Publishers. 6 (2): 193–205. doi:10.1023/A:1009921407000. eISSN 1573-1340. ISSN 1385-3139. OCLC 67306353. S2CID 10716391. Archived (PDF) from the original on 2021-08-23. Retrieved 2021-08-23. (14 pages)