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

Index Nav: | [Date Index] [Subject Index] [Author Index] [Thread Index] | |
---|---|---|

Message Nav: | [Date Prev] [Date Next] | [Thread Prev] [Thread Next] |

*To*: gcc at gcc dot gnu dot org*Subject*: Re: Change definition of complex::norm*From*: Brad Lucier <lucier at math dot purdue dot edu>*Date*: Wed, 31 Oct 2001 21:45:37 -0500 (EST)*Cc*: hjstein at bloomberg dot com, lucier at math dot purdue dot edu, nbecker at fred dot net, gdr at codesourcery dot com, bkoz at redhat dot com

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

**Follow-Ups**:**Re: Change definition of complex::norm***From:*Benjamin Kosnik

**Re: Change definition of complex::norm***From:*Gabriel Dos Reis

Index Nav: | [Date Index] [Subject Index] [Author Index] [Thread Index] | |
---|---|---|

Message Nav: | [Date Prev] [Date Next] | [Thread Prev] [Thread Next] |