Subject: Re: WANTED: SOURCE CODE FOR EXP, LN, ARCTAN, ARCSIN, ARCCOS From: Mike van Scherrenburg Date: 1997/11/11 Message-ID: <3468A42B.3FE1@ix.netcom.com> Newsgroups: sci.physics [More Headers] Someone put you on the track of Newton-Raphson. Well well, that's good textbook stuff, but that's not the only way to do it. If you don't need too high an accuracy to the functions, you should use some variant of the CORDIC technique. Several (hardware) graphics processors, some flight controllers, and at least the FPU for the 80386 used this technique. I'm not up on if the latest FPUs still use this technique or not. In the typical usage, one presumes the representation of the result as a binary number. The technique has an inner loop that adds about 1 bit of accuracy to the result per iteration of a loop. The loop consists of a few additions and (binary) shifts. There is also typically one multiplication or division by a constant after the iterations to scale the result properly. The CORDIC technique was designed to simplify the inner loop. There's nothing inherently "binary" in the technique. Perhaps you could re-do the math presuming a factorial representation, in which case, if the iteration converges, it ought to give a correct answer using addition and (factorial) shifts. This might help a lot, as I don't think doing binary shifts in factorial number bases is any faster than division by an arbitrary constant. Without "factorial" shifts, CORDIC would supply no advantage, as Newton-Raphson techniques converge more quickly. The reference: Graphics Gems, edited by Andrew S. Glassner, Academic Press, (C)1990. This book (and its four successors) is a collection of short articles and C code. The article of interest is: Fixed-Point Trigonometry with CORDIC Iterations, pp494-497. C code on pages 773-774. The author is Ken Turkowski, than at Apple in Cupertino, so you might be able to chase him down. Hope this helps. ------------------------------------ Subject: Re: Exponential function architecture From: Ray Andraka Date: 1997/04/20 Message-ID: <335AC759.73A0@ids.net> Newsgroups: comp.arch.arithmetic,comp.lsi,comp.dsp,comp.arch.fpga Mourad Aberbour wrote: > > Hi there, > Did anybody had any experience with designing an architecture for > the computation of the exponential function? > I would appreciate any kind of help. > > regards, > -- You can compute the exponential using the CORDIC algorithm, or more efficiently with a CORDIC-like algorithm that works on an incremental expression of the exponential. The CORDIC algorithms are iterative shift-add algorithms useful for computation of trig, inverse trig, hyperbolic and inverse hyperbolic, square roots and other functions. Generally speaking, the results improve by one bit for each iteration (the exponent is an exception). I have extensive experience using the CORDIC algorithm, including design of a log/exponential. I've got one paper on my website that discusses a bit serial construction of a CORDIC vector magnitude. The structure for the exponential is similar. The first method, described by Walther in his paper on elementary functions (ca 1971) derives the exponential from hyperbolic CORDIC functions. The more efficient method uses a small look-up (one entry per iteration) to compute the Log: The incremental expression for an exponential, bx+dx = bx * bdx, can be restated as bx+dx =bx + bx * 2-i if bdx = 1 + 2-i. For this to be useful, the control element must be incremented by logb(1 + 2-i) for each iteration that the decision function is true. Since the value of the log element is not used by the exponent element, the same results are obtained if the control element is initialized with x and decremented to zero. The resulting cordic equations are: Fi+1 = Fi + di*Fi*2-i and Xi+1 = Xi - di*Logb(1 + 2-i). where di = 1 if Xi+1 >= 0, 0 otherwise, Logb(1 + 2-i) are read from a constants table, X0 is the input value (x), and F0=1. F will contain bx after the iterations. Initializing F with some other value than 1 has the effect of scaling the results by that value. For natural base, a scale factor of 1/4 is required to achieve valid outputs for any input 0 < x < 1. Of course, if you don't need much precision, a look-up table may be a better choice. ------------------------ Subject: Re: Problem: Calcing log from simple math From: Ray Andraka Date: 1995/11/15 Message-ID: <48bet1$9ga@paperboy.ids.net> Newsgroups: sci.math,comp.arch.arithmetic,comp.graphics.algorithms,comp.infosystems.www.authoring.cgi sjledet@ledet.com (Sterling Ledet) wrote: > > This is a math exercise with some real world application. The new > LiveScript programming language built into Netscape 2.0 has just the > following arithmetic operators: > > Addition + > Subtraction - > Multiplication * > Division / > Modulus > Bitwise AND & > Bitwise OR | > Bitwise XOR ^ > Left Shift (<<) > Right Shift (>>) > Double right shift (>>>) > > Is it possible to come up with a formula or subroutine that can > calculate logarithms using just the above operations? > Logs can be computed using shifts and adds plus a few constants (one per iteration). The algorithm is a CORDIC style algorithm, which is simply an iterative approach. The equations needed to compute a log in any base are: F[i+1] = F[i] + d[i]*F[i]*2^(-i) and X[i+1] = X[i] + d[i]*Log(1+2^(-i)) where d[i] = 1 if F[i+1] Date: 1997/05/05 Message-ID: <336E8913.1895@ids.net> Newsgroups: comp.arch.arithmetic,sci.math.num-analysis Bernard Danloy wrote: > > The problem of computing Sinus & Cosinus has been recently considered > in a thread entitled " !!! Fastest Sinus & Cosinus algorithms !!! " > and devoted to the construction of some " trigonometric table " > > I would like to get some more informations and/or comments related > to the computation of one single random sinus or cosinus > Bernard, I wrote about this earlier in the thread, but perhaps you missed it. There is a shift-add algorithm known as CORDIC that can be used to simultaneously compute the sine and cosine of a random angle. The CORDIC algorithm works by rotating a unit vector through the desired angle using a series of fixed rotation angles. The x and y components of the rotated vector are interpreted as cosine and sine respectively. The incremental rotation angles are chosen so that their arctangents are negative powers of two. Rearranging the Givens rotation yields: x(i+1)=cos(v)[x(i)-y(n)*2^(-i)*d(i)] y(i+1)=cos(v)[y(i)+x(n)*2^(-i)*d(i)] where d(i) is +/-1 and controls the direction of the rotation at that iteration At every iteration we rotate either +/- the rotation angle for that iteration. That makes the cos terms constant, so they can be dropped out and treated as a system gain. For Sin and Cos, the system gain can be compensated for by diminishing the unit vector (which starts on the x axis) by the gain factor (approx 1.6) Each iteration of this algorithm (which uses a progressively smaller rotation angle) improves the accuracy of the result by one bit. This is the algorithm used by many (especially early ones) scientific calculators as well as by the intel 8087 math coprocessor. I've got a paper on my website that goes into more detail on the CORDIC algorithm (it is a paper on high speed bit serial techniques; the case study is a CORDIC vector magnitude processor). My website is http://www.ids.net/~randraka ------------------------ Subject: New book on elementary functions From: jmmuller@ens-lyon.fr (Jean-Michel Muller) Date: 1997/05/12 Message-ID: <5l78dv$9bb@cri.ens-lyon.fr> Newsgroups: comp.arch.arithmetic Dear colleagues, This is to announce the imminent outcome of my book, "Elementary functions, algorithms and Implementation", published by Birkhauser Boston. Birkhauser's web page for the book is: http://www.birkhauser.com/cgi-win/ISBN/0-8176-3990-X it includes a description of the book, ordering information, and will include additional resources and selected code. My own web page for the book is: http://www.ens-lyon.fr/~jmmuller/book_functions.html Contents: 1. Introduction 2. Computer Arithmetic 2.1 Floating-Point Arithmetic 2.1.1 Floating-Point formats 2.1.2 Rounding modes 2.1.3 Subnormal numbers and exceptions 2.1.4 ULPs 2.1.5 Testing your computational environment 2.2 Redundant Number Systems 2.2.1 Signed-digit number systems 2.2.2 Radix-2 redundant number systems I. Algorithms Based on Polynomial Approximation and/or Table Lookup 3. Polynomial Approximations 3.1 Least Squares Polynomial Approximations 3.1.1 Legendre polynomials 3.1.2 Chebyshev polynomials 3.1.3 Jacobi polynomials 3.2 Least Maximum Approximations 3.3 Speed of Convergence 3.4 Rational Approximations 3.5 Actual Computation 3.6 Example: the Cyrix FastMath Processor 3.7 Algorithms and Architectures 3.7.1 The E-Method 3.7.2 Estrin's Method 3.8 Miscellaneous 4. Table-Based Methods 4.1 Introduction 4.2 Table-Driven Algorithms 4.2.1 Tang's algorithm for exp(x) in IEEE floating-point arithmetic 4.2.2 ln(x) on [1,2] 4.2.3 sin(x) on [0, pi/4] 4.3 Gal's Accurate Tables Method 4.4 Methods Requiring Specialized Hardware 4.4.1 Wong and Goto, logarithm 4.4.2 Wong and Goto, exponential II. Shift-and-Add Algorithms 5. Shift-and-Add algorithms 5.1 The Restoring and Nonrestoring Algorithms 5.2 Simple Algorithms for Exponentials and Logarithms 5.2.1 The restoring algorithm for exponentials 5.2.2 The restoring algorithm for logarithms 5.3 Faster Algorithms 5.3.1 Faster computation of exponentials 5.3.2 Faster computation of logarithms 5.4 Baker's Predictive Algorithm 5.5 Bibliographic notes 6. The CORDIC Algorithm 6.1 Introduction 6.2 The Conventional Iteration 6.3 Scale Factor Compensation 6.4 CORDIC With Redundant Number Systems 6.4.1 Signed-digit implementation 6.4.2 Carry-save implementation 6.4.3 The variable scale factor problem 6.5 The Double Rotation Method 6.6 Branching CORDIC 6.7 Differential CORDIC 6.8 Computation of cos-1 and sin-1 6.9 Variations on CORDIC 7. Other Shift-and-Add Algorithms 7.1 High-Radix Algorithms 7.1.1 Ercegovac's radix-16 algorithms 7.2 The BKM Algorithm 7.2.1 The BKM iteration 7.2.2 Computation of the exponential function (E-mode) 7.2.3 Computation of the logarithm function (L-mode) 7.2.4 Application to the computation of elementary functions III. Range Reduction, Final Rounding and Exceptions 8. Range Reduction 8.1 Introduction 8.2 Cody and Waite's Method for Range Reduction 8.3 Worst Cases for Range Reduction 8.3.1 A few basic notions on continued fractions 8.3.2 Finding worst cases using continued fractions 8.4 The Payne and Hanek Algorithm 8.5 The Modular Algorithm 8.5.1 Fixed-point reduction 8.5.2 Floating-point reduction 8.5.3 Architectures for Modular Reduction 9. Final Rounding 9.1 Introduction 9.2 Monotonicity 9.3 Exact Rounding: Presentation of the Problem 9.4 Some Experiments 9.5 A "Probabilistic" Approach 9.6 Upper Bounds on m 9.6.1 Frequency of Failures 9.6.2 Computing with one million bits 10. Miscellaneous 10.1 Exceptions 10.1.1 NaNs 10.1.2 Exact Results 10.2 Notes on xy 10.3 Multiple Precision --------------------------------------------------------------------- Jean-Michel Muller, CNRS, Lab. LIP, Ecole Normale Superieure de Lyon 46 Allee d'Italie, 69364 Lyon Cedex 07, FRANCE. Tel. +33 4 72 72 82 29 Fax. +33 4 72 72 80 80 from abroad 04 72 72 82 29 04 72 72 80 80 from France http://www.ens-lyon.fr/~jmmuller email: jmmuller@lip.ens-lyon.fr ------------------------------ Subject: CORDIC: A fast algorithm for trigonometric calculations. From: Ori Berger Date: 1997/04/28 Message-ID: <33646509.49D@iol.co.il> Newsgroups: comp.graphics.algorithms,comp.dsp Due to public demand, here is a short description of the cordic algorithem: We begin with the standard two dimensional rotation matrix (abberviating cos theta to 'c' and sin theta to 's'): [[ c s] [-s c]] Take out the cosine as a scalar out of the matrix, (abbreviating tan theta to t), to get: c * [[ 1 t] [-t 1]] Now, the idea is as follows: To achieve a rotation by some degree, we can use a fixed, predetermined set of (monotonically decreasing) rotations iteratively, each one taken either at sign + or at sign -, until we reach the correct rotation. Now for the computational savings: First, we select the preset rotations to be of angles having the property that tan theta = t = 2^k (for some k, positive or negative); Second, we fix the number of rotations in advance. The computational saving gained by the first suggestion, is that all matrix multiplications amount to additions, subtractions and shifts. The second suggestion makes sure that multiplication by all cos theta (= c) terms above can be put into one constant, because cos theta = cos -theta, so that it does not matter if we take it by a plus sign or a minus sign. Thus, this algorithem in its general form, rotates a given (x,y) pair by a given degree and multiplies the length of the vector by a constant - This may be a useful thing by itself; However, to calculate the sine and cosine of an angle (both at once), you would rotate x=1/(constant due to cosines), y=0 and you would get as a result x=cos theta and y=sin theta. Thus the algorithem generally looks like this (in p-code, types omitted, and untested so it may contain bugs and typos): input: (x,y) - cartesian pair to rotate (theta) - degree to rotate by output (x,y) - rotated (theta) - difference between requested and achieved degree. shift[] = {2,1,0,-1,-2,...}; rotation[] = {atan(2^2), atan(2^1), atan(2^0), atan(2^-1), atan(2^-2), ....}; factor = cos(atan(2^2))*cos(atan(2^1))* cos(atan(2^0))*cos(atan(2^-1))*cos(atan(2^-2))*...; for i = 1 to length(shifts) { /* Target is theta=0 */ if (theta >= 0) { /* Rotate right by matrix multiplication */ d_theta = -rotation[i]; newx = x - (y >> shift[i]); newy = y + (x >> shift[i]); } else { /* Rotate left by matrix multiplication */ d_theta = rotation[i]; newx = x + (y >> shift[i]); newy = y - (x >> shift[i]); } theta = theta + d_theta; x = newx; y = newy; } /* Apply all constant factors at once */ x = x * factor; y = y * factor; /* x and y now contain rotated values */ A few notes are in order: 1. The shift array can usually be disposed of, as it can be the loop index; Note that RIGHT shifts by a negative amount should be interpreted as LEFT shifts by the same absolute value. The rotation array and the scalar factor are of course "real" values, as are x, y and theta; They are usually realised as *fixed* point values - It makes little sense to use this algorithem with floating point values. 2. In fact, the entire code can be unrolled; for example, if we optimize for the Pentium, it may be profitable to replace some of the shifts with effective address calculations. 3. Unless having range problems, the cosine factor should only be applied at the *end* of the calculation, because it is always less than unity (although not significantly so), and would impair fixed point range if used at the beginning. 4. As previously mentioned, in order to calculate both the sine and the cosine, remove the factor multiplication at the end and instead, initialize x to be 1/factor and y to be 0. 5. At the end, theta contains how much more you would need to rotate to get exact results; This information may or may not be useful, but keep in mind that the number of elements in the table generally determines the precision that can be attained, approximately one entry for one bit of precision. With a slight modification, the same algorithem can be used to translate (x,y) coordinates into (r,theta): Initialize theta=0 at the beginning, and instead of selecting a left or right rotation according to theta's sign, select it according to y's sign; When we get to y=0 (or as close as can be reached using the prescribed precision), x will contain the modulos and theta will contain the angle. Also, it can be used to calculate exponential functions: Just use hyperbolic functions instead of trigonometric functions (And of course, the appropriate rotation matrix, which is: [[ ch sh [ sh ch ]] if I am not mistaken), and recall that exp(x) = cosh(x)+sinh(x). Also, logarithems may be taken by noting that if y = exp(x) than: cosh(x)=(y+1/y)/2; sinh(x)=(y-1/y)/2; and using the same slight modification described above to invert the algorithem (although the size of the table needed here might be much larger). I claimed that I have a copy of a Dr. Dobbs magazine that contains this info, but was unable to find it; Maybe someone out there with the CD-ROM copy of Dr. Dobbs can search for the word "CORDIC" and tell us all where it could be found. Anyway, another reference is: "Graphics Gems I; CORDIC vector rotations, pages 494-497". This is taken from the cumulative index at the end of Graphics Gems V, so don't blame me. If my memory serves me right, the article is by Ken Turkowski, and there is some example code on the attached diskette (The contents of which are available in some FTP sites). Remarks anyone? Ori Berger (orib@writeme.com) ---------------------------------------- Subject: Re: CORDIC: A fast algorithm for trigonometric calculations. From: tbryan@mathworks.com (Tom Bryan) Date: 1997/04/28 Message-ID: Newsgroups: comp.graphics.algorithms,comp.dsp Thanks for the CORDIC description. Here is the original reference to the CORDIC algorithm. Jack E. Volder, The CORDIC trigonometric computing technique, IRE Transactions on Electronic Computers (now called IEEE Transactions on Computers), Volume 8, September 1959, pages 330-334. Also, the following textbook has a good description of CORDIC. John P. Hayes, Computer Architecture and Organization, McGraw-Hill, New York, second edition, 1988, p. 282 Best wishes, Tom +=== Tom Bryan ========================== tbryan@mathworks.com ===+ | The MathWorks, Inc. info@mathworks.com | | 24 Prime Park Way http://www.mathworks.com | | Natick, MA 01760-1500 ftp.mathworks.com | +=== Tel: 508-647-7669 ==== Fax: 508-647-7002 ====================+ ------------------------------------------ bject: Re: Algorithms for computing logarithms From: israel@math.ubc.ca (Robert Israel) Date: 1996/12/18 Message-ID: <597iil$fqt@nntp.ucs.ubc.ca> Newsgroups: sci.math In article <32b72865.835574@news.utdallas.edu>, jmathis@utdallas.edu (James Mathis) writes: |> I would appreciate any information anyone can provide on algorithms |> for computing logarithms for arbitrary numbers. What algorithms do |> pocket calculators use? I find it hard to believe that they access a |> logarithm table. Most use a CORDIC algorithm, not only for log but also for many other functions (sqrt, exp, trig, hyperbolic, inverse trig and inverse hyperbolic). Everything is done with additions and shifts. Some use polynomial approximations. A base-2 CORDIC algorithm for ln(x) might proceed as follows. First scale (using ln(2^n x) = n ln(2) + ln(x)) so 1 <= x < 2. Now x can be written as a product of factors of the form 1/(1 - 2^(-i)) (each appearing 0, 1 or 2 times) for i from 2 to infinity (but only take it up to a finite n because of limited machine accuracy), and ln(x) is the sum of the logarithms of those factors (precomputed and stored in a table). The algorithm uses a loop like this: L:= 0; X:= x; for i from 2 to n do while X - 2^(-i)*X >= 1 do X:= X - 2^(-i)*X; L:= L + ln(1/(1-2^(-i)); For more CORDIC info, see, for example, http://ds0.internic.net/pub/graph-ti/calc-apps/info/cordic.txt For an extensive bibliography of CORDIC algorithms, see http://devil.ece.utexas.edu Robert Israel israel@math.ubc.ca Department of Mathematics (604) 822-3629 University of British Columbia fax 822-6074 Vancouver, BC, Canada V6T 1Y4 ---------------------------- Subject: Re: CORDIC Algorithm Question From: boytim@ (Matt Boytim) Date: 1996/03/18 Message-ID: <4ijll0$915@kocrsv08.delcoelect.com> Newsgroups: sci.math The CORDIC algorithm is easy to understand with the proper geometric interpretation. Two good references: Volder, J.E., "The CORDIC Trigonometric computing technique", IRE Trans. on Electronic Computers, Sept. 1959 Hu, Y.H., "CORDIC-Based VLSI architectures for digital signal processing", IEEE Signal Processing Magazine, July 1992 Regards, Matt ----- Subject: Re: Cordics Alogorithms From: Jens Olav Nygaard Date: 1996/06/08 Message-ID: <31B9E526.167E@math.uio.no> Newsgroups: sci.math.num-analysis [More Headers] Here are some references for the CORDIC (COordinate Rotation DIgital Computer) algorithm(s): Charles W. Schelin, Calculator Function Approximation, American Mathematics Monthly, 1983, p. 317-325. Jack E. Volder, The CORDIC Computing Technique, IRE Transactions on Electronic Computers, September 1959, p. 330-334. J. S. Walther, A Unified Algorithm for Elementary Functions, Joint Computer Conference Proceedings, Spring 1971, p. 379-385. (Note that the second is from 1959!) I'm not sure if these methods are very much in use today... Jens Olav Nygaard ----------------------- Subject: Meggitt's algorithms From: baykov@cslab.felk.cvut.cz (Baykov Vladimir) Date: 1995/05/10 Message-ID: Newsgroups: comp.lsi I suggest to pay attention of the people in Computer Arithmetic field to the paper of J.E.Meggitt, published in 1962 in IBM J. Res &Dev. Dr. Meggitt was the first who suggested the algorithms similar to CORDIC, but using normalised constants and normalised variables reducible to zero. It allowed him to improve significantly the accuracy of the iterative al- gorithms. In addition relative errors of the fixed-point representation numbers in Meggitt's algorithms in contrast with Volder's algorithms practicaly don't depend on the word size, like in the floating-point number representation. In other words, Dr.Meggitt suggested not only "pseudodivision and pseudo- multiplication processes" (as indicated in the title of his paper), but also "pseudofloating- point number representation" for CORDIC method. All who are interested in details I can suggest to find in "Computer and Control Abstracts", Author Indexes, beginning from 1971, last name: BAIKOV V. ----------------------------- Subject: errors of CORDIC From: baykov@cslab.felk.cvut.cz (Baykov Vladimir) Date: 1995/05/06 Message-ID: Newsgroups: comp.arch.arithmetic The errors of Volder"s algorithm (CORDIC) are very serious. I began to study them in 1967. The major part of my PhD thesis (1972 year) was related just to these problems. It is necessary to differ maximum and root-mean-square errors of CORDIC. I showed that root-mean-square error (counted in number of not-accuracy of right-most bits) is propotional to the square root of the number of iterations. For example, for function sin x if the word size equals to 16 then r-m-s error is 3-4 least significant bits, and for 32 word size r-m-s error is 5-6 l.s.b. The maximum errors in 3-4 times greater. All these formulas and also the results of statistical experiments you can find in my books , published in 1975 and 1985. To decrease errors of computation, normalising of the iterative variables is very effective. It was suggested by J.Meggit in 1962 . I checked it very thouroughly on a computer in double precision representation. It is realy so. Vladimir Baykov ------------------------------------- Subject: Re: Hardware special functions (CORDIC?) From: jmmuller@ens-lyon.fr (Jean-Michel Muller) Date: 1995/05/04 Message-ID: <3oa0rf$s5p@cri.ens-lyon.fr> Newsgroups: comp.arch.arithmetic In article <3o8a95$7or@lyra.csx.cam.ac.uk>, nmm1@cus.cam.ac.uk (Nick Maclaren) writes: > >I am extremely impressed by the accuracy of most hardware special >function implementations, especially as the software ones of the >1980s and earlier were pretty bad. A previous poster mentioned that >these use mainly CORDIC methods. > >Is this true, and can anyone describe a suitable reference for me >to look them up in? I am a software person not a hardware one, >in case that is relevant, but am not scared of mathematics. > There are two classes of algorithms that are popular for special function implementations. The first one is a mix-up of table-lookup and polynomial interpolation (see for instance P.T.P. Tang, Table-Lookup Algorithms for Special Functions and Their Error Analysis, Proc. of 10th IEEE Symposium on Computer Arithmetic, Grenoble, France, June 1991. I've heard that this algorithm is implemented in the Pentium, but I did not check this information), and the second one is the class of the shift-and-add CORDIC-like methods (the first reference is Volder, "The CORDIC computing technique", IRE Trans Comp., Sept. 1959, but many improvements have been done by Walther, "A unified algorithm for elementary functions", Joint Comp. Conf., Vol. 38, 1971. Some recent contributions can be found in Takagi et al., IEEE Trans. Comp., Sept. 1991, or Duprat et al., IEEE Trans. Comp., Feb. 1993). To my knowledge, CORDIC-like methods were implemented in the 8087 and 68881, and in many pocket computers. Hope this helps -- ---------------------------------------------------------------------- Jean-Michel Muller, CNRS, Lab. LIP, Ecole Normale Superieure de Lyon 46 Allee d'Italie, 69364 Lyon Cedex 07, FRANCE Tel. +33 72 72 82 29 Fax. +33 72 72 80 80 jmmuller@lip.ens-lyon.fr ------------------------------ Subject: Re: Problem: Calcing log from simple math From: jmmuller@ens-lyon.fr (Jean-Michel Muller) Date: 1995/12/01 Message-ID: <49mcqd$gal@cri.ens-lyon.fr> Newsgroups: sci.math,comp.arch.arithmetic,comp.graphics.algorithms [More Headers] In article <49l2ni$2tq@cnn.Princeton.EDU>, ah@dynastar.cs.princeton.edu (Alejo Hausner) writes: >In a recent article, djb@silverton.berkeley.edu (D. J. Bernstein) writes: >> Ray Andraka wrote: >> > F[i+1] = F[i] + F[i]*d[i]*2^(-i) >> > X[i+1] = X[i] + Log(1+2^(-i))*d[i] >> >> ..... > >Actually, I saw this algorithm in a 100-year old table of logarithms. > Of course, it was base 10, not base 2. This class of algorithms is much older: to my knowledge, the first who used a similar algorithm was Henry Briggs (1561-1631), a contem- porary of Neper. Using his algorithm, he was able to publish 15-digit accurate tables of logarithms in his "Arithmetica Logarithmica" (1624). To my opinion, the only algorithm that it is accurate to refer to as a Cordic algorithm is Volder's algorithm: X[n+1] = X[n] -d[n]Y[n]2^(-n) Y[n+1] = Y[n] +d[n]X[n]2^(-n) Z[n+1] = Z[n] - d[n] arctan 2^(-n) or Walther's generalization of it. But even that algorithm can be viewed as a particular case (using complex numbers) of Briggs' algorithm: if you write T[n] = X[n] + iY[n] (with i^2 = -1), then you will get: T[n+1] = T[n] (1+id[n]2[-n]) and Z[n] will be the imaginary part of the sequence L[n+1] = L[n] - Log (1+id[n]2[-n]). ---------------------------------------------------------------------- Jean-Michel Muller, CNRS, Lab. LIP, Ecole Normale Superieure de Lyon 46 Allee d'Italie, 69364 Lyon Cedex 07, FRANCE Tel. +33 72 72 82 29 Fax. +33 72 72 80 80 jmmuller@lip.ens-lyon.fr ------------------------------- Subject: Re: Problem: Calcing log from simple math From: Ray Andraka Date: 1995/11/15 Message-ID: <48bet1$9ga@paperboy.ids.net> Newsgroups: sci.math,comp.arch.arithmetic,comp.graphics.algorithms,comp.infosystems.www.authoring.cgi [More Headers] sjledet@ledet.com (Sterling Ledet) wrote: > > This is a math exercise with some real world application. The new > LiveScript programming language built into Netscape 2.0 has just the > following arithmetic operators: > > Addition + > Subtraction - > Multiplication * > Division / > Modulus > Bitwise AND & > Bitwise OR | > Bitwise XOR ^ > Left Shift (<<) > Right Shift (>>) > Double right shift (>>>) > > Is it possible to come up with a formula or subroutine that can > calculate logarithms using just the above operations? > Logs can be computed using shifts and adds plus a few constants (one per iteration). The algorithm is a CORDIC style algorithm, which is simply an iterative approach. The equations needed to compute a log in any base are: F[i+1] = F[i] + d[i]*F[i]*2^(-i) and X[i+1] = X[i] + d[i]*Log(1+2^(-i)) where d[i] = 1 if F[i+1]