Re: Change definition of complex::norm

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

nbecker@fred.net writes:

 >   template<typename _Tp>
 >     inline _Tp
 >     norm(const complex<_Tp>& __z)
 >     {
 >       _Tp __res = abs(__z);
 >       return __res * __res;
 >     }
 > 
 > where abs is defined as:
 > 
 >   template<typename _Tp>
 >     inline _Tp
 >     abs(const complex<_Tp>& __z)
 >     {
 >       _Tp __x = __z.real();
 >       _Tp __y = __z.imag();
 >       const _Tp __s = abs(__x) + abs(__y);
 >       if (__s == _Tp())  // well ...
 >         return __s;
 >       __x /= __s; 
 >       __y /= __s;
 >       return __s * sqrt(__x * __x + __y * __y);
 >     }
 > 
 > I claim:
 > 
 > 1) the current definition is not as useful, since it only works where
 >    sqrt is defined
 > 
 > 2) the current definition is much less efficient
 > 
 > 3) the current definition is possibly less accurate

I agree on all 3 points.  x^2 + y^2 is definitely preferable to
abs(z)^2.  In fact, it would be better to do it the other way around -
define abs(x) as norm_factor*norm(z/norm_factor)^2.

Aside from doing that, the above abs() can be improved.  Dividing by
__s = abs(__x) + abs(__y) isn't the best thing in the world.  Dividing
by __s reduces accuracy.  Computing __s can overflow when abs(__z)
wouldn't, and dividing by it can cause the computation of __x*__x to
underflow.

Better would be to use frexp, ldexp & modf to put together a more
appropriate normalizer than the sum of the absolute values.  This
would improve accuracy.

Even just using the larger of abs(__x) & abs(__y) as the normalizer
(as is done in the xlispstat code) instead of the sum of the two would
be an improvement.  Then you'd have (doing away with annoying
underscores, and before hand optimizing):

   x = abs(z.real());
   y = abs(z.imag());

   big    = max(x,y);
   little = min(x,y);

   normalized_y = little/big;

   return(big * sqrt(1 + normalized_y * normalized_y));

which saves a multiplication & a division at the cost of an additional
test, and won't overflow & underflow as often.

Using frexp, ldexp & modf would make multiplying & dividing by the
normalization factor a simpler operation than normal multiplication &
division, and would further improve accuracy.

-- 
Harvey Stein
Bloomberg LP
hjstein@bloomberg.com