Re: Change definition of complex::norm

[HJSTEIN at bloomberg dot com (Harvey J. Stein)]

Brad Lucier <lucier@math.purdue.edu> writes:

 > double
 > abs_1(double x, double y) {
 > 
 >   double abs_x = (x > 0. ? x : -x);
 >   double abs_y = (y > 0. ? y : -y);
 >   double s     = (abs_x > abs_y ? abs_x : abs_y);
 > 
 >   if (s == 0.)
 >     return s;
 > 
 >   x /= s;
 >   y /= s;
 > 
 >   return s * sqrt( x * x + y * y);
 > }

Wouldn't this be a little better:

double
abs_3(double x, double y) {

  double abs_x = (x > 0. ? x : -x);
  double abs_y = (y > 0. ? y : -y);
  double s, min;

   if (abs_x > abs_y) {
      s   = abs_x;
      min = y;
   } else {
      s   = abs_y;
      min = x;
  }
  if (s == 0.)
    return s;

  min /= s;

  return s * sqrt( 1 + min * min);
}

It avoids computing (max(abs_x, abs_y)/s)^2, which has to be 1, thus
saving a division & a multiplication.

As for abs(z)^2 doing better as norm(z) than real(z)^2 + imag(z)^2, I
think it's due more to the normalization done by abs(z) than anything
else.  In general, using the sqrt(x)^2 *should* be a little less
accurate, and you could add normalization to the real(z)^2 + imag(z)^2
if it really helps...

-- 
Harvey Stein
Bloomberg LP
hjstein@bloomberg.com