Path: utzoo!utgpu!news-server.csri.toronto.edu!mailrus!ames!haven!decuac!e2big.mko.dec.com!bacchus.pa.dec.com!decwrl!wuarchive!usc!orion.oac.uci.edu!ucivax!milne
From: milne@ics.uci.edu (Alastair Milne)
Newsgroups: comp.lang.pascal
Subject: Re: TP: Exponentiation?
Message-ID: <26C3B699.8294@ics.uci.edu>
Date: 11 Aug 90 07:41:13 GMT
References: <950036@hpclapd.HP.COM> <2041@bimacs.BITNET> <17410@dime.cs.umass.edu>
Organization: UC Irvine Department of ICS
Lines: 45

In <17410@dime.cs.umass.edu> eli@smectos.gang.umass.edu (Eli Brandt) writes:

>In article <2041@bimacs.BITNET> ehrlich@bimacs.biu.ac.il.UUCP (Gideon Ehrlich) writes:
>>
>>Very easy. Either use    x**y = exp(y*ln(x))
>>           or write your own function ,using Tailor sequence.

    I was taught you shouldn't use exp(y*ln(x)) unless you really have
    to, because both of them are bound to be implemented in your library
    as approximations, and often not very good ones, so the result will
    get back a certain amount of error and multiply it by error.  Not good.

>you should consider these modifications:
>   1) Check for negatives and 0^0, obviously.
>   2) A significant fraction of calls will probably involve integer powers.  
>     The way to handle this is to use a simple for loop for low powers 
>     and a repeated squaring loop for higher powers.  The crossover point 
>     depends on the machine and compiler but is usually around 5 to 10.

   I was shown one that handles any non-negative integer exponent, with limited
   complexity (linear, I think, though I'm too tired to remember just now).
   You may remember the shift-and-add implementation of multiplication.
   This is a shift-and-multiply for exponentiation.  Start your power at 1.
   If the low bit of the exponent is 1, multiply the base into the power.
   Either way, square the base and shift the exponent right 1.  
   Repeat until the exponent is 0 and return the accumulated power.  
   At most, 32 multiplications (assuming a 16-bit word), usually much less.  
   And if the base is also an integer, 16 of the multiplications become 
   left shifts.

   For floating-point bases, this reduces the error a lot (unless your
   squaring operation is poor).

   This should be checked, since I'm doing it from memory, but I believe
   it's correct.

   As for floating-point exponents, all I can suggest is finding your 
   own very good approximations of ln and exp.  (I believe there's
   a continued fraction for either exp(x) or x*exp(x) that behaves 
   very well and has amazingly few operations.)

   That, however, plumbs the depths of what I remember from numeric
   analysis.  Hope it helps.

   Alastair Milne