Re: Change definition of complex::norm

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

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