00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
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
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
00055
00056
00057
00058
00059
00060 inline double additive_error (const complex<double>& z) {
00061 return ldexp (abs (Re (z)) + abs (Im (z)), -49); }
00062
00063 }
00064 #endif // __MMX_COMPLEX_DOUBLE_HPP