C-XSC - A C++ Class Library for Extended Scientific Computing  2.5.4
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: complex.cpp,v 1.29 2014/01/30 17:23:44 cxsc Exp $ */
00025 /* New versions of product(...), quotient(...); Blomquist, 26.10.02; */
00026 
00027 #include "complex.hpp"
00028 #include "dot.hpp"
00029 #include "cdot.hpp"
00030 #include "real.hpp"
00031 #include "rmath.hpp"
00032 #include "idot.hpp"
00033 #include "interval.hpp"
00034 
00035 #include "cinterval.hpp"
00036 
00037 namespace cxsc {
00038 
00039 complex::complex(const cdotprecision & a) throw()
00040 {
00041    *this=rnd(a);
00042 }
00043 
00044 complex & complex::operator =(const cdotprecision & a) throw()
00045 {
00046    return *this=rnd(a);
00047 }
00048 
00049 bool operator== (const complex & a, const dotprecision & b)    throw() { return !a.im && a.re==b; }
00050 bool operator== (const dotprecision & a, const complex & b)    throw() { return !b.im && a==b.re; }
00051 bool operator!= (const complex & a, const dotprecision & b)    throw() { return !!a.im || a.re!=b; }
00052 bool operator!= (const dotprecision & a, const complex & b)    throw() { return !!b.im || a!=b.re; }
00053 
00054 
00055 
00056 // ---- Ausgabefunkt. ---------------------------------------
00057 
00058 std::ostream & operator << (std::ostream &s, const complex& a) throw()
00059 {
00060    s << '('          
00061      << a.re << ','  
00062      << a.im       
00063      << ')';
00064    return s;
00065 }
00066 std::string & operator << (std::string &s, const complex &a) throw()
00067 {
00068    s+='(';
00069    s << a.re;
00070    s+=',';
00071    s << a.im; 
00072    s+=')';
00073    return s;
00074 }
00075 
00076 std::istream & operator >> (std::istream &s, complex &a) throw()
00077 {
00078    char c;
00079 
00080    skipeolnflag = inpdotflag = true;
00081    c = skipwhitespacessinglechar (s, '(');
00082    if (inpdotflag) 
00083       s.putback(c);
00084 
00085    s >> a.re;
00086 
00087    skipeolnflag = inpdotflag = true;
00088    c = skipwhitespacessinglechar (s, ',');
00089    if (inpdotflag) s.putback(c);
00090 
00091    s >> a.im >> RestoreOpt;
00092 
00093    if (!waseolnflag) 
00094    {
00095       skipeolnflag = false, inpdotflag = true;
00096       c = skipwhitespaces (s);
00097       if (inpdotflag && c != ')') 
00098          s.putback(c);
00099    }
00100 
00101    /*if (a.re > a.im) {
00102       errmon (ERR_INTERVAL(EMPTY));
00103    } */         
00104    return s;
00105 }
00106 
00107 std::string & operator >> (std::string &s, complex &a) throw()
00108 {
00109    s = skipwhitespacessinglechar (s, '(');
00110    s >> SaveOpt >> RndDown >> a.re;
00111    s = skipwhitespacessinglechar (s, ',');
00112    s >> RndUp >> a.im >> RestoreOpt;
00113    s = skipwhitespaces (s);
00114 
00115    if (s[0] == ')') 
00116       s.erase(0,1);
00117 
00118     /*if (a.re > a.im) {
00119       errmon (ERR_INTERVAL(EMPTY));
00120     } */                         
00121    return s;
00122 }
00123 
00124 void operator >>(const std::string &s,complex &a) throw()
00125 {
00126    std::string r(s);
00127    r>>a;
00128 }
00129 void operator >>(const char *s,complex &a) throw()
00130 {
00131    std::string r(s);
00132    r>>a;
00133 }
00134 
00135 complex exp(const complex& x) throw() { return mid(exp(cinterval(x))); }
00136 complex expm1(const complex& x) throw() { return mid(expm1(cinterval(x))); }
00137 complex exp2(const complex& x) throw() { return mid(exp2(cinterval(x))); }
00138 complex exp10(const complex& x) throw() { return mid(exp10(cinterval(x))); }
00139 complex cos(const complex& x) throw() { return mid(cos(cinterval(x))); }
00140 complex sin(const complex& x) throw() { return mid(sin(cinterval(x))); }
00141 complex cosh(const complex& x) throw() { return mid(cosh(cinterval(x))); }
00142 complex sinh(const complex& x) throw() { return mid(sinh(cinterval(x))); }
00143 
00144 complex tan(const complex& x) throw() { return mid(tan(cinterval(x))); }
00145 complex cot(const complex& x) throw() { return mid(cot(cinterval(x))); }
00146 complex tanh(const complex& x) throw() { return mid(tanh(cinterval(x))); }
00147 complex coth(const complex& x) throw() { return mid(coth(cinterval(x))); }
00148 
00149 real arg(const complex& x) throw() { return mid(arg(cinterval(x))); }
00150 real Arg(const complex& x) throw() { return mid(Arg(cinterval(x))); }
00151 
00152 complex ln(const complex& x) throw() { return mid(ln(cinterval(x))); }
00153 complex lnp1(const complex& x) throw() { return mid(lnp1(cinterval(x))); }
00154 complex log2(const complex& x) throw() { return mid(log2(cinterval(x))); }
00155 complex log10(const complex& x) throw() { return mid(log10(cinterval(x))); }
00156 
00157 complex sqr(const complex& x) throw() { return mid(sqr(cinterval(x))); }
00158 
00159 complex sqrt(const complex& x) throw() { return mid(sqrt(cinterval(x))); }
00160 complex sqrtp1m1(const complex& x) throw() { return mid(sqrtp1m1(cinterval(x))); }
00161 complex sqrt1px2(const complex& x) throw() { return mid(sqrt1px2(cinterval(x))); }
00162 complex sqrtx2m1(const complex& x) throw() { return mid(sqrtx2m1(cinterval(x))); }
00163 complex sqrt1mx2(const complex& x) throw() { return mid(sqrt1mx2(cinterval(x))); }
00164 complex sqrt(const complex& x, int d) throw() 
00165 { return mid(sqrt(cinterval(x),d)); }
00166 
00167 std::list<complex> sqrt_all( const complex& c )
00168 { 
00169         complex z;
00170         z = sqrt(c);
00171                 
00172         std::list<complex> res;
00173         res.push_back(  z );
00174         res.push_back( -z );
00175 
00176         return res;
00177 }       // end sqrt_all
00178 
00179 std::list<complex> sqrt_all( const complex& z, int n )
00180                         //
00181 //  sqrt_all(z,n) computes a list of n values approximating all n-th roots of z
00182                         //
00183 //  For n >=3, computing the optimal approximations is very expensive
00184 //  and thus not considered cost-effective.
00185                         //
00186 //  Hence, the polar form is used to calculate sqrt_all(z,n)
00187                         //
00188 {
00189         std::list<complex> res;
00190 
00191         if( n == 0 )
00192         {
00193                 res.push_back( complex(1,0) );
00194                 return res;
00195         }
00196         else if( n == 1 )
00197         {
00198                 res.push_back(z);
00199                 return res;
00200         }
00201         else if( n == 2 ) return sqrt_all( z );
00202         else
00203         {
00204                 real
00205                                 arg_z = arg( z ), root_abs_z = sqrt( abs( z ), n );
00206 
00207                 for(int k = 0; k < n; k++)
00208                 {
00209                         real arg_k = ( arg_z + 2 * k * mid(Pi_interval) ) / n;
00210 
00211                         res.push_back( complex( root_abs_z * cos( arg_k ),
00212                                                                          root_abs_z * sin( arg_k ) ) );
00213                 }
00214                 return res;
00215         }
00216 }
00217 //-- end sqrt_all -------------------------------------------------------------
00218         
00219 
00220 complex power_fast(const complex& x,int d) throw() 
00221 { return mid(power_fast(cinterval(x),d)); }
00222 complex power(const complex& x,int d) throw() 
00223 { return mid(power(cinterval(x),d)); }
00224 complex pow(const complex& x, const real& r) throw() 
00225 { return mid(pow(cinterval(x),interval(r))); }
00226 complex pow(const complex& x, const complex& y) throw() 
00227 { return mid(pow(cinterval(x),cinterval(y))); }
00228 
00229 complex asin(const complex& x) throw() { return mid(asin(cinterval(x))); }
00230 complex acos(const complex& x) throw() { return mid(acos(cinterval(x))); }
00231 complex asinh(const complex& x) throw() { return mid(asinh(cinterval(x))); }
00232 complex acosh(const complex& x) throw() { return mid(acosh(cinterval(x))); }
00233 complex atan(const complex& x) throw() { return mid(atan(cinterval(x))); }
00234 complex acot(const complex& x) throw() { return mid(acot(cinterval(x))); }
00235 complex atanh(const complex& x) throw() { return mid(atanh(cinterval(x))); }
00236 complex acoth(const complex& x) throw() { return mid(acoth(cinterval(x))); }
00237 
00238 } // namespace cxsc
00239