Relay-Version: version B 2.10 5/3/83; site utzoo.UUCP
Posting-Version: version B 2.10.2 9/18/84; site brl-tgr.ARPA
Path: utzoo!linus!decvax!genrad!panda!talcott!harvard!seismo!brl-tgr!tgr!milne@icse.uci.edu
From: milne@icse.uci.edu (Alastair Milne)
Newsgroups: net.micro
Subject: Re: Arctan approximation for generating pi
Message-ID: <1122@brl-tgr.ARPA>
Date: Tue, 31-Dec-85 04:14:16 EST
Article-I.D.: brl-tgr.1122
Posted: Tue Dec 31 04:14:16 1985
Date-Received: Wed, 1-Jan-86 04:57:03 EST
Sender: news@brl-tgr.ARPA
Lines: 66


>  Of course there are many better types of approximations than power series:
>  ....
>  about a point very close to zero, as some people suggested.

    I quite agree.  I did not, however, want to clog the net with discussions
    of numeric analysis considerations: it is clogged enough already
    with other things.  The question at hand was how to approximate
    arctan, and the suggestions I had seen so far, though correct,
    would have been costly both to their implementors and their users.
    I simply wanted to point out why, and suggest the best alternative
    I know of, specifically for arctan.  Certainly you could (for
    instance) apply Chebyshev economisation to the power series, and
    with x as small as 1/5 and 1/239, the economised series should
    converge very rapidly.  But it still won't compete with the
    continued fraction.  (Notice that I didn't even mention Horner's
    method.)

    I grant I hadn't thought about dynamically changing accuracy, which
    is one advantage that a plain power series does have.  I don't
    recall just now whether that can be done in any way with continued 
    fractions.

>  It is true that a**x will be computed as exp(x*ln a), x real, [well known
>  formula] but a**n will usually be implemented as a*a*...*a (n-1 multiplies)
>  n integer, at least that's how IBM Fortran compilers do it. But then
>  in generating the next term in a power series the previous term is
>  multiplied by x and, usually, divided by some constant so the exponentiation
>  problem should never arise in practice. Incidently, I once heard a
>  great horror story how a**x = exp(x*ln a) almost sunk an engineering
>  Ph.D. thesis here at UofT. It happened nearly 20 years ago before the
>  builtin functions we all rely on got as good as they are now.

     I'm not sure whether the practices of IBM Fortran compilers
     constitute recommendations.  Many would say just the opposite.
     Either way, I'm not sure increasingly long series of
     floating-point multiplications are much to be preferred: greater
     accuracy, but still expensive to compute.  And, though the next
     term of a power series SHOULD be generated the way you say, the
     average programmer might not implement it that way on seeing its
     arithmetic description.  Furthermore, I'm sure none of us wants to
     encourage programming techniques based on knowledge of how our
     favourite compiler implements feature x.

     Horror stories in numeric analysis!  Yes, I imagine we could start a
     whole new distribution list for them alone, if it hasn't actually been
     done already.  There was a 100x100 sparse matrix that had to be inverted
     ... well, never mind.  Let's hope that we're generating them a little
     less frequently these days.  

>  Wirth probably didn't include exponentiation in Pascal because of
>  the weird arithmetic on the CDC6600 computer, ...
>  multiply a number by one and change the last few bits!

   Wirth's design considerations for Pascal were more oriented around
   human than machine (though you have made it much clearer to me why
   the routines PACK and UNPACK were included).  Specifically,
   concerning exponentiation, Wirth said that it was misused in most of
   the places where it occurred, and that he was not going to include
   something in his new language which he was then going to have to
   teach his students not to use.  Teaching was, after all, his primary
   aim in creating Pascal.


Thanks for the reply,
Alastair Milne