C-XSC - A C++ Class Library for Extended Scientific Computing  2.5.4
l_cmath.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: l_cmath.cpp,v 1.11 2014/01/30 17:23:46 cxsc Exp $ */
00025 
00026 #include "l_cmath.hpp"
00027 
00028 namespace cxsc {
00029         
00030         l_complex sqrt(const l_complex& z) throw()
00031 // Computation of sqrt(z); stagprec <= stagmax=19 determines the maximum 
00032 // accuracy of about 16*19 = 304 decimal digits.
00033 // The branch cut of this sqrt-function is the negative imaginary axis.
00034 // stagmax=19 is predetermined by the internal used function 
00035 // l_real sqrt(const l_real&), which has its own maximum precision of 
00036 // stagprec=19;
00037         {
00038                 l_real x=Re(z), y=Im(z), w;
00039                 if ( zero_(x) && zero_(y) ) return l_complex(0,0);
00040     // Now z != 0
00041                 int stagsave = stagprec, stagmax = 19;
00042                 if (stagprec > stagmax) stagprec = stagmax;
00043     // stagprec>19 makes no sense, because stagmax(sqrt(...))=19;
00044                 int ex = expo(x[1]), exy = expo(y[1]);
00045                 if (exy > ex) ex = exy; // ex: maximum of the exponents;
00046                 ex = 400-ex; // ex: optimal scaling factor
00047                 if (ex%2) ex--;  // ex is now even;
00048                 times2pown(x,ex);  times2pown(y,ex); // scaling with 2^ex
00049                 bool neg = sign(x[1]) < 0;
00050                 if (neg) x = -x;
00051                 w = abs( l_complex(x,y) ) + x;
00052                 times2pown(w,1);
00053                 w = sqrt(w);
00054                 if (neg)
00055                 {
00056                         x = abs(y)/w;
00057                         y = sign(y[1]) < 0 ? -w : w;
00058                         times2pown(y,-1); 
00059                 } else
00060                 {
00061                         x = w;
00062                         times2pown(x,-1); 
00063                         y /= w;
00064                 }
00065     ex /= 2; // Backscaling of the current result with 2^(-ex):
00066          times2pown(x,-ex);  times2pown(y,-ex);
00067          stagprec = stagsave; // restore old value of the stagprec variable
00068          return l_complex(x,y);
00069         } // sqrt
00070         
00071         l_complex sqrtp1m1(const l_complex& z) throw()
00072         { return mid(sqrtp1m1(l_cinterval(z))); }
00073         
00074         l_complex sqrt1px2(const l_complex& z) throw()
00075         { return mid(sqrt1px2(l_cinterval(z))); }
00076         
00077         l_complex sqrtx2m1(const l_complex& z) throw()
00078         { return mid(sqrtx2m1(l_cinterval(z))); }
00079         
00080         l_complex sqrt1mx2(const l_complex& z) throw()
00081         { return mid(sqrt1mx2(l_cinterval(z))); }
00082 
00083         l_complex exp(const l_complex& z) throw()
00084         { return mid(exp(l_cinterval(z))); }
00085         
00086         l_complex expm1(const l_complex& z) throw()
00087         { return mid(expm1(l_cinterval(z))); }
00088         
00089         l_complex exp2(const l_complex& z) throw()
00090         { return mid(exp2(l_cinterval(z))); }
00091         
00092         l_complex exp10(const l_complex& z) throw()
00093         { return mid(exp10(l_cinterval(z))); }
00094 
00095         l_complex sin(const l_complex& z) throw()
00096         { return mid(sin(l_cinterval(z))); }
00097 
00098         l_complex cos(const l_complex& z) throw()
00099         { return mid(cos(l_cinterval(z))); }
00100 
00101         l_complex tan(const l_complex& z) throw()
00102         { return mid(tan(l_cinterval(z))); }
00103 
00104         l_complex cot(const l_complex& z) throw()
00105         { return mid(cot(l_cinterval(z))); }    
00106         
00107         l_complex asin(const l_complex& z) throw()
00108         { return mid(asin(l_cinterval(z))); }
00109 
00110         l_complex acos(const l_complex& z) throw()
00111         { return mid(acos(l_cinterval(z))); }
00112 
00113         l_complex atan(const l_complex& z) throw()
00114         { return mid(atan(l_cinterval(z))); }
00115 
00116         l_complex acot(const l_complex& z) throw()
00117         { return mid(acot(l_cinterval(z))); }   
00118         
00119         l_complex sinh(const l_complex& z) throw()
00120         { return mid(sinh(l_cinterval(z))); }
00121 
00122         l_complex cosh(const l_complex& z) throw()
00123         { return mid(cosh(l_cinterval(z))); }
00124 
00125         l_complex tanh(const l_complex& z) throw()
00126         { return mid(tanh(l_cinterval(z))); }
00127 
00128         l_complex coth(const l_complex& z) throw()
00129         { return mid(coth(l_cinterval(z))); }
00130         
00131         l_complex asinh(const l_complex& z) throw()
00132         { return mid(asinh(l_cinterval(z))); }
00133 
00134         l_complex acosh(const l_complex& z) throw()
00135         { return mid(acosh(l_cinterval(z))); }
00136 
00137         l_complex atanh(const l_complex& z) throw()
00138         { return mid(atanh(l_cinterval(z))); }
00139 
00140         l_complex acoth(const l_complex& z) throw()
00141         { return mid(acoth(l_cinterval(z))); }
00142 
00143    // sqrt_all(c) computes a list of 2 values for all square roots of c
00144         std::list<l_complex> sqrt_all( const l_complex& c )
00145         { 
00146                 l_complex lc;
00147                 lc = sqrt(c);
00148                 
00149                 std::list<l_complex> res;
00150                 res.push_back(  lc );
00151                 res.push_back( -lc );
00152 
00153                 return res;
00154         }       // end sqrt_all
00155         
00156         l_complex sqrt(const l_complex& z, int n) throw()
00157         { return mid(sqrt(l_cinterval(z),n)); }         
00158                         
00159         l_real arg(const l_complex& z) throw()
00160         { return mid(arg(l_cinterval(z))); }
00161         
00162         l_real Arg(const l_complex& z) throw()
00163         { return mid(Arg(l_cinterval(z))); }
00164         
00165         std::list<l_complex> sqrt_all( const l_complex& z, int n )
00166 //
00167 //  sqrt_all(z,n) computes a list of n values approximating all n-th roots of z
00168 //
00169 //  For n >=3, computing the optimal approximations is very expensive
00170 //  and thus not considered cost-effective.
00171 //
00172 //  Hence, the polar form is used to calculate sqrt_all(z,n)
00173 //
00174         {
00175                 std::list<l_complex> res;
00176 
00177                 if( n == 0 )
00178                 {
00179                         res.push_back( l_complex(1,0) );
00180                         return res;
00181                 }
00182                 else if( n == 1 )
00183                 {
00184                         res.push_back(z);
00185                         return res;
00186                 }
00187                 else if( n == 2 ) return sqrt_all( z );
00188                 else
00189                 {
00190                         l_real
00191                                 arg_z = arg( z ), root_abs_z = sqrt( abs( z ), n );
00192 
00193                         for(int k = 0; k < n; k++)
00194                         {
00195                                 l_real arg_k = ( arg_z + 2 * k * mid(Pi_l_interval()) ) / n;
00196 
00197                                 res.push_back( l_complex( root_abs_z * cos( arg_k ),
00198                                                                                    root_abs_z * sin( arg_k ) ) );
00199                         }
00200                         return res;
00201                 }
00202         }
00203         //
00204 //-- end sqrt_all -------------------------------------------------------------
00205         
00206         l_complex ln(const l_complex& z) throw()
00207         { return mid(ln(l_cinterval(z))); }
00208         
00209         l_complex lnp1(const l_complex& z) throw()
00210         { return mid(lnp1(l_cinterval(z))); }
00211         
00212         l_complex power_fast(const l_complex& z, int n) throw()
00213         {
00214                 if( n == 0 )
00215                         return l_complex(1,0);
00216            else 
00217                         if( n == 1 ) return z;
00218               else 
00219                                 if( n == -1 ) return 1 / z;
00220                  else 
00221                                         if( n == 2 ) return sqr(z); 
00222         else
00223         {
00224                 l_real abs_z = abs(z);
00225 
00226                 if( n < 0 && abs_z == 0.0 )
00227            //  z contains 0
00228                         cxscthrow (STD_FKT_OUT_OF_DEF
00229                                 ("l_complex power_fast(const l_complex& z, int n ); z == 0."));
00230                 if( abs_z == 0.0 )
00231                         return l_complex(0,0);
00232                 else
00233                 {
00234                         l_real arg_z = arg(z);
00235                         l_real abs_z_n = exp( n * ln( abs_z ) );
00236 
00237                         return l_complex( abs_z_n * cos( n * arg_z ),
00238                                                                         abs_z_n * sin( n * arg_z ) );
00239                 }
00240         }
00241 }
00242 
00243 l_complex power(const l_complex& z, int n) throw()
00244 { return mid( power(l_cinterval(z),n) ); }
00245         
00246 l_complex log2(const l_complex& z) throw()
00247 { return mid(log2(l_cinterval(z))); }
00248 
00249 l_complex log10(const l_complex& z) throw()
00250 { return mid(log10(l_cinterval(z))); }
00251         
00252 l_complex pow(const l_complex& z, const l_real& p) throw()
00253 { return mid( pow( l_cinterval(z) , l_interval(p) ) ); }
00254         
00255 l_complex pow(const l_complex& z, const l_complex& p) throw()   
00256 { return mid( pow( l_cinterval(z) , l_cinterval(p) ) ); }
00257 
00258 
00259 
00260 } // namespace cxsc