C-XSC - A C++ Class Library for Extended Scientific Computing  2.5.4
lx_complex.cpp
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