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

Per induction, the iteration

x_{n+1}=x_{n}\cdot (1+d_{n}\cdot 2^{-n})

with $x_{0}=1$ and $d_{n}\in \{0,1\}$ is identical to

x_{n+1}=\prod _{i=0}^{n}(1+2^{-i})^{d_{i}}

If all $d_{n}=0$ , then all $x_{n}=1$ . If all $d_{n}=1$ then $x_{\infty }\approx 4{,}768$ .^{[nb 1]} Choosing appropriate $d_{n}$ any real number $x$ in the range $1\leqq x\lessapprox 4{,}768$ can be shown as limit of a sequence.

Further

y_{n+1}=y_{n}+d_{n}\cdot \ln(1+2^{-n})

with $y_{0}=0$ or equivalent

y_{n+1}=\sum _{i=0}^{n}d_{i}\cdot \ln(1+2^{-i})=\ln \left(\prod _{i=0}^{n}(1+2^{-i})^{d_{i}}\right)

.

For numerical calculations $A_{n}=\ln(1+2^{-n})$ is realized through a pre-calculated table.

It can be derived that $y_{n}=\ln(x_{n})$ for all $n$ . Similarly, for the logarithm: $0\leqq y=\ln(x)\lessapprox 1{,}562$ .

Logarithm function

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

Exponential function

To calculate the exponential function (E-mode), the algorithm in each iteration tests if $y_{n}+\ln(1+2^{-n})\leq y$ . If so, it calculates $x_{n+1}$ and $y_{n+1}$ . After $N$ iterations the value of the function is known with an error of $\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;
 }

