Re: Change definition of complex::norm

[Brad Lucier <lucier at math dot purdue dot edu>]

> | There is a small group of arguments for which this
> | formula gives zero when the current formula gives a nonzero result,
> | but this now seems to me to be irrelevant---there are many nonzero z
> | for which both formulas give zero for the norm, and the fact that
> | there are a few more of these cases with the simpler formula does not
> | seem important to me.
> 
> I presume we disagree on exactly where the line should be drwan. From
> my point of view, I consider norm() should deliver an answer as close
> as possible to the real value; that means, for non-subnormal a value
> acceptable with respect to relative error criterion, and for
> subnormal, a value acceptable with respect to absolute error criterion.
> Which is why I prefer to retain the current implementation for norm().

Perhaps I wasn't clear in my previous post.

My point (which was in the part of my post that you don't quote) is
that the formula

norm(z) = real(z)^2 + imag(z)^2

does give good relative error for normal results, and good absolute
error for denormal results.

I should point out that all this is simpler on x86 *if* one can
assume that all floating-point register spills are spilt to extended
precision (which we can't assume with the current gcc).  Then the two formulas

abs(z) = sqrt(real(z)^2 + imag(z)^2)

and

norm(z) = real(z)^2 + imag(z)^2

provide full accuracy and safety for both double- and single-precision
types _Tp, and one would need the more complicated formulas only for
extended-precision arguments.

Brad