Change definition of complex::norm

Change definition of complex::norm

[nbecker]

We talked about this before, but it appears that in 3.0.2 this still
isn't changed.

There are clear disadvantages to the current definition of
complex::norm in terms of abs, not the least of which is that it won't work for
complex<int>, for example.   I can't imagine any drawback to defining
norm as

real(z) * real(z) + imag(z) * imag(z).

Can we please change it?
-----BEGIN PGP PUBLIC KEY BLOCK-----
Version: GnuPG v1.0.6 (GNU/Linux)
Comment: Processed by Mailcrypt 3.5.6 and Gnu Privacy Guard < http://www.gnupg.org/ >

mQGiBDZa9ZcRBACYGMoAmHIBUR19lDLZNJhgxGtqVchV7OiwniGIE0UpwRj08fDX
/KO7/cXgXDZqFEgHF98e6Gbm4efyyC7seP4Ye8Av3n8h8PMv307lQieJd5qQVvwx
vwJWGHsX1EOsv/Suzb2ZcYllU4dgrdBIkRLQ5tsJPiWtxjsfBsONGqWmIwCghmQA
GayzNTFUUy0JkGP8SEJRycED/0GvchcxrSnN0FDe5HqM2YzNOnQYGEasAgRSNoG7
O87uudA3j+Hh4GQSD7VgleYArCXqfaNd8pj+EY0ykGjcTJk07aAl+Ib8UrQ8eNk/
RON0+/ZRN6QGqte1lokR969AgVFDQaHV0IctElZdpRg+JbKUiBn3iYaY7xfYYr1z
M6l/A/4v7HkRTfoMsEae+vhuatmekXpV7rrcmhAjLdaUWbamNrkp7N6fnDMQcRjJ
DA/9VBV8qjokGu2PEj+HQGZb52y1+/S+wmUbKlS/EkYMME9gEDuUBFhHlC6xbYg1
akcddicTFhNHtwNQ9GFliIaJzU1Mt7LumB03/Cy0A9PouNUhv7QhTmVhbCBELiBC
ZWNrZXIgPG5iZWNrZXJAZnJlZC5uZXQ+iFcEExECABcFAjZa9ZcDCwQDBRUDAgYB
AxYCAQIXgAAKCRCtdGDCLVoO090GAJsFFd/nUF315R0Snt97iV39JP/OTQCeNAaU
5MsmAJHGFcXXj9AkMRoguzu5AQ0ENlr1pRAEAKpFYKuYC++L4RuzreeuKObO15SS
LXgUo0A/q9Hm3VFQw/FaWShBilVKjw6C7lUFnajx0uzy3EhczjitdcHewXyOH/9T
1WyqtiJG9CJTRgQkA1vDSgLBqLQ8so4saOF0bT/66iaiBE9Rbl1yRvjJh5lIULJr
BG2WhHfh/xWl2KS/AAMFBACQ/DrlJe9ooOQAuuUFK8P1A1t4zN5u9gvVMLhpxnr+
KYFa4+GdP3939lRTb7smtVxh9gote66kTmH776sqx7Sn8/Vsx5DOEKpikTlQ9IPR
mXu8Oe9skh+rJcOrjSOH7fSsYqqH7O1GArw0l82bBwA6Xz86vWfyHj/Slo0YXxey
QohGBBgRAgAGBQI2WvWlAAoJEK10YMItWg7TDiEAn3kIiU3z9lbtF4UexjL8zWIv
QszbAJ4om+wo1penO8/y9uI0UOgJQZUtJg==
=Q5Ab
-----END PGP PUBLIC KEY BLOCK-----

Re: Change definition of complex::norm

[Harvey J. Stein]

nbecker@fred.net writes:

 > We talked about this before, but it appears that in 3.0.2 this still
 > isn't changed.
 > 
 > There are clear disadvantages to the current definition of
 > complex::norm in terms of abs, not the least of which is that it won't work for
 > complex<int>, for example.   I can't imagine any drawback to defining
 > norm as
 > 
 > real(z) * real(z) + imag(z) * imag(z).

What's the current definition?

-- 
Harvey Stein
Bloomberg LP
hjstein@bloomberg.com

Re: Change definition of complex::norm

[nbecker]

  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
-----BEGIN PGP PUBLIC KEY BLOCK-----
Version: GnuPG v1.0.6 (GNU/Linux)
Comment: Processed by Mailcrypt 3.5.6 and Gnu Privacy Guard < http://www.gnupg.org/ >

