Matthieu Schaller wrote:
Dear all,
I have recently been playing a lot on numerical analysis problems containing complex numbers and functions. I have been using the SL std::complex<> library and could observe that some of the algorithms weren't optimized in cases where complex numbers have to be multiplied by pure imaginary numbers (i.e. x*sqrt(-1)). I have implemented a class to represent these numbers and overloaded all relevant operators in a more optimized fashion thanks to the knowledge that the real part of these complex numbers is 0. I could observe an important speed-up in some calculations. As an example, this piece of code (1D Schroedinger equation using finite differences):
for(size_t i(0); i<nbPoints; ++i) { Psi_new[i] = Psi_old[i] + ii * hbar * deltaT / (12.*m*deltaX*deltaX) * (-Psi_cur[i-2] + 16.*Psi_cur[i-1] - 30.*Psi_cur[i] + 16.*Psi_cur[i+1] - Psi_cur[i+2]) - 2. * ii * V(x[i]) / hbar * deltaT * Psi_cur[i];
}
is more than 2 times faster (g++ 4.6.1 core i7 2820QM) when the imaginary unit (constant) variable ii is defined to be
Imaginary<double> ii(1.);
than when it is defined as
std::complex<double> ii(0.,1.);
Some other expressions also benefit from this improvement and can run faster. Expressions including square roots or exponentials of pure imaginary numbers are improved a lot.
The range of application of such a class is quite narrow but it may be helpful to speed-up some simple numerical problems. The source code as well as a working example is available here: http://code.google.com/p/cpp-imaginary-numbers/
Hi, I'm sure that if the performances can be improved, which seems logic, the Boost community will be interested by your library. Some comments related to the interface. I will replace T& imag(); const T& imag() const; by T imag() const; As there is no reason to provide non-const access to the internal representation. If these operations are needed by the non-member operators, maybe you can declare them private and move non-member operators as members, or declare them friends. Comparing imaginary and reals seems confusing to me ///Returns true if x equals y template<typename T> inline bool operator==(const Imaginary<T>& x, const T& y); template<typename T> inline bool operator==(const T& x, const Imaginary<T>& y); The following divide assign operator is missing: Imaginary<T>& operator/=(const Imaginary<T>& rhs) I would also expect an implicit conversion to complex as an imaginary is also a complex. and a kind of imaginary downcast from complex. template <typename T> imaginary<T> imaginary_cast(complex<T> const & rhs); My advice would be to boostify the code, write the documentation including as much performances tests as you consider could help to consider this library a should have. Good luck, Vicente -- View this message in context: http://boost.2283326.n4.nabble.com/Determining-interest-Pure-imaginary-numbe... Sent from the Boost - Dev mailing list archive at Nabble.com.