Message-ID: <32E56283.122@pobox.oleane.com>
Date: Wed, 22 Jan 1997 01:42:43 +0100
From: Francois Charton <deef AT pobox DOT oleane DOT com>
Organization: CCMSA
MIME-Version: 1.0
To: Orlando Andico <orly AT gibson DOT eee DOT upd DOT edu DOT ph>
CC: djgpp AT delorie DOT com
Subject: Re: floating-point speed
References: <Pine DOT SGI DOT 3 DOT 93 DOT 970121155211 DOT 2597A-100000 AT gibson DOT eee DOT upd DOT edu DOT ph>
Content-Type: text/plain; charset=us-ascii
Content-Transfer-Encoding: 7bit

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