/* 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.henorm.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/c/math/enorm.c.html © Copyright 2001 by John Halleck, All Rights Reserved. This snapshot was last modified on May 17th, 2001 And the underlying file was last modified on May 11th, 2000