Re: Square root

Gregory Toomey wrote:


Of course not. The Newton-Raphson iteration is much faster.

However, in some cases, the use of a square root indicates that the
actual thing to be computed is a Pythagorean sum. That can be calculated
better.

/**
  * Compute the Pythagorean Sum (I.E. what is the distance given Delta X,
  * and Delta Y)
  *
  * Traditionally, this is done as sqrt(x*x+y*y), but this has the
  * problem that the squared quantity can overflow long before the result
  * would have.
  *
  * Cleve Moler and Donald Morrison show (as documented in their article
  * in the IBM Journal of Research and Development, VOL27 No 6 NOV 83)
  * that this can 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 numerically stable. It
  * is very fast.
  *
  * This code started life in MASM Assembley Language on a Univac. That
  * code was then ported to Fortran IV, on a Univac. That code was then
  * ported to C. That code was then ported to Ada. That code was then
  * used as the base of this port. SO.. it has a lot of inherited uglies.
  */

public class PythagoreanSum {

     /**
      * Calculate a Pythagorean sum
      *
      * @param x
      * @param y
      * @return sqrt(x*x + y*y), calculated directly
      */
     static public double sum(final double x, final double y) {
         final double delta_x = Math.abs(x);
         final double delta_y = Math.abs(y);
         double guess, errval;
         if (delta_x > delta_y) {
             guess = delta_x;
             errval = delta_y;
         } else {
             guess = delta_y;
             errval = delta_x;
         }
         // The following code is the loop unrolling of:
         // temp = error / guess;
         // r = temp * temp;
         // s = r / (r + 4.0);
         // guess += 2.0 * s * guess;
    // error *= s;
         // The loop is unrolled because checking for convergence takes
         // more time than the computation.
         // 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 the
         // Length with each iteration dramaticly shrinking ERROR.
         // After this we have one significant figure in the result
         double temp = errval / guess;
         double r = temp * temp;
         double s = r / (r + 4.0);
         guess += 2.0 * s * guess;
         // After this we have 6 significant digits in the result
         errval *= s;
         temp = errval / guess;
         r = temp * temp;
         s = r / (r + 4.0);
         guess += 2.0 * s * guess;
         // After this we have 20 significant digits in the result
         errval *= s;
         temp = errval / guess;
         r = temp * temp;
         s = r / (r + 4.0);
         guess += 2.0 * s * guess;
         return guess;
     }

     /**
      * Tester
      * @param args
      */
     public static void main(String[] args) {
         for (double i = 3.0, j = 4.0, k = 4.0;
             i < 1000.0; i += 2.0, j += 4.0, k += j) {
             System.out.println (i + '\t' + k + '\t' + sum(i, k));
         }
     }

}

--
John W. Kennedy
"The blind rulers of Logres
Nourished the land on a fallacy of rational virtue."
   -- Charles Williams. "Taliessin through Logres: Prelude"