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