Re: [boost] Determining interest: Pure imaginary number library

Subject: Re: [boost] Determining interest: Pure imaginary number library
From: Matthias Schabel (boost_at_[hidden])
Date: 2011-11-29 12:09:02

>>> I also think that it will be worth defining a boost::complex class on the style of the rejected standard proposal that will integrate better with the imaginary class.
>>> The standard complex class has a constraint that a boost::complex class could avoid, it only accepts the builtin double and float types. This complex class could accept any type conforming to the expected Concept. I'm sure that others are expecting a complex class that can be used with specific classes, as arbitrary precision, ...
> This can be done. The boost::complex class could simply be an implementation of the n1869 draft by Thorsten Ottosen then:
> 
> http://www.open-std.org/JTC1/SC22/WG21/docs/papers/2005/n1869.html
> 
> With the addition of the TR1 reciprocal trigonometric functions and the small updates brought by C++11.
> On top of this, some mathematical operators could be specialized for float and double if some performance improvements can be obtained in those cases.
> I will work on this.
>> 
>> Before we all get too carried away, there is an issue with "mycomplex<my_big_num>" which is that if "my_big_num" implements expression templates to optimize arithmetic operations, then these expression template end up in the wrong place with "mycomplex<my_big_num>" : which is to say inside mycomplex's implementation, instead of on mycomplex's own external operators (which is where you really need the expression templates given that each mycomplex temporary will contain 2 values).
> 
> I may be wrong but I don't see how complex numbers could benefit from template expressions. The two "components" get mixed in almost every expression and cannot be decoupled. The complex<> class always ends up using the operators of its template argument which in the case of a GMP-like number means that the appropriate template operators are called.
> The complex number operators always end up using real number arithmetic. So mycomplex<my_big_num> would use your "improved" operators.
> 
>> 
>> And finally, there's already a small amount of complex number code in Boost under Boost.Math (just adds some extra TR1 functions), so if adding imaginary support is small enough, we could adopt it under Boost.Math I guess - or vice versa if the new code is large enough to warrant having it's own library.  Just thinking out loud here that it would be better to have all the complex number code in one place if possible (and it makes sense to do so).
>> 
> The TR1 operators could be moved into (next to) this new complex<> class in order to group everything. I would anyway have to implement them if I want the class to be as close as possible to the SL-one as these mathematical functions are now part of the standards.
> 
> Thanks for your input.

We over in Boost.Units land are also interested in a boost::complex implementation that plays better with units. Here's a snippet of code from Janek Kozicki that doesn't quite work the way it should:

Hi, I wrote a short function that integrates Huyghens spherical wave
for diffraction. (This is a single slit experiment btw).

I have some questions, marked with (*)

double integral(
          const quantity<length> x	// x - position where light intensity is calculated
        , const quantity<length> R	// R - distance from slit to screen
        , const quantity<length> D	// D - slit width
        , const quantity<length> lambda	// lambda - wavelength
        , double steps			// steps - number of integration steps // (*)
)
{
   complex<double> i(0,1);
   complex<double> a(0);
   static const double pi = 3.14159265358979323846;        // (**)
   quantity<wavenumber,complex<double> > k=2*pi/lambda;
   quantity<length> d=-0.5*D;   // d - position on the slit
   quantity<length> dd=D/steps; // dd - width of integral step
   // integrate from -d to d: exp(i*k*r)/r
   for( ; d<= 0.5*D ; d+=dd )
   {
       quantity<length> r=sqrt(pow<2>(R)+pow<2>(x-d));
       a += (dd*exp(i*k*r)/r).value();                    // (***) -> shouldn't need .value() here...
   }
   return pow<2>(abs(a));
}

We have a skeleton of a complex class that partially works with units in lib/units/examples/complex.cpp, but it is not fleshed out. Seems like it would be silly to have multiple competing implementations...

Matthias