## "Fossies" - the Fresh Open Source Software Archive

As a special service "Fossies" has tried to format the requested text file into HTML format (style: standard) with prefixed line numbers. Alternatively you can here view or download the uninterpreted source code file.

```    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
```