This is the mail archive of the gcc@gcc.gnu.org mailing list for the GCC project.

Index Nav: Message Nav: [Date Index] [Subject Index] [Author Index] [Thread Index] [Date Prev] [Date Next] [Thread Prev] [Thread Next]

# Re: Change definition of complex::norm

• To: gdr at codesourcery dot com (Gabriel Dos Reis)
• Subject: Re: Change definition of complex::norm
• From: Brad Lucier <lucier at math dot purdue dot edu>
• Date: Thu, 1 Nov 2001 08:58:43 -0500 (EST)
• Cc: lucier at math dot purdue dot edu (Brad Lucier), gcc at gcc dot gnu dot org, hjstein at bloomberg dot com, nbecker at fred dot net, gdr at codesourcery dot com, bkoz at redhat dot com

```> | norm_2 uses the definition in std_complex.h (with the fixed abs, i.e.,
> | abs_1). norm_1 uses the simpler, faster, algorithm for norm proposed
> | by nbecker.  Here, the simpler algorithm gives an anwer that loses
> | all precision.  On the other hand, I can't judge how important it
> | is that a simpler, faster, algorithm gives 0.0 as the answer instead
> | of 4.940656e-324.
>
> In that case, I prefer to retain the current implementation and
> document the option.  After all, nbecker can just explicitly
> specialize for his type.

I'm sorry---I posted before working out all the issues (don't we all ;-).

The relevant error tolerance when calculating a quantity x is roughly

\max (|x| * 2^{-53}, 2^{-53-1024})

The first quantity is the relative tolerance, the second is the absolute
tolerance for numbers in the denormal range.  One would like the
computation of any quantity x to be within a small multiple of this
tolerance.  There can also be qualitative issues---returning a finite
answer as often as possible, etc.

Under this error criteria, the formula

\real(z)^2 + \imag(z)^2

is as valid as what is currently used now.  It is also much faster
and simpler, and works for a larger range of base types, as nbecker
points out.  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.  (The complex /= uses norm, so /= will fail
very slightly more often with the new formula than with the current
formula, but this problem seems trivial compared with much more serious
problems in the current naive definition of /=, and fixing /= seems
nearly impossible without assuming IEEE 754 arithmetic, perhaps a
full implementation.)

So I recommend that std_complex.h adopt

\real(z)^2 + \imag(z)^2

as the definition of norm.