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: 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