mQGiBDZa9ZcRBACYGMoAmHIBUR19lDLZNJhgxGtqVchV7OiwniGIE0UpwRj08fDX
/KO7/cXgXDZqFEgHF98e6Gbm4efyyC7seP4Ye8Av3n8h8PMv307lQieJd5qQVvwx
vwJWGHsX1EOsv/Suzb2ZcYllU4dgrdBIkRLQ5tsJPiWtxjsfBsONGqWmIwCghmQA
GayzNTFUUy0JkGP8SEJRycED/0GvchcxrSnN0FDe5HqM2YzNOnQYGEasAgRSNoG7
O87uudA3j+Hh4GQSD7VgleYArCXqfaNd8pj+EY0ykGjcTJk07aAl+Ib8UrQ8eNk/
RON0+/ZRN6QGqte1lokR969AgVFDQaHV0IctElZdpRg+JbKUiBn3iYaY7xfYYr1z
M6l/A/4v7HkRTfoMsEae+vhuatmekXpV7rrcmhAjLdaUWbamNrkp7N6fnDMQcRjJ
DA/9VBV8qjokGu2PEj+HQGZb52y1+/S+wmUbKlS/EkYMME9gEDuUBFhHlC6xbYg1
akcddicTFhNHtwNQ9GFliIaJzU1Mt7LumB03/Cy0A9PouNUhv7QhTmVhbCBELiBC
ZWNrZXIgPG5iZWNrZXJAZnJlZC5uZXQ+iFcEExECABcFAjZa9ZcDCwQDBRUDAgYB
AxYCAQIXgAAKCRCtdGDCLVoO090GAJsFFd/nUF315R0Snt97iV39JP/OTQCeNAaU
5MsmAJHGFcXXj9AkMRoguzu5AQ0ENlr1pRAEAKpFYKuYC++L4RuzreeuKObO15SS
LXgUo0A/q9Hm3VFQw/FaWShBilVKjw6C7lUFnajx0uzy3EhczjitdcHewXyOH/9T
1WyqtiJG9CJTRgQkA1vDSgLBqLQ8so4saOF0bT/66iaiBE9Rbl1yRvjJh5lIULJr
BG2WhHfh/xWl2KS/AAMFBACQ/DrlJe9ooOQAuuUFK8P1A1t4zN5u9gvVMLhpxnr+
KYFa4+GdP3939lRTb7smtVxh9gote66kTmH776sqx7Sn8/Vsx5DOEKpikTlQ9IPR
mXu8Oe9skh+rJcOrjSOH7fSsYqqH7O1GArw0l82bBwA6Xz86vWfyHj/Slo0YXxey
QohGBBgRAgAGBQI2WvWlAAoJEK10YMItWg7TDiEAn3kIiU3z9lbtF4UexjL8zWIv
QszbAJ4om+wo1penO8/y9uI0UOgJQZUtJg==
=Q5Ab
-----END PGP PUBLIC KEY BLOCK-----

Re: Change definition of complex::norm

