BKM algorithm

Last updated April 26, 2023

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

the BKM algorithm takes advantage of a basic property of logarithms

\ln(ab)=\ln(a)+\ln(b)

Using Pi notation, this identity generalizes to

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

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

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

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

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

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

Then the output is computed iteratively by the definition

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\to \infty }y_{k}=y

for any $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$ is the known value and $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

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

to

\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}\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):

doublelog_e(constdoubleArgument,constintBits=53)// 1 <= Argument <= 4.768462058{doublex=1.0,y=0.0,s=1.0;for(intk=0;k<Bits;k++){doubleconstz=x+x*s;if(z<=Argument){x=z;y+=A_e[k];}s*=0.5;}returny;}

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

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

doublelog_2(constdoubleArgument,constintBits=53)// 1 <= Argument <= 4.768462058{doublex=1.0,y=0.0,s=1.0;for(intk=0;k<Bits;k++){doubleconstz=x+x*s;if(z<=Argument){x=z;y+=A_2[k];}s*=0.5;}returny;}

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})\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):

doubleexp(constdoubleArgument,constintBits=54)// 0 <= Argument <= 1.5620238332{doublex=1.0,y=0.0,s=1.0;for(intk=0;k<Bits;k++){doubleconstz=y+A_e[k];if(z<=Argument){y=z;x=x+x*s;}s*=0.5;}returnx;}

Tables for examples

static const double A_e [] = ...

staticconstdoubleA_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 [] = ...

staticconstdoubleA_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

Related Research Articles

<span class="mw-page-title-main">Exponential function</span> Mathematical function, denoted exp(x) or e^x

The exponential function is a mathematical function denoted by $or . Unless otherwise specified, the term generally refers to the positive-valued function of a real variable, although it can be extended to the complex numbers or generalized to other mathematical objects like matrices or Lie algebras. The exponential function originated from the notion of exponentiation, but modern definitions allow it to be rigorously extended to all real arguments, including irrational numbers. Its ubiquitous occurrence in pure and applied mathematics led mathematician Walter Rudin to opine that the exponential function is "the most important function in mathematics".$

In information theory, the entropy of a random variable is the average level of "information", "surprise", or "uncertainty" inherent to the variable's possible outcomes. Given a discrete random variable $, which takes values in the alphabet and is distributed according to :$

In mathematics, the logarithm is the inverse function to exponentiation. That means the logarithm of a number $x$ to the base $b$ is the exponent to which $b$ must be raised, to produce $x$ . For example, since $1000 = 10 3$ , the logarithm base 10 of $1000$ is $3$ , or $log 10 (1000) = 3$ . The logarithm of $x$ to base $b$ is denoted as $log b (x)$ , or without parentheses, $log b x$ , or even without the explicit base, $log x$ , when no confusion is possible, or when the base does not matter such as in big O notation.

The natural logarithm of a number is its logarithm to the base of the mathematical constant $e$ , which is an irrational and transcendental number approximately equal to $2.718 281828459$ . The natural logarithm of $x$ is generally written as $ln x$ , $log e x$ , or sometimes, if the base $e$ is implicit, simply $log x$ . Parentheses are sometimes added for clarity, giving $ln(x)$ , $log e (x)$ , or $log(x)$ . This is done particularly when the argument to the logarithm is not a single symbol, so as to prevent ambiguity.

In computability theory, a primitive recursive function is, roughly speaking, a function that can be computed by a computer program whose loops are all "for" loops. Primitive recursive functions form a strict subset of those general recursive functions that are also total functions.

In mathematics, the Lambert $W$ function, also called the omega function or product logarithm, is a multivalued function, namely the branches of the converse relation of the function $f (w) = we w$ , where $w$ is any complex number and $e w$ is the exponential function.

<span class="mw-page-title-main">Exponentiation</span> Mathematical operation

In mathematics, exponentiation is an operation involving two numbers, the base and the exponent or power. Exponentiation is written as $b n$ , where $b$ is the base and $n$ is the power; this is pronounced as " $b$ (raised) to the $n$ ". When $n$ is a positive integer, exponentiation corresponds to repeated multiplication of the base: that is, $b n$ is the product of multiplying $n$ bases:

In mathematics, the binary logarithm is the power to which the number $2$ must be raised to obtain the value $n$ . That is, for any real number $x$ ,

