I've seen and used all the old tricks for doing ln, log2, 2^x, exp(x). I'm using an old fixed-point processor (DSP56XXX) and need to do a fast x^y approximation (to within a few percent) where x, y < 1. I'd like to do it directly -- i.e., no long sequence of log, add, mul, and exp operations because I'm tight on RAM/ROM (and curiously, not too tight for CPU cycles). I can only spare a handful of words but more than a handful of cycles. Any ideas? Glenn Zelniker
Computing x^y, sort-of quick, somewhat dirty
Started by ●December 22, 2005
Reply by ●December 23, 20052005-12-23
Hi Glenn, nice to have you back! Glenn Z wrote:> I've seen and used all the old tricks for doing ln, log2, > 2^x, exp(x). I'm using an old fixed-point processor > (DSP56XXX) and need to do a fast x^y approximation (to > within a few percent) where x, y < 1. I'd like to do it > directly -- i.e., no long sequence of log, add, mul, and exp > operations because I'm tight on RAM/ROM (and curiously, not > too tight for CPU cycles). I can only spare a handful of > words but more than a handful of cycles.The way I would do it is by computing 2^(y log2(x)). Both exp and log can be approximated with a polynomials for (x,y) in the given domain - you can use the same code for the polynomial evaluation, just change the coeffcients and the input arguments to save code space. log2(x) ~= p1(x) , 0 < x < 1 2^x ~= p2(x) , 0 < x < 1 => x^y = 2^(y log2(x)) = p2(y p1(x)). log2(x) has unbounded output range for 0 < x < 1, this will inherently be a problem in any implementation (specially fixed-point). Regards, Andor
Reply by ●December 23, 20052005-12-23
Glenn Zelniker wrote:> I've seen and used all the old tricks for doing ln, log2, 2^x, exp(x). > I'm using an old fixed-point processor (DSP56XXX) and need to do a fast > x^y approximation (to within a few percent) where x, y < 1. I'd like to > do it directly -- i.e., no long sequence of log, add, mul, and exp > operations because I'm tight on RAM/ROM (and curiously, not too tight > for CPU cycles). I can only spare a handful of words but more than a > handful of cycles. > > Any ideas?Well, if all else fails, you could always decide what order you want for the polynomial, and then obtain a "best-fitting" polynomial for a set of values. The fitting could even be optimized following the PDF of your variables, as to reduce the error around the points that occur more often. For instance, you could take a grid of 101 by 101 values (x going from 0 to 1 in steps of 0.01, and same thing for y); then, compute the *actual* value of x^y, and then write one equation for each of the 101*101 points. You end up with a bit over 10 thousand *linear* equations for your few unknowns (the unknowns are the polynomial coefficients) -- well, or less; you could do a grid of 10 by 10, or 20 by 20, etc. You could then try different polynomials, different orders, and play with it until the approximation that you obtain gives you a worst-case error that is acceptable. It may not be as crazy as it sounds -- honest! :-) HTH, Carlos --
Reply by ●December 23, 20052005-12-23
On Thu, 22 Dec 2005 22:46:38 -0500, Glenn Zelniker wrote:> I've seen and used all the old tricks for doing ln, log2, > 2^x, exp(x). I'm using an old fixed-point processor > (DSP56XXX) and need to do a fast x^y approximation (to > within a few percent) where x, y < 1. I'd like to do it > directly -- i.e., no long sequence of log, add, mul, and exp > operations because I'm tight on RAM/ROM (and curiously, not > too tight for CPU cycles). I can only spare a handful of > words but more than a handful of cycles. >How about the generalized binomial theorem? It'll be crappy as x -> 0, but if you can live with that, it is simple to compute. The radius of convergence (as shown) for non-natural n is 0 < x < 2. double binomial(double x, double n, int terms) { double y = 1.0; double z = 1.0; double i = 1.0; x = 1.0 - x; // It really computes (1 + y)^n while (terms-- > 0) { z *= -x*n / i; n -= 1.0; i += 1.0; y += z; } return y; } -- Regards, Bob Monsen "A poor portrayal is about as effective as a good one for most people. They don't see the defects; they see a symbol which inspires their deepest emotions; it recalls to them the Agony and Sacrifice of God." "Jubal, I thought you weren't a Christian?" "Does that make me blind to human emotion? I was saying that the crummiest painted plaster crucifix can evoke emotions in the human heart so strong that many have died for them. The artistry with which such a symbol is wrought is irrelevant.
Reply by ●December 24, 20052005-12-24
Glenn Zelniker wrote:> I've seen and used all the old tricks for doing ln, log2, > 2^x, exp(x). I'm using an old fixed-point processor > (DSP56XXX) and need to do a fast x^y approximation (to > within a few percent) where x, y < 1. I'd like to do it > directly -- i.e., no long sequence of log, add, mul, and exp > operations because I'm tight on RAM/ROM (and curiously, not > too tight for CPU cycles). I can only spare a handful of > words but more than a handful of cycles. > > Any ideas? > > Glenn ZelnikerHi Glen. Some time back I did an audio spectrum analyzer for the 56002 devaluation kit. The dB scale conversion used a table lookup/interpolate scheme. While you were actually looking for X^Y, you should be able to simply change the calculation line in the table generating macro: lg2 set @LOG(tmp)/@LOG(2) ;calc log2 of X to be able to achieve 2^X or other nonlinear functions. The table covers a 2:1 numeric range. The upper significant bits of A are used to index into the table, the LSB's are used to interpolate. For more entries, the error introduced by linear interpolation is reduced. Note that it has only been tested for correctness using a limited set of values. You may have to determine how small the table can be while preserving sufficient accuracy for your application. It is reasonably quick & has constant exec time for in-range numbers. You are welcome to the rest of the spec analyzer code if you wish. The only thing that might reduce its appeal is that it requires the use of a CMOS analog switch chip (4051 or similar) to multiplex the DC voltage on several pot wipers into codec channel B input for instrument control purposes such as gain, span, peak hold, video averaging etc. Output is from codec channel A to an oscilloscope in the form of a pseudo video signal with sync pulses and calibration markers. Jim Adamthwaite ;-------------------------------------------------------------------- ;"log2mac.ASM" - logarithm (base 2) macro ; Uses lookup table technique ; Input: A:R0 = (x) normalised f.p., A range = +0.5..+0.9999999 ; Output: A:R0 = log2(A) ; Destroys: lots LOG2 macro move #$ffff,m1 ;ensure linear modulus tst a #>(1.0/32768.0),x1 ;test (x), prep 15-bit shift const jne _log2a ;(x) = zero? move #>$003f,r0 ;yes, prep rslt = -infinity move #>$800000,a jmp _log2x ;and quit _log2a move a,x0 ;shift (X), split into A1:A0 mpy x0,x1,a #>(Log2Bas-$80),x1 ;prep X1 = table ptr add x1,a a0,b ;A1 = adrs of first point lsr b a1,r1 ;b = US fract, R1 = &tbl(i) move b1,x0 ;X0 = US fract move x:(r1)+,a ;A = tbl(i), point R1 to next move x:(r1),b ;B = tbl(i+1) sub a,b ;B = tbl(i+1) - tbl(i) move b,y0 ;slope to Y0 mac x0,y0,a #0,b ;A = tbl(i) + (slope * fract) _log2b move r0,b2 ;fetch FP exponent of (x) asr b ;align B to suit natural binary point add b,a #0,r0 ;A2:A1:A0 = log2(x) do #24,_log2x ;cnvt (restore) to FP norm r0,a ;use DO to be interruptible nop nop _log2x endm ;-------------------------------------------------------------------- ; Table of logarithms (to base 2) for LOG2 function. ; Indexing values of X start at 0.5 & end at +1.0 LOG2TBL macro log2base org x:log2base tmp set 0.5 ;start value of X DUP 128+1 lg2 set @LOG(tmp)/@LOG(2) ;calc log2 of X dc lg2 tmp set tmp+(1.0/256.0) ;adjust X ENDM ENDM ;--------------------------------------------------------------------
Reply by ●December 24, 20052005-12-24
Jim Adamthwaite wrote: ...> the 56002 devaluation kit.... What did the price go down to? (Spell checkers! Fah!) Jerry -- Engineering is the art of making what you want from things you can get. �����������������������������������������������������������������������
Reply by ●December 24, 20052005-12-24
Jerry Avins <jya@ieee.org> writes:> > What did the price go down to? (Spell checkers! Fah!) >Taken from http://en.wikipedia.org/wiki/Spell_checker: Eye halve a spelling chequer, It came with my pea sea, It plainly marques four my revue Miss steaks eye kin knot sea. Eye strike a key and type a word And weight four it two say Weather eye am wrong oar write It shows me strait a weigh. As soon as a mist ache is maid It nose bee fore two long And eye can put the error rite Its rarely ever wrong. Eye have run this poem threw it I'm shore your pleased two no Its letter perfect in it's weigh, My chequer tolled me sew.
Reply by ●December 25, 20052005-12-25
Jerry Avins wrote:> Jim Adamthwaite wrote: > > ... > >> the 56002 devaluation kit. > > ... > > What did the price go down to? (Spell checkers! Fah!) > > JerryIt was intentional. Jim A.
Reply by ●December 26, 20052005-12-26
Bob Monsen wrote:> > How about the generalized binomial theorem? It'll be crappy as x > -> 0, but if you can live with that, it is simple to compute. The radius > of convergence (as shown) for non-natural n is 0 < x < 2. > > double binomial(double x, double n, int terms) > { > double y = 1.0; > double z = 1.0; > double i = 1.0; > > x = 1.0 - x; // It really computes (1 + y)^n > > while (terms-- > 0) > { > z *= -x*n / i; > n -= 1.0; > i += 1.0; > > y += z; > } > > return y; > }This is exactly the kind of thing I was looking for. Thank you, thank you. It'll only take a few minutes to code up when I get back to the lab, but I thought I'd jump the gun and ask a leading question. For x not too close to zero, how many iterations is typical to get, say, ten or so bits of accuracy in the result? GZ
Reply by ●December 26, 20052005-12-26
Andor wrote:> Hi Glenn, nice to have you back!Nice to be back, Andor. Are you still at Weiss?> The way I would do it is by computing 2^(y log2(x)). Both exp and log > can be approximated with a polynomials for (x,y) in the given domain - > you can use the same code for the polynomial evaluation, just change > the coeffcients and the input arguments to save code space. > > log2(x) ~= p1(x) , 0 < x < 1 > 2^x ~= p2(x) , 0 < x < 1 > > => x^y = 2^(y log2(x)) = p2(y p1(x)).Believe it or not, I don't even have room for the coefficients! I've already coded it up this way and the whole thing takes too many lines of code once you see it written down. Thanks for taking the time to respons, though. Best, GZ