html/doxygen/complex__double_8hpp_source.html

00001
00002 /******************************************************************************
00003 * MODULE     : complex_double.hpp
00004 * DESCRIPTION: Complexified doubles
00005 * COPYRIGHT  : (C) 2008  Joris van der Hoeven
00006 *******************************************************************************
00007 * This software falls under the GNU general public license and comes WITHOUT
00008 * ANY WARRANTY WHATSOEVER. See the file $TEXMACS_PATH/LICENSE for more details.
00009 * If you don't have this file, write to the Free Software Foundation, Inc.,
00010 * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
00011 ******************************************************************************/
00012
00013 #ifndef __MMX_COMPLEX_DOUBLE_HPP
00014 #define __MMX_COMPLEX_DOUBLE_HPP
00015 #include <basix/double.hpp>
00016 #include <numerix/complex.hpp>
00017 #include <numerix/rounded.hpp>
00018 namespace mmx {
00019
00020 /******************************************************************************
00021 * Avoid overflows for complex<double>
00022 ******************************************************************************/
00023
00024 STMPL inline complex<double>
00025 invert (const complex<double>& z) {
00026   if (abs (Re (z)) >= abs (Im (z))) {
00027     double f= Im (z) / Re (z);
00028     double c= 1.0 / (Re (z) + Im (z) * f);
00029     return complex<double> (c, -c * f);
00030   }
00031   else {
00032     double f= Re (z) / Im (z);
00033     double c= 1.0 / (Re (z) * f + Im (z));
00034     return complex<double> (c * f, -c);
00035   }
00036 }
00037
00038 STMPL inline complex<double>
00039 operator / (const complex<double>& z1, const complex<double>& z2) {
00040   if (abs (Re (z2)) >= abs (Im (z2))) {
00041     double f= Im (z2) / Re (z2);
00042     double c= 1.0 / (Re (z2) + Im (z2) * f);
00043     return complex<double> (c * (Re (z1) + Im (z1) * f),
00044                             c * (Im (z1) - Re (z1) * f));
00045   }
00046   else {
00047     double f= Re (z2) / Im (z2);
00048     double c= 1.0 / (Re (z2) * f + Im (z2));
00049     return complex<double> (c * (Re (z1) * f + Im (z1)),
00050                             c * (Im (z1) * f - Re (z1)));
00051   }
00052 }
00053
00054 // TODO: abs, sqrt and in-place operations
00055
00056 /******************************************************************************
00057 * Floating point related operations
00058 ******************************************************************************/
00059
00060 inline double additive_error (const complex<double>& z) {
00061   return ldexp (abs (Re (z)) + abs (Im (z)), -49); }
00062
00063 } // namespace mmx
00064 #endif // __MMX_COMPLEX_DOUBLE_HPP