/* enorm.c compute Euclidean distance functions
* © Copyright 1999 by John Halleck
* All rights reserved
*
* Copyright applies to this implementation, not to the underlying
* algorithm.
*/
/* Version of August 18th, 1999 */
#include "enorm.h"
double enorm2 (double X, double Y) {
/*
*
* Compute the Euclidean vector norm for a 2 dimension vector.
* (I.E. what is the distance given Delta X, DeltaY)
* **** THIS IS A NORMALIZED RESULT.
* Traditionally, this is done as sqrt(x*x+y*y).
* Cleve Moler and Donald Morrison discovered (as documented in their article
* in the IBM Journal of Research and Development, VOL27 No 6 NOV 83)
* that this could actually computed more directly.
*
* Advantages of this algorithm:
* It never squares X or Y, so it never overflows unless the
* answer does.
* It is very numericly stable.
* It is very very very fast.
*/
/* This code was ported from Fortran, so it shows a lot of Fortranisms. */
double GUESS;
/* Our current guess */
double ERROR;
/* The amount our error is wrong by. */
double temp;
/* Just a place for the intermediate values. */
double DX, DY;
/* Delta X and Y */
double R, S;
/* Intermediate values */
/* See the original paper for the complete details */
/* Set up first guess (approximation) */
DX = X>=0 ? X : -X;
DY = Y>=0 ? Y : -Y;
GUESS = DX > DY ? DX : DY;
ERROR = DX < DY ? DX : DY;
/* The following code is the loop unrolling of the following:
*
* temp = ERROR/GUESS;
* R = temp*temp;
* S = R / (4.0+R);
* GUESS += 2.0 * S * GUESS ;
* ERROR *= S;
*
* 1 iteration gives 1 significant digits in the result.
* 2 iterations give 5 significant digits,
* 3 iterations give 20 significant digits,
* 4 iterations give 62 significant digits.
* Etc. (The number of significant decimal digits roughly triples
* with each iteration.)
*
* The loop invariant is that SQRT(GUESS**2+ERROR**2) is ENORM2,
* with each iteration dramaticly shrinking ERROR.
*/
/* After this section we have one significant figure in the result */
temp = ERROR/GUESS;
R = temp*temp;
S = R / (4.0+R);
GUESS += 2.0 * S * GUESS ;
ERROR *= S;
/* After this section we have 5 significant digits in the result */
temp = ERROR/GUESS;
R = temp*temp;
S = R / (4.0+R);
GUESS += 2.0 * S * GUESS ;
ERROR *= S;
/* After this section we have 20 significant digits in the result. */
temp = ERROR/GUESS;
R = temp*temp;
S = R / (4.0+R);
GUESS += 2.0 * S * GUESS ;
ERROR *= S;
/* If we were to do it again, we'd have 62 significant digits... */
return GUESS;
}
This page is http://www.cc.utah.edu/~nahaj/cave/survey/code/c/enorm.c.html
© Copyright 2000 by John Halleck, All Rights Reserved.
This snapshot was last modified on August 23rd, 2000
And the underlying file was last modified on May 11th, 2000