complex number policy for beast

[Tim Janik]

hi Stefan.

after the recent hassle we had with varying types of complex numbers:
- C++ complex numbers from <complex>
- C99 complex numbers from <complex.h> in C mode (including complex.h
   breaks for C++ sources)
- our own old comlex numbers GslComplex
- our own partially ported numbers BseComplex
- ellf/bsefilter-ellf.c own complex numbers EllfComplex

and the resulting confusion in source files about which complex types
and statements to use, i've decided that we'll follow this policy:

- <complex> and <complex.h> are never included in header files
   (this keeps our headers C++ and C conform)
- C sources *only* use BseComplex and assorted functions
- C++ source files can use BseComplex but may opt to include <complex>
   and use that. (they have to convert back-and-forth between std::complex
   and BseComplex for API though).

comitted as:

Sat Nov  4 17:38:28 2006  Tim Janik  <timj@...>

         * bsemath.h: don't include neither of the C++ <complex> or the C99
         <complex.h> header files. removed compat typedefs. we'll simply resort
         to using BseComplex in headers and use <complex> in C++ files.


---
ciaoTJ
_______________________________________________
beast mailing list
beast@...
http://mail.gnome.org/mailman/listinfo/beast

[Tim Janik]

On Sat, 4 Nov 2006, Tim Janik wrote:

> - C++ source files can use BseComplex but may opt to include <complex>
>   and use that. (they have to convert back-and-forth between std::complex
>   and BseComplex for API though).

after reading some more C++ headers, i have to add that using C++ complex
numbers can not be recommended from an accuracy point of view, from
/usr/include/c++/3.4/complex:

   // 26.2.5/13
   // XXX: This is a grammar school implementation.
   template<typename _Tp>
     template<typename _Up>
     complex<_Tp>&
     complex<_Tp>::operator*=(const complex<_Up>& __z)
     {
       const _Tp __r = _M_real * __z.real() - _M_imag * __z.imag();
       _M_imag = _M_real * __z.imag() + _M_imag * __z.real();
       _M_real = __r;
       return *this;
     }

this has many cancellation/etc problems, described e.g. in Numerical Recipes
and is usually worked around by in math packages, or e.g. the _Cdiv() from the
C9x specification.
also our bse_complex_div() from bsemath.[hc]* has none of these problems
(similar things hold for our other complex number functions), so if accuracy
matters at all, we better use BseComplex over std::complex<double>.

---
ciaoTJ
_______________________________________________
beast mailing list
beast@...
http://mail.gnome.org/mailman/listinfo/beast

[Stefan Westerfeld]

   Hi!

On Sat, Nov 04, 2006 at 06:41:58PM +0100, Tim Janik wrote:

> On Sat, 4 Nov 2006, Tim Janik wrote:
>
> >- C++ source files can use BseComplex but may opt to include <complex>
> >  and use that. (they have to convert back-and-forth between std::complex
> >  and BseComplex for API though).
>
> after reading some more C++ headers, i have to add that using C++ complex
> numbers can not be recommended from an accuracy point of view, from
> /usr/include/c++/3.4/complex:
>
>   // 26.2.5/13
>   // XXX: This is a grammar school implementation.
>   template<typename _Tp>
>     template<typename _Up>
>     complex<_Tp>&
>     complex<_Tp>::operator*=(const complex<_Up>& __z)
>     {
>       const _Tp __r = _M_real * __z.real() - _M_imag * __z.imag();
>       _M_imag = _M_real * __z.imag() + _M_imag * __z.real();
>       _M_real = __r;
>       return *this;
>     }
>
> this has many cancellation/etc problems, described e.g. in Numerical Recipes
> and is usually worked around by in math packages, or e.g. the _Cdiv() from
> the
> C9x specification.
> also our bse_complex_div() from bsemath.[hc]* has none of these problems
> (similar things hold for our other complex number functions), so if accuracy
> matters at all, we better use BseComplex over std::complex<double>.

This problem has been brought up on the gcc mailing lists, here is a
link:

   http://gcc.gnu.org/ml/gcc/2003-04/msg00814.html

The relevant information is that while this is the generic
implementation, there are template specializations for the
std::complex<float>, std::complex<double> and std::complex<long double>
case, which use compiler builtin __complex__ types / operations. So
as I understand it there should be not more numeric problems than in
gccs C9x implementation.

   Cu... Stefan
--
Stefan Westerfeld, Hamburg/Germany, http://space.twc.de/~stefan
_______________________________________________
beast mailing list
beast@...
http://mail.gnome.org/mailman/listinfo/beast