C-XSC - A C++ Class Library for Extended Scientific Computing
2.5.4
|
00001 /* 00002 ** CXSC is a C++ library for eXtended Scientific Computing (V 2.5.4) 00003 ** 00004 ** Copyright (C) 1990-2000 Institut fuer Angewandte Mathematik, 00005 ** Universitaet Karlsruhe, Germany 00006 ** (C) 2000-2014 Wiss. Rechnen/Softwaretechnologie 00007 ** Universitaet Wuppertal, Germany 00008 ** 00009 ** This library is free software; you can redistribute it and/or 00010 ** modify it under the terms of the GNU Library General Public 00011 ** License as published by the Free Software Foundation; either 00012 ** version 2 of the License, or (at your option) any later version. 00013 ** 00014 ** This library is distributed in the hope that it will be useful, 00015 ** but WITHOUT ANY WARRANTY; without even the implied warranty of 00016 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00017 ** Library General Public License for more details. 00018 ** 00019 ** You should have received a copy of the GNU Library General Public 00020 ** License along with this library; if not, write to the Free 00021 ** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00022 */ 00023 00024 /* CVS $Id: lx_complex.cpp,v 1.9 2014/01/30 17:23:47 cxsc Exp $ */ 00025 00026 /* 00027 ** F. Blomquist, University of Wuppertal, 19.09.2007; 00028 */ 00029 00030 #include "lx_complex.hpp" 00031 #include "lx_cinterval.hpp" 00032 00033 namespace cxsc { 00034 // ---------------------------------------------------------------------------- 00035 // ----- Functions and operators related to type lx_complex ------------------- 00036 // ---------------------------------------------------------------------------- 00037 00038 l_complex & l_complex::operator = (const lx_complex& a) throw() 00039 { 00040 l_real u,v; 00041 u = Re(a); v = Im(a); 00042 return *this = l_complex(u,v); 00043 } 00044 00045 complex & complex::operator = (const lx_complex& a) throw() 00046 { 00047 l_complex x; 00048 complex y; 00049 00050 x = a; 00051 y = x; 00052 return *this = y; 00053 } 00054 00055 std::string & operator >> (std::string& s, lx_complex& a) throw() 00056 // Writes string s to variable a of type lx_complex; 00057 // and returns an empty string s; 00058 // Example: s = "({2,3} , {4,5})" delivers a value a 00059 // with: a ~ 10^2 * 3.0 + i*10^4 * 5.0; 00060 { 00061 string su; 00062 int i; 00063 std::cout << "Halo 1" << std::endl; 00064 s = skipwhitespacessinglechar (s, '('); 00065 std::cout << "s = " << s << std::endl; 00066 i = s.find("}"); 00067 std::cout << "i = " << i << std::endl; 00068 su = s.substr(0,i+1); 00069 std::cout << "su = " << su << std::endl; 00070 su >> a.re; 00071 00072 s.erase(0,i+1); 00073 s = skipwhitespacessinglechar (s, ','); 00074 std::cout << "s = " << s << std::endl; 00075 s >> a.im; 00076 00077 s = ""; 00078 return s; 00079 } 00080 00081 void operator >> (const std::string &s, lx_complex &a) throw() 00082 { 00083 // Writes string s to variable a of type lx_real; 00084 std::string r(s); 00085 r >> a; 00086 } 00087 00088 void operator >> (const char *s, lx_complex& a) throw() 00089 { 00090 std::string r(s); 00091 r >> a; 00092 } 00093 00094 lx_real abs (const lx_complex& a) throw() 00095 { return sqrtx2y2(a.re,a.im); } 00096 lx_real abs2 (const lx_complex& a) throw() 00097 { return a.re*a.re + a.im*a.im; } 00098 00099 // ----------------------------------------------------------------------------- 00100 // --------- Elementary functions related to type lx_complex ------------------- 00101 // ----------------------------------------------------------------------------- 00102 00103 lx_complex sqr(const lx_complex& z) throw() { return (z*z); } 00104 lx_complex sqrt(const lx_complex& z) throw() 00105 { return mid(sqrt(lx_cinterval(z))); } 00106 00107 lx_complex sqrt(const lx_complex& z,int n) throw() 00108 { return mid( sqrt(lx_cinterval(z),n) ); } 00109 00110 lx_complex exp(const lx_complex& z) throw() 00111 { return mid(exp(lx_cinterval(z))); } 00112 00113 lx_complex exp2(const lx_complex& z) throw() 00114 { return mid(exp2(lx_cinterval(z))); } 00115 00116 lx_complex exp10(const lx_complex& z) throw() 00117 { return mid(exp10(lx_cinterval(z))); } 00118 00119 lx_complex sin(const lx_complex& z) throw() 00120 { return mid(sin(lx_cinterval(z))); } 00121 00122 lx_complex cos(const lx_complex& z) throw() 00123 { return mid(cos(lx_cinterval(z))); } 00124 00125 lx_complex tan(const lx_complex& z) throw() 00126 { return mid(tan(lx_cinterval(z))); } 00127 00128 lx_complex cot(const lx_complex& z) throw() 00129 { return mid(cot(lx_cinterval(z))); } 00130 00131 lx_complex asin(const lx_complex& z) throw() 00132 { return mid(asin(lx_cinterval(z))); } 00133 00134 lx_complex acos(const lx_complex& z) throw() 00135 { return mid(acos(lx_cinterval(z))); } 00136 00137 lx_complex atan(const lx_complex& z) throw() 00138 { return mid(atan(lx_cinterval(z))); } 00139 00140 lx_complex acot(const lx_complex& z) throw() 00141 { return mid(acot(lx_cinterval(z))); } 00142 00143 lx_complex sinh(const lx_complex& z) throw() 00144 { return mid(sinh(lx_cinterval(z))); } 00145 00146 lx_complex cosh(const lx_complex& z) throw() 00147 { return mid(cosh(lx_cinterval(z))); } 00148 00149 lx_complex tanh(const lx_complex& z) throw() 00150 { return mid(tanh(lx_cinterval(z))); } 00151 00152 lx_complex coth(const lx_complex& z) throw() 00153 { return mid(coth(lx_cinterval(z))); } 00154 00155 lx_complex asinh(const lx_complex& z) throw() 00156 { return mid(asinh(lx_cinterval(z))); } 00157 00158 lx_complex acosh(const lx_complex& z) throw() 00159 { return mid(acosh(lx_cinterval(z))); } 00160 00161 lx_complex atanh(const lx_complex& z) throw() 00162 { return mid(atanh(lx_cinterval(z))); } 00163 00164 lx_complex acoth(const lx_complex& z) throw() 00165 { return mid(acoth(lx_cinterval(z))); } 00166 00167 // sqrt_all(c) computes a list of 2 values for all square roots of c 00168 std::list<lx_complex> sqrt_all( const lx_complex& c ) 00169 { 00170 lx_complex lc; 00171 lc = sqrt(c); 00172 00173 std::list<lx_complex> res; 00174 res.push_back( lc ); 00175 res.push_back( -lc ); 00176 00177 return res; 00178 } // end sqrt_all 00179 00180 lx_real arg(const lx_complex& z) throw() 00181 { return mid(arg(lx_cinterval(z))); } 00182 00183 lx_real Arg(const lx_complex& z) throw() 00184 { return mid(Arg(lx_cinterval(z))); } 00185 00186 std::list<lx_complex> sqrt_all( const lx_complex& z, int n ) 00187 // 00188 // sqrt_all(z,n) computes a list of n values approximating all n-th roots of z 00189 // 00190 // For n >=3, computing the optimal approximations is very expensive 00191 // and thus not considered cost-effective. 00192 // 00193 // Hence, the polar form is used to calculate sqrt_all(z,n) 00194 // 00195 { 00196 std::list<lx_complex> res; 00197 00198 if( n == 0 ) 00199 { 00200 res.push_back( lx_complex(l_real(1),l_real(0)) ); 00201 return res; 00202 } 00203 else if( n == 1 ) 00204 { 00205 res.push_back(z); 00206 return res; 00207 } 00208 else if( n == 2 ) return sqrt_all( z ); 00209 else 00210 { 00211 lx_real 00212 arg_z = arg( z ), root_abs_z = sqrt( abs( z ), n ); 00213 00214 for(int k = 0; k < n; k++) 00215 { 00216 lx_real arg_k = ( arg_z + 2 * k * Pi_lx_real() ) / n; 00217 res.push_back( lx_complex( root_abs_z * cos( arg_k ), 00218 root_abs_z * sin( arg_k ) ) ); 00219 } 00220 return res; 00221 } 00222 } 00223 // 00224 //-- end sqrt_all ------------------------------------------------------------- 00225 00226 lx_complex ln(const lx_complex& z) throw() 00227 { return mid(ln(lx_cinterval(z))); } 00228 00229 lx_complex log2(const lx_complex& z) throw() 00230 { return mid(log2(lx_cinterval(z))); } 00231 lx_complex log10(const lx_complex& z) throw() 00232 { return mid(log10(lx_cinterval(z))); } 00233 00234 lx_complex power_fast(const lx_complex& z, const real& n) throw() 00235 { 00236 if( n == 0 ) 00237 return lx_complex(l_real(1),l_real(0)); 00238 else 00239 if( n == 1 ) return z; 00240 else 00241 if( n == -1 ) return 1 / z; 00242 else 00243 if( n == 2 ) return sqr(z); 00244 else 00245 { 00246 lx_real abs_z = abs(z); 00247 if( ((n < 0) && (abs_z == 0.0)) || !(Is_Integer(n))) 00248 // z contains 0 00249 cxscthrow (STD_FKT_OUT_OF_DEF( 00250 "lx_complex power_fast(const lx_complex& z, const real& n ); z = 0 or n is not integer.")); 00251 if( abs_z == 0.0 ) 00252 return lx_complex(l_real(0),l_real(0)); 00253 else 00254 { 00255 lx_real arg_z = arg(z); 00256 lx_real abs_z_n = exp( n * ln( abs_z ) ); 00257 00258 return lx_complex( abs_z_n * cos( n * arg_z ), 00259 abs_z_n * sin( n * arg_z ) ); 00260 } 00261 } 00262 } // End: power_fast 00263 00264 lx_complex power(const lx_complex& x, const real& n) throw() 00265 { 00266 if( !(Is_Integer(n)) ) 00267 // n is not an integer 00268 cxscthrow(STD_FKT_OUT_OF_DEF( 00269 "lx_complex power(const lx_complex& z, const real& n); n is not integer.")); 00270 00271 real zhi(2.0), N(n), r; 00272 lx_complex one = lx_complex(l_real(1),l_real(0)); 00273 double dbl; 00274 lx_complex y, neu, X(x); 00275 00276 if (x == one) y = x; 00277 else 00278 if (N == 0.0) y = one; 00279 else 00280 { 00281 if (N == 1) y = x; 00282 else 00283 if (N == 2) y = sqr(x); 00284 else 00285 { 00286 if (N < 0) 00287 { 00288 X = 1.0/X; 00289 N = -N; 00290 } 00291 // Initialisierung 00292 if ( !Is_Integer(N/2) ) 00293 y = X; 00294 else y = one; 00295 neu = sqr(X); // neu = X*X; 00296 do 00297 { 00298 dbl = _double(N/zhi); 00299 dbl = floor(dbl); 00300 r = (real) dbl; 00301 if ( !Is_Integer( r/2 ) ) 00302 y *= neu; 00303 zhi += zhi; 00304 if (zhi <= N) 00305 neu = sqr(neu); // neu = neu * neu; 00306 } while (zhi <= N); 00307 } 00308 } 00309 return y; 00310 } // end power(z,n) 00311 00312 lx_complex pow(const lx_complex& z, const lx_real& p) throw() 00313 { return mid( pow( lx_cinterval(z) , lx_interval(p) ) ); } 00314 00315 lx_complex pow(const lx_complex& z, const lx_complex& p) throw() 00316 { return mid( pow( lx_cinterval(z) , lx_cinterval(p) ) ); } 00317 00318 lx_complex sqrt1px2(const lx_complex& z) throw() 00319 { return mid(sqrt1px2(lx_cinterval(z))); } 00320 00321 lx_complex sqrt1mx2(const lx_complex& z) throw() 00322 { return mid(sqrt1mx2(lx_cinterval(z))); } 00323 00324 lx_complex sqrtx2m1(const lx_complex& z) throw() 00325 { return mid(sqrtx2m1(lx_cinterval(z))); } 00326 00327 lx_complex sqrtp1m1(const lx_complex& z) throw() 00328 { return mid(sqrtp1m1(lx_cinterval(z))); } 00329 00330 lx_complex expm1(const lx_complex& z) throw() 00331 { return mid(expm1(lx_cinterval(z))); } 00332 00333 lx_complex lnp1(const lx_complex& z) throw() 00334 { return mid(lnp1(lx_cinterval(z))); } 00335 00336 } // namespace cxsc