[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

Re: Change definition of complex::norm

[Gabriel Dos Reis]

nbecker@fred.net writes:

| We talked about this before, but it appears that in 3.0.2 this still
| isn't changed.

Yes, we discussed it and it didn't changed because I didn't took the
time to make the change. Sorry.

-- Gaby
CodeSourcery, LLC                       http://www.codesourcery.com

Re: Change definition of complex::norm

[nbecker]

>>>>> "Gabriel" == Gabriel Dos Reis <gdr@codesourcery.com> writes:

    Gabriel> nbecker@fred.net writes:
    Gabriel> | We talked about this before, but it appears that in 3.0.2 this still
    Gabriel> | isn't changed.

    Gabriel> Yes, we discussed it and it didn't changed because I didn't took the
    Gabriel> time to make the change. Sorry.

OK, so do I need to take any further action, or will you generate the
patch?  I don't currently have a copy of the source, so I'd prefer
someone else make the patch.

-----BEGIN PGP PUBLIC KEY BLOCK-----
Version: GnuPG v1.0.6 (GNU/Linux)
Comment: Processed by Mailcrypt 3.5.6 and Gnu Privacy Guard < http://www.gnupg.org/ >

mQGiBDZa9ZcRBACYGMoAmHIBUR19lDLZNJhgxGtqVchV7OiwniGIE0UpwRj08fDX
/KO7/cXgXDZqFEgHF98e6Gbm4efyyC7seP4Ye8Av3n8h8PMv307lQieJd5qQVvwx
vwJWGHsX1EOsv/Suzb2ZcYllU4dgrdBIkRLQ5tsJPiWtxjsfBsONGqWmIwCghmQA
GayzNTFUUy0JkGP8SEJRycED/0GvchcxrSnN0FDe5HqM2YzNOnQYGEasAgRSNoG7
O87uudA3j+Hh4GQSD7VgleYArCXqfaNd8pj+EY0ykGjcTJk07aAl+Ib8UrQ8eNk/
RON0+/ZRN6QGqte1lokR969AgVFDQaHV0IctElZdpRg+JbKUiBn3iYaY7xfYYr1z
M6l/A/4v7HkRTfoMsEae+vhuatmekXpV7rrcmhAjLdaUWbamNrkp7N6fnDMQcRjJ
DA/9VBV8qjokGu2PEj+HQGZb52y1+/S+wmUbKlS/EkYMME9gEDuUBFhHlC6xbYg1
akcddicTFhNHtwNQ9GFliIaJzU1Mt7LumB03/Cy0A9PouNUhv7QhTmVhbCBELiBC
ZWNrZXIgPG5iZWNrZXJAZnJlZC5uZXQ+iFcEExECABcFAjZa9ZcDCwQDBRUDAgYB
AxYCAQIXgAAKCRCtdGDCLVoO090GAJsFFd/nUF315R0Snt97iV39JP/OTQCeNAaU
5MsmAJHGFcXXj9AkMRoguzu5AQ0ENlr1pRAEAKpFYKuYC++L4RuzreeuKObO15SS
LXgUo0A/q9Hm3VFQw/FaWShBilVKjw6C7lUFnajx0uzy3EhczjitdcHewXyOH/9T
1WyqtiJG9CJTRgQkA1vDSgLBqLQ8so4saOF0bT/66iaiBE9Rbl1yRvjJh5lIULJr
BG2WhHfh/xWl2KS/AAMFBACQ/DrlJe9ooOQAuuUFK8P1A1t4zN5u9gvVMLhpxnr+
KYFa4+GdP3939lRTb7smtVxh9gote66kTmH776sqx7Sn8/Vsx5DOEKpikTlQ9IPR
mXu8Oe9skh+rJcOrjSOH7fSsYqqH7O1GArw0l82bBwA6Xz86vWfyHj/Slo0YXxey
QohGBBgRAgAGBQI2WvWlAAoJEK10YMItWg7TDiEAn3kIiU3z9lbtF4UexjL8zWIv
QszbAJ4om+wo1penO8/y9uI0UOgJQZUtJg==
=Q5Ab
-----END PGP PUBLIC KEY BLOCK-----

Re: Change definition of complex::norm

[Gabriel Dos Reis]

HJSTEIN@bloomberg.com (Harvey J. Stein) writes:

[...]

| 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.

Oh, we can certainly make a PhD dissertation about what is the best
thing in the world.

|  Dividing
| by __s reduces accuracy.  Computing __s can overflow when abs(__z)
| wouldn't, 

Really?

| and dividing by it can cause the computation of __x*__x to
| underflow.

But that is harmless (see f.e. by G.W Stewart in Matrix Algorithms, vol 1).

[...]

| 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.

G.W Stewart reported an interesting real world experience (Matrix
Algorithms, vol 1, page 144):

   If the scale factor in the algorithm Euclid [the above mentionned
   algorithm] is replaced by max{|a|, |b|}, the results may be less
   accurate on a hexadecimal machine.  The reason is that the number

       \sqrt{(a/s)^2 + (b/s)^2}

   is a little bit greater than one so tht the leading three bits in its
   representation are zero.  I discovered this fact after two days of
   trying to figure out why an algorithm I had coded consistently
   returned answers about a decimal digit less accurate than the
   algorithm it was meant to replace.  Such are the minutiae of
   computer arithmetic.


[...]

| Using frexp, ldexp & modf

Those are not acceptable in the generic algorithm abs() or norm..

-- Gaby

Re: Change definition of complex::norm

[Gabriel Dos Reis]

Benjamin Kosnik <bkoz@redhat.com> writes:

| Thanks for this mail: it was very clear and detailed. 
| 
| > > 2 The effect of instantiating the template complex for any type other than 
| > > float, double or long double is unspecified.
| > 
| > Point (2) seems to turn the issue of implementation of <complex> templates
| > and operations for, e.g., int or long, into a QOI issue.  
| 
| More than that. The standard specifies float, double, and long double 
| specializations, so there is definite room for optimizations for floating 
| point types, which will be the most used anyway.
| 
| For user-defined types, I think the generic complex template will fall 
| down. In this case, I think the smart thing to do is allow the code to 
| compile, but have the library get out of the way by then having the used 
| member functions in the given translation unit be undefined at 
| link time. This allows users to define their own specializations, if they 
| really want to do this.

They already can do so with the current sitiuation.

[...]

| Great. The libstdc++-v3 numerics testsuite is pretty anemic at this 
| point: perhaps this could be added?

Yes, but the qanswer may vary from one machine to another.  We'll need
to add more machinery which is most of the important point why the
numeric testsuite is so anemic.

[...]

| > 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.
| 
| Good question. I doubt there is any precision in this number, but who 
| knows.

How can one know?

| Physicists? Gaby? 

In doubt, I prefer to be conservative and deliver the best that be
computed within the machine precision, i.e. retain the current
definition. 

-- Gaby

Re: Change definition of complex::norm

[Gabriel Dos Reis]

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

| I've done a bit of investigation, and can offer some information 
| and suggestions, but someone else will need to generate the patch,
| as I know very little about this part of C++.
| 
| First, I find the following statement in the standard:
| 
| > 26.2 Complex numbers [lib.complex.numbers]
| > ...
| > 2 The effect of instantiating the template complex for any type other than 
| > float, double or long double is unspecified.
| 
| Point (2) seems to turn the issue of implementation of <complex> templates
| and operations for, e.g., int or long, into a QOI issue. 

Yes, this was already discussed.  And I agreed that we should make it
possible for type other than floa, double and long double to be used
where possible.

[...]

| The results are the same on sparc-sun-solaris28 and i686-pc-linux-gnu:
| 
| [lucier@curie ~]$ ./test
| abs_1: 1.694881e+308
| abs_2: nan
| norm_1: 0.000000e+00
| norm_2: 4.940656e-324
| 
| abs_2 uses the algorithm in std_complex.h; abs_1 uses the max of
| the absolute values of x and y as the divisor.  Getting a NaN
| when the proper answer is finite is not a "Good Thing", so I think
| that the definition in std_complex.h should be changed.

OK.

| 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.

-- Gaby

Re: Change definition of complex::norm

[Brad Lucier]

> > 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.
> 
> Good question. I doubt there is any precision in this number, but who 
> knows. Physicists? Gaby? 

The input values were chosen so that the "correct" answer was in the low
denormal range.  In this range, relative precision makes little sense,
but the computed answer is close to the best answer in absolute precision.

Brad

Re: Change definition of complex::norm

[Benjamin Kosnik]

Thanks for this mail: it was very clear and detailed. 

> > 2 The effect of instantiating the template complex for any type other than 
> > float, double or long double is unspecified.
> 
> Point (2) seems to turn the issue of implementation of <complex> templates
> and operations for, e.g., int or long, into a QOI issue.  

More than that. The standard specifies float, double, and long double 
specializations, so there is definite room for optimizations for floating 
point types, which will be the most used anyway.

For user-defined types, I think the generic complex template will fall 
down. In this case, I think the smart thing to do is allow the code to 
compile, but have the library get out of the way by then having the used 
member functions in the given translation unit be undefined at 
link time. This allows users to define their own specializations, if they 
really want to do this.

(Think about what happens if generic definitions are provided, and are 
wrong for the user-defined type. The user has to resort to link trickery.)

This is related to my point about facets, from my checkin earlier today.

I'm curious as to what other C++ users think is the best policy. Should 
these be treated on a case-by-case basis, or should the entire library 
try to conform to some general "instantion policy?"

Hmmmmm.

> Here's a C test program compiled on x86 and sparc with

Great. The libstdc++-v3 numerics testsuite is pretty anemic at this 
point: perhaps this could be added?

> The results are the same on sparc-sun-solaris28 and i686-pc-linux-gnu:
> 
> [lucier@curie ~]$ ./test
> abs_1: 1.694881e+308
> abs_2: nan
> norm_1: 0.000000e+00
> norm_2: 4.940656e-324
> 
> abs_2 uses the algorithm in std_complex.h; abs_1 uses the max of
> the absolute values of x and y as the divisor.  Getting a NaN
> when the proper answer is finite is not a "Good Thing", so I think
> that the definition in std_complex.h should be changed.

Seems like pretty concrete proof to me.

> 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.

Good question. I doubt there is any precision in this number, but who 
knows. Physicists? Gaby? 

-benjamin

Re: Change definition of complex::norm

[Brad Lucier]

I've done a bit of investigation, and can offer some information 
and suggestions, but someone else will need to generate the patch,
as I know very little about this part of C++.

First, I find the following statement in the standard:

> 26.2 Complex numbers [lib.complex.numbers]
> ...
> 2 The effect of instantiating the template complex for any type other than 
> float, double or long double is unspecified.

Point (2) seems to turn the issue of implementation of <complex> templates
and operations for, e.g., int or long, into a QOI issue.  So, the point of
nbecker@fred.net:

 > 1) the current definition is not as useful, since it only works where
 > sqrt is defined

