/* enorm.c compute Euclidean distance functions
 * (C) 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;
}