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*: nbecker at fred dot net*Subject*: Re: Change definition of complex::norm*From*: HJSTEIN at bloomberg dot com (Harvey J. Stein)*Date*: 31 Oct 2001 11:07:51 -0500*Cc*: gcc at gcc dot gnu dot org*References*: <x88snbzncba.fsf@adglinux1.hns.com> <kiwsnbzamtv.fsf@blinky.bloomberg.com> <x88wv1blv5r.fsf@adglinux1.hns.com>

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

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

**References**:**Change definition of complex::norm***From:*nbecker

**Re: Change definition of complex::norm***From:*Harvey J. Stein

**Re: Change definition of complex::norm***From:*nbecker

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

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