In mathematics, the inverse trigonometric functions are the inverse functions of the trigonometric functions. Specifically, they are the inverses of the sine, cosine, tangent, cotangent, secant, and cosecant functions, and are used to obtain an angle from any of the angle's trigonometric ratios. Inverse trigonometric functions are widely used in engineering, navigation, physics, and geometry.

<span class="mw-page-title-main">Tetration</span> Repeated exponentiation

In mathematics, tetration is an operation based on iterated, or repeated, exponentiation. There is no standard notation for tetration, though $and the left-exponent x b are common.$

CORDIC, also known as Volder's algorithm, or: Digit-by-digit methodCircular CORDIC, Linear CORDIC, Hyperbolic CORDIC, and Generalized Hyperbolic CORDIC, is a simple and efficient algorithm to calculate trigonometric functions, hyperbolic functions, square roots, multiplications, divisions, and exponentials and logarithms with arbitrary base, typically converging with one digit per iteration. CORDIC is therefore also an example of digit-by-digit algorithms. CORDIC and closely related methods known as pseudo-multiplication and pseudo-division or factor combining are commonly used when no hardware multiplier is available, as the only operations it requires are additions, subtractions, bitshift and lookup tables. As such, they all belong to the class of shift-and-add algorithms. In computer science, CORDIC is often used to implement floating-point arithmetic when the target platform lacks hardware multiply for cost or space reasons.

<span class="mw-page-title-main">Iterated logarithm</span> Inverse function to a tower of powers

In computer science, the iterated logarithm of $, written log *, is the number of times the logarithm function must be iteratively applied before the result is less than or equal to . The simplest formal definition is the result of this recurrence relation:$

Methods of computing square roots are numerical analysis algorithms for approximating the principal, or non-negative, square root of a real number. Arithmetically, it means given $, a procedure for finding a number which when multiplied by itself, yields; algebraically, it means a procedure for finding the non-negative root of the equation; geometrically, it means given two line segments, a procedure for constructing their geometric mean.$

The Anderson–Darling test is a statistical test of whether a given sample of data is drawn from a given probability distribution. In its basic form, the test assumes that there are no parameters to be estimated in the distribution being tested, in which case the test and its set of critical values is distribution-free. However, the test is most often used in contexts where a family of distributions is being tested, in which case the parameters of that family need to be estimated and account must be taken of this in adjusting either the test-statistic or its critical values. When applied to testing whether a normal distribution adequately describes a set of data, it is one of the most powerful statistical tools for detecting most departures from normality. K-sample Anderson–Darling tests are available for testing whether several collections of observations can be modelled as coming from a single population, where the distribution function does not have to be specified.

In mathematics, a complex logarithm is a generalization of the natural logarithm to nonzero complex numbers. The term refers to one of the following, which are strongly related:

This is a summary of differentiation rules, that is, rules for computing the derivative of a function in calculus.

In mathematics, the super-logarithm is one of the two inverse functions of tetration. Just as exponentiation has two inverse functions, roots and logarithms, tetration has two inverse functions, super-roots and super-logarithms. There are several ways of interpreting super-logarithms:

<span class="mw-page-title-main">Fast inverse square root</span> Root-finding algorithm

Fast inverse square root, sometimes referred to as Fast InvSqrt or by the hexadecimal constant 0x5F3759DF, is an algorithm that estimates $, the reciprocal of the square root of a 32-bit floating-point number in IEEE 754 floating-point format. This operation is used in digital signal processing to normalize a vector, such as scaling it to length 1. For example, computer graphics programs use inverse square roots to compute angles of incidence and reflection for lighting and shading. Predated by similar video game algorithms, this one is best known for its implementation in 1999 in Quake III Arena, a first-person shooter video game heavily based on 3D graphics. The algorithm only started appearing on public forums between 2002 and 2003. Computation of square roots usually depends upon many division operations, which for floating point numbers are computationally expensive. The fast inverse square generates a good approximation with only one division step.$

In cryptography, Very Smooth Hash (VSH) is a provably secure cryptographic hash function invented in 2005 by Scott Contini, Arjen Lenstra and Ron Steinfeld. Provably secure means that finding collisions is as difficult as some known hard mathematical problem. Unlike other provably secure collision-resistant hashes, VSH is efficient and usable in practice. Asymptotically, it only requires a single multiplication per log(n) message-bits and uses RSA-type arithmetic. Therefore, VSH can be useful in embedded environments where code space is limited.

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)
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. 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.