NetNews Usenet Archive 1992 #16/NN_1992_16.iso/spool/comp/programm/2052

Text File | 1992-07-22 | 7.7 KB | 180 lines

Newsgroups: comp.programming
Path: sparky!uunet!zaphod.mps.ohio-state.edu!wupost!gumby!yale!yale.edu!ira.uka.de!rz.uni-karlsruhe.de!stepsun.uni-kl.de!uklirb!posthorn!vier!neuhaus
From: [e-mail] (Stephan Neuhaus (HiWi Mattern))
Subject: Re: Calculating a cube root
Message-ID: <neuhaus.711800570@vier>
Sender: [e-mail] (News system account)
Nntp-Posting-Host: vier.informatik.uni-kl.de
Organization: University of Kaiserslautern, Germany
References: <[e-mail]> <[e-mail]> <[e-mail]> <[e-mail]>
Date: Wed, 22 Jul 1992 10:22:50 GMT
Lines: 167

[e-mail] (Dave Schaumann) writes:

>In article <[e-mail]>, jeroen@cs (Jeroen Fokker) writes:
>>
>>[e-mail] (Randi Joseph) asked:
>>
>>> Can anyone present an example routine in C
>>> to calculate a cube root of a number?  Thanks!
>>
>>[e-mail] (Frank Bennett) replied:
>>
>>> y = exp(log(x)/3.0 );
>>> you didn't say it had to be fast!
>>
>>A faster function uses Newton's method, using iteration:

>I don't think so, at least for most values.  Library routines like exp(),
>log(), and pow() are generally defined using a polynomial expansion, which
>require perhaps 5 floating point multiplies, and yield a reasonably accurate
>answer.

This is precisely the problem with this solution.  It is
possible---and quite cheaply at that---to make the answer correct to
maybe two ulps using Newton's method.  Depending on the magnitude of
x, this may not be the case for the exp (log (x)/3.0) method, for if
you look at the relative error, you will see that it is less than

    x^{(2/3)\epsilon} - 1 + \epsilon x^{(2/3)\epsilon}

(Send email if you want the details.)  \epsilon here is the machine
precision, which is usually defined to be half an ulp.  For IEEE
single precision, this is about 1e-8, for double it is 1e-16,
approximately.  Since the second \epsilon term is smaller than the
first by a factor of 1e8, we can discard it and look at the remaining
terms.

This function of x behaves like a polynomial, and is always greater
than or equal to -1.  Moreover, it is strictly increasing, and will
surpass every threshold sooner or later.  Let's look at the threshold
for x where the relative error becomes larger than two ulps (=
4\epsilon).

    x^{(2/3)\epsilon} - 1 >= 4\epsilon
<==>    x^{(2/3)\epsilon} >= 1 + 4\epsilon
<==>    x >= (1 + 4\epsilon)^{3/(2\epsilon)}    (since x, \epsilon > 0)

For \epsilon = 1e-8, we get x >= 404, and for \epsilon = 1e-16, we get
x >= 782.  Therefore, large arguments tend to make the exp (log
(x)/3.0) solution more inaccurate than need be.

>Looking at your answer, we see

>>    #define EPSILON 0.00001

>which of course limits the accuracy quite severly.  Of course, one could
>write the algorithm so that it terminates when the value stops converging.

No.  That is not possible.  A theorem of Numerical Analysis states
that if you implement a theoretically at least linearly convergent
series on a machine, then either the limit is exactly a machine
number, in which case the process will halt, or it isn't, in which
case there exists a region of indeterminacy around the limit in which
your values will get stuck.  Your sequence will become periodic.  The
precise a-priori formula is (in TeX notation)

|\hat{x}_n - x| < \Delta/(1-q) + {q^n \over 1-q}|\hat{x}_1 - \hat{x}_0|

Where \Delta is the worst-case error that you make when you calculate
the next iterate from the most recent, q is the constant of
convergence, x is the exact limit, and \hat{x}_k are the iterates.  As
you can see, the formula does not tell the whole story, since it only
says something about the maximum error, but then, the maximum error
will not drop below \Delta/(1-q), plus I can give you the complete
theorem with proof if you want.

