The Math Forum



Search All of the Math Forum:

Views expressed in these public forums are not endorsed by NCTM or The Math Forum.


Math Forum » Discussions » sci.math.* » sci.math

Notice: We are no longer accepting new posts, but the forums will continue to be readable.

Topic: Fast exponent and logarithm, given initial estimate
Replies: 29   Last Post: Nov 8, 2004 2:31 AM

Advanced Search

Back to Topic List Back to Topic List Jump to Tree View Jump to Tree View   Messages: [ Previous | Next ]
Glen Low

Posts: 41
Registered: 12/13/04
Re: Fast exponent and logarithm, given initial estimate
Posted: Nov 8, 2004 2:31 AM
  Click to see the message monospaced in plain text Plain Text   Click to reply to this topic Reply


> I need a fast exponent and logarithm routine for floats. It should
> preferably use only adds, subtracts and multiplies; divides and square
> roots are permissible if necessary. It should not use large lookup
> tables. An additional wrinkle (or help) is that there's a exponent and
> a logarithm estimate available. (In case you're wondering, this is the
> situation with the Altivec SIMD instruction set.)

Just thought I'd provide a quick summary of our findings in case it's
helpful for others.

Exponent is fairly straightforward.

Suppose x = a + b
=> exp(x) = exp(a) * exp(b)

We have accurate 2^n for any integer n, so if b = n * ln(2), we have

exp(x) = exp(a) * 2^n

So choose b = floor (x / ln(2)) * ln(2)
and a = x - b, restricts a to [0, ln(2)]

Performing a 6th degree polynomial minimax approximation yields a
relative error < 10e-8 on a as Gert noted.

Now err(exp(x)) = err(exp(a)*exp(b)) = a*err(exp(b)) + b*err(exp(a))
err(exp(x))/exp(x) = err(exp(b))/exp(b) + err(exp(a))/exp(a)

So if b is accurate, then the relative error of x = relative error in
a, so we're done.

Logarithm is more complicated.

Suppose x = a * b
=> ln(x) = ln(a) + ln(b)

Now if b = 2^n, then ln(b) = n * ln(2).

However look at the error analysis:

err(lg2(x)) = err(lg2(a) + lg2(b)) = err(lg2(a)) + err(lg2(b))
Suppose b is accurate, so we have
err(lg2(x)) = err(lg2(a))
Two cases: either b = 1 or b != 1.
For b = 1,
err(lg2(x))/lg2(x) = err(lg2(a))/lg2(a) i.e the relative error is
equal
For b != 1
err(lg2(x))/lg2(x)
= err(lg2(a))/(lg2(a) + lg2(b))
= {err(lg2(a))/lg2(a)} / {1 + lg2(b)/lg2(a)}, a != 1
= relative error in a / (1 + lg2(b)/lg2(a) }

So the relative error is sensitive to the choice of a and b. In
particular, some of the proposed intervals for a were [1,2),
[1,2^i/8), [1/4, 1/2) all suffer from the problem that lg2(b)/lg2(a)
can be close or equal to -1, thus causing the relative error to
balloon out around x = 1.

Instead I used the interval [.75, 1.5) and if the relative error is
well bounded in that interval, then ln(b)/ln(a) where b != 1 cannot be
close to -1.

Using the standard IEEE 754 representation we know we can express
x = a * 2^n, where a in [1.0, 2.0), except for x = 0.

So if a in [0.75, 2.0) simply subtract 1 from n to get [.75, 1.0) and
we get our interval [.75, 1.5).

Then I performed a 3,4 rational polynomial minimax approximation on
that range on ln(x)/(x-1) to get rid of the zero at x = 1 with
relative error of almost 5e-9. Nicely enough, even though numerator
and denominator are evaluated using Horner's method, the scheduler
should be able to interleave them in the CPU pipeline since they are
independent.

Both exponent and logarithm have about 8 multiply-add pipeline depth
and about 5 ulps difference from the scalar version. Better, both are
about 6x - 7x faster than calls to the scalar library version.

Any comments, critiques would be appreciated.

The algorithms will eventually make their way into my macstl library,
version 0.2 will be dual licensed under RPL and a proprietary license.
I'll make the announcement here when I'm ready.

<a href="http://www.pixelglow.com/macstl/">http://www.pixelglow.com/macstl/</a>

Cheers,
Glen Low, Pixelglow Software
www.pixelglow.com

P.S. I used the 30 day trial version of Mathematica and an online
crash course in rational approximation to get the results. Much thanks
to Gert and other contributors in this thread.



Date Subject Author
10/18/04
Read Fast exponent and logarithm, given initial estimate
Glen Low
10/18/04
Read Re: Fast exponent and logarithm, given initial estimate
Jeremy Watts
10/19/04
Read Re: Fast exponent and logarithm, given initial estimate
Peter Spellucci
10/19/04
Read Re: Fast exponent and logarithm, given initial estimate
Glen Low
10/18/04
Read Re: Fast exponent and logarithm, given initial estimate
bv
10/19/04
Read Re: Fast exponent and logarithm, given initial estimate
Glen Low
10/19/04
Read Re: Fast exponent and logarithm, given initial estimate
George Russell
10/19/04
Read Re: Fast exponent and logarithm, given initial estimate
Glen Low
10/20/04
Read Re: Fast exponent and logarithm, given initial estimate
George Russell
10/20/04
Read Re: Fast exponent and logarithm, given initial estimate
Glen Low
10/21/04
Read Re: Fast exponent and logarithm, given initial estimate
Christer Ericson
10/21/04
Read Re: Fast exponent and logarithm, given initial estimate
Glen Low
10/22/04
Read Re: Fast exponent and logarithm, given initial estimate
Christer Ericson
10/19/04
Read Re: Fast exponent and logarithm, given initial estimate
Martin Brown
10/19/04
Read Re: Fast exponent and logarithm, given initial estimate
Glen Low
10/19/04
Read Re: Fast exponent and logarithm, given initial estimate
Richard Mathar
10/19/04
Read Re: Fast exponent and logarithm, given initial estimate
Glen Low
10/20/04
Read Re: Fast exponent and logarithm, given initial estimate
Gert Van den Eynde
10/20/04
Read Re: Fast exponent and logarithm, given initial estimate
Glen Low
10/20/04
Read Re: Fast exponent and logarithm, given initial estimate
Richard Mathar
10/21/04
Read Re: Fast exponent and logarithm, given initial estimate
Gert Van den Eynde
10/21/04
Read Re: Fast exponent and logarithm, given initial estimate
bv
10/22/04
Read Re: Fast exponent and logarithm, given initial estimate
Glen Low
10/22/04
Read Re: Fast exponent and logarithm, given initial estimate
Peter Spellucci
10/22/04
Read Re: Fast exponent and logarithm, given initial estimate
Glen Low
10/23/04
Read Re: Fast exponent and logarithm, given initial estimate
bv
10/24/04
Read Re: Fast exponent and logarithm, given initial estimate
Gert Van den Eynde
10/25/04
Read Re: Fast exponent and logarithm, given initial estimate
Peter Spellucci
10/20/04
Read Re: Fast exponent and logarithm, given initial estimate
Gert Van den Eynde
11/8/04
Read Re: Fast exponent and logarithm, given initial estimate
Glen Low

Point your RSS reader here for a feed of the latest messages in this topic.

[Privacy Policy] [Terms of Use]

© The Math Forum at NCTM 1994-2018. All Rights Reserved.