| Author |
Message |
Glenn Zelniker
Guest
|
 Posted:
Fri Dec 23, 2005 9:15 am Post subject:
Computing x^y, sort-of quick, somewhat dirty |
|
|
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 |
|
| Back to top |
|
 |
Andor
Guest
|
| Posted:
Fri Dec 23, 2005 3:33 pm Post subject:
Re: Computing x^y, sort-of quick, somewhat dirty |
|
|
Hi Glenn, nice to have you back!
Glenn Z wrote:
| Quote: | 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 |
|
| Back to top |
|
|
Carlos Moreno
Guest
|
| Posted:
Fri Dec 23, 2005 5:16 pm Post subject:
Re: Computing x^y, sort-of quick, somewhat dirty |
|
|
Glenn Zelniker wrote:
| Quote: | 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
-- |
|
| Back to top |
|
|
Bob Monsen
Guest
|
| Posted:
Sat Dec 24, 2005 1:15 am Post subject:
Re: Computing x^y, sort-of quick, somewhat dirty |
|
|
On Thu, 22 Dec 2005 22:46:38 -0500, Glenn Zelniker wrote:
| Quote: | 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. |
|
| Back to top |
|
|
Jim Adamthwaite
Guest
|
| Posted:
Sat Dec 24, 2005 9:15 am Post subject:
Re: Computing x^y, sort-of quick, somewhat dirty |
|
|
Glenn Zelniker wrote:
| Quote: | 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
|
Hi 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
;-------------------------------------------------------------------- |
|
| Back to top |
|
|
Peter K.
Guest
|
| Posted:
Sat Dec 24, 2005 5:16 pm Post subject:
Re: Computing x^y, sort-of quick, somewhat dirty |
|
|
Jerry Avins <jya@ieee.org> writes:
| Quote: |
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. |
|
| Back to top |
|
|
Jerry Avins
Guest
|
| Posted:
Sat Dec 24, 2005 5:16 pm Post subject:
Re: Computing x^y, sort-of quick, somewhat dirty |
|
|
Jim Adamthwaite wrote:
....
| Quote: | 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.
ŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻŻ |
|
| Back to top |
|
|
Jim Adamthwaite
Guest
|
| Posted:
Sun Dec 25, 2005 9:15 am Post subject:
Re: Computing x^y, sort-of quick, somewhat dirty |
|
|
Jerry Avins wrote:
| Quote: | Jim Adamthwaite wrote:
...
the 56002 devaluation kit.
...
What did the price go down to? (Spell checkers! Fah!)
Jerry
|
It was intentional.
Jim A. |
|
| Back to top |
|
|
Glenn Zelniker
Guest
|
| Posted:
Mon Dec 26, 2005 9:15 am Post subject:
Re: Computing x^y, sort-of quick, somewhat dirty |
|
|
Bob Monsen wrote:
| Quote: |
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 |
|
| Back to top |
|
|
Glenn Zelniker
Guest
|
| Posted:
Mon Dec 26, 2005 9:15 am Post subject:
Re: Computing x^y, sort-of quick, somewhat dirty |
|
|
Andor wrote:
| Quote: | Hi Glenn, nice to have you back!
|
Nice to be back, Andor. Are you still at Weiss?
| Quote: | 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 |
|
| Back to top |
|
|
Bob Monsen
Guest
|
| Posted:
Tue Dec 27, 2005 1:16 am Post subject:
Re: Computing x^y, sort-of quick, somewhat dirty |
|
|
On Mon, 26 Dec 2005 01:11:39 -0500, Glenn Zelniker wrote:
| Quote: | 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?
|
I don't know. However, you can easily run a test program on a PC to find
out, at least with double precision on a PC (which is what I did when I
blasted this out the other night). In fact, here is the output from my
trivial test program, massaged to keep error, defined as (x^n -
binomal())/(x^n) = (e) under 0.001. c is the count of iterations required
to do this (at least when compared against the linux pow() function).
0.01^0=1 c=0 e=0
0.01^0.1=0.631582 c=304 e=0.000624887
0.01^0.2=0.398502 c=340 e=0.000395107
0.01^0.3=0.251437 c=352 e=0.00024846
0.01^0.4=0.158647 c=352 e=0.000157706
0.01^0.5=0.1001 c=345 e=9.95747e-05
0.01^0.6=0.0631585 c=332 e=6.27259e-05
0.01^0.7=0.03985 c=313 e=3.92514e-05
0.01^0.8=0.0251439 c=284 e=2.50121e-05
0.01^0.9=0.0158648 c=238 e=1.58481e-05
0.11^0=1 c=0 e=0
0.11^0.1=0.802706 c=26 e=0.00077064
0.11^0.2=0.643734 c=29 e=0.000633494
0.11^0.3=0.516237 c=30 e=0.000512505
0.11^0.4=0.413934 c=31 e=0.000356675
0.11^0.5=0.331968 c=30 e=0.000305913
0.11^0.6=0.266213 c=29 e=0.0002414
0.11^0.7=0.213497 c=27 e=0.00020497
0.11^0.8=0.171199 c=25 e=0.000152972
0.11^0.9=0.137292 c=21 e=0.00012396
0.21^0=1 c=0 e=0
0.21^0.1=0.856235 c=13 e=0.000731476
0.21^0.2=0.73244 c=15 e=0.000552894
0.21^0.3=0.626689 c=15 e=0.000557423
0.21^0.4=0.536144 c=15 e=0.000485893
0.21^0.5=0.458641 c=15 e=0.000383062
0.21^0.6=0.392424 c=14 e=0.000382881
0.21^0.7=0.335643 c=14 e=0.000250614
0.21^0.8=0.287215 c=12 e=0.000285108
0.21^0.9=0.245642 c=11 e=0.00017262
0.31^0=1 c=0 e=0
0.31^0.1=0.890271 c=8 e=0.000790612
0.31^0.2=0.791887 c=9 e=0.000712535
0.31^0.3=0.704201 c=10 e=0.000466796
0.31^0.4=0.626382 c=10 e=0.000424319
0.31^0.5=0.557125 c=10 e=0.000348943
0.31^0.6=0.495678 c=9 e=0.000436995
0.31^0.7=0.440807 c=9 e=0.000299584
0.31^0.8=0.392128 c=8 e=0.000305788
0.31^0.9=0.348756 c=7 e=0.000238144
0.41^0=1 c=0 e=0
0.41^0.1=0.915271 c=6 e=0.00057166
0.41^0.2=0.837114 c=7 e=0.00043884
0.41^0.3=0.765784 c=7 e=0.000478193
0.41^0.4=0.700476 c=7 e=0.000450814
0.41^0.5=0.640697 c=7 e=0.00038464
0.41^0.6=0.585993 c=7 e=0.000300093
0.41^0.7=0.536178 c=6 e=0.000444724
0.41^0.8=0.490309 c=6 e=0.000274357
0.41^0.9=0.448513 c=5 e=0.000278416
0.51^0=1 c=0 e=0
0.51^0.1=0.935651 c=4 e=0.000768879
0.51^0.2=0.874487 c=5 e=0.000481308
0.51^0.3=0.817635 c=5 e=0.000542395
0.51^0.4=0.764414 c=5 e=0.000529082
0.51^0.5=0.71461 c=5 e=0.000467326
0.51^0.6=0.668017 c=5 e=0.000377657
0.51^0.7=0.62444 c=5 e=0.000275912
0.51^0.8=0.584013 c=4 e=0.000492096
0.51^0.9=0.545755 c=4 e=0.000231799
0.61^0=1 c=0 e=0
0.61^0.1=0.952465 c=3 e=0.000692779
0.61^0.2=0.906207 c=4 e=0.000337185
0.61^0.3=0.862571 c=4 e=0.000388886
0.61^0.4=0.820989 c=4 e=0.000388396
0.61^0.5=0.781376 c=4 e=0.000351407
0.61^0.6=0.743649 c=4 e=0.000291023
0.61^0.7=0.707725 c=4 e=0.000217998
0.61^0.8=0.673934 c=3 e=0.000548125
0.61^0.9=0.641177 c=3 e=0.000267026
0.71^0=1 c=0 e=0
0.71^0.1=0.967216 c=2 e=0.000884672
0.71^0.2=0.934101 c=3 e=0.000306058
0.71^0.3=0.902718 c=3 e=0.000363198
0.71^0.4=0.872347 c=3 e=0.000373498
0.71^0.5=0.842963 c=3 e=0.00034821
0.71^0.6=0.814542 c=3 e=0.000297387
0.71^0.7=0.78706 c=3 e=0.00022992
0.71^0.8=0.760492 c=3 e=0.000153582
0.71^0.9=0.735216 c=2 e=0.00047748
0.81^0=1 c=0 e=0
0.81^0.1=0.979375 c=2 e=0.000227138
0.81^0.2=0.959112 c=2 e=0.000380484
0.81^0.3=0.939209 c=2 e=0.000469107
0.81^0.4=0.919668 c=2 e=0.000501881
0.81^0.5=0.900487 c=2 e=0.0004875
0.81^0.6=0.881668 c=2 e=0.000434474
0.81^0.7=0.863209 c=2 e=0.000351136
0.81^0.8=0.845112 c=2 e=0.000245646
0.81^0.9=0.827375 c=2 e=0.000125993
0.91^0=1 c=0 e=0
0.91^0.1=0.991 c=1 e=0.000386735
0.91^0.2=0.982 c=1 e=0.000685359
0.91^0.3=0.973 c=1 e=0.000896699
0.91^0.4=0.963028 c=2 e=4.95754e-05
0.91^0.5=0.953987 c=2 e=4.82986e-05
0.91^0.6=0.945028 c=2 e=4.3173e-05
0.91^0.7=0.937 c=1 e=0.000885495
0.91^0.8=0.928 c=1 e=0.000672554
0.91^0.9=0.919 c=1 e=0.000377131
--
Regards,
Bob Monsen
Our minds are finite, and yet even in those circumstances of finitude, we
are surrounded by possibilities that are infinite, and the purpose of human
life is to grasp as much as we can out of that infinitude.
- Alfred North Whitehead |
|
| Back to top |
|
|
Glenn Zelniker
Guest
|
| Posted:
Wed Dec 28, 2005 1:16 am Post subject:
Re: Computing x^y, sort-of quick, somewhat dirty |
|
|
Bob,
Do you have a reference for the iterative gen. binomial algorithm? I
need to look into the convergence properties a bit more carefully. I
couldn't find anything by googling. Also, could you post or email a
snippet of the code you used for the table you generated? I'm not sure
about the pseudocode you posted the other day.
GZ |
|
| Back to top |
|
|
Bob Monsen
Guest
|
| Posted:
Wed Dec 28, 2005 9:15 am Post subject:
Re: Computing x^y, sort-of quick, somewhat dirty |
|
|
On Tue, 27 Dec 2005 15:01:13 -0800, Glenn Zelniker wrote:
| Quote: | Bob,
Do you have a reference for the iterative gen. binomial algorithm? I
need to look into the convergence properties a bit more carefully. I
couldn't find anything by googling. Also, could you post or email a
snippet of the code you used for the table you generated? I'm not sure
about the pseudocode you posted the other day.
GZ
|
The binomial theorem was Fermat and Pascal, I think. Newton was the first
to use it with negative exponents.
http://mathworld.wolfram.com/BinomialTheorem.html
Here is the example:
----------------------------
#include <stdio.h>
#include <math.h>
double binomial_power(double x, double n, int * cnt)
{
double y = 1.0; /* result */
double z = 1.0; /* running term */
double i = 1.0; /* for n! */
double p = pow(x,n);
x = 1.0 - x; /* computes (1+z)^n */
while (fabs((p-y)/p) > 0.001)
{
z *= -x*n / i;
n -= 1.0;
i += 1.0;
y += z;
*cnt = *cnt + 1;
}
return y;
}
test_binomial()
{
double x = 0.1;
double n = 0.1;
int i, j;
x = .01;
for (i = 0; i < 10; i++)
{
n = 0.0;
for (j = 0; j < 10; j++)
{
int cnt = 0;
double b;
b = binomial_power(x,n,&cnt);
printf("%lg^%lg=%lg\tc=%d\te=%lg\n", x, n, b, cnt, fabs(pow(x,n)-b));
n += 0.1;
}
x += 0.1;
}
}
int main()
{
test_binomial();
}
--
Regards,
Bob Monsen
One cannot inquire into the foundations and nature of mathematics without
delving into the question of the operations by which the mathematical
activity of the mind is conducted. If one failed to take that into account,
then one would be left studying only the language in which mathematics is
represented rather than the essence of mathematics.
- Luitzen Brouwer |
|
| Back to top |
|
|
|
|
|
|
FAQ
Memberlist
Register
Profile
Log in to check your private messages
Log in