is true, but is perhaps not relevant.

Here's a C test program compiled on x86 and sparc with

gcc -W -Wall -O2 -o test test.c -lm -ffloat-store

#include <math.h>
#include <stdio.h>

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);
}

double
abs_2(double x, double y) {

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

  if (s == 0.)
    return s;

  x /= s;
  y /= s;

  return s * sqrt( x * x + y * y);
}

double
norm_1 (double x, double y) {

  double x_sq = x * x;
  double y_sq = y * y;

  return x_sq + y_sq;

}

double
norm_2 (double x, double y) {

  double s = abs_1 (x, y);

  return s * s;

}



int main() {

  printf ("abs_1: %e\n", abs_1 (1.79769313486231570e+308/1.5, 1.79769313486231570e+308/1.5));

  printf ("abs_2: %e\n", abs_2 (1.79769313486231570e+308/1.5, 1.79769313486231570e+308/1.5));

  printf ("norm_1: %e\n", norm_1 (sqrt(4.940656458412465441765687928682e-324)/2.0, sqrt(4.940656458412465441765687928682e-324)/2.0));

  printf ("norm_2: %e\n", norm_2 (sqrt(4.940656458412465441765687928682e-324)/2.0, sqrt(4.940656458412465441765687928682e-324)/2.0));
 
  return 0;
}