Tables for examples

 static const double A_e [] = // A_e[k] = ln (1 + 0.5^k)
 {
   0.693147180559945297099404706000, 0.405465108108164392935428259000, 0.223143551314209769962616701000,
   0.117783035656383447138088388000, 0.060624621816434840186291518000, 0.030771658666753686222134530000,
   0.015504186535965253358272343000, 0.007782140442054949100825041000, 0.003898640415657322636221046000,
   0.001951220131261749216850870000, 0.000976085973055458892686544000, 0.000488162079501351186957460000,
   0.000244110827527362687853374000, 0.000122062862525677363338881000, 0.000061033293680638525913091000,
   0.000030517112473186377476993000, 0.000015258672648362398138404000, 0.000007629365427567572417821000,
   0.000003814689989685889381171000, 0.000001907346813825409407938000, 0.000000953673861659188260005000,
   0.000000476837044516323000000000, 0.000000238418550679858000000000, 0.000000119209282445354000000000,
   0.000000059604642999033900000000, 0.000000029802321943606100000000, 0.000000014901161082825400000000,
   0.000000007450580569168250000000, 0.000000003725290291523020000000, 0.000000001862645147496230000000,
   0.000000000931322574181798000000, 0.000000000465661287199319000000, 0.000000000232830643626765000000,
   0.000000000116415321820159000000, 0.000000000058207660911773300000, 0.000000000029103830456310200000,
   0.000000000014551915228261000000, 0.000000000007275957614156960000, 0.000000000003637978807085100000,
   0.000000000001818989403544200000, 0.000000000000909494701772515000, 0.000000000000454747350886361000,
   0.000000000000227373675443206000, 0.000000000000113686837721610000, 0.000000000000056843418860806400,
   0.000000000000028421709430403600, 0.000000000000014210854715201900, 0.000000000000007105427357600980,
   0.000000000000003552713678800490, 0.000000000000001776356839400250, 0.000000000000000888178419700125,
   0.000000000000000444089209850063, 0.000000000000000222044604925031, 0.000000000000000111022302462516,
   0.000000000000000055511151231258, 0.000000000000000027755575615629, 0.000000000000000013877787807815,
   0.000000000000000006938893903907, 0.000000000000000003469446951954, 0.000000000000000001734723475977,
   0.000000000000000000867361737988, 0.000000000000000000433680868994, 0.000000000000000000216840434497,
   0.000000000000000000108420217249, 0.000000000000000000054210108624, 0.000000000000000000027105054312,
 };

 static const double A_2 [] = // A_2[k] = log_2 (1 + 0.5^k)
 {
   1.0000000000000000000000000000000000000000000000000000000000000000000000000000,
   0.5849625007211561814537389439478165087598144076924810604557526545410982276485,
   0.3219280948873623478703194294893901758648313930245806120547563958159347765589,
   0.1699250014423123629074778878956330175196288153849621209115053090821964552970,
   0.0874628412503394082540660108104043540112672823448206881266090643866965081686,
   0.0443941193584534376531019906736094674630459333742491317685543002674288465967,
   0.0223678130284545082671320837460849094932677948156179815932199216587899627785,
   0.0112272554232541203378805844158839407281095943600297940811823651462712311786,
   0.0056245491938781069198591026740666017211096815383520359072957784732489771013,
   0.0028150156070540381547362547502839489729507927389771959487826944878598909400,
   0.0014081943928083889066101665016890524233311715793462235597709051792834906001,
   0.0007042690112466432585379340422201964456668872087249334581924550139514213168,
   0.0003521774803010272377989609925281744988670304302127133979341729842842377649,
   0.0001760994864425060348637509459678580940163670081839283659942864068257522373,
   0.0000880524301221769086378699983597183301490534085738474534831071719854721939,
   0.0000440268868273167176441087067175806394819146645511899503059774914593663365,
   0.0000220136113603404964890728830697555571275493801909791504158295359319433723,
   0.0000110068476674814423006223021573490183469930819844945565597452748333526464,
   0.0000055034343306486037230640321058826431606183125807276574241540303833251704,
   0.0000027517197895612831123023958331509538486493412831626219340570294203116559,
   0.0000013758605508411382010566802834037147561973553922354232704569052932922954,
   0.0000006879304394358496786728937442939160483304056131990916985043387874690617,
   0.0000003439652607217645360118314743718005315334062644619363447395987584138324,
   0.0000001719826406118446361936972479533123619972434705828085978955697643547921,
   0.0000000859913228686632156462565208266682841603921494181830811515318381744650,
   0.0000000429956620750168703982940244684787907148132725669106053076409624949917,
   0.0000000214978311976797556164155504126645192380395989504741781512309853438587,
   0.0000000107489156388827085092095702361647949603617203979413516082280717515504,
   0.0000000053744578294520620044408178949217773318785601260677517784797554422804,
   0.0000000026872289172287079490026152352638891824761667284401180026908031182361,
   0.0000000013436144592400232123622589569799954658536700992739887706412976115422,
   0.0000000006718072297764289157920422846078078155859484240808550018085324187007,
   0.0000000003359036149273187853169587152657145221968468364663464125722491530858,
   0.0000000001679518074734354745159899223037458278711244127245990591908996412262,
   0.0000000000839759037391617577226571237484864917411614198675604731728132152582,
   0.0000000000419879518701918839775296677020135040214077417929807824842667285938,
   0.0000000000209939759352486932678195559552767641474249812845414125580747434389,
   0.0000000000104969879676625344536740142096218372850561859495065136990936290929,
   0.0000000000052484939838408141817781356260462777942148580518406975851213868092,
   0.0000000000026242469919227938296243586262369156865545638305682553644113887909,
   0.0000000000013121234959619935994960031017850191710121890821178731821983105443,
   0.0000000000006560617479811459709189576337295395590603644549624717910616347038,
   0.0000000000003280308739906102782522178545328259781415615142931952662153623493,
   0.0000000000001640154369953144623242936888032768768777422997704541618141646683,
   0.0000000000000820077184976595619616930350508356401599552034612281802599177300,
   0.0000000000000410038592488303636807330652208397742314215159774270270147020117,
   0.0000000000000205019296244153275153381695384157073687186580546938331088730952,
   0.0000000000000102509648122077001764119940017243502120046885379813510430378661,
   0.0000000000000051254824061038591928917243090559919209628584150482483994782302,
   0.0000000000000025627412030519318726172939815845367496027046030028595094737777,
   0.0000000000000012813706015259665053515049475574143952543145124550608158430592,
   0.0000000000000006406853007629833949364669629701200556369782295210193569318434,
   0.0000000000000003203426503814917330334121037829290364330169106716787999052925,
   0.0000000000000001601713251907458754080007074659337446341494733882570243497196,
   0.0000000000000000800856625953729399268240176265844257044861248416330071223615,
   0.0000000000000000400428312976864705191179247866966320469710511619971334577509,
   0.0000000000000000200214156488432353984854413866994246781519154793320684126179,
   0.0000000000000000100107078244216177339743404416874899847406043033792202127070,
   0.0000000000000000050053539122108088756700751579281894640362199287591340285355,
   0.0000000000000000025026769561054044400057638132352058574658089256646014899499,
   0.0000000000000000012513384780527022205455634651853807110362316427807660551208,
   0.0000000000000000006256692390263511104084521222346348012116229213309001913762,
   0.0000000000000000003128346195131755552381436585278035120438976487697544916191,
   0.0000000000000000001564173097565877776275512286165232838833090480508502328437,
   0.0000000000000000000782086548782938888158954641464170239072244145219054734086,
   0.0000000000000000000391043274391469444084776945327473574450334092075712154016,
   0.0000000000000000000195521637195734722043713378812583900953755962557525252782,
   0.0000000000000000000097760818597867361022187915943503728909029699365320287407,
   0.0000000000000000000048880409298933680511176764606054809062553340323879609794,
   0.0000000000000000000024440204649466840255609083961603140683286362962192177597,
   0.0000000000000000000012220102324733420127809717395445504379645613448652614939,
   0.0000000000000000000006110051162366710063906152551383735699323415812152114058,
   0.0000000000000000000003055025581183355031953399739107113727036860315024588989,
   0.0000000000000000000001527512790591677515976780735407368332862218276873443537,
   0.0000000000000000000000763756395295838757988410584167137033767056170417508383,
   0.0000000000000000000000381878197647919378994210346199431733717514843471513618,
   0.0000000000000000000000190939098823959689497106436628681671067254111334889005,
   0.0000000000000000000000095469549411979844748553534196582286585751228071408728,
   0.0000000000000000000000047734774705989922374276846068851506055906657137209047,
   0.0000000000000000000000023867387352994961187138442777065843718711089344045782,
   0.0000000000000000000000011933693676497480593569226324192944532044984865894525,
   0.0000000000000000000000005966846838248740296784614396011477934194852481410926,
   0.0000000000000000000000002983423419124370148392307506484490384140516252814304,
   0.0000000000000000000000001491711709562185074196153830361933046331030629430117,
   0.0000000000000000000000000745855854781092537098076934460888486730708440475045,
   0.0000000000000000000000000372927927390546268549038472050424734256652501673274,
   0.0000000000000000000000000186463963695273134274519237230207489851150821191330,
   0.0000000000000000000000000093231981847636567137259618916352525606281553180093,
   0.0000000000000000000000000046615990923818283568629809533488457973317312233323,
   0.0000000000000000000000000023307995461909141784314904785572277779202790023236,
   0.0000000000000000000000000011653997730954570892157452397493151087737428485431,
   0.0000000000000000000000000005826998865477285446078726199923328593402722606924,
   0.0000000000000000000000000002913499432738642723039363100255852559084863397344,
   0.0000000000000000000000000001456749716369321361519681550201473345138307215067,
   0.0000000000000000000000000000728374858184660680759840775119123438968122488047,
   0.0000000000000000000000000000364187429092330340379920387564158411083803465567,
   0.0000000000000000000000000000182093714546165170189960193783228378441837282509,
   0.0000000000000000000000000000091046857273082585094980096891901482445902524441,
   0.0000000000000000000000000000045523428636541292547490048446022564529197237262,
   0.0000000000000000000000000000022761714318270646273745024223029238091160103901,
 };

Notes

^ For more digits see (sequence A081845 in the OEIS).

References

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) [1][2]
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. 3807. Society of Photo-Optical Instrumentation Engineers (SPIE): 2–9. Bibcode:1999SPIE.3807....2B. doi:10.1117/12.367631. Archived (PDF) from the original on 2020-06-09. Retrieved 2020-06-09. (8 pages) [3]
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). 25 (2). Kluwer Academic Publishers / Institut national de recherche en informatique et en automatique (INRIA): 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) [4][5][6]
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. [7] (1+11 pages)
Didier, Laurent-Stéphane; Rico, Fabien (2004-12-01). "High Radix BKM Algorithm". Numerical Algorithms. SCAN'2002 International Conference. 37 (1–4 [4]). Springer Science+Business Media, LLC: 113–125. doi:10.1023/B:NUMA.0000049459.69390.ff. eISSN 1572-9265. ISSN 1017-1398. S2CID 2761452. [8]
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.