"Fossies" - the Fresh Open Source Software Archive

Member "brlcad-7.32.4/doc/notes/hypot.txt" (29 Jul 2021, 3070 Bytes) of package /linux/misc/brlcad-7.32.4.tar.bz2:

    1 From ken@turtlevax.UUCP (Ken Turkowski) Sun Nov 18 15:41:29 1984
    2 Date: 18 Nov 84 20:41:29 GMT
    3 
    4 > >> 	It would also be nice if sin( x ) and cos( x ) could be
    5 > >> 	computed simultaneously with reduced cost.  I doubt if this is
    6 > >> 	possible but would like to know if it is.
    7 >
    8 > > 	sine = sin(x);
    9 > > 	cosine = sqrt(1.0 - sine*sine);
   10 >
   11 > I wish people would check things out before posting them.  On our
   12 > 4.2BSD system, the sqrt() approach is actually SLOWER than using
   13 > cos().  (Under UNIX System V it is very slightly faster.)  Now,
   14 > if I wanted cos(x)^2, I would certainly use 1 - sin(x)^2, which
   15 > I already knew about (I too went to high school).
   16 >
   17 > What I had in mind was more along the lines of a CORDIC algorithm
   18 > or some other obscure but useful approach.
   19 
   20 In Tom Duff's SIGGRAPH '84 tutorial on Numerical Methods for Computer
   21 Graphics, he writes (quote until end of article):
   22 
   23 Calculating Euclidean distance probably accounts for 90% of the square
   24 roots expended in computer graphics applications.  Aside from the
   25 obvious geometric uses, most illumination models require that the unit
   26 normal to each surface be computed at each pixel.
   27 
   28 [Moler-Morrison, "Replacing Square Roots by Pythagorean Sums", IBM
   29 Journal of Research and Development", vol. 27, no. 6, pp. 577-581, Nov.
   30 1983] describes an algorithm for computing sqrt(a^2+b^2) which is fast,
   31 robust and portable.  The naive approach to this problem (just writing
   32 sqrt(a*a+b*b) ) has the problem that for many cases when the result is
   33 a representable floating point number, the intermediate results may
   34 overflow or underflow.  This seriously restricts the usable range of a
   35 and b.  Moler and Morrison's method never causes harmful overflow or
   36 underflow when the result is in range.  The algorithm has _____cubic
   37 convergence, and depending on implementation details may be even faster
   38 than sqrt(a*a+b*b).  Here is a C function that implements their
   39 algorithm:
   40 
   41 double hypot(p,q)
   42 double p,q;
   43 {
   44     double r, s;
   45     if (p < 0) p = -p;
   46     if (q < 0) q = -q;
   47     if (p < q) { r = p; p = q; q = r; }
   48     if (p == 0) return 0;
   49     for (;;) {
   50 	r = q / p;
   51 	r *= r;
   52 	s = r + 4;
   53 	if (s == 4)
   54 	    return p;
   55 	r /= s;
   56 	p += 2 * r * p;
   57 	q *= r;
   58     }
   59 }
   60 
   61 This routine produces a result good to 6.5 digits after 2 iterations,
   62 to 20 digits after 3 iterations, and to 62 digits after 4 iterations.
   63 In normal use you would probably unwind the loop as many times as you
   64 need and omit the test altogether.  Moler and Morrison's paper
   65 describes the generalization of this algorithm to n dimensions, and
   66 exhibits a related square root algorithm which also has cubic
   67 convergence.  [Dubrelle, "A Class of Numerical Methods for the
   68 Computation of Pythagorean Sums", IBM Journal of Research and
   69 Development, vol. 27, no. 6, pp. 582-589, Nov. 1983] gives a geometric
   70 justification of the algorithm and presents a set of generalizations
   71 that have arbitrarily large order of convergence.
   72 --
   73 Ken Turkowski @ CADLINC, Palo Alto (soon Menlo Park), CA
   74 UUCP: {amd,decwrl,flairvax,nsc}!turtlevax!ken
   75 ARPA: turtlevax!ken@DECWRL.ARPA