The results are the same on sparc-sun-solaris28 and i686-pc-linux-gnu:

[lucier@curie ~]$ ./test
abs_1: 1.694881e+308
abs_2: nan
norm_1: 0.000000e+00
norm_2: 4.940656e-324

abs_2 uses the algorithm in std_complex.h; abs_1 uses the max of
the absolute values of x and y as the divisor.  Getting a NaN
when the proper answer is finite is not a "Good Thing", so I think
that the definition in std_complex.h should be changed.

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.

Brad Lucier

Re: Change definition of complex::norm

[dewar]

<< If the scale factor in the algorithm Euclid [the above mentionned
 algorithm] is replaced by max{|a|, |b|}, the results may be less
 accurate on a hexadecimal machine.  The reason is that the number
>>

Yes, well hex machines are a numerical abomination, and this is a good
example of why :-)

Re: Change definition of complex::norm

[lucier]

Some comments on various parts of this thread, with no attempt
to get the attributions right.

Re:

> In fact, it would be better to do it the other way around -
> define abs(x) as norm_factor*norm(z/norm_factor)^2.

Surely, this formula is a typo/thinko.

> |  Dividing
> | by __s reduces accuracy.  Computing __s can overflow when abs(__z)
> | wouldn't, 
> 
> Really?

Yes, if __x and __y are between MAX_DOUBLE/2.0 and sqrt(2.)*MAX_DOUBLE/2.0, then
abs(__z) does not overflow, but __s = abs(__x) + abs(__y) does.

> G.W Stewart reported an interesting real world experience (Matrix
> Algorithms, vol 1, page 144):
> 
>    If the scale factor in the algorithm Euclid [the above mentionned
>    algorithm] is replaced by max{|a|, |b|}, the results may be less
>    accurate on a hexadecimal machine.  The reason is that the number
> 
>        \sqrt{(a/s)^2 + (b/s)^2}
> 
>    is a little bit greater than one so tht the leading three bits in its
>    representation are zero.  I discovered this fact after two days of
>    trying to figure out why an algorithm I had coded consistently
>    returned answers about a decimal digit less accurate than the
>    algorithm it was meant to replace.  Such are the minutiae of
>    computer arithmetic.

Yes, non--base 2 arithmetics are fun.  20 years ago (how time flies)
I wrote an elementary function library for UCSD Pascal on a base 100
machine, the TI 990.  Floating-point format was 7 base-100 digits for
the fraction part in 7 bytes, one byte for the sign and the biased
exponent.  Thank god for IEEE 754.

