Re: My Homework: how to calculate natural log, e?

Rich Townsend wrote:
....
> And likewise, with x = 0.01 with Intel ifort:
>
> Naturally, e = 1.01005016708417 1.01005015761443
> x?
>
> So GT's claim of 1ulp precision are demonstrably false. Now there's
> a surprise!

Not only is Gerry's solution wrong, it's *_S_L_ O_W_*!!!
Just the loop:

do i = 0, kmax
pwr = 2 ** real(i,8)
Y(0,i) = (1.d0 + x/pwr) ** pwr
enddo

is slower than the professional versions of EXP() in
their entirety. This actually explains why Gerry's
solution is not accurate as well. 1.0d0+X/pwr (even
with pwr equal to a power of two) throws away a lot
of the low order bits of X. Raising that to a large
power (even when that power is itself a power of
two) amplifies that error a lot. Y, and the calculations
using and defining it, should all be in extended, or even
quad precision. Gerry's solution does well if X is a
small integer (no low order bits to get lost).

The actual algorithms used in hardware and run-time
support libraries that work on hardware without the
function built-in are much different. They usually work
by reducing the exponent to a narrower range (like -2^(-64)
through 2^(-64)), evaluating the EXP of that smaller
argument using a low-order polynomial or rational
polynomial interpolation, and multiplying the result by
some large power of 2 to undo the consequences of the
argument reduction. These can be done a lot quicker
and with a lot less silicon than Gerry's turkey. Here's
one that's still not as good as the real professionals, but
it's a lot better than Gerry's:

integer :: i, j
real(kind(0.0d0)) :: x, s, xsave, ans, lib_ans
real(kind(0.0d0)),parameter :: Lg2 = log(2.0d0), overLg2 = 1.0d0/log(2.0d0)

read(*,*)xsave
x = xsave*overLg2
i = x
x = xsave-i*Lg2
s = 1.0d0
do j = 17,1,-1
s = 1+x*s/j
end do
ans=s*2.0d0**i ! the answer
lib_ans = exp(xsave) ! to compare
! output the generated answer, the library's answer and the
! approximate relative error in ULPs
write(*,'(g25.19)') ans, lib_ans, abs(lib_ans-ans)/lib_ans/epsilon(xsave)

This should be within an ULP or so everywhere if all
the intermediates (including Lg2 and overLg2) were in
extended precision instead of just double. As it is,
the more the input magnitude exceeds Lg2, the larger the
typical error (almost in direct proportion). Gets to be a
couple of hundred ULPs when the input is around 700.

--
J. Giles

"I conclude that there are two ways of constructing a software
design: One way is to make it so simple that there are obviously
no deficiencies and the other way is to make it so complicated
that there are no obvious deficiencies." -- C. A. R. Hoare


.

Relevant Pages

  • Re: So far OT its over the horizon: Anyone into 3D modelling?
    ... > screws, 2.5 metre ball screws have a certain minimum cost, even if you ... > cut out ALL the middle men you'll be looking at IRO 500 quid per axis ... you could use precision threaded bar and a couple of nuts to ... > eliminate backlash, but it's gonna wear quickly, take more power to ...
    (uk.rec.motorcycles)
  • Re: So far OT its over the horizon: Anyone into 3D modelling?
    ... screws, 2.5 metre ball screws have a certain minimum cost, even if you ... cut out ALL the middle men you'll be looking at IRO 500 quid per axis ... you could use precision threaded bar and a couple of nuts to ... eliminate backlash, but it's gonna wear quickly, take more power to ...
    (uk.rec.motorcycles)
  • Re: For those who missed it -another PC milestone
    ... The processing power of the GPU is now available for general purpose ... do you mean "single and double precision" or "single or double ...
    (uk.politics.misc)
  • Re: What sports car are YOU?
    ... > You are a Ferrari 360 Modena! ... power, passion, precision, and style. ...
    (rec.sport.football.college)
Loading