>But this still leaves the problem of speed.  A quick check shows that your
>algorithm iterates as much as 15 times for values less than 1000. Clearly,
>unless you have a /very/ inefficient exp() and log(), using them will be a
>speed win as well.

Typical exp and log routines will spend some time looking at the
argument and performing argument reduction.  The inner loops in a
Newton step can be heavily optimized by a good compiler.  Furthermore,
a good initial guess for the cube root will cut down the number of
iterations enormously.

>There's also another advantage to using library routines.  Presumably, someone
>was paid to make sure they are correct.

Yes, but only to so many ulps.  There is (IMHO) no excuse for making a
numerical routine less stable and accurate than need be.  The Newton
solution wins over the exp/log solution for x > 1000, approximately.

>>        return( a-b<EPSILON && b-a<EPSILON );
>                 ^^^^^^^^^^^^^^^^^^^^^^^^^^

>This expression will always yield 0.

As some follow-uppers have already remarked, this is not true.  Worse
still is the fact that no one seems to have recognized the danger of
using the absolute error instead of the relative error.  Using the
absolute error means that you claim that x^{1/3} behaves everywhere
like the constant 1.  It's far better to root out the case x = 0 first
and then use

    fabs (a - b) < EPSILON*min (a, b)

I use min here to get a worst-case bound on the relative error.  Of
course, you will have to make sure that a and b are both positive, but
this is not a problem, since Newton converges monotonically (I think).

For huge machine values of a and b, the absolute error cannot be less
than EPSILON, unless a and b are exactly equal.  For very small values
of a and b, the criterion may be satisfied immediately, even though
your values might not even be close to the cube root of x.

Another approach that has been used in other routines, e.g., square
root, is the following:  Write x as

    x = 2^e * f    1.0 <= f < 2.0

Then,

    x^{1/3} = 2^{e/3} * f^{1/3}.

Algorithm:

    1. Write e as 3*e' + e''    (0 <= e'' < 3)
    2. Compute f' = 2^{e''} * f    (1.0 <= f' < 8.0)
    3. Calculate r = f'^{1/3} by Newton's method (however, see below)
    4. Set result = 2^{e'} * r

To get an initial guess for f', it usually pays to approximate the
function x^{1/3} over [1,8) by a polynomial of degree 3 (say) with an
approximately equal-ripple error curve.  This requires some modest
amount of work but will in general provide you with excellent initial
guesses.  For example, a 3rd degree polynomial for sqrt(x) over [1,2)
will give you three to four correct decimals, making only three Newton
steps necessary for 24 decimals (more than IEEE double has).

Careful argument reduction pays both in terms of speed (because you
can do an a-priori error analysis to find out how many Newton steps
you need) and accuracy (because of the fixed, rather small interval).

Everyone who is interested in these topics should read:

    F.S. Acton: _Numerical Methods that (Usually) Work_, Harper &
        Row, 1990 (?)

    P.J. Plauger: _The Standard C Library_, Prentice-Hall, 1992

    W. Cody and W. Waite: _Software Manual for the Elementary
        Functions_, (Publisher unknown, I can look it up.  You
        can find the reference in Plauger's book, too.)

Sorry this took so long.

--
+------- No .sig? OH, COME ON! --------+-"The derivative snuggles  close to-+
| email:  [e-mail] | the function---whatever to snuggle |
| snail mail: Stephan Neuhaus/Hilgard- | means; I'm too old for that"       |
| ring 32/6750 Kaiserslautern/Germany  |                       -- Math Prof |
--
+------- No .sig? OH, COME ON! --------+-"The derivative snuggles  close to-+
| email:  [e-mail] | the function---whatever to snuggle |
| snail mail: Stephan Neuhaus/Hilgard- | means; I'm too old for that"       |
| ring 32/6750 Kaiserslautern/Germany  |                       -- Math Prof |