Brad Lucier

Re: Change definition of complex::norm

[Benjamin Kosnik]

If there is a problem, let's fix it, and add a test case.

You'll need: 

1) ChangeLog entry
2) patch
3) testcase submission, probably to testsuite/26_numerics/complex_value.cc

In particular, please show in C++ code that:

- the current definition is not as useful, since it only works where
   sqrt is defined
- the current definition is much less efficient
- the current definition is possibly less accurate

thanks,
benjamin

Re: Change definition of complex::norm

[dewar]

Robert said
>> real(z) * real(z) + imag(z) * imag(z).
> It looks like this will generate incorrect overflow in some cases ...

but he was wrong :-)
I agree this is a reasonable suggestion (I don't know what bothered me, but
I see no objection).

Re: Change definition of complex::norm

[nbecker]

>>>>> "dewar" == dewar  <dewar@gnat.com> writes:

    >>> real(z) * real(z) + imag(z) * imag(z).
    dewar> It looks like this will generate incorrect overflow in some cases ...

I don't see how that is possible.  Can you give an example?
-----BEGIN PGP PUBLIC KEY BLOCK-----
Version: GnuPG v1.0.6 (GNU/Linux)
Comment: Processed by Mailcrypt 3.5.6 and Gnu Privacy Guard < http://www.gnupg.org/ >

mQGiBDZa9ZcRBACYGMoAmHIBUR19lDLZNJhgxGtqVchV7OiwniGIE0UpwRj08fDX
/KO7/cXgXDZqFEgHF98e6Gbm4efyyC7seP4Ye8Av3n8h8PMv307lQieJd5qQVvwx
vwJWGHsX1EOsv/Suzb2ZcYllU4dgrdBIkRLQ5tsJPiWtxjsfBsONGqWmIwCghmQA
GayzNTFUUy0JkGP8SEJRycED/0GvchcxrSnN0FDe5HqM2YzNOnQYGEasAgRSNoG7
O87uudA3j+Hh4GQSD7VgleYArCXqfaNd8pj+EY0ykGjcTJk07aAl+Ib8UrQ8eNk/
RON0+/ZRN6QGqte1lokR969AgVFDQaHV0IctElZdpRg+JbKUiBn3iYaY7xfYYr1z
M6l/A/4v7HkRTfoMsEae+vhuatmekXpV7rrcmhAjLdaUWbamNrkp7N6fnDMQcRjJ
DA/9VBV8qjokGu2PEj+HQGZb52y1+/S+wmUbKlS/EkYMME9gEDuUBFhHlC6xbYg1
akcddicTFhNHtwNQ9GFliIaJzU1Mt7LumB03/Cy0A9PouNUhv7QhTmVhbCBELiBC
ZWNrZXIgPG5iZWNrZXJAZnJlZC5uZXQ+iFcEExECABcFAjZa9ZcDCwQDBRUDAgYB
AxYCAQIXgAAKCRCtdGDCLVoO090GAJsFFd/nUF315R0Snt97iV39JP/OTQCeNAaU
5MsmAJHGFcXXj9AkMRoguzu5AQ0ENlr1pRAEAKpFYKuYC++L4RuzreeuKObO15SS
LXgUo0A/q9Hm3VFQw/FaWShBilVKjw6C7lUFnajx0uzy3EhczjitdcHewXyOH/9T
1WyqtiJG9CJTRgQkA1vDSgLBqLQ8so4saOF0bT/66iaiBE9Rbl1yRvjJh5lIULJr
BG2WhHfh/xWl2KS/AAMFBACQ/DrlJe9ooOQAuuUFK8P1A1t4zN5u9gvVMLhpxnr+
KYFa4+GdP3939lRTb7smtVxh9gote66kTmH776sqx7Sn8/Vsx5DOEKpikTlQ9IPR
mXu8Oe9skh+rJcOrjSOH7fSsYqqH7O1GArw0l82bBwA6Xz86vWfyHj/Slo0YXxey
QohGBBgRAgAGBQI2WvWlAAoJEK10YMItWg7TDiEAn3kIiU3z9lbtF4UexjL8zWIv
QszbAJ4om+wo1penO8/y9uI0UOgJQZUtJg==
=Q5Ab
-----END PGP PUBLIC KEY BLOCK-----

Re: Change definition of complex::norm

[dewar]

>>real(z) * real(z) + imag(z) * imag(z).

It looks like this will generate incorrect overflow in some cases ...