Mail Archives: djgpp/1997/01/22/04:42:10
Orlando Andico wrote:
>
> i've been reading the past few posts about FP speed.. i would like to
> know if it's possible/profitable to have a modified pow() function. is the
> pow() done by the FP unit? i'm working on an MPEG audio codec and i found
> out that 60% of the time is spent in pow()..
>
60% is a lot : there are two ways to spend time in math functions :
having slow functions and calling them too much... can't you reduce the
number of calls?
1/ by simplifying the formulae,
2/ by saving on repeated calculations of pow(x,y) for the same (x,y),
3/ by using multiplies instead of pow() when y is an integer,
pow(x,2)=x*x
pow(x,y+1)=x*pow(x,y) (this can save one pow() when you need both
pow(x,y) and pow(x,y+1)...)
It is hard to say more with little information about what you calculate,
could you post the formula?
This done, you can try to speed pow()up:
There are two pow() in DJGPP, the one in libc (linked by default if you
don't specify -lm) is in assembly, it uses the FP instructions. The one
in libm uses a formula (I suppose it also uses the FP instruction, when
compiled on any FP capable machine...).
The two versions do not run at the same speed, on my machine the libc
pow() is faster than the libm one. If your program spends so much time in
pow(), it might be worth experimenting.
Writing a replacement for pow() is feasible, but will only be profitable
if
1/ the x and y in pow(x,y) are limited to a restricted range of values
2/ you need less precision on your replacement function (say single
precision instead of double, or even less).
Here is a general method for calculating approximations to pow(x,y), I
think the libm pow() is built on this scheme.
pow(x,y)=exp(y*log(x)),
decompose y*log(x) as
y*log(x)= k*log(2) - b, where b is a positive number , between 0 and
log(2), then
pow(x,y)=2^k exp(-b)
and use some approximation to the exponential function (or lookup
table), to compute exp(-b)
Here are a few simple one (from Handbook of Mathematical Functions, by
Abramowitz/Stegun):
exp(-x)=1-.9664*x+.3536x*x
(error under 3e-3)
exp(-x)=1-.9998684x+.4982926*x^2-.1595332*x^3+.0293642x^4
(error under 3e-5)
Once you have k and exp(-b), you can rebuild the floating point result
(4 byte float), by multiplying exp(-b) by 2^24 (this yield the mantissa),
and adding 128 to k to get the exponent.
(Hope this is understandable... sorry if it is not)
Hope this helps,
Francois
- Raw text -