- Generated by
1.7.5.1
|
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_cinterval.cpp,v 1.12 2014/01/30 17:23:47 cxsc Exp $ */ 00025 00026 #include "lx_cinterval.hpp" 00027 00028 namespace cxsc { 00029 00030 // ----------------------------------------------------------------------- 00031 // --------- Functions and operators related to type lx_cinterval -------- 00032 // ----------------------------------------------------------------------- 00033 00034 l_cinterval & l_cinterval::operator = (const lx_cinterval &a) throw() 00035 // This operator is declared in l_cinterval.hpp !! 00036 { 00037 l_interval lre,lim; 00038 lre = Re(a); 00039 lim = Im(a); 00040 l_cinterval lc(lre,lim); 00041 00042 return (*this) = lc; 00043 } 00044 00045 cinterval & cinterval::operator = (const lx_cinterval &a) throw() 00046 // This operator is declared in cinterval.hpp !! 00047 { 00048 l_cinterval lc; 00049 cinterval z; 00050 00051 lc = a; 00052 z = lc; 00053 00054 return (*this) = z; 00055 } 00056 00057 // ------------------------ Input -------------------------------------- 00058 /* 00059 std::string & operator >> (std::string &s, lx_cinterval &a) throw() 00060 // Funktionierende Version mit Zweier-Exponenten 00061 { 00062 lx_interval Lar, Lai; 00063 00064 s = skipwhitespacessinglechar (s, '('); 00065 s = s >> SaveOpt >> Lar; 00066 s = skipwhitespacessinglechar (s, ','); 00067 s = s >> Lai >> RestoreOpt; 00068 s = skipwhitespacessinglechar (s, ')'); 00069 s = skipwhitespaces (s); // erasing whitespaces 00070 a = lx_cinterval(Lar,Lai); 00071 00072 return s; 00073 } 00074 */ 00075 std::string & operator >> (std::string &s, lx_cinterval &a) throw() 00076 // Funktionierende Version mit Zehner-Exponenten, z.B. 00077 // s = ( {+1, [4,4.01]} , {-2, [7,7]} ) 00078 { 00079 lx_interval Lar, Lai; 00080 string str(s); 00081 int i; 00082 00083 str = skipwhitespacessinglechar (str, '('); 00084 i = str.find("}"); 00085 str.erase(i+1); 00086 str >> SaveOpt >> Lar; 00087 00088 i = s.find("}"); 00089 s.erase(0,i+1); 00090 s = skipwhitespacessinglechar (s, ','); 00091 s >> Lai >> RestoreOpt; 00092 s = ""; 00093 00094 a = lx_cinterval(Lar,Lai); 00095 00096 return s; 00097 } 00098 00099 void operator >> (const std::string &s, lx_cinterval &a) throw() 00100 { 00101 // Writes strings s to variable a of type lx_cinterval; 00102 std::string r(s); 00103 r >> a; 00104 } 00105 00106 void operator >> (const char *s, lx_cinterval &a) throw() 00107 { 00108 std::string r(s); 00109 r >> a; 00110 } 00111 00112 std::istream & operator >> (std::istream &s, lx_cinterval &a) throw() 00113 { 00114 lx_interval Lar, Lai; 00115 char c; 00116 00117 std::cerr << "Real part: {Exponent to base 10, [a,b]} = ?" 00118 << std::endl; 00119 s >> Lar; 00120 std::cerr << "Img. part: {Exponent to base 10, [a,b]} = ?" 00121 << std::endl; 00122 s >> Lai >> RestoreOpt; 00123 a = lx_cinterval(Lar,Lai); 00124 00125 if (!waseolnflag) 00126 { 00127 skipeolnflag = false; inpdotflag = true; 00128 c = skipwhitespaces (s); 00129 if (inpdotflag && c != ')') 00130 s.putback(c); 00131 } 00132 return s; 00133 } 00134 00135 00136 // ----------------------------------------------------------------------------- 00137 // ---------- Elementary functions related to type lx_cinterval ---------------- 00138 // ----------------------------------------------------------------------------- 00139 00140 // -------------------------------------------------------------------------- 00141 // ------------------------- Help functions --------------------------------- 00142 // -------------------------------------------------------------------------- 00143 00144 int max(int a, int b) 00145 { 00146 int res(a); 00147 if (b>a) res = b; 00148 return res; 00149 } 00150 00151 00152 // -------------------------------------------------------------------------- 00153 // --------------------------- sqr function --------------------------------- 00154 // -------------------------------------------------------------------------- 00155 00156 lx_cinterval sqr(const lx_cinterval &z) throw() 00157 { 00158 lx_interval rez(Re(z)), reza(abs(rez)), 00159 imz(Im(z)), imza(abs(imz)); 00160 lx_real 00161 irez = Inf(reza), 00162 srez = Sup(reza), 00163 iimz = Inf(imza), 00164 simz = Sup(imza); 00165 lx_interval 00166 hxl(irez), hxu(srez), hyl(iimz), hyu(simz); 00167 lx_real resxl, resxu; 00168 resxl = Inf( (hxl-hyu)*(hxl+hyu) ); 00169 resxu = Sup( (hxu-hyl)*(hxu+hyl) ); 00170 00171 hxl = rez * imz; 00172 times2pown(hxl,1); 00173 00174 return lx_cinterval( lx_interval(resxl,resxu), hxl ); 00175 } // sqr(...) 00176 00177 // ------------------------------------------------------------------------- 00178 // ---------------------- Begin square root -------------------------------- 00179 // ------------------------------------------------------------------------- 00180 00181 lx_interval Sqrt_zpx_m2( const lx_interval &x, const lx_interval &y ) 00182 // z = x + i*y; 00183 // Calculating res = sqrt( 2*(|z| + |x|) ); 00184 { 00185 lx_interval res; 00186 res = sqrtx2y2(x,y) + abs(x); // New: 22.05.2008; 00187 times2pown(res,1); 00188 return sqrt(res); 00189 } 00190 00191 lx_interval Sqrt_zpx_d2( const lx_interval &x, const lx_interval &y ) 00192 // z = x + i*y; 00193 // Calculating res = sqrt( (|z| + |x|)/2 ); 00194 { 00195 lx_interval res; 00196 res = sqrtx2y2(x,y) + abs(x); // New: 22.05.2008; 00197 times2pown(res,-1); 00198 return sqrt(res); 00199 } 00200 00201 lx_interval Re_Sqrt_Point(const lx_interval &rez, const lx_interval &imz) 00202 // Real part of analytic square root of a POINT INTERVAL only. 00203 // Do not use this as a general function 00204 // - it is only a subroutine for sqrt(...) 00205 // Blomquist, 03.05.2007; 00206 { 00207 lx_interval res; 00208 lx_real irez = Inf( rez ), iimz = Inf( imz ); 00209 00210 if (eq_zero(iimz)) 00211 res = (ge_zero(irez))? sqrt(rez) : lx_interval(0); 00212 else // iimz <> 0 00213 res = (ge_zero(irez))? Sqrt_zpx_d2(rez,imz) : 00214 abs(iimz)/Sqrt_zpx_m2(rez,imz); 00215 return res; 00216 } 00217 00218 lx_interval Im_Sqrt_Point(const lx_interval &rez, const lx_interval &imz) 00219 // Imaginary part of analytic square root of a POINT INTERVAL only. 00220 // Do not use this as a general function 00221 // - it is only a subroutine for sqrt(...) 00222 // Blomquist, 03.05.2007; 00223 { 00224 lx_interval res; 00225 lx_real irez = Inf( rez ), iimz = Inf( imz ); 00226 00227 if (eq_zero(iimz)) 00228 res = (ge_zero(irez))? lx_interval(0) : sqrt(-rez); 00229 else 00230 if (ge_zero(irez)) res = iimz/Sqrt_zpx_m2(rez,imz); 00231 else 00232 { 00233 res = Sqrt_zpx_d2(rez,imz); 00234 if (se_zero(iimz)) res = -res; 00235 } 00236 return res; 00237 } 00238 00239 lx_cinterval sqrt(const lx_cinterval& z) throw() 00240 { 00241 lx_cinterval y; 00242 lx_real 00243 irez = Inf( Re(z) ), 00244 srez = Sup( Re(z) ), 00245 iimz = Inf( Im(z) ), 00246 simz = Sup( Im(z) ); 00247 lx_interval 00248 hxl( irez ), hxu( srez ), hyl( iimz ), hyu( simz ); 00249 lx_real 00250 resxl, resxu, resyl, resyu; 00251 00252 if ( sm_zero(irez) && sm_zero(iimz) && ge_zero(simz) ) 00253 cxscthrow(STD_FKT_OUT_OF_DEF( 00254 "lx_cinterval sqrt(const lx_cinterval& z); z not in principal branch.")); 00255 00256 if (ge_zero(iimz)) 00257 { 00258 resxl = Inf( Re_Sqrt_Point( hxl,hyl ) ); 00259 resxu = Sup( Re_Sqrt_Point( hxu,hyu ) ); 00260 00261 resyl = Inf( Im_Sqrt_Point( hxu,hyl ) ); 00262 resyu = Sup( Im_Sqrt_Point( hxl,hyu ) ); 00263 } 00264 else 00265 if (se_zero(simz)) 00266 { 00267 resxl = Inf( Re_Sqrt_Point( hxl,hyu ) ); 00268 resxu = Sup( Re_Sqrt_Point( hxu,hyl ) ); 00269 00270 resyl = Inf( Im_Sqrt_Point( hxl,hyl ) ); 00271 resyu = Sup( Im_Sqrt_Point( hxu,hyu ) ); 00272 } 00273 else 00274 { 00275 resxl = Inf( sqrt(hxl) ); 00276 resxu = ( -iimz>simz ? Sup( Re_Sqrt_Point( hxu,hyl ) ) 00277 : Sup( Re_Sqrt_Point( hxu,hyu ) ) ); 00278 00279 resyl = Inf( Im_Sqrt_Point( hxl,hyl ) ); 00280 resyu = Sup( Im_Sqrt_Point( hxl,hyu ) ); 00281 } 00282 y = lx_cinterval( lx_interval(resxl,resxu), lx_interval(resyl,resyu) ); 00283 00284 return y; 00285 } // sqrt(...) 00286 00287 // sqrt_all(z) computes a list of 2 intervals containing all square roots of z 00288 // 00289 std::list<lx_cinterval> sqrt_all( const lx_cinterval& z ) throw() 00290 { 00291 lx_real 00292 irez = Inf( Re(z) ), 00293 srez = Sup( Re(z) ), 00294 iimz = Inf( Im(z) ), 00295 simz = Sup( Im(z) ); 00296 lx_interval hxl( irez ), hxu( srez ), hyl( iimz ), hyu( simz ); 00297 lx_real resxl, resxu, resyl, resyu; 00298 lx_cinterval w; 00299 00300 if( irez < 0.0 && iimz <= 0.0 && simz >= 0.0 ) 00301 // z contains negative real values 00302 { 00303 if( iimz == 0.0 ) 00304 // z in upper half plane 00305 // principal values can be used 00306 { 00307 // min( Re ( sqrt( z ) ) ) in lower left corner 00308 // max( Re ( sqrt( z ) ) ) in upper right corner 00309 resxl = Inf( Re_Sqrt_Point( hxl, hyl ) ); 00310 resxu = Sup( Re_Sqrt_Point( hxu, hyu ) ); 00311 // min( Im ( sqrt( z ) ) ) in lower right corner 00312 // max( Im ( sqrt( z ) ) ) in upper left corner 00313 resyl = Inf( Im_Sqrt_Point( hxu, hyl ) ); 00314 resyu = Sup( Im_Sqrt_Point( hxl, hyu ) ); 00315 } 00316 else 00317 { 00318 if( simz == 0.0 ) 00319 // z in lower half plane 00320 // principal values can be used in lower corners 00321 { 00322 // z in lower half plane 00323 // min( Re ( sqrt( z ) ) ) in upper left corner 00324 // max( Re ( sqrt( z ) ) ) in lower right corner 00325 resxl = 0; 00326 resxu = Sup( Re_Sqrt_Point( hxu, hyl ) ); 00327 // min( Im ( sqrt( z ) ) ) in lower left corner 00328 // max( Im ( sqrt( z ) ) ) in upper right corner 00329 resyl = Inf( Im_Sqrt_Point( hxl, hyl ) ); 00330 resyu = ( srez > 0.0 ? lx_real(0.0,l_real(0)) : -Inf(sqrt(-hxu) ) ); 00331 } 00332 else 00333 // 0 is interior point of Im( z ) 00334 { 00335 if( srez <= 0.0 ) 00336 { 00337 // 0 is no interior point of Re (z ) 00338 // in quadrant III, Re( z ) = Im_Sqrt_Point(|x|,y), 00339 // Im( z ) = Re_Sqrt_Point(|x|,y) 00340 // min( Re ( sqrt( z ) ) ) in lower right corner = quadrant II/III 00341 // max( Re ( sqrt( z ) ) ) in upper right corner = quadrant II 00342 resxl = Inf( Im_Sqrt_Point(-hxu, hyl ) ); 00343 resxu = Sup( Re_Sqrt_Point( hxu, hyu ) ); 00344 // min( Im ( sqrt( z ) ) ) on real line 00345 // max( Im ( sqrt( z ) ) ) in lower or upper left corner 00346 resyl = Inf( sqrt( -hxu ) ); 00347 resyu = ( - iimz > simz ? Sup( Re_Sqrt_Point( -hxl, hyl ) ) : 00348 Sup( Im_Sqrt_Point( hxl, hyu ) ) ); 00349 } 00350 else 00351 // 0 is interior point of z 00352 // here, the principal values apply in all corner points 00353 { 00354 resxl = 0; 00355 // max( Re ( sqrt( z ) ) ) in lower or upper right corner 00356 resxu = ( - iimz > simz ? Sup( Re_Sqrt_Point( hxu, hyl ) ) : 00357 Sup( Re_Sqrt_Point( hxu, hyu ) ) ); 00358 // min( Im ( sqrt( z ) ) ) in lower left corner 00359 // max( Im ( sqrt( z ) ) ) in upper left corner 00360 resyl = Inf( Im_Sqrt_Point( hxl, hyl ) ); 00361 resyu = Sup( Im_Sqrt_Point( hxl, hyu ) ); 00362 } 00363 } 00364 } 00365 w = lx_cinterval( lx_interval(resxl,resxu), lx_interval(resyl,resyu ) ); 00366 } 00367 else 00368 // sqrt( z ) is well-defined 00369 w = sqrt( z ); 00370 00371 std::list<lx_cinterval> res; 00372 res.push_back( w ); 00373 res.push_back( -w ); 00374 00375 return res; 00376 } 00377 //-- end sqrt_all ------------------------------------------------------------- 00378 // ---------------------- End square root ----------------------------------- 00379 00380 // *************************************************************************** 00381 // *************************************************************************** 00382 // *** Single-valued functions **** 00383 // *************************************************************************** 00384 // *************************************************************************** 00385 00386 00387 // *************************************************************************** 00388 // *** Power operator pow is not listed here, since it relies on the **** 00389 // *** (multi-valued) logarithm **** 00390 // *************************************************************************** 00391 00392 00393 // *************************************************************************** 00394 // *** The hyperbolic functions exp, sin, cos, sinh, cosh are separable: **** 00395 // *** Their real and imaginary parts are products of real functions **** 00396 // *************************************************************************** 00397 // *** With Re(z)=x, Im(z)=y : **** 00398 // *** **** 00399 // *** exp : Re(exp(z)) = exp(x) * cos(y) **** 00400 // *** Im(exp(z)) = exp(x) * sin(y) **** 00401 // *** **** 00402 // *** sin : Re(sin(z)) = sin(x) * cosh(y) **** 00403 // *** Im(sin(x)) = cos(x) * sinh(y) **** 00404 // *** **** 00405 // *** cos : Re(cos(z)) = cos(x) * cosh(y) **** 00406 // *** Im(sin(x)) = -sin(x) * sinh(y) **** 00407 // *** **** 00408 // *** sinh : Re(sinh(z)) = sinh(x) * cos(y) **** 00409 // *** Im(sinh(z)) = cosh(x) * sin(y) **** 00410 // *** **** 00411 // *** cosh : Re(cosh(z)) = cosh(x) * cos(y) **** 00412 // *** Im(cosh(z)) = sinh(x) * sin(y) **** 00413 // *** **** 00414 // *************************************************************************** 00415 00416 // -------------------------- exp(...) -------------------------------------- 00417 00418 lx_cinterval exp(const lx_cinterval& z) throw() 00419 { 00420 int stagsave = stagprec, 00421 stagmax = 39; 00422 if (stagprec > stagmax) 00423 stagprec = stagmax; 00424 00425 lx_interval lreal(Re(z)), limg(Im(z)); 00426 lx_cinterval y; 00427 lx_interval 00428 A( exp(lreal) ), 00429 B( limg ); 00430 y = lx_cinterval( A*cos( B ) , A*sin( B ) ); 00431 00432 stagprec = stagsave; 00433 y = adjust(y); 00434 00435 return y; 00436 } 00437 00438 lx_cinterval exp2(const lx_cinterval& z) throw() 00439 { 00440 return exp(z*Ln2_lx_interval()); 00441 } 00442 00443 lx_cinterval exp10(const lx_cinterval& z) throw() 00444 { 00445 return exp(z*Ln10_lx_interval()); 00446 } 00447 00448 // -------------------------- cos(...) -------------------------------------- 00449 00450 lx_cinterval cos(const lx_cinterval& z) throw() 00451 { 00452 int stagsave = stagprec, 00453 stagmax = 39; 00454 if (stagprec > stagmax) 00455 stagprec = stagmax; 00456 00457 lx_interval lreal(Re(z)),limg(Im(z)); 00458 lx_cinterval y; 00459 lx_interval 00460 A( lreal ), 00461 B( limg ); 00462 y = lx_cinterval( cos( A )*cosh( B ) , -sin( A )*sinh( B ) ); 00463 00464 stagprec = stagsave; 00465 y = adjust(y); 00466 00467 return y; 00468 } 00469 00470 // -------------------------- sin(...) -------------------------------------- 00471 00472 lx_cinterval sin(const lx_cinterval& z) throw() 00473 { 00474 int stagsave = stagprec, 00475 stagmax = 39; 00476 if (stagprec > stagmax) 00477 stagprec = stagmax; 00478 00479 lx_interval lreal(Re(z)),limg(Im(z)); 00480 lx_cinterval y; 00481 lx_interval 00482 A( lreal ), 00483 B( limg ); 00484 y = lx_cinterval( sin( A )*cosh( B ) , cos( A )*sinh( B ) ); 00485 00486 stagprec = stagsave; 00487 y = adjust(y); 00488 00489 return y; 00490 } 00491 00492 // -------------------------- cosh(...) ------------------------------------- 00493 00494 lx_cinterval cosh(const lx_cinterval& z) throw() 00495 { 00496 int stagsave = stagprec, 00497 stagmax = 39; 00498 if (stagprec > stagmax) 00499 stagprec = stagmax; 00500 00501 lx_interval lreal(Re(z)),limg(Im(z)); 00502 lx_cinterval y; 00503 lx_interval 00504 A( lreal ), 00505 B( limg ); 00506 y = lx_cinterval( cos( B )*cosh( A ) , sin( B )*sinh( A ) ); 00507 00508 stagprec = stagsave; 00509 y = adjust(y); 00510 00511 return y; 00512 } 00513 00514 // -------------------------- sinh(...) ------------------------------------- 00515 00516 lx_cinterval sinh(const lx_cinterval& z) throw() 00517 { 00518 int stagsave = stagprec, 00519 stagmax = 39; 00520 if (stagprec > stagmax) 00521 stagprec = stagmax; 00522 00523 lx_interval lreal(Re(z)),limg(Im(z)); 00524 lx_cinterval y; 00525 lx_interval 00526 A( lreal ), 00527 B( limg ); 00528 y = lx_cinterval( cos( B )*sinh( A ) , sin( B )*cosh( A ) ); 00529 00530 stagprec = stagsave; 00531 y = adjust(y); 00532 00533 return y; 00534 } 00535 00536 //----------------------------------------------------------------------------- 00537 // 00538 // Part I: Multi-valued functions 00539 // 00540 //----------------------------------------------------------------------------- 00541 // 00542 // 00543 00544 lx_interval Atan(const lx_interval& y,const lx_interval& x) 00545 // Calculating an inclusion of atan(y/x) with x<>[0,0]; 00546 // This help function must only be used for POINT intervals y,x !! 00547 // This function avoids internal overflow by calculating y/x. 00548 // Blomquist, 25.05.2008; 00549 { 00550 lx_interval res(0),u; 00551 l_interval yl(li_part(y)); 00552 int ex_yl(expo_gr(yl)), 00553 signy(sign(Inf(y))), signx(sign(Inf(x))); 00554 real ex_y(expo(Inf(y))), ex_x(expo(Inf(x))); 00555 bool neg; 00556 00557 if (ex_yl > -1000000) 00558 { // y != 0; 00559 neg = signy * signx == -1; 00560 if (ex_y > ex_x+4197) 00561 { // We now assume y>0 and x>0: 00562 res = Pi_lx_interval(); 00563 times2pown(res,-1); 00564 if (neg) res = -res; 00565 } 00566 else 00567 { 00568 if (ex_x >= 9007199254738946.0) 00569 // Now y/x can generate an error! 00570 { 00571 if (ex_y<-5217) 00572 { 00573 res = lx_interval( lx_real(0,l_real(0)), 00574 lx_real(-Max_Int_R,l_real(minreal)) ); 00575 if (neg) res = -res; 00576 } 00577 else // ex_y >= -5217 00578 { 00579 res = x; 00580 times2pown(res,-2045); 00581 u = y; 00582 times2pown(u,-2045); 00583 res = atan(u/res); // Division now without an error! 00584 } 00585 } 00586 else 00587 if (ex_x <= -9007199254738891.0) 00588 { 00589 res = x; 00590 times2pown(res,2101); 00591 u = y; 00592 times2pown(u,2101); 00593 res = atan(u/res); 00594 } 00595 else 00596 res = atan(y/x); 00597 } 00598 } 00599 00600 return res; 00601 } 00602 00603 lx_interval Atan(const lx_interval& y, const lx_real& x) 00604 { 00605 return Atan(y,lx_interval(x)); 00606 } 00607 00608 //----------------------------------------------------------------------------- 00609 // 00610 // Section 1: Argument functions 00611 // 00612 //----------------------------------------------------------------------------- 00613 00614 //-- Arg: analytic argument function -------------------------------- 040916 -- 00615 // 00616 // (i) Arg(Z) subset (-pi,pi). 00617 // (ii) Arg([0,0]) = 0. 00618 // (iii) Arg(Z) is undefined if Z contains negative real numbers. 00619 // (iv) Otherwise, Arg(Z) is the interval hull of { Arg(z) | z in Z, z<>0 }. 00620 // 00621 // atan is the inverse function of tan(t), t in (-pi/2,pi/2). 00622 // 00623 lx_interval Arg(const lx_cinterval& z) throw() 00624 { 00625 lx_real 00626 srez = Sup( Re(z) ), 00627 irez = Inf( Re(z) ), 00628 simz = Sup( Im(z) ), 00629 iimz = Inf( Im(z) ); 00630 00631 lx_interval 00632 hxl(irez), hxu(srez), hyl(iimz), hyu(simz), Pid2; // Point intervals 00633 00634 lx_real resl, resu; 00635 00636 Pid2 = Pid2_lx_interval(); 00637 00638 if( iimz > 0.0 ) 00639 // case I: Im(z) > 0 00640 { 00641 resl = ( srez > 0.0 ? Inf( Atan( hyl,hxu ) ) : 00642 ( srez < 0.0 ? Inf( Atan( hyu,hxu ) + Pi_lx_interval() ) : 00643 Inf( Pid2 ) ) ); 00644 resu = ( irez > 0.0 ? Sup( Atan( hyu,hxl ) ) : 00645 ( irez < 0.0 ? Sup( Atan( hyl,hxl ) + Pi_lx_interval() ) : 00646 Sup( Pid2 ) ) ); 00647 return lx_interval( resl, resu ); 00648 } 00649 else 00650 { 00651 if( simz < 0.0 ) 00652 // case II: Im(z) < 0 00653 { 00654 resl = ( irez < 0.0 ? Inf( Atan( hyu,hxl ) - Pi_lx_interval() ) : 00655 ( irez > 0.0 ? Inf( Atan( hyl,hxl ) ) : -Sup( Pid2 ) ) ); 00656 resu = ( srez < 0.0 ? Sup( Atan( hyl,hxu ) - Pi_lx_interval() ) : 00657 ( srez > 0.0 ? Sup( Atan( hyu,hxu ) ) : -Inf(Pid2) ) ); 00658 return lx_interval( resl, resu ); 00659 } 00660 else 00661 // 0 in Im(z) 00662 { 00663 if( irez > 0.0 ) 00664 // case III: Re(z) > 0 00665 // z contains positive real values 00666 { 00667 resl = iimz < 0.0 ? Inf( Atan( hyl,hxl ) ) : lx_real(0.0); 00668 return lx_interval( resl, Sup( Atan( hyu,hxl ) ) ); 00669 } 00670 else 00671 // z contains nonpositive real numbers 00672 { 00673 if( irez < 0.0 ) 00674 { 00675 // case IV: z contains negative real numbers 00676 cxscthrow (STD_FKT_OUT_OF_DEF("lx_interval Arg(const lx_cinterval& z); z contains negative real numbers")); 00677 return lx_interval(0.0); 00678 } 00679 else 00680 // case V: 0 in z, but z doesn't contain negative real numbers 00681 { 00682 if( srez > 0.0 ) 00683 // diam( Re(z) > 0.0 ) 00684 { 00685 resl = iimz < 0.0 ? -Sup(Pid2) : lx_real(0.0); 00686 resu = simz > 0.0 ? Sup(Pid2) : lx_real(0.0); 00687 return lx_interval( resl, resu ); 00688 } 00689 else 00690 // Re(z) == 0.0 00691 { 00692 if( eq_zero(iimz) && eq_zero(simz) ) 00693 // Z == 0 00694 return lx_interval(0.0); 00695 else 00696 { 00697 resl = ( iimz < 0.0 ? - Sup(Pid2) : Inf(Pid2) ); 00698 resu = ( simz > 0.0 ? Sup(Pid2) : -Inf(Pid2) ); 00699 return lx_interval( resl, resu ); 00700 } 00701 } 00702 } 00703 } 00704 } 00705 } 00706 } 00707 // 00708 //-- end Arg ------------------------------------------------------------------ 00709 00710 //-- arg: non-analytic argument function -------------------------------------- 00711 // 00712 // (i) arg(Z) is defined for all Z in IC. 00713 // (ii) arg([0,0]) = 0. 00714 // (iii) arg(Z) subset [-pi,3*pi/2). 00715 // (iv) arg(Z) == Arg(Z) if Arg(Z) is well-defined. 00716 // 00717 00718 lx_interval arg( const lx_cinterval& z ) throw() 00719 { 00720 lx_real 00721 srez = Sup( Re(z) ), 00722 irez = Inf( Re(z) ), 00723 simz = Sup( Im(z) ), 00724 iimz = Inf( Im(z) ); 00725 00726 lx_real resl, resu; 00727 lx_interval Pid2; 00728 Pid2 = Pid2_lx_interval(); 00729 00730 if( sm_zero(irez) && se_zero(iimz) && ge_zero(simz) ) 00731 // z contains negative real values 00732 { 00733 if( gr_zero(srez) ) 00734 // 0 in z and 0 interior point of Re(z) 00735 { 00736 resl = ( sm_zero(iimz)? - Sup( Pi_lx_interval() ) : lx_real(0.0) ); 00737 resu = ( ( sm_zero(iimz) && eq_zero(simz) ) ? lx_real(0.0) : 00738 Sup( Pi_lx_interval() ) ); 00739 return lx_interval( resl, resu ); 00740 } 00741 else 00742 { // srez <= 0.0 00743 if( iimz == simz ) 00744 // z is real interval containing no positive values 00745 return Pi_lx_interval(); 00746 else 00747 // sup( Re(z) ) <= 0, diam( Im(z) ) > 0 00748 { 00749 if( eq_zero(srez) ) 00750 { 00751 resl = ( gr_zero(simz)? Inf(Pid2) : 00752 -Sup( Pi_lx_interval() ) ); 00753 resu = ( sm_zero(iimz)? 00754 ( gr_zero(simz)? Sup(3*Pid2) : -Inf(Pid2) ) : 00755 Sup( Pi_lx_interval() ) ); 00756 return lx_interval( resl, resu ); 00757 } 00758 else 00759 // sup( Re(z) ) < 0, diam( Im(z) ) > 0 00760 { 00761 lx_interval hyl(iimz), hyu(simz); 00762 resl = ( simz > 0.0 ? Inf( Atan( hyu,srez ) + Pi_lx_interval() ) : 00763 -Sup( Pi_lx_interval() ) ); 00764 resu = ( iimz < 0.0 ? ( simz > 0.0 ? Sup( Atan( hyl,srez ) + 00765 Pi_lx_interval() ) : Sup( Atan( hyl,srez ) - 00766 Pi_lx_interval() ) ) : Sup( Pi_lx_interval() ) ); 00767 return lx_interval( resl, resu ); 00768 } 00769 } 00770 } 00771 } 00772 else 00773 // Arg(z) is well-defined 00774 return Arg( z ); 00775 } 00776 // 00777 //-- end arg ------------------------------------------------------------------ 00778 00779 //-------------------- sqrt(z,n): analytic n-th root -------------------------- 00780 // 00781 lx_interval Re_Sqrt_point(const lx_interval& rez, const lx_interval& imz, 00782 int n ) 00783 // before: unsigned int n ---------- 040624 -- 00784 // 00785 // Real part of analytic n-th root. 00786 // Do not use this as a general function 00787 // - it's only a subroutine for sqrt(z,n) and sqrt_all(z,n). 00788 // The calculation is validated but largely overestimating 00789 // if (rez+I*imz) is not a complex number. 00790 // Blomquist, 25.05.2008; 00791 { 00792 lx_interval a = sqrtx2y2(rez,imz); 00793 if( eq_zero(Sup(a)) ) 00794 // a == 0 00795 return lx_interval(0); 00796 else 00797 return sqrt(a, n ) * 00798 cos( Arg( lx_cinterval( rez, imz ) ) / n ); 00799 } 00800 00801 lx_interval Im_Sqrt_point(const lx_interval& rez, const lx_interval& imz, 00802 int n ) // before: unsigned int n --- 040624 -- 00803 // 00804 // Imaginary part of analytic n-th root. 00805 // Do not use this as a general function 00806 // - it's only a subroutine for sqrt(z,n) and sqrt_all(z,n). 00807 // The calculation is validated but largely overestimating 00808 // if (rez+I*imz) is not a complex number. 00809 // Blomquist, 25.05.2008; 00810 { 00811 lx_interval a = sqrtx2y2(rez,imz); 00812 if( eq_zero(Sup(a)) ) 00813 // a == 0 00814 return lx_interval(0); 00815 else 00816 return sqrt(a, n ) * 00817 sin( Arg( lx_cinterval( rez, imz ) ) / n ); 00818 } 00819 00820 lx_cinterval sqrt(const lx_cinterval& z, int n ) throw() 00821 // 00822 // Analytic n-th root function 00823 // sqrt(z,n) is undefined if z contains negative real numbers. 00824 // 00825 { 00826 if( n == 0 ) return lx_cinterval(lx_interval(1)); 00827 if( n == 1 ) return z; 00828 if( n == 2 ) return sqrt(z); 00829 else 00830 { 00831 lx_real 00832 irez = Inf( Re(z) ), 00833 srez = Sup( Re(z) ), 00834 iimz = Inf( Im(z) ), 00835 simz = Sup( Im(z) ); 00836 lx_interval hxl( irez ), hxu( srez ), hyl( iimz ), hyu( simz ); 00837 lx_real resxl, resxu, resyl, resyu; 00838 00839 if( sm_zero(irez) && se_zero(iimz) && ge_zero(simz) ) 00840 { 00841 // z contains negative real values 00842 cxscthrow(STD_FKT_OUT_OF_DEF("lx_cinterval sqrt(const lx_cinterval& z, int n ); z contains negative real values.")); 00843 return z; 00844 } 00845 else 00846 { 00847 if( sm_zero(simz) ) 00848 { 00849 // z in lower half plane 00850 lx_cinterval hres = sqrt( lx_cinterval( Re(z), -Im(z) ), n ); 00851 return lx_cinterval( Re(hres), -Im(hres) ); 00852 } 00853 else 00854 { 00855 if( iimz > 0.0 ) 00856 { 00857 // z in upper half plane 00858 lx_interval tangle = tan( ( Pi_lx_interval() * n ) 00859 / ( 2 * ( n-1 ) ) ); 00860 lx_real tanglel = Inf( tangle ), tangleu = Sup( tangle ); 00861 // 00862 // min( Re( Root( z ) ) ) 00863 // 00864 if ( ge_zero(irez) || Sup( hyl / irez ) <= tanglel ) 00865 // lower boundary right of phi = n*Pi/(2n-2) 00866 // min( Re( Root( z ) ) ) in lower left corner 00867 resxl = Inf( Re_Sqrt_point( hxl, hyl, n ) ); 00868 else 00869 { 00870 if( sm_zero(srez) && Inf( hyl / srez ) >= tangleu ) 00871 // lower boundary left of phi = n*Pi/(2n-2) 00872 // min( Re( Root( z ) ) ) in lower right corner 00873 resxl = Inf( Re_Sqrt_point( hxu, hyl, n ) ); 00874 else 00875 // lower boundary intersects phi = n*Pi/(2n-2) 00876 // min( Re( Root( z ) ) ) on intersection 00877 resxl = Inf( Re_Sqrt_point( iimz / tangle , hyl, n ) ); 00878 } 00879 // 00880 // max( Re( Root( z ) ) ) 00881 // 00882 if ( ge_zero(irez) || Sup( hyu / irez ) <= tanglel ) 00883 // upper boundary right of phi = n*Pi/(2n-2) 00884 // max( Re( Root( z ) ) ) in upper right corner 00885 resxu = Sup( Re_Sqrt_point( lx_interval(srez), 00886 lx_interval(simz), n ) ); 00887 else 00888 { 00889 if ( sm_zero(srez) && Inf( hyu / srez ) >= tangleu ) 00890 // upper boundary left of phi = n*Pi/(2n-2) 00891 // max( Re( Root( z ) ) ) in upper left corner 00892 resxu = Sup( Re_Sqrt_point( hxl, hyu, n ) ); 00893 else 00894 // upper boundary intersects phi = n*Pi/(2n-2) 00895 // max( Re(Root( z )) ) on upper left or right corner 00896 resxu = max( Sup( Re_Sqrt_point( hxl, hyu, n ) ), 00897 Sup( Re_Sqrt_point( hxu, hyu, n ) ) ); 00898 } 00899 // 00900 // min( Im( Root( z ) ) ) 00901 // 00902 if ( ge_zero(srez) || Sup( hyl / srez ) <= tanglel ) 00903 // right boundary right of phi = n*Pi/(2n-2) 00904 // min( Im( Root( z ) ) ) in lower right corner 00905 resyl = Inf( Im_Sqrt_point( hxu, hyl, n ) ); 00906 else 00907 { 00908 if( Inf( hyu / srez ) >= tangleu ) 00909 // right boundary left of phi = n*Pi/(2n-2) 00910 // min( Im( Root( z ) ) ) in upper right corner 00911 resyl = Inf( Im_Sqrt_point( hxu, hyu, n ) ); 00912 else 00913 // right boundary intersects phi = n*Pi/(2n-2) 00914 // min( Im( Root( z ) ) ) on intersection 00915 resyl = Inf(Im_Sqrt_point( hxu, srez * tangle, n )); 00916 } 00917 // 00918 // max( Im( Root( z ) ) ) 00919 // 00920 if( ge_zero(irez) || Sup( hyl / irez ) <= tanglel ) 00921 // left boundary right of phi = n*Pi/(2n-2) 00922 // max( Im( Root( z ) ) ) in upper left corner 00923 resyu = Sup( Im_Sqrt_point( hxl, hyu, n ) ); 00924 else 00925 { 00926 if( Inf( hyu / irez ) >= tangleu ) 00927 // left boundary left of phi = n*Pi/(2n-2) 00928 // max( Im( Root( z ) ) ) in lower left corner 00929 resyu = Sup( Im_Sqrt_point( hxl, hyl, n ) ); 00930 else 00931 // left boundary intersects phi = n*Pi/(2n-2) 00932 // max( Im(Root( z )) ) on lower or upper left corner 00933 resyu = max( Sup( Im_Sqrt_point( hxl, hyl, n ) ), 00934 Sup( Im_Sqrt_point( hxl, hyu, n ) ) ); 00935 } 00936 } 00937 else 00938 { 00939 // z intersects positive real axis 00940 // min( Re( Root( z ) ) ) on real line 00941 // max( Re( Root( z ) ) ) in lower or upper right corner 00942 resxl = ( eq_zero(irez)? lx_real(0.0) : Inf( sqrt( hxl, n ) ) ); 00943 resxu = ( - iimz > simz ? Sup( Re_Sqrt_point( hxu, hyl, n ) ) : 00944 Sup( Re_Sqrt_point( hxu, hyu, n ) ) ); 00945 // min( Im ( sqrt( z ) ) ) in lower left corner 00946 // max( Im ( sqrt( z ) ) ) in upper left corner 00947 resyl = Inf( Im_Sqrt_point( hxl, hyl, n ) ); 00948 resyu = Sup( Im_Sqrt_point( hxl, hyu, n ) ); 00949 } 00950 return lx_cinterval( lx_interval( resxl, resxu ), 00951 lx_interval( resyl, resyu ) ); 00952 } 00953 } 00954 } 00955 } 00956 // 00957 //-- end sqrt(z,n) ----------------------------------------------------------- 00958 00959 //-- sqrt_all ------------------------------------------------------- 040628 -- 00960 // 00961 std::list<lx_cinterval> sqrt_all( const lx_cinterval& z, int n ) throw() 00962 // 00963 // sqrt_all(z,n) computes a list of n intervals containing all n-th roots of z 00964 // 00965 // For n >=3, computing the optimal interval bounds is very expensive 00966 // and thus not considered cost-effective. 00967 // 00968 // Hence, the polar form is used to calculate sqrt_all(z,n) 00969 // (observed overestimation of Re() and Im() in test cases: 5-15%). 00970 // 00971 // z is enclosed into an interval in polar coordinates 00972 // (i.e. a sector of an annulus), from which the roots 00973 // (also polar intervals) are computed error-free 00974 // (save roundoff, which is enclosed). 00975 // 00976 // The the final inclusion of the roots into a rectangular complex interval 00977 // involves some additional overestimation. 00978 // 00979 { 00980 std::list<lx_cinterval> res; 00981 00982 if( n == 0 ) 00983 { 00984 res.push_back( lx_cinterval( lx_interval(0,l_interval(1)), 00985 lx_interval(0,l_interval(0)) ) ); 00986 return res; 00987 } 00988 else if( n == 1 ) 00989 { 00990 res.push_back(z); 00991 return res; 00992 } 00993 else if( n == 2 ) return sqrt_all( z ); 00994 else 00995 { 00996 lx_interval 00997 arg_z = arg( z ), root_abs_z = sqrt( abs(z), n ); 00998 00999 for(int k = 0; k < n; k++) 01000 { 01001 lx_interval arg_k = ( arg_z + 2 * k * Pi_l_interval() ) / n; 01002 01003 res.push_back( lx_cinterval( root_abs_z * cos( arg_k ), 01004 root_abs_z * sin( arg_k ) ) ); 01005 } 01006 return res; 01007 } 01008 } 01009 // 01010 //-- end sqrt_all ------------------------------------------------------------- 01011 01012 01013 //----------------------------------------------------------------------------- 01014 // 01015 // Section 2: Logarithms 01016 // 01017 //----------------------------------------------------------------------------- 01018 01019 // ------------------- Ln: analytic natural logarithm ------------------------- 01020 // 01021 // Ln(z) is undefined if z contains zero; z must not touch the negative real 01022 // axis from below; 01023 // 01024 lx_cinterval Ln(const lx_cinterval& z) throw() 01025 { // Blomquist, 24.09.2007; 01026 int stagsave = stagprec, 01027 stagmax = 30; 01028 if (stagprec > stagmax) stagprec = stagmax; 01029 01030 lx_cinterval y; 01031 lx_real 01032 srez = Sup( Re(z) ), 01033 simz = Sup( Im(z) ), 01034 iimz = Inf( Im(z) ); 01035 lx_interval a1( abs(Re(z)) ), a2( abs(Im(z)) ); 01036 if ( eq_zero(Inf(a1)) && eq_zero(Inf(a2)) ) 01037 // z contains 0 01038 cxscthrow(STD_FKT_OUT_OF_DEF( 01039 "lx_cinterval Ln(const lx_cinterval& z); z contains 0")); 01040 if ( sm_zero(srez) && sm_zero(iimz) && ge_zero(simz) ) 01041 cxscthrow(STD_FKT_OUT_OF_DEF( 01042 "lx_cinterval Ln(const lx_cinterval& z); z not allowed")); 01043 01044 y = lx_cinterval( ln_sqrtx2y2(Re(z),Im(z)), arg(z) ); 01045 01046 stagprec = stagsave; 01047 y = adjust(y); 01048 01049 return y; 01050 } 01051 // 01052 //-- end Ln ------------------------------------------------------------------- 01053 01054 //-- ln: non-analytic natural logarithm ----------------------------- 040923 -- 01055 // 01056 // ln(z) is undefined if z contains zero. 01057 // 01058 lx_cinterval ln(const lx_cinterval& z) throw() 01059 { 01060 int stagsave = stagprec, 01061 stagmax = 30; 01062 if (stagprec > stagmax) stagprec = stagmax; 01063 01064 lx_cinterval y; 01065 lx_interval a1( abs(Re(z)) ), a2( abs(Im(z)) ); 01066 if ( eq_zero(Inf(a1)) && eq_zero(Inf(a2)) ) 01067 // z contains 0 01068 cxscthrow(STD_FKT_OUT_OF_DEF 01069 ("lx_cinterval ln(const lx_cinterval& z); z contains 0")); 01070 y = lx_cinterval( ln_sqrtx2y2(Re(z),Im(z)), arg(z) ); 01071 01072 stagprec = stagsave; 01073 y = adjust(y); 01074 01075 return y; 01076 } 01077 // 01078 //-- end ln ------------------------------------------------------------------- 01079 01080 // ---------------------- log2, log10: analytic functions --------------------- 01081 // 01082 // log2(z),log10 are undefined if z contains zero; z must not touch the 01083 // negative real axis from below; 01084 // 01085 lx_cinterval log2(const lx_cinterval& z) throw() 01086 { // Blomquist, 30.11.2008; 01087 return Ln(z) / Ln2_lx_interval(); 01088 } 01089 01090 lx_cinterval log10(const lx_cinterval& z) throw() 01091 { // Blomquist, 30.11.2008; 01092 return Ln(z) / Ln10_lx_interval(); 01093 } 01094 01095 //----------------------------------------------------------------------------- 01096 // 01097 // Section 3: Power functions 01098 // 01099 //----------------------------------------------------------------------------- 01100 01101 //-- power_fast --------------------------------------------------------------- 01102 // 01103 // Fast, validated power function for integer powers, based on Moivre's formula. 01104 // If n is not an integer value, an error message is generated. 01105 // Medium amount of overestimation. 01106 // 01107 lx_cinterval power_fast(const lx_cinterval& z, const real& n) throw() 01108 { 01109 if( n == 0 ) return lx_cinterval(lx_interval(1)); 01110 else if( n == 1 ) return z; 01111 else if( n == -1 ) return 1 / z; 01112 else if( n == 2 ) return sqr(z); 01113 else 01114 // n >= 3 or n <= -2 01115 // If z is a large interval, z^n is a distorted set, for which the 01116 // inclusion into a complex rectangle contains large overestimation. 01117 // For example, if n * Arg( z ) intersects the cartesian axes 01118 // more than once then 0 is contained in the rectangular enclosure 01119 // of z^n. 01120 // For n <= -2, also inversion of z or z^n is required; 01121 // both inversions lead to large overestimation of the resulting interval. 01122 // 01123 // In power_fast(z,n), z is enclosed into a sector S of 01124 // an annulus. S^n is again some sector S' of a (different) annulus. 01125 // S' is calculated exactly (apart from roundoff), and then enclosed 01126 // into a rectangle. There is a certain amount of overestimation 01127 // but the calculations are quit cheap. 01128 // The implementation is base on Moivre's formula. 01129 { 01130 lx_interval abs_z = abs(z); 01131 01132 if( ((n < 0) && eq_zero(Inf(abs_z))) || !(Is_Integer(n))) 01133 // z contains 0 01134 cxscthrow (STD_FKT_OUT_OF_DEF("lx_cinterval power_fast(const lx_cinterval& z, const real& n); z contains 0 or n is not integer.")); 01135 01136 if( eq_zero(Sup(abs_z)) ) return lx_cinterval( lx_interval(0)); 01137 else 01138 { 01139 lx_interval arg_z = arg(z); 01140 lx_interval abs_z_n = power(abs_z,n); // Legendre algorithm 01141 01142 return lx_cinterval( abs_z_n * cos( n * arg_z ), 01143 abs_z_n * sin( n * arg_z ) ); 01144 } 01145 } 01146 } 01147 // 01148 //-- end power_fast ----------------------------------------------------------- 01149 01150 //------------------------------------ power ---------------------------------- 01151 // 01152 // power function based on the Legendre algotizkm with an integer valued 01153 // exponent n of type real. 01154 // 01155 01156 lx_cinterval power( const lx_cinterval& x, const real& n ) throw() 01157 { 01158 if( !(Is_Integer(n)) ) 01159 // n is not an integer 01160 cxscthrow(STD_FKT_OUT_OF_DEF("lx_cinterval power(const lx_cinterval& z, const real& n); n is not integer.")); 01161 01162 real one(1.0), zhi(2.0), N(n), r; 01163 double dbl; 01164 lx_cinterval y,neu,X(x); 01165 01166 if (x == one) y = x; 01167 else 01168 if (N == 0.0) y = one; 01169 else 01170 { 01171 if (N == 1) y = x; 01172 else 01173 if (N == 2) y = sqr(x); 01174 else 01175 { 01176 if (N < 0) 01177 { 01178 X = 1.0/X; 01179 N = -N; 01180 } 01181 // Initialisierung 01182 if ( !Is_Integer(N/2) ) 01183 y = X; 01184 else y = one; 01185 neu = sqr(X); // neu = X*X; 01186 do { 01187 dbl = _double(N/zhi); 01188 dbl = floor(dbl); 01189 r = (real) dbl; 01190 if ( !Is_Integer( r/2 ) ) 01191 y *= neu; 01192 zhi += zhi; 01193 if (zhi <= N) 01194 neu = sqr(neu); // neu = neu * neu; 01195 } while (zhi <= N); 01196 } 01197 } 01198 return y; 01199 } 01200 // 01201 //-- end power ---------------------------------------------------------------- 01202 01203 //------------------------------------ pow ------------------------------------ 01204 // 01205 // Analytic power function for a real interval exponent, based on Ln. 01206 // 01207 01208 lx_cinterval pow( const lx_cinterval& z, const lx_interval& p ) throw() 01209 { 01210 return exp( p*Ln(z) ); 01211 } 01212 01213 // 01214 //-- end pow ------------------------------------------------------------------ 01215 01216 //------------------------------------ pow ------------------------------------ 01217 // 01218 // Analytic power function for a complex interval exponent, based on Ln. 01219 // 01220 01221 lx_cinterval pow(const lx_cinterval& z, const lx_cinterval& p) throw() 01222 { 01223 return exp( p*Ln(z) ); 01224 } 01225 01226 // 01227 //-- end pow ------------------------------------------------------------------ 01228 01229 //----------------------------------------------------------------------------- 01230 // 01231 // Section 2: tan, cot, tanh, coth 01232 // 01233 // The implementation of cot, tanh, and coth is based on tan: 01234 // 01235 // cot( z ) = tan( pi/2 - z ) 01236 // tanh( z ) = transp( i * tan( transp( i * z ) ) 01237 // coth( z ) = i * cot( i * z ) = i * tan( pi/2 - i * z ) 01238 // 01239 //----------------------------------------------------------------------------- 01240 01241 //---------------------------------- tan -------------------------------------- 01242 // 01243 // Complex tangent function 01244 // 01245 01246 void horizontal_check( const lx_interval& hy, const lx_interval& cosh_2y, 01247 lx_real irez, lx_real srez, const lx_interval& hxl, 01248 const lx_interval& hxu, lx_real& resxl, lx_real& resxu ) 01249 // 01250 // Subroutine of tangent function. 01251 // Check intersections with extremal curves of tan on a horizontal boundary. 01252 // This subroutine is only called if an intersection occurs. 01253 // In this case, sinh( 2 * hy ) <> 0.0 (poles are handled before). 01254 // 01255 // There may be 1 or 2 intersections. 01256 // If intersections lie next to a boundary of rez, then it is impossible to 01257 // decide if there are 1 or 2 intersections. 01258 // In this case, 2 intersections are assumed 01259 // (valid enclosure, at the expense of a potential slight overestimation). 01260 // 01261 { 01262 bool both = false, left = false, right = false; 01263 lx_interval hlp, hlp1; 01264 01265 hlp = Pid2_lx_interval(); 01266 hlp1 = lx_interval(srez) - lx_interval(irez); 01267 if (Inf(hlp1) > Inf(hlp)) 01268 // 2 intersections 01269 both = true; 01270 else 01271 { 01272 lx_interval 01273 res_l = cos( 2 * hxl ) - cosh_2y, 01274 res_u = cos( 2 * hxu ) - cosh_2y; 01275 01276 if( Inf( res_l * res_u ) > 0.0 ) 01277 // 2 intersections 01278 both = true; 01279 else 01280 { 01281 if( Sup( res_l * res_u ) < 0.0 ) 01282 { 01283 // 1 intersection (3 intersections are PI() apart) 01284 // neither of the intervals res_l and res_u contains zero 01285 if( Inf( res_l ) > 0.0 ) 01286 // "left" intersection 01287 left = true; 01288 else 01289 // "right" intersection 01290 right = true; 01291 } 01292 else 01293 // 01294 // 1 (or both) intersections lie next to a boundary point 01295 // here, the task is to decide if 2 intersections occurs 01296 // if only 1 intersection takes place, then which one? 01297 // 01298 { 01299 lx_interval 01300 sin_2xl = sin( 2 * hxl ), 01301 sin_2xu = sin( 2 * hxu ); 01302 01303 if( !Disjoint( lx_interval(0), res_l ) ) 01304 // intersection on the left boundary 01305 { 01306 if( ge_zero(Inf(sin_2xl)) ) 01307 // "left" intersection 01308 { 01309 left = true; 01310 // remove the intersection by changing res_l! 01311 res_l = -lx_interval(1); 01312 } 01313 else 01314 { 01315 if( se_zero(Sup(sin_2xl)) ) 01316 // "right" intersection 01317 { 01318 right = true; 01319 // remove the intersection by changing res_l! 01320 res_l = lx_interval(1); 01321 } 01322 else 01323 // zero is interior point of sin_2xl 01324 // if the real sine function has optimal precision, 01325 // this case should never happen 01326 both = true; 01327 } 01328 } 01329 01330 if( !Disjoint( lx_interval(0), res_u ) ) 01331 // intersection on the right boundary 01332 { 01333 if( ge_zero(Inf(sin_2xu)) ) 01334 // "left" intersection 01335 { 01336 left = true; 01337 // remove the intersection by changing res_u! 01338 res_u = lx_interval(1); 01339 } 01340 else 01341 { 01342 if( se_zero(Sup(sin_2xu)) ) 01343 // "right" intersection 01344 { 01345 right = true; 01346 // remove the intersection by changing res_u! 01347 res_u = -lx_interval(1); 01348 } 01349 else 01350 // zero is interior point of sin_2xu 01351 // if the real sine function has optimal precision, 01352 // this case should never happen 01353 both = true; 01354 } 01355 } 01356 // 01357 // finally, check if there is a second intersection 01358 // 01359 if( Inf( res_l * res_u ) < 0.0 ) 01360 both = true; 01361 } 01362 } 01363 } 01364 // 01365 // Calculate extremal values 01366 // 01367 lx_interval re_tan = 1 / sinh( 2 * abs( hy ) ); 01368 01369 // "left" intersection, sin( 2x ) > 0 01370 if( left || both ) 01371 { 01372 resxl = min( resxl, Inf( re_tan ) ); 01373 resxu = max( resxu, Sup( re_tan ) ); 01374 } 01375 01376 // "right" intersection, sin( 2x ) < 0 01377 if( right || both ) 01378 { 01379 resxl = min( resxl, -Sup( re_tan ) ); 01380 resxu = max( resxu, -Inf( re_tan ) ); 01381 } 01382 } // end horizontal_check 01383 01384 01385 lx_cinterval Tan(const lx_cinterval& z) throw() 01386 { 01387 lx_cinterval y; 01388 lx_interval 01389 rez = Re(z), // rez = z.re(), 01390 imz = Im(z); // imz = z.im(); 01391 01392 lx_real 01393 irez = Inf(rez), 01394 srez = Sup(rez), 01395 iimz = Inf(imz), 01396 simz = Sup(imz); 01397 01398 lx_interval hxl(irez), hxu(srez), hyl(iimz), hyu(simz); 01399 01400 lx_real resxl, resxu, resyl, resyu; 01401 // 01402 // 1st: check for poles 01403 // 01404 if( ( !Disjoint( lx_interval(0), imz ) ) && 01405 ( !Disjoint( lx_interval(0), cos( rez ) ) ) ) 01406 cxscthrow (STD_FKT_OUT_OF_DEF( 01407 "lx_cinterval tan(const lx_cinterval& z); Pole(s) in z")); 01408 // 01409 // 2nd: real part 01410 // 01411 // evaluate tan on vertical boundaries 01412 // 01413 lx_interval 01414 cos_2rez = cos( 2 * rez ), sinh_imz_2 = sqr( sinh( imz ) ); 01415 01416 lx_interval 01417 re_tan_l=sin( 2*hxl ) / ( 2*( sqr( cos( hxl ) ) + sinh_imz_2 ) ), 01418 re_tan_u=sin( 2*hxu ) / ( 2*( sqr( cos( hxu ) ) + sinh_imz_2 ) ); 01419 01420 resxl = min( Inf( re_tan_l ), Inf( re_tan_u ) ); 01421 resxu = max( Sup( re_tan_l ), Sup( re_tan_u ) ); 01422 01423 // 01424 // look for extremal values on horizontal boundaries 01425 // if one of the horizontal boundaries is the x-axes, 01426 // then the complex tangent is the real tangent on this 01427 // boundary, and due to monotonicity, its range 01428 // is already included in the ranges of the vertical boundaries 01429 // 01430 if( irez < srez ) 01431 { 01432 lx_interval 01433 cosh_2yl = - 1 / cosh( 2 * hyl ), 01434 cosh_2yu = - 1 / cosh( 2 * hyu ); 01435 01436 if( !Disjoint( cos_2rez, cosh_2yl ) && iimz != 0.0 ) 01437 //extremal curve intersects lower boundary 01438 horizontal_check(hyl,cosh_2yl,irez,srez,hxl,hxu,resxl,resxu); 01439 01440 if( !Disjoint( cos_2rez, cosh_2yu ) && simz != 0.0 ) 01441 //extremal curve intersects upper boundary 01442 horizontal_check(hyu,cosh_2yu,irez,srez,hxl,hxu,resxl,resxu); 01443 } 01444 // 01445 // 3rd: imaginary part 01446 // 01447 // evaluate tan on horizontal boundaries 01448 // 01449 lx_interval cos_rez_2 = sqr( cos( rez ) ); 01450 01451 lx_interval 01452 im_tan_l = sinh( 2*hyl ) / (2*(cos_rez_2 + sqr( sinh( hyl ) ))), 01453 im_tan_u = sinh( 2*hyu ) / (2*(cos_rez_2 + sqr( sinh( hyu ) ))); 01454 01455 resyl = min( Inf( im_tan_l ), Inf( im_tan_u ) ); 01456 resyu = max( Sup( im_tan_l ), Sup( im_tan_u ) ); 01457 01458 // 01459 // look for extremal values on vertical boundaries 01460 // here, the situation is simpler than for the real part 01461 // if 2 intersections with extremal curves appear , 01462 // one lies in the lower half plane, the other in the upper half plane 01463 // 01464 lx_interval 01465 cos_2xl = cos( 2 * hxl ), cos_2xu = cos( 2 * hxu ); 01466 lx_interval im_tan; 01467 01468 if( iimz < 0.0 ) 01469 // intersection in lower half plane? 01470 { 01471 lx_interval 01472 imz_h = lx_interval( iimz, min( simz, lx_real(0.0) ) ), 01473 cosh_2imz = - 1 / cosh( 2 * imz_h ); 01474 01475 if( ( !Disjoint( cosh_2imz, cos_2xl ) ) ) 01476 //extremal curve intersects left boundary 01477 //in this case, sin( 2 * xl ) <> 0.0 (no poles here!) 01478 { 01479 im_tan = - 1 / abs( sin( 2 * hxl ) ); 01480 resyl = min( resyl, Inf( im_tan ) ); 01481 resyu = max( resyu, Sup( im_tan ) ); 01482 } 01483 if( ( !Disjoint( cosh_2imz, cos_2xu ) ) ) 01484 //extremal curve intersects right boundary 01485 //in this case, sin( 2 * xu ) <> 0.0 (no poles here!) 01486 { 01487 im_tan = - 1 / abs( sin( 2 * hxu ) ); 01488 resyl = min( resyl, Inf( im_tan ) ); 01489 resyu = max( resyu, Sup( im_tan ) ); 01490 } 01491 } 01492 if( simz > 0.0 ) 01493 // intersection in upper half plane? 01494 { 01495 lx_interval 01496 imz_h = lx_interval( max( iimz, lx_real(0.0) ), simz ), 01497 cosh_2imz = - 1 / cosh( 2 * imz_h ); 01498 01499 if( ( !Disjoint( cosh_2imz, cos_2xl ) ) ) 01500 //extremal curve intersects left boundary 01501 //in this case, sin( 2 * xl ) <> 0.0 (no poles here!) 01502 { 01503 im_tan = + 1 / abs( sin( 2 * hxl ) ); 01504 resyl = min( resyl, Inf( im_tan ) ); 01505 resyu = max( resyu, Sup( im_tan ) ); 01506 } 01507 if( ( !Disjoint( cosh_2imz, cos_2xu ) ) ) 01508 //extremal curve intersects right boundary 01509 //in this case, sin( 2 * xu ) <> 0.0 (no poles here!) 01510 { 01511 im_tan = + 1 / abs( sin( 2 * hxu ) ); 01512 resyl = min( resyl, Inf( im_tan ) ); 01513 resyu = max( resyu, Sup( im_tan ) ); 01514 } 01515 } 01516 01517 y = lx_cinterval( lx_interval(resxl,resxu ),lx_interval(resyl,resyu ) ); 01518 01519 return y; 01520 } // Tan 01521 01522 lx_cinterval tan(const lx_cinterval& z) throw() 01523 { 01524 // tan(z) has the poles z_k = pi*(1+2k)/2; |k| in {0,1,2,3,...}. 01525 // z = z_k + eps = pi*(1+2k)/2 + eps; With |eps|<<1 k can be calculated 01526 // by: k = Re(z)/pi - 0.5; With this k the complex value eps is given 01527 // by: eps = z - pi*(1+2k)/2; pi = 3.1415926... ; 01528 // It holds: 01529 // tan(z) = tan(z_k+eps) = tan[pi*(1+2k)/2 + eps] 01530 // = tan[pi/2 + pi*k + eps] = tan[pi/2 + eps] = -1 / tan(eps); 01531 // Definitions: u = Re(eps); u = abs(u); v = Im(eps); v = abs(v); 01532 // if (Sup(u)<S && Sup(v)<S) tan(z) = -1 / Tan(eps); 01533 // else tan(z) = Tan(z); S = 1e-15; 01534 // Thus, near the poles tan(z) is calculated in higher accuracy with 01535 // -1 / Tan(eps); 01536 // Blomquist, 27.09.2007; 01537 const real S = 1e-15; 01538 int stagsave = stagprec, 01539 stagmax = 39; 01540 if (stagprec > stagmax) 01541 stagprec = stagmax; 01542 01543 lx_cinterval y,eps; 01544 l_interval rezl,ul,vl; 01545 01546 lx_interval 01547 rez = Re(z), Pih; // rez = z.re(), 01548 rezl = rez; 01549 interval re,u_,v_; 01550 re = rezl; 01551 real x(mid(re)),k, pi(3.14159265358979323); 01552 lx_interval u,v; 01553 01554 int s; 01555 x = x/pi - 0.5; 01556 s = sign(x); 01557 k = (s>=0)? cutint(x+0.5) : cutint(x-0.5); 01558 if (k == 9007199254740992.0) 01559 cxscthrow (STD_FKT_OUT_OF_DEF( 01560 "lx_cinterval tan(const lx_cinterval& z); z out of range")); 01561 Pih = Pi_lx_interval(); 01562 rez = k*Pih; 01563 times2pown(Pih,-1); // Pih = pi/2 01564 rez = rez + Pih; // rez: inclusion of k*pi + pi/2 01565 eps = z - rez; 01566 01567 u = Re(eps); u = abs(u); 01568 v = Im(eps); v = abs(v); 01569 ul = u; vl = v; 01570 u_ = ul; v_ = vl; 01571 y = (Sup(u_)<S && Sup(v_)<S)? -lx_cinterval(1) / Tan(eps) : Tan(z); 01572 01573 stagprec = stagsave; 01574 y = adjust(y); 01575 01576 return y; 01577 } // tan() 01578 01579 lx_cinterval cot(const lx_cinterval& z) throw() 01580 { 01581 // cot(z) has the poles z_k = k*pi; |k| in {0,1,2,3,...}. 01582 // z = z_k + eps = k*pi + eps; With |eps|<<1 k can be calculated 01583 // by: k = Re(z)/pi; With this k the komplex value eps is given 01584 // by: eps = z - k*pi; pi = 3.1415926... ; 01585 // It holds: 01586 // cot(z) = cot(z_k+eps) = cot(k*pi + eps) 01587 // = cot(eps) = 1 / tan(eps); 01588 // Definitions: u = Re(eps); u = abs(u); v = Im(eps); v = abs(v); 01589 // if (Sup(u)<S && Sup(v)<S) cot(z) = 1 / tan(eps); 01590 // else cot(z) = tan(pi/2 - z); S = 1e-15; 01591 // Thus, near the poles cot(z) is calculated in higher accuracy with 01592 // 1 / tan(eps); 01593 // Blomquist, 28.09.2007; 01594 const real S = 1e-15; 01595 int stagsave = stagprec, 01596 stagmax = 39; 01597 if (stagprec > stagmax) 01598 stagprec = stagmax; 01599 01600 lx_cinterval y,eps; 01601 l_interval rezl,ul,vl; 01602 01603 lx_interval rez = Re(z); 01604 rezl = rez; 01605 interval re,u_,v_; 01606 re = rezl; 01607 real x(mid(re)), k, pi(3.14159265358979323); 01608 lx_interval u,v; 01609 01610 int s; 01611 x = x/pi; 01612 s = sign(x); 01613 k = (s>=0)? cutint(x+0.5) : cutint(x-0.5); 01614 if (k == 9007199254740992.0) 01615 cxscthrow (STD_FKT_OUT_OF_DEF( 01616 "lx_cinterval cot(const lx_cinterval& z); z out of range")); 01617 eps = z - Pi_lx_interval()*k; 01618 u = Re(eps); u = abs(u); 01619 v = Im(eps); v = abs(v); 01620 ul = u; vl = v; 01621 u_ = ul; v_ = vl; 01622 if (Sup(u)<S && Sup(v)<S) 01623 y = 1 / Tan(eps); 01624 else 01625 { 01626 u = Pi_lx_interval(); 01627 times2pown(u,-1); // u = pi/2 01628 y = Tan(u - z); 01629 } 01630 01631 stagprec = stagsave; 01632 y = adjust(y); 01633 01634 return y; 01635 } // cot 01636 01637 //---------------------------------- tanh ----------------------------------- 01638 // 01639 // tanh( z ) = transp( i * tan( transp( i * z ) ) 01640 // 01641 lx_cinterval tanh(const lx_cinterval& z) throw() 01642 { 01643 int stagsave = stagprec, 01644 stagmax = 39; 01645 if (stagprec > stagmax) stagprec = stagmax; 01646 01647 lx_cinterval res = tan( lx_cinterval( Im(z), Re(z) ) ),y; 01648 y = lx_cinterval( Im(res), Re(res) ); 01649 01650 stagprec = stagsave; 01651 y = adjust(y); 01652 01653 return y; 01654 } 01655 // 01656 //-- end tanh ------------------------------------------------------- 01657 01658 //-- coth ----------------------------------------------------------- 01659 // 01660 // coth( z ) = i * cot( i * z ); 01661 // 01662 01663 lx_cinterval coth(const lx_cinterval& z) throw() 01664 { // coth( z ) = i * cot( i * z ); 01665 int stagsave = stagprec, 01666 stagmax = 39; 01667 if (stagprec > stagmax) stagprec = stagmax; 01668 01669 lx_cinterval zh = lx_cinterval( -Im(z), Re(z) ); // zh = i*z; 01670 lx_cinterval res = cot(zh); 01671 zh = lx_cinterval( -Im(res), Re(res) ); 01672 01673 stagprec = stagsave; 01674 zh = adjust(zh); 01675 01676 return zh; 01677 } 01678 // 01679 //-- end coth ----------------------------------------------------------------- 01680 01681 //----------------------------------------------------------------------------- 01682 //----------------------------------------------------------------------------- 01683 // 01684 // Section 5: asin, acos, asinh, acosh 01685 // 01686 // The implementation of acos, asinh, and acosh is based on asin: 01687 // 01688 // acos( z ) = pi / 2 - asin( z ) 01689 // asinh( z ) = i * asin( - i * z ) 01690 // acosh( z ) = i * acos( z ) = +/- i * ( pi / 2 - asin( z ) ) 01691 // 01692 // Only principal values are computed. 01693 // 01694 //----------------------------------------------------------------------------- 01695 01696 //------------------------------------ asin ----------------------------------- 01697 // 01698 // Analytic inverse sine function 01699 // 01700 lx_interval f_aux_asin(const lx_interval& x, const lx_interval& y) 01701 // 01702 // auxiliary function for evaluating the imaginary part of asin(z); 01703 // f_aux_asin(z) = ( |z+1| + |z-1| ) / 2; z = x + i*y; 01704 // = ( sqrt[(x+1)^2 + y^2] + sqrt[(x-1)^2 + y^2] )/2; 01705 // Blomquist, 01.10.2007; 01706 { 01707 lx_interval res; 01708 // interval sqr_y = sqr( y ); 01709 // interval res = ( sqrt( sqr( x + 1.0 ) + sqr_y ) + 01710 // sqrt( sqr( x - 1.0 ) + sqr_y ) ) / 2.0; 01711 res = abs(x); 01712 if (y != 0.0 || Inf(res) < 1.0) 01713 { 01714 res = sqrtx2y2(x+1.0,y) + sqrtx2y2(x-1.0,y); // Blomquist; 01715 times2pown(res,-1); 01716 } 01717 01718 // Now we correct a possible overestimation of the lower bound 01719 // of res. 01720 // It holds: 01721 // (sqrt[(x+1)^2 + y^2] + sqrt[(x-1)^2 + y^2])/2 >= Max(1,|x|); 01722 lx_real hlb = max( lx_real(1.0), abs( Sup(x) ) ); 01723 if( Inf(res) < hlb ) // this is an invalid overestimation! 01724 res = lx_interval( hlb, Sup(res) ); 01725 return res; 01726 } 01727 01728 lx_interval ACOSH_f_aux(const lx_interval& x, const lx_interval& y) 01729 // Function for calculating the imaginary part of asin(z), 01730 // z = x + i*y; 01731 // Notations: 01732 // f_aux_asin(x,y) = alpha(x,y); 01733 // alpha(x,y) := ( sqrt[(x+1)^2 + y^2] + sqrt[(x-1)^2 + y^2] )/2; 01734 // return values: 01735 // 1. res = acosh( f_aux_asin(x,y) ), if |x|>2 or |y|>2; 01736 // 2. res = ln[ alpha + sqrt(alpha^2 - 1) ]; d = alpha(x,y) -1 01737 // = ln[1 + sqrt(d)*(sqrt(d) + sqrt(2+d))], if it holds: 01738 // |x|<=2 und |y|<=2; x,y are point intervals only! 01739 // Blomquist, 02.06.2008; 01740 01741 { 01742 lx_interval res, xa(abs(x)), ya(abs(y)), S1, S2, S3; 01743 lx_real rx((Inf(xa))), ry(Inf(ya)); 01744 01745 if (rx>2.0 || ry>2.0) res = acosh( f_aux_asin(x,y) ); 01746 else 01747 { // Now with improvements: 01748 if (rx == 1.0) 01749 { 01750 res = ya/2; S2 = res/2; 01751 S1 = ya*(0.5 + S2/(sqrt1px2(res)+1)); // S1 = delta 01752 res = lnp1(S1 + sqrt(S1*(2+S1))); 01753 } 01754 else 01755 if (rx<1.0) // 0 <= x < +1; 01756 { 01757 S1 = 1+xa; 01758 S1 = 1/(sqrtx2y2(S1,ya) + S1); 01759 S2 = 1-xa; 01760 S2 = 1/(sqrtx2y2(S2,ya) + S2); 01761 S1 = S1 + S2; 01762 times2pown(S1,-1); // S1 = {...}/2 01763 S2 = ya * sqrt(S1); // S2 = sqrt(delta) 01764 res = lnp1( S2*(S2 + sqrt(2+sqr(ya)*S1)) ); 01765 } 01766 else // 1 < x <= 2; 01767 { 01768 if (ya == 0) 01769 { 01770 S1 = xa-1; 01771 if (eq_zero(Inf(S1))) 01772 S1 = lx_interval(Inf(lx_interval(-2097,l_interval(1))), 01773 Sup(S1)); 01774 } 01775 else // 1 < x <= 2 and 0 < |y| <= 2 01776 { 01777 S1 = xa-1; 01778 if (eq_zero(Inf(S1))) 01779 S1 = lx_interval(Inf(lx_interval(-2097,l_interval(1))), 01780 Sup(S1)); 01781 S2 = ya; 01782 times2pown(S2,1048); 01783 if (Inf(S1)>Inf(S2)) 01784 S1 *= lx_interval( lx_real(1),Sup(One_p_lx_interval()) ); 01785 else 01786 { 01787 res = sqr(ya); 01788 S3 = res / (sqrtx2y2(S1,ya) + S1); // S3 = v; 01789 S2 = 1+xa; 01790 S2 = res / (sqrtx2y2(S2,ya) + S2); // S2 = u; 01791 times2pown(S1,1); // S1 = w=2*(x-1) 01792 S1 = (S3 + S2) + S1; // S1 = u+v+w 01793 times2pown(S1,-1); // S1 = delta 01794 } 01795 } 01796 res = lnp1( S1+sqrt(S1*(2+S1)) ); 01797 } 01798 } 01799 return res; 01800 } // ACOSH_f_aux 01801 01802 lx_interval Asin_beta(const lx_interval& x, const lx_interval& y ) 01803 // Calculating the improved real part of asin(z); Blomquist 22.01.2007; 01804 // Re(asin(z)) = asin[ 2x/(sqrt[(x+1)^2+y^2] + sqrt[(x-1)^2+y^2]) ]=asin[beta]; 01805 // Improvements for beta --> +1 and beta --> -1 are necessary because of 01806 // nearly vertical tangents of the real asin(t)-function for |t|-->+1; 01807 // This function should only be used for POINT intervals x,y! z = x + i*y; 01808 // Blomquist, 30.05.2008; 01809 { 01810 const real c1 = 0.75; 01811 bool neg_x, fertig(false); 01812 lx_real Infxa(Inf(x)), L_y(Inf(y)); 01813 int ex_xl(expo_gr(lr_part(Infxa))), 01814 ex_yl(expo_gr(lr_part(L_y))); 01815 real ex_x(expo(Infxa)), ex_y(expo(L_y)); 01816 lx_interval res,beta,abs_beta,delta,w_delta,xa,Ne,Kla; 01817 01818 if (ex_xl<-1000000) 01819 res = 0; 01820 else 01821 { // x != 0 01822 if ( ex_yl>-1000000 && ex_y>=9007199254739972.0-ex_yl && 01823 ex_x <= 2097-ex_xl ) 01824 { // It holds: |x/y| <= 2^(-9007199254737871) and |y| != 0 01825 // res in [0,eps] or res in [-eps,0], with 01826 // eps = 2^(-9007199254737870); 01827 xa = x; 01828 neg_x = Inf(x)<0; 01829 if (neg_x) xa = -xa; // Inf(xa) >0 : 01830 res = lx_interval( lx_real(0), lx_real(-9007199254737870.0,l_real(1)) ); 01831 if (neg_x) 01832 res = -res; 01833 } 01834 else 01835 { 01836 Kla = sqrtx2y2(1-x,y); 01837 Ne = sqrtx2y2(1+x,y) + Kla; 01838 beta = x / ( Ne/2 ); 01839 if (Inf(beta)<-1) Inf(beta) = -1; 01840 if (Sup(beta)> 1) Sup(beta) = 1; 01841 abs_beta = abs(beta); 01842 if (Inf(abs_beta) < c1) 01843 res = asin(beta); // Normal calculation 01844 else 01845 { // Inf(abs_beta)>=c1; Calculation now with improvements: 01846 Ne = Ne + 2.0; 01847 // Ne = 2 + sqrt[(1+x)^2 + y^2] + sqrt[(1-x)^2 + y^2]; 01848 xa = x; 01849 neg_x = Inf(x)<0; 01850 if (neg_x) xa = -xa; // Inf(xa) >0 : 01851 Infxa = Inf(xa); 01852 if (Infxa > 1.0) 01853 { // x > +1; 01854 if (y == 0.0) 01855 fertig = true; // due to delta == 0; 01856 else 01857 { 01858 beta = xa - 1; 01859 delta = sqrt(beta); 01860 delta = lx_interval(lx_real(-2100,7.4699)) * delta; 01861 // delta = 7.4699*(2^-2100)*sqrt(x-1); 01862 if (Sup(abs(y)) < Inf(delta)) 01863 fertig = true; 01864 else 01865 { 01866 delta = sqr(y/beta); // delta = (y/(x-1))^2 01867 delta = beta * sqrtp1m1(delta); 01868 times2pown(Ne,-1); 01869 delta /= Ne; 01870 } 01871 } 01872 } 01873 else 01874 if (Infxa == 1.0) 01875 { // x = 1; 01876 delta = abs(y); 01877 times2pown(delta,1); 01878 delta /= Ne; 01879 } 01880 else 01881 { // 0.75 <= x < 1; 01882 if (y == 0.0) 01883 delta = 1 - xa; 01884 else // 0.75 <= x < 1 and y != 0: 01885 { 01886 beta = 1 - xa; // beta > 0; 01887 delta = Kla + beta; 01888 times2pown(delta,1); 01889 delta /= (2+Ne); 01890 } 01891 } 01892 res = Pi_lx_interval(); 01893 times2pown(res,-1); // res = pi/2; 01894 if (!fertig) 01895 res -= asin( sqrt(delta*(2-delta)) ); 01896 if (neg_x) res = -res; 01897 } 01898 } 01899 } // x != 0 01900 01901 return res; 01902 } // Asin_beta(...) 01903 01904 lx_cinterval asin(const lx_cinterval& z) throw() 01905 { 01906 int stagsave = stagprec, 01907 stagmax = 30; 01908 if (stagprec>stagmax) stagprec = stagmax; 01909 01910 lx_cinterval res; 01911 lx_interval rez = Re(z), imz = Im(z); 01912 01913 lx_real 01914 irez = Inf(rez), 01915 srez = Sup(rez), 01916 iimz = Inf(imz), 01917 simz = Sup(imz); 01918 01919 lx_interval hxl(irez), hxu(srez), hyl(iimz), hyu(simz); 01920 01921 lx_real resxl, resxu, resyl, resyu; 01922 01923 bool bl = (sm_zero(iimz)) && (gr_zero(simz)), 01924 raxis = (eq_zero(iimz)) && (eq_zero(simz)); 01925 01926 // 01927 // 1st: check for singularities 01928 // 01929 if ( (irez<-1 && (bl || (sm_zero(iimz) && eq_zero(simz)))) || 01930 (srez >1 && (bl || (eq_zero(iimz) && gr_zero(simz)))) ) 01931 cxscthrow(STD_FKT_OUT_OF_DEF( 01932 "lx_cinterval asin(const lx_cinterval& z); z contains singularities.")); 01933 01934 // 01935 // 2nd: real part 01936 // 01937 if (sm_zero(iimz) && gr_zero(simz)) 01938 // z intersects [-1,1] 01939 { 01940 if (se_zero(irez)) 01941 resxl = Inf( asin(hxl) ); 01942 else 01943 resxl = Inf(Asin_beta(hxl,lx_interval( max(-iimz,simz) )) ); 01944 if (sm_zero(srez)) 01945 resxu = Sup(Asin_beta(hxu,lx_interval( max(-iimz,simz) )) ); 01946 else 01947 resxu = Sup( asin(hxu) ); 01948 } 01949 else 01950 { 01951 if ( (ge_zero(iimz) && ge_zero(irez)) || 01952 (se_zero(simz) && se_zero(irez)) ) 01953 // left boundary in quadrants I or III 01954 // min( Re( z ) ) in upper left corner 01955 resxl = Inf( Asin_beta(hxl,hyu) ); // Blomquist, 19.06.2005; 01956 else 01957 // left boundary in quadrants II or IV 01958 // min( Re( z ) ) in lower left corner 01959 resxl = Inf( Asin_beta(hxl,hyl) ); // Blomquist, 19.06.2005; 01960 if ( (ge_zero(iimz) && ge_zero(srez)) || (se_zero(simz) && se_zero(srez) ) ) 01961 // right boundary in quadrants I or III 01962 // max( Re( z ) ) in lower right corner 01963 resxu = Sup( Asin_beta(hxu,hyl) ); // Blomquist, 19.06.2005; 01964 else 01965 // right boundary in quadrants II or IV 01966 // max( Re( z ) ) in upper right corner 01967 resxu = Sup( Asin_beta(hxu,hyu) ); // Blomquist, 19.06.2005; 01968 } 01969 01970 // 01971 // 3rd: imaginary part 01972 // 01973 if (raxis) 01974 { // Interval argument is now a subset of the real axis. 01975 // Blomquist, 16.06.2005; 01976 if (sm_zero(srez)) resyl = Inf( ACOSH_f_aux( hxu, hyu )); 01977 else resyl = -Sup( ACOSH_f_aux( hxu, hyu )); 01978 if (gr_zero(irez)) resyu = -Inf( ACOSH_f_aux( hxl, hyu )); 01979 else resyu = Sup( ACOSH_f_aux( hxl, hyu )); 01980 } 01981 else 01982 if (se_zero(simz)) 01983 // z in lower half plane 01984 // min( Im( z ) ) in point with max |z| 01985 // max( Im( z ) ) in point with min |z| 01986 { 01987 if (irez < -srez) 01988 // most of z in quadrant III 01989 { 01990 resyl = -Sup( ACOSH_f_aux( hxl, hyl ) ); 01991 if (sm_zero(srez)) 01992 resyu = -Inf( ACOSH_f_aux( hxu, hyu ) ); 01993 else 01994 resyu = -Inf( ACOSH_f_aux( lx_interval(0),hyu ) ); 01995 } 01996 else 01997 // most of z in quadrant IV 01998 { 01999 resyl = -Sup( ACOSH_f_aux( hxu, hyl ) ); 02000 if (gr_zero(irez)) 02001 resyu = -Inf( ACOSH_f_aux( hxl, hyu ) ); 02002 else 02003 resyu = -Inf( ACOSH_f_aux( lx_interval(0),hyu ) ); 02004 } 02005 } 02006 else 02007 if (ge_zero(iimz)) 02008 // z in upper half plane 02009 // min( Im( z ) ) in point with min |z| 02010 // max( Im( z ) ) in point with max |z| 02011 { 02012 if( irez < -srez ) // if( irez + srez < 0.0 ) 02013 // most of z in quadrant II 02014 { 02015 resyu = Sup( ACOSH_f_aux( hxl, hyu ) ); 02016 if (sm_zero(srez)) 02017 resyl = Inf( ACOSH_f_aux( hxu, hyl ) ); 02018 else 02019 resyl = Inf( ACOSH_f_aux( lx_interval(0), hyl ) ); 02020 } 02021 else 02022 // most of z in quadrant I 02023 { 02024 resyu = Sup( ACOSH_f_aux( hxu, hyu ) ); 02025 if (gr_zero(irez)) 02026 resyl = Inf( ACOSH_f_aux( hxl, hyl ) ); 02027 else 02028 resyl = Inf( ACOSH_f_aux( lx_interval(0), hyl ) ); 02029 } 02030 } 02031 else 02032 // z intersects imaginary axes 02033 // min( Im( z ) ) in point in lower half plane with max |z| 02034 // max( Im( z ) ) in point in upper half plane with max |z| 02035 { 02036 if( irez < -srez ) // if( irez + srez < 0.0 ) 02037 // most of z in quadrants II and IV 02038 { 02039 resyl = - Sup( ACOSH_f_aux( hxl, hyl ) ); 02040 resyu = Sup( ACOSH_f_aux( hxl, hyu ) ); 02041 } 02042 else 02043 { 02044 resyl = - Sup( ACOSH_f_aux( hxu, hyl ) ); 02045 resyu = Sup( ACOSH_f_aux( hxu, hyu ) ); 02046 } 02047 } 02048 02049 res = lx_cinterval( lx_interval(resxl,resxu),lx_interval(resyl,resyu) ); 02050 stagprec = stagsave; 02051 res = adjust(res); 02052 return res; 02053 } 02054 // 02055 //-- end asin ----------------------------------------------------------------- 02056 02057 //----------------------------------------------------------------------------- 02058 //----------------------------- acos ------------------------------------------ 02059 02060 lx_interval BETA_xy(const lx_interval& x, const lx_interval& y) 02061 // Calculation of beta = x / ( sqrt[(x+1)^2 + y^2] + sqrt[(x-1)^2 + y^2]/2 ); 02062 // Only for the internal use for calculating acos(beta); 02063 { 02064 lx_interval res,xa; 02065 l_interval xl,yl; 02066 real exx,exy; 02067 int exxl,exyl; 02068 interval z; 02069 bool neg; 02070 02071 xa = x; 02072 xl = li_part(xa); yl = li_part(y); 02073 neg = Inf(xl)<0; 02074 if (neg) xa = -xa; 02075 exx = expo(xa); exxl = expo_gr(Inf(xl)); 02076 exy = expo(y); exyl = expo_gr(Inf(yl)); 02077 02078 if (exyl<-100000) // y = 0; 02079 if (Inf(xa)<1) res = xa; 02080 else res = 1.0; 02081 else // y != 0; 02082 if (exxl<-100000) // x = 0 02083 res = 0; 02084 else 02085 { 02086 z = interval(exy) + (exyl-exxl-2103); 02087 if (exx<Inf(z)) res = 0; 02088 else 02089 res = xa / f_aux_asin(xa,y); 02090 } 02091 if (Sup(res)>1) res = 1.0; 02092 02093 if (neg) res = -res; 02094 return res; 02095 } // BETA_xy(...) 02096 02097 lx_interval Asin_arg(const lx_interval& x, const lx_interval& y ) 02098 // Asin_arg calculats for point intervals x,y with 02099 // beta := 2*|x| / ( sqrt((x+1)^2+y^2) + sqrt((x-1)^2+y^2) ) 02100 // and delta := 1 - |beta| an inclusion of: 02101 // arcsin[ sqrt(delta)*sqrt(2-delta) ]. 02102 // The point interval x may be negative! 02103 // Blomquist, 06.06.2008; 02104 { 02105 const real c2 = -9007199254739968.0; 02106 lx_interval res,xa,F_xa,S1,S2,xm1,xp1,w_delta,T; 02107 lx_real Infxa; 02108 bool neg_x; 02109 xa = x; 02110 neg_x = Inf(x)<0; 02111 if (neg_x) xa = -xa; // Inf(xa) > 0 : 02112 Infxa = Inf(xa); 02113 02114 if (Infxa > 1.0) 02115 { // x > +1; 02116 if (y == 0.0) 02117 res = 0.0; 02118 else // y != 0; 02119 { 02120 F_xa = xa; 02121 times2pown(F_xa,c2); 02122 02123 if (Inf(abs(y)) > Inf(F_xa)) 02124 // y > x * 2^(-9007199254739968): 02125 // Now the calculation of w_delta according to (1) is 02126 // possible and additional the calculation of the total 02127 // argument of the asin-function: 02128 { 02129 xm1 = xa-1.0; 02130 if (Inf(xm1)<=0) 02131 xm1 = lx_interval(Inf( lx_interval(-2097,l_interval(1)) ), 02132 Sup(xm1)); 02133 xp1 = xa+1.0; 02134 S1 = sqrtx2y2(xm1,y); 02135 S2 = sqrtx2y2(xp1,y); 02136 w_delta = sqrt(2 + S1 + S2); 02137 w_delta = w_delta * sqrt(S1+xm1); 02138 w_delta = Sqrt2_lx_interval()*abs(y) / w_delta; 02139 T = (Sqrt2_lx_interval()-w_delta)*(Sqrt2_lx_interval()+w_delta); 02140 T = w_delta * sqrt(T); 02141 res = asin( T ); 02142 } 02143 else 02144 { 02145 // y <= x * 2^(-9007199254739968) 02146 if (Inf(xa) >= 2) 02147 { 02148 real exx,exxl,exy,exyl; 02149 interval z; 02150 double dbl; 02151 02152 exx = expo(xa); exy = expo(y); 02153 exxl = expo_gr(li_part(xa)); 02154 exyl = expo_gr(li_part(y)); 02155 02156 z = interval(exy)-interval(exx) + (exyl-exxl+1079); 02157 // 1079 = 1074 + 5; see documentation 02158 exx = Sup(z); 02159 if (exx<-Max_Int_R) 02160 exyl = -Max_Int_R; 02161 else // Sup(z) >= -9007199254740991 02162 { 02163 if (diam(z)==0) exyl = Sup(z); 02164 else 02165 { 02166 dbl = floor(_double(exx)); 02167 exyl = real(dbl) + 1; // exyl is integer! 02168 } 02169 } 02170 02171 res = lx_interval(lx_real(0.0), 02172 lx_real(exyl,l_real(minreal))); 02173 } 02174 else // 1 < Inf(xa) < 2; 02175 { 02176 T = sqrt(xa*(xa-1.0)); 02177 if (Inf(T) <= 0) 02178 { 02179 times2pown(xa,-2097); 02180 xa = sqrt(xa); 02181 times2pown(xa,3000); 02182 res = abs(y); 02183 times2pown(res,3001); 02184 res = res / xa; 02185 } 02186 else 02187 { 02188 res = abs(y); 02189 times2pown(res,1); 02190 res = res/T; 02191 } 02192 res = lx_interval(lx_real(0.0),Sup(res)); 02193 } 02194 } 02195 } 02196 } 02197 else 02198 if (Infxa == 1.0) 02199 { // x = 1; 02200 T = abs(y); 02201 res = 2 + T + sqrtx2y2(lx_interval(0,l_interval(2)),T); 02202 times2pown(T,1); 02203 T = T / res; 02204 res = asin(sqrt(T*(2-T))); 02205 } 02206 else 02207 { // 0.75 <= x < 1; 02208 if (y == 0.0) 02209 { 02210 T = sqrt1mx2(xa); 02211 res = asin(T); 02212 } 02213 else 02214 { 02215 T = 1 - xa; 02216 S1 = sqrtx2y2(T,y); 02217 S2 = sqrtx2y2(x+1,y); 02218 res = S1 + T; 02219 times2pown(res,1); 02220 T = 2 + S1 + S2; 02221 T = res / T; 02222 T = sqrt( T*(2-T) ); 02223 res = asin(T); 02224 } 02225 } 02226 return res; 02227 02228 } // Asin_arg 02229 02230 lx_interval Acos_beta(const lx_interval& x, const lx_interval& y) 02231 // Calculating the improved real part of acos(z); Blomquist 05.06.2005; 02232 // Re(acos(z)) = acos[ 2x / (sqrt[(x+1)^2+y^2] + sqrt[(x-1)^2+y^2]) ] 02233 { 02234 const real c1 = 0.75; 02235 lx_interval res(0),beta; 02236 beta = BETA_xy(x,y); 02237 if (Sup(beta)<c1) 02238 if (Sup(beta)<-c1) 02239 { // Improvement for beta --> -1 02240 res = Pi_lx_interval() - Asin_arg(x,y); 02241 } 02242 else res = acos(beta); // Normal calculation 02243 else // Sup(beta)>=c1 02244 res = Asin_arg(x,y); 02245 return res; 02246 } // Acos_beta(...) 02247 02248 // 02249 lx_cinterval acos(const lx_cinterval& z) throw() 02250 { 02251 lx_interval 02252 rez = Re(z), 02253 imz = Im(z); 02254 02255 lx_real 02256 irez = Inf(rez), 02257 srez = Sup(rez), 02258 iimz = Inf(imz), 02259 simz = Sup(imz); 02260 02261 lx_interval 02262 hxl(irez), hxu(srez), hyl(iimz), hyu(simz); 02263 02264 bool bl = iimz< 0.0 && simz>0.0, 02265 raxis = iimz==0.0 && simz==0; 02266 lx_real resxl, resxu, resyl, resyu; 02267 // 02268 // 1st: check for singularities 02269 // 02270 if( (irez<-1 && (bl || (iimz<0 && simz==0))) || 02271 (srez >1 && (bl || (iimz==0 && simz>0))) ) 02272 cxscthrow(STD_FKT_OUT_OF_DEF("lx_cinterval acos(const lx_cinterval& z); z contains singularities.")); 02273 // 02274 // 2nd: real part 02275 // 02276 // Blomquist, 05.02.2007; 02277 02278 if( iimz < 0.0 && simz > 0.0 ) 02279 // z intersects [-1,1] on the x-axis 02280 { 02281 if( irez <= 0.0 ) resxu = Sup( acos( hxl ) ); 02282 else resxu = Sup( Acos_beta(hxl,lx_interval( max(-iimz,simz) )) ); 02283 02284 if( srez < 0.0 ) 02285 resxl = Inf( Acos_beta(hxu,lx_interval(max(-iimz,simz))) ); 02286 else resxl = Inf( acos( hxu ) ); 02287 } 02288 else 02289 { 02290 if (irez<0 && srez>0) 02291 // z intersects the posizive or negative y-axis 02292 if (ge_zero(iimz)) 02293 { 02294 resxl = Inf( Acos_beta(hxu,hyl) ); 02295 resxu = Sup( Acos_beta(hxl,hyl) ); 02296 } 02297 else 02298 { 02299 resxl = Inf( Acos_beta(hxu,hyu) ); 02300 resxu = Sup( Acos_beta(hxl,hyu) ); 02301 } 02302 else 02303 { 02304 if ( (ge_zero(iimz) && ge_zero(irez)) || 02305 (se_zero(simz) && sm_zero(irez)) ) 02306 // left boundary in quadrants I or III 02307 // min( Re( z ) ) in lower right corner 02308 resxl = Inf( Acos_beta(hxu,hyl) ); 02309 else 02310 // left boundary in quadrants II or IV 02311 // min( Re( z ) ) in upper right corner 02312 resxl = Inf( Acos_beta(hxu,hyu) ); 02313 02314 if ( (ge_zero(iimz) && gr_zero(srez)) || 02315 (se_zero(simz) && se_zero(srez)) ) 02316 // right boundary in quadrants I or III 02317 // max( Re( z ) ) in upper left corner 02318 resxu = Sup( Acos_beta(hxl,hyu) ); 02319 else 02320 // right boundary in quadrants II or IV 02321 // max( Re( z ) ) in lower left corner 02322 resxu = Sup( Acos_beta(hxl,hyl) ); 02323 } 02324 } 02325 02326 // 02327 // 3rd: imaginary part 02328 // 02329 if (raxis) 02330 { // Interval argument is now a subset of the real axis. 02331 // Blomquist, 16.06.2005; 02332 if (srez<0.0) resyl = Inf( ACOSH_f_aux( hxu, hyu )); 02333 else resyl = -Sup( ACOSH_f_aux( hxu, hyu )); 02334 if (irez>0.0) resyu = -Inf( ACOSH_f_aux( hxl, hyu )); 02335 else resyu = Sup( ACOSH_f_aux( hxl, hyu )); 02336 } 02337 else 02338 if( simz <= 0.0 ) 02339 // z in lower half plane 02340 // min( Im( z ) ) in point with max |z| 02341 // max( Im( z ) ) in point with min |z| 02342 { 02343 if (irez < -srez) 02344 // most of z in quadrant III 02345 { 02346 resyl = -Sup( ACOSH_f_aux( hxl, hyl ) ); 02347 if( srez < 0.0 ) 02348 resyu = -Inf( ACOSH_f_aux( hxu, hyu ) ); 02349 else 02350 resyu = -Inf( ACOSH_f_aux( lx_interval(0), hyu ) ); 02351 } 02352 else 02353 // most of z in quadrant IV 02354 { 02355 resyl = -Sup( ACOSH_f_aux( hxu, hyl ) ); 02356 if( irez > 0.0 ) 02357 resyu = -Inf( ACOSH_f_aux( hxl, hyu ) ); 02358 else 02359 resyu = -Inf( ACOSH_f_aux( lx_interval(0), hyu ) ); 02360 } 02361 } 02362 else 02363 if( ge_zero(iimz) ) 02364 // z in upper half plane 02365 // min( Im( z ) ) in point with min |z| 02366 // max( Im( z ) ) in point with max |z| 02367 { 02368 if( irez < -srez ) // if( irez + srez < 0.0 ) 02369 // most of z in quadrant II 02370 { 02371 resyu = Sup( ACOSH_f_aux( hxl, hyu ) ); 02372 if( sm_zero(srez) ) 02373 resyl = Inf( ACOSH_f_aux( hxu, hyl ) ); 02374 else 02375 resyl = Inf( ACOSH_f_aux( lx_interval(0), hyl ) ); 02376 } 02377 else 02378 // most of z in quadrant I 02379 { 02380 resyu = Sup( ACOSH_f_aux( hxu, hyu ) ); 02381 if( irez > 0.0 ) 02382 resyl = Inf( ACOSH_f_aux( hxl, hyl ) ); 02383 else 02384 resyl = Inf( ACOSH_f_aux( lx_interval(0), hyl ) ); 02385 } 02386 } 02387 else 02388 // z intersects imaginary axes 02389 // min( Im( z ) ) in point in lower half plane with max |z| 02390 // max( Im( z ) ) in point in upper half plane with max |z| 02391 { 02392 if( irez < -srez ) // if( irez + srez < 0.0 ) 02393 // most of z in quadrants II and IV 02394 { 02395 resyl = -Sup( ACOSH_f_aux( hxl, hyl ) ); 02396 resyu = Sup( ACOSH_f_aux( hxl, hyu ) ); 02397 } 02398 else 02399 { 02400 resyl = -Sup( ACOSH_f_aux( hxu, hyl ) ); 02401 resyu = Sup( ACOSH_f_aux( hxu, hyu ) ); 02402 } 02403 } 02404 02405 return lx_cinterval( lx_interval( resxl, resxu ),-lx_interval( resyl, resyu ) ); 02406 02407 } // acos(...) 02408 02409 //-- end acos ----------------------------------------------------------------- 02410 02411 //-- asinh -------------------------------------------------------------------- 02412 // 02413 lx_cinterval asinh( const lx_cinterval& z ) throw() 02414 // 02415 // asinh( z ) = i * asin( -i * z ) 02416 // 02417 { 02418 lx_cinterval res = asin( lx_cinterval( Im(z), -Re(z) ) ); 02419 return lx_cinterval( -Im(res), Re(res) ); 02420 } 02421 // 02422 //-- end asinh ---------------------------------------------------------------- 02423 02424 //-- acosh -------------------------------------------------------------------- 02425 // 02426 lx_cinterval acosh( const lx_cinterval& z ) throw() 02427 // 02428 // acosh( z ) = i * acos( z ) = +/- i * ( pi / 2 - asin( z ) ) 02429 // 02430 { 02431 lx_interval 02432 rez = Re(z), 02433 imz = Im(z); 02434 02435 lx_real 02436 irez = Inf(rez), 02437 srez = Sup(rez), 02438 iimz = Inf(imz), 02439 simz = Sup(imz); 02440 02441 lx_interval hxl(irez), hxu(srez), hyl(iimz), hyu(simz); 02442 02443 lx_real resxl, resxu, resyl, resyu; 02444 02445 // 02446 // 1st: check for singularities 02447 // 02448 if( ( se_zero(iimz) && ge_zero(simz) ) && ( irez < 1.0 ) ) 02449 cxscthrow(STD_FKT_OUT_OF_DEF( 02450 "lx_cinterval acosh( const lx_cinterval& z ); z contains singularities.")); 02451 // With this restriction the complex interval argument and the real axis 02452 // must not have any common point, if irez < +1; 02453 // So for example the negative real axis must not be touched from above if 02454 // irez<1, although this should be possible if the principal branch is 02455 // considered! So the above restriction is too widely in 02456 // some cases; Blomquist, 21.06.2005; 02457 // 02458 // 2nd: z in upper half plane (or on the real axis) 02459 // acosh( z ) = + i * ( pi / 2 - asin( z ) ) 02460 // 02461 if( gr_zero(iimz) ) 02462 { 02463 lx_cinterval res = acos(z); 02464 return lx_cinterval( -Im(res),Re(res) ); 02465 } 02466 // 02467 // 3rd: z in lower half plane 02468 // acosh( z ) = - i * ( pi / 2 - asin( z ) ) 02469 // 02470 if( sm_zero(simz) ) 02471 { 02472 // cinterval res = HALFPI() - asin( z ); 02473 lx_cinterval res = acos(z); // Blomquist, 14.06.2005 02474 return lx_cinterval( Im(res), -Re(res) ); 02475 } 02476 // 02477 // z intersects [1,infinity) 02478 // 02479 // real part 02480 // minimum on the left on real axes, maximum in lower or upper right corner 02481 // 02482 resxl = Inf( acosh( hxl ) ); 02483 lx_interval ytilde( max( -iimz, simz ) ); 02484 // resxu = Sup( acosh( f_aux_asin( hxu, ytilde ) ) ); 02485 resxu = Sup( ACOSH_f_aux(hxu,ytilde) ); // Blomquist, 14.06.2005; 02486 // 02487 // imaginary part 02488 // minimum in lower left corner, maximum in upper left corner 02489 // 02490 // resyl = -Sup( acos( hxl / f_aux_asin( hxl, hyl ) ) ); 02491 resyl = -Sup( Acos_beta(hxl,hyl) ); // Blomquist, 14.06.2005; 02492 // resyu = Sup( acos( hxl / f_aux_asin( hxl, hyu ) ) ); 02493 resyu = Sup( Acos_beta(hxl,hyu) ); 02494 02495 return lx_cinterval(lx_interval( resxl, resxu ),lx_interval( resyl, resyu )); 02496 } 02497 // 02498 //-- end acosh ---------------------------------------------------------------- 02499 02500 02501 //-- atan --------------------------------------------------------------------- 02502 // 02503 02504 bool sign_test(const lx_interval& x, int s_org) 02505 // Only for internal use by function re_atan(...) 02506 { 02507 bool bl1,bl2,alter; 02508 if (diam(li_part(x))>0) 02509 { 02510 bl1 = eq_zero(Sup(x)) && (s_org==1); 02511 bl2 = eq_zero(Inf(x)) && (s_org==-1); 02512 alter = bl1 || bl2; 02513 } 02514 else 02515 alter = sign(Sup(li_part(x))) != s_org; 02516 return alter; 02517 } // sign_test(...) 02518 02519 void re_atan( const lx_interval& y, lx_interval& x, lx_interval& res ) 02520 // Calculating an inclusion res of 1 - y^2 - x^2; 02521 // x is always a point interval and y =[y1,y2], with y1 <= y2; 02522 // Blomquist, 10.06.2008; 02523 { 02524 const real c1 = 4503599627369982.0; 02525 lx_interval ya(abs(y)), One(lx_interval(0,l_interval(1))), xa; 02526 lx_real R,U; 02527 real ex_x,ex_xl,ex_y,ex_yl,s; 02528 int signx; 02529 bool alter; 02530 02531 ex_x = expo(x); 02532 ex_xl = expo_gr(Inf(li_part(x))); 02533 ex_y = expo(y); 02534 ex_yl = expo_gr(Sup(li_part(ya))); 02535 ya = y; 02536 signx = sign(Inf(li_part(x))); 02537 02538 if (ex_xl < -1000000) 02539 res = 1.0; 02540 else // x <> 0: 02541 if (ex_yl < -1000000) // y == [0,0] and x <> [0,0]; 02542 if(ex_x > c1) 02543 { 02544 res = 1/x - x; 02545 x = -1; // Eigentlich sollte hier x = +1 stehen ??? 02546 } 02547 else // y == [0,0] and x^2 without overflow 02548 { 02549 res = (1-x)*(1+x); 02550 if (Sup(res)>1) // to avoid an overestimation 02551 res = lx_interval(Inf(res),lx_real(1.0)); 02552 } 02553 else // y <> [0,0] and x <> [0,0]; 02554 if (ex_x>c1) // |x| --> infty; 02555 if (ex_y>c1) // |x| --> infty and |y| --> infty; 02556 { 02557 s = c1 - max(ex_x,ex_y); // s < 0; and s is an integer value! 02558 times2pown_neg(One,s); // 2^s 02559 times2pown_neg(x,s); // x*2^s 02560 times2pown(ya,s); // y*2^s 02561 res = (One-x)*(One+x) - sqr(ya); 02562 times2pown_neg(x,s); // x*2^2s 02563 } 02564 else // |x| --> infty and |y| not--> infty; 02565 { 02566 res = sqr(y); xa = abs(x); 02567 if (res == 1.0) 02568 res = xa; 02569 else 02570 { 02571 res = res - 1.0; // res = y^2 - 1; 02572 R = Sup(abs(res)); 02573 ex_y = expo(R); 02574 ex_yl = expo_gr(lr_part(R)); 02575 ex_xl = ex_xl - 1051; 02576 if (ex_y+ex_yl < 2*(ex_x+ex_xl)) 02577 res = xa * 02578 lx_interval(Inf(One_m_lx_interval()), 02579 Sup(One_p_lx_interval())); 02580 else 02581 res = res/xa + xa; 02582 } 02583 x = (Sup(x)>0)? -1.0 : 1.0; 02584 } 02585 else // |x| not to infty 02586 if (ex_y>c1) // |y| --> infty; 02587 { 02588 s = c1 - ex_y; // 02589 res = (x-1)*(x+1); 02590 R = Sup( abs(res) ); 02591 times2pown_neg(x,s); // x*2^s 02592 times2pown(ya,s); // y*2^s 02593 if (eq_zero(R)) 02594 res = sqr(ya); 02595 else 02596 { // R > 0; 02597 ya = sqr(ya); // ya = (y*2^s)^2; 02598 times2pown_neg(res,s); 02599 times2pown_neg(res,s); 02600 // res = {(x-1)*(x+1)}*2^2s; 02601 res = res + ya; 02602 } 02603 times2pown_neg(x,s); // (x*2^s)*2^s = x*2^(2s) 02604 x = -x; 02605 } 02606 else // |x|,|y| not to +infty 02607 if (ex_y<-c1) // |y| ---> 0, y<>[0,0]; 02608 if (abs(x)==1) 02609 { 02610 times2pown(x,9007199254738894.0); 02611 x = -x; 02612 times2pown(ya,4503599627369447.0); 02613 res = sqr(ya); 02614 } 02615 else 02616 { 02617 res = (x-1)*(x+1) + sqr(ya); 02618 x = -x; 02619 } 02620 else // |x|,|y| not to +infty and |y| not to 0 02621 if (ex_x < -c1) // now: |x| --> 0 02622 { 02623 res = (y-1)*(y+1); 02624 if (res == 0.0) // alpha = -1/x; 02625 { 02626 x = -1; 02627 res = x; 02628 } 02629 else 02630 { 02631 res += sqr(x); 02632 x = -x; 02633 } 02634 } 02635 else // x^2 and y2 can now be evaluated without any 02636 // integer overflow! 02637 res = (1-x)*(1+x) - sqr(y); 02638 alter = sign_test(x,signx); 02639 if (alter) 02640 { 02641 x = -x; 02642 res = -res; 02643 } 02644 } // re_atan 02645 02646 void re_vert( const lx_real& x, const lx_interval& hx, 02647 const lx_real& rew_inf, const lx_real& rew_sup, 02648 lx_real& resxl, lx_real& resxu ) 02649 // 02650 // Subroutine of analytic inverse tangent function. 02651 // Evaluate real part on a vertical boundary. 02652 // 02653 { 02654 if( eq_zero(x) ) 02655 // singularities have been handled before, hence Re( w ) > 0 02656 { 02657 resxl = 0.0; 02658 resxu = 0.0; 02659 } 02660 else 02661 { 02662 lx_interval hx2(hx),Pid2,Pid4; 02663 Pid4 = Pid4_lx_interval(); 02664 times2pown(hx2,1); 02665 if( gr_zero(x) ) 02666 // w in quadrants I and/or II 02667 // atan is the inverse function of tan(t), t in (-pi/2,pi/2). 02668 { 02669 resxl = gr_zero(rew_sup)? Inf( Atan( hx2,rew_sup )/2.0 ) 02670 : ( sm_zero(rew_sup)? Inf( (Atan( hx2,rew_sup ) + Pi_lx_interval() )/2.0 ) 02671 : Inf( Pid4 ) ); 02672 02673 resxu = gr_zero(rew_inf)? Sup( Atan( hx2,rew_inf )/2.0 ) 02674 : ( sm_zero(rew_inf)? Sup( (Atan( hx2,rew_inf ) + Pi_lx_interval())/2.0 ) 02675 : Sup( Pid4 ) ); 02676 } 02677 else 02678 // w in quadrants III and/or IV 02679 { 02680 resxl = sm_zero(rew_inf)? Inf( (Atan( hx2,rew_inf ) - Pi_lx_interval())/2.0 ) 02681 : ( gr_zero(rew_inf)? Inf( Atan( hx2,rew_inf )/2.0 ) 02682 : -Sup( Pid4 ) ); 02683 resxu = sm_zero(rew_sup)? Sup( (Atan( hx2,rew_sup ) - Pi_lx_interval())/2.0 ) 02684 : ( gr_zero(rew_sup)? Sup( Atan( hx2,rew_sup )/2.0 ) 02685 : -Inf( Pid4 ) ); 02686 } 02687 } 02688 } // re_vert 02689 02690 lx_interval T_atan(const lx_real& x) 02691 // Calculating an inclusion of 02692 // ln[ 1+2/(sqrt(1+x^2)-1) ]. 02693 // x will be handeld as a point interval [x]=[x,x], with x>0. 02694 // Blomquist, 10.06.2008; 02695 { 02696 const real c1 = 4503599627367859.0; 02697 lx_interval res, 02698 ix(x); // ix is point interval with x>0; 02699 real ex_x(expo(ix)); 02700 02701 if (ex_x<-c1) 02702 { 02703 res = sqrt1px2(ix) + 1; 02704 times2pown(res,1); 02705 res += sqr(ix); // res = 2(1+sqrt(1+x^2)) + x^2; 02706 ix = ln(ix); 02707 times2pown(ix,1); // ix = 2*ln(x); 02708 res = ln(res) - ix; 02709 } 02710 else // -c1 <= ex_x 02711 if (ex_x<1150) // -c1 <= ex_x < 1150 02712 res = lnp1(2/sqrtp1m1(sqr(ix))); 02713 else res = lnp1(2/(sqrt1px2(ix)-1)); 02714 02715 return res; 02716 } // T_atan 02717 02718 lx_interval Q_atan(const lx_interval& x, const lx_interval& y) 02719 { 02720 // x: abs(Re(z)); So x is a real interval x = [x1,x2], with 0<=x1<x2. 02721 // y: Inf(Im(z)); So y is a point interval, with y >= 0. 02722 // Q_atan returns an inclusion of ln[1 + 4y/(x^2+(1-y)^2)] 02723 // Tested in detail; Blomquist, 13.06.2008; 02724 02725 const real c1 = 4503599627367859.0, 02726 c2 = 9007199254740990.0, 02727 c3 = 4503599627370495.0; 02728 const lx_real S = lx_real(-c2,MinReal); 02729 02730 double Dbl; 02731 int r; 02732 lx_real R; 02733 lx_interval res(0.0); 02734 real ex_y(expo(y)), ex_yl(expo_gr(li_part(y))), 02735 ex_x(expo(Sup(x))),ex_xl,exx,exxl,dbl,s,up,exy,exyl; 02736 interval dbli,z; 02737 02738 ex_xl = expo_gr(lr_part(Sup(x))); 02739 if (ex_xl<-100000) ex_x = 0; 02740 if (ex_yl>-100000) // y = [y,y], y>0; 02741 { 02742 if (y==1.0) 02743 { // y = 1; 02744 if (ex_x < -c1) 02745 { 02746 res = ln(x); 02747 times2pown(res,1); // res = 2*ln(x) 02748 res = ln(4+sqr(x)) - res; 02749 } 02750 else 02751 if (ex_x > c1) 02752 { 02753 R = Inf(x); 02754 exx = expo(R); 02755 exxl = expo_gr(lr_part(R)); 02756 if (514-exxl < -c3+exx) 02757 res = lx_interval(lx_real(0),S); 02758 else 02759 if (3 - exxl >= -c3 + exx) // dbl >= -c2 02760 { 02761 dbl = -2*(exx+exxl-3); // without any rounding! 02762 res = lx_interval(lx_real(0),lx_real(dbl,l_real(1))); 02763 } 02764 else // dbl < -c2 02765 { 02766 s = c2-exx; s = (s - exx) -2*exxl + 6; 02767 r = (int) _double(s); 02768 res = lx_interval(lx_real(0),lx_real(-c2,comp(0.5,r+1))); 02769 } 02770 } 02771 else // x^2 now without integer overflow 02772 res = lnp1(sqr(2/x)); 02773 } 02774 else // Now: y>0 and y != [1,1] and 0<=x1<=x2; 02775 { 02776 if (ex_x>c1 || ex_y>c1 ) 02777 { 02778 R = Inf(x); 02779 exx = expo(R); exxl = expo_gr(lr_part(R)); 02780 R = Inf(y-1.0); // R != 0: 02781 exy = expo(R); exyl = expo_gr(lr_part(R)); 02782 02783 if (exxl<-100000) // x1 = 0; 02784 dbli = 2*( interval(exy) + (exyl-2) ); 02785 else 02786 { 02787 z = 2*( interval(exx) + (exxl-2) ); 02788 dbli = 2*( interval(exy) + (exyl-2) ); 02789 if (Sup(z) > Sup(dbli)) dbli = z; 02790 } 02791 dbli = (interval(exy) - dbli) + (exyl + 3); 02792 dbl = Sup(dbli); 02793 // Now it holds: Sup{4y/(x^2+(1-y)^2)} < 2^dbl; 02794 // and this upper bound 2^dbl is guaranteed, 02795 // even if dbl is not an integer value! 02796 up = Inf( interval(-c2) - 1022); 02797 if (dbl<up) 02798 res = lx_interval(lx_real(0),S); 02799 else // dbl >= up 02800 if (dbl>=-c2) 02801 { 02802 if (!(Is_Integer(dbl))) 02803 { // rounding to the next greater integer value. 02804 Dbl = floor(_double(dbl)) + 1; // Dbl: integer value 02805 dbl = (real) Dbl; 02806 } 02807 res = lx_interval(lx_real(0),lx_real(dbl,l_real(1))); 02808 } 02809 else 02810 { 02811 s = Sup(interval(c2) + dbl); 02812 if ( !(Is_Integer(s)) ) 02813 { // rounding to the next greater integer value. 02814 Dbl = floor(_double(s)) + 1; 02815 r = (int) Dbl; 02816 } 02817 else r = (int) _double(s); 02818 res = lx_interval(lx_real(0),lx_real(-c2,comp(0.5,r+1))); 02819 } 02820 } 02821 else // Now it holds: y>0 and y!=1 and 0<=x1<=x2; 02822 // and: x^2, (1-y)^2 produce no overflow for 02823 // x2, |1-y| --> +infty; 02824 // furthermore |1-y| produces no integer overflow 02825 // for y --> +1. 02826 { 02827 res = y; 02828 times2pown(res,2); // res = 4*y; 02829 res = res/(sqr(x) + sqr(1-y)); 02830 res = lnp1(res); // res: inclusion of Q(x,y); 02831 } 02832 } 02833 } 02834 return res; 02835 } // Q_atan 02836 02837 lx_cinterval atan( const lx_cinterval& z ) throw() 02838 { 02839 int stagsave = stagprec, 02840 stagmax = 30; 02841 if (stagprec>stagmax) stagprec = stagmax; 02842 02843 lx_cinterval res; 02844 lx_interval 02845 rez = Re(z), 02846 imz = Im(z); 02847 02848 lx_real 02849 irez = Inf(rez), 02850 srez = Sup(rez), 02851 iimz = Inf(imz), 02852 simz = Sup(imz); 02853 02854 lx_interval // all theses intervals are point intervals! 02855 hxl(irez), hxu(srez), hyl(iimz), hyu(simz); 02856 02857 lx_real 02858 resxl, resxu, resyl, resyu; 02859 // 02860 // 1st: check for singularities 02861 // 02862 if( (se_zero(irez) && ge_zero(srez)) && ( iimz <= -1.0 || simz >= 1.0 ) ) 02863 cxscthrow(STD_FKT_OUT_OF_DEF("lx_cinterval atan( const lx_cinterval& z ); points of the branch cuts are not allowed in z.")); 02864 // 02865 // 2nd: real part 02866 // Re( atan( z ) ) = Arg( w ) / 2, where w = 1 - x^2 - y^2 + i * 2x ) 02867 // 02868 // evaluate atan on vertical boundaries 02869 // 02870 lx_interval 02871 // y_sqr = sqr( imz ), 02872 // rew_l = (1 - y_sqr) - sqr( hxl ), 02873 // Blomquist; before: rew_l = 1 - sqr(hxl) - y_sqr, 02874 // rew_u = (1 - y_sqr) - sqr( hxu ); 02875 // Blomquist; before: rew_u = 1 - sqr(hxu) - y_sqr; 02876 rew_l, rew_u; 02877 // ------------------------------ Blomquist --------------------------------- 02878 // ---------- Improvements for Im(z) = [1,1] or Im(z) = [-1,-1] 02879 bool sqrImz_1 = (iimz==simz) && (iimz==1.0 || iimz==-1.0); 02880 // Test for Im(z) = [1,1] or [-1,-1] 02881 02882 if (sqrImz_1) 02883 { 02884 rew_l = -abs(hxl); hxl = lx_interval(0,l_interval(sign(irez))); 02885 rew_u = -abs(hxu); hxu = lx_interval(0,l_interval(sign(srez))); 02886 } 02887 else 02888 { 02889 // rew_l = (1 - sqr( imz )) - sqr( hxl ); 02890 re_atan(imz,hxl,rew_l); 02891 // rew_u = (1 - sqr( imz )) - sqr( hxu ); 02892 re_atan(imz,hxu,rew_u); 02893 } 02894 // ------------------------------ Blomquist; 22.02.05; -------------------- 02895 02896 // 02897 // left boundary 02898 // 02899 lx_real rew_inf = Inf( rew_l ); 02900 lx_real rew_sup = Sup( rew_l ); 02901 re_vert( irez, hxl, rew_inf, rew_sup, resxl, resxu ); 02902 02903 // 02904 // right boundary 02905 // 02906 rew_inf = Inf( rew_u ); 02907 rew_sup = Sup( rew_u ); 02908 lx_real res_l, res_u; 02909 re_vert( srez, hxu, rew_inf, rew_sup, res_l, res_u ); 02910 02911 if (res_l<resxl) resxl = res_l; 02912 if (res_u>resxu) resxu = res_u; 02913 02914 // 02915 // look for extremal values on horizontal boundaries 02916 // since atan( x+iy ) = atan( x-iy ), 02917 // intersections can be considered in the upper half plane 02918 // 02919 lx_real abs_y_min = Inf( abs( imz ) ); 02920 if( abs_y_min > 1.0 ) 02921 { 02922 lx_interval 02923 abs_hyl = lx_interval( abs_y_min ), 02924 // abs_hxl = sqrt( sqr( abs_hyl ) - 1.0 ); 02925 abs_hxl = sqrtx2m1(abs_hyl); // Blomquist; 02926 02927 if( Sup( abs_hxl ) > irez && Inf( abs_hxl ) < srez ) 02928 // extremal curve intersects lower boundary of x+i|y| in quadrant I 02929 // intersection in Q I or Q IV: update minimum 02930 // resxl = inf( atan( abs_y_min / abs_hxl ) ) / 2.0; 02931 resxl = Inf( (Pi_lx_interval() - atan( 1.0 / abs_hxl ))/2.0 ); 02932 else 02933 if ( -Inf( abs_hxl ) > irez && -Sup( abs_hxl ) < srez ) 02934 // extremal curve intersects lower boundary of x+i|y| in quadrant II 02935 // intersection in Q II or Q III: update maximum 02936 resxu = Sup( (atan( 1.0 / abs_hxl ) - Pi_lx_interval())/2.0 ); 02937 } 02938 // 3rd: imaginary part 02939 // Im( atan( z ) ) = +/- Ln( 1 +/- 4y/( x^2 + (1 -/+ y)^2 ) ) / 4 02940 // 02941 // evaluate atan on horizontal boundaries 02942 lx_interval 02943 abs_rez = abs(rez), 02944 im_atan_l, im_atan_u; 02945 02946 if ( sm_zero(iimz) ) 02947 // im_atan_l = -ln( 1 - 4 * hyl / ( x_sqr + sqr( 1 + hyl ) ) ) / 4.0; 02948 // im_atan_l = -lnp1(-4 * hyl / ( x_sqr + sqr( 1 + hyl ) )) / 4.0; // Blomquist 02949 im_atan_l = -Q_atan(abs_rez,-hyl); // Blomquist 02950 else 02951 // im_atan_l = ln( 1 + 4 * hyl / ( x_sqr + sqr( 1 - hyl ) ) ) / 4.0; 02952 im_atan_l = Q_atan(abs_rez,hyl); // Blomquist 02953 times2pown(im_atan_l,-2); // Division by 4 02954 if ( sm_zero(simz) ) 02955 // im_atan_u = -ln( 1 - 4 * hyu / ( x_sqr + sqr( 1 + hyu ) ) ) / 4.0; 02956 // im_atan_u = -lnp1(-4 * hyu / ( x_sqr + sqr( 1 + hyu ) ) ) / 4.0; // Blomquist 02957 im_atan_u = -Q_atan(abs_rez,-hyu); // Blomquist 02958 else 02959 // im_atan_u = ln( 1 + 4 * hyu / ( x_sqr + sqr( 1 - hyu ) ) ) / 4.0; 02960 im_atan_u = Q_atan(abs_rez,hyu); // Blomquist 02961 times2pown(im_atan_u,-2); // Division by 4 02962 resyl = min( Inf( im_atan_l ), Inf( im_atan_u ) ); 02963 resyu = max( Sup( im_atan_l ), Sup( im_atan_u ) ); 02964 // 02965 // look for extremal values on vertical boundaries, 02966 // if vertical boundaries intersect extremal curves 02967 // 02968 lx_real abs_x_min = Inf( abs( rez ) ); 02969 lx_interval 02970 x_extr = lx_interval( abs_x_min ), 02971 // y_extr = sqrt( 1.0 + sqr( x_extr ) ); 02972 y_extr = sqrt1px2(x_extr); // Blomquist; 02973 02974 if( Inf( y_extr ) < simz && Sup( y_extr ) > iimz ) 02975 // extremal curve intersects left boundary of |x|+iy in quadrant I 02976 // update maximum 02977 // resyu = Sup( ln( 1 + 4 * y_extr / ( sqr( x_extr ) + sqr( 1 - y_extr ) ) ) ) / 4.0; 02978 // resyu = Sup( T_atan(abs_x_min)/4.0 ); // Blomquist 02979 { 02980 rez = T_atan(abs_x_min); 02981 times2pown(rez,-2); 02982 resyu = Sup(rez); 02983 } 02984 02985 if( -Sup(y_extr) < simz && -Inf(y_extr) > iimz ) 02986 // extremal curve intersects left boundary of |x|+iy in quadrant IV 02987 // update minimum 02988 // resyl = -Sup( ln( 1 + 4 * y_extr / ( sqr( x_extr ) + sqr( 1 - y_extr ) ) ) ) / 4.0; 02989 // resyl = -Sup( T_atan(abs_x_min)/4.0 ); // Blomquist 02990 { 02991 rez = T_atan(abs_x_min); 02992 times2pown(rez,-2); 02993 resyl = -Sup(rez); 02994 } 02995 02996 res = lx_cinterval( lx_interval(resxl,resxu),lx_interval(resyl,resyu) ); 02997 stagprec = stagsave; 02998 res = adjust(res); 02999 return res; 03000 } 03001 // 03002 //-- end atan ----------------------------------------------------------------- 03003 03004 lx_cinterval acot( const lx_cinterval& z ) throw() 03005 { 03006 int stagsave = stagprec, 03007 stagmax = 30; 03008 if (stagprec>stagmax) stagprec = stagmax; 03009 03010 lx_cinterval res; 03011 lx_interval 03012 rez = Re(z), 03013 imz = Im(z); 03014 03015 lx_real 03016 irez = Inf(rez), 03017 srez = Sup(rez), 03018 iimz = Inf(imz), 03019 simz = Sup(imz); 03020 03021 lx_interval // all theses intervals are point intervals! 03022 hxl(irez), hxu(srez), hyl(iimz), hyu(simz); 03023 03024 lx_real 03025 resxl, resxu, resyl, resyu; 03026 // 03027 // 1st: check for singularities 03028 // 03029 if( (se_zero(irez) && ge_zero(srez)) && ( iimz <= 1.0 || simz >= -1.0 ) ) 03030 cxscthrow(STD_FKT_OUT_OF_DEF("lx_cinterval acot( const lx_cinterval& z ); points of the branch cuts are not allowed in z.")); 03031 // 03032 // 2nd: real part 03033 // Re( atan( z ) ) = Arg( w ) / 2, where w = 1 - x^2 - y^2 + i * 2x ) 03034 // 03035 // evaluate atan on vertical boundaries 03036 // 03037 lx_interval 03038 // y_sqr = sqr( imz ), 03039 // rew_l = (1 - y_sqr) - sqr( hxl ), // Blomquist; before: rew_l = 1 - sqr(hxl) - y_sqr, 03040 // rew_u = (1 - y_sqr) - sqr( hxu ); // Blomquist; before: rew_u = 1 - sqr(hxu) - y_sqr; 03041 rew_l, rew_u; 03042 03043 // ------------------------------ Blomquist --------------------------------- 03044 // ---------- Improvements for Im(z) = [1,1] or Im(z) = [-1,-1] 03045 bool sqrImz_1 = (iimz==simz) && (iimz==1.0 || iimz==-1.0); 03046 // Test for Im(z) = [1,1] or [-1,-1] 03047 03048 if (sqrImz_1) 03049 { 03050 rew_l = abs(hxl); hxl = lx_interval(0,l_interval(sign(irez))); 03051 rew_u = abs(hxu); hxu = lx_interval(0,l_interval(sign(srez))); 03052 } 03053 else 03054 { 03055 // rew_l = (1 - sqr( imz )) - sqr( hxl ); 03056 re_atan(imz,hxl,rew_l); 03057 rew_l = -rew_l; 03058 // rew_u = (1 - sqr( imz )) - sqr( hxu ); 03059 re_atan(imz,hxu,rew_u); 03060 rew_u = -rew_u; 03061 } 03062 // ------------------------------ Blomquist; 22.02.05; -------------------- 03063 03064 // 03065 // left boundary 03066 // 03067 lx_real rew_inf = Inf( rew_l ); 03068 lx_real rew_sup = Sup( rew_l ); 03069 re_vert( irez, hxl, rew_inf, rew_sup, resxl, resxu ); 03070 03071 // 03072 // right boundary 03073 // 03074 rew_inf = Inf( rew_u ); 03075 rew_sup = Sup( rew_u ); 03076 lx_real res_l, res_u; 03077 re_vert( srez, hxu, rew_inf, rew_sup, res_l, res_u ); 03078 03079 if (res_l<resxl) resxl = res_l; 03080 if (res_u>resxu) resxu = res_u; 03081 03082 // 03083 // look for extremal values on horizontal boundaries 03084 // since atan( x+iy ) = atan( x-iy ), 03085 // intersections can be considered in the upper half plane 03086 // 03087 lx_real abs_y_min = Inf( abs( imz ) ); 03088 03089 if( abs_y_min > 1.0 ) 03090 { 03091 lx_interval 03092 abs_hyl = lx_interval( abs_y_min ), 03093 // abs_hxl = sqrt( sqr( abs_hyl ) - 1.0 ); 03094 abs_hxl = sqrtx2m1(abs_hyl); // Blomquist; 03095 03096 if( Sup( abs_hxl ) > irez && Inf( abs_hxl ) < srez ) 03097 // extremal curve intersects lower boundary of x+i|y| in quadrant I 03098 // intersection in Q I or Q IV: update maximum 03099 // resxl = inf( atan( abs_y_min / abs_hxl ) ) / 2.0; 03100 resxu = Sup( atan( 1.0 / abs_hxl ) / 2.0 ); 03101 if ( -Inf( abs_hxl ) > irez && -Sup( abs_hxl ) < srez ) 03102 // extremal curve intersects lower boundary of x+i|y| in quadrant II 03103 // intersection in Q II or Q III: update minimum 03104 resxl = -Sup( atan( 1.0 / abs_hxl ) /2.0 ); 03105 } 03106 03107 // 3rd: imaginary part 03108 // Im( atan( z ) ) = +/- Ln( 1 +/- 4y/( x^2 + (1 -/+ y)^2 ) ) / 4 03109 // 03110 // evaluate atan on horizontal boundaries 03111 lx_interval 03112 abs_rez = abs(rez), 03113 im_atan_l, im_atan_u; 03114 03115 if ( sm_zero(iimz) ) 03116 // im_atan_l = -ln( 1 - 4 * hyl / ( x_sqr + sqr( 1 + hyl ) ) ) / 4.0; 03117 // im_atan_l = -lnp1(-4 * hyl / ( x_sqr + sqr( 1 + hyl ) )) / 4.0; // Blomquist 03118 im_atan_l = -Q_atan(abs_rez,-hyl); // Blomquist 03119 else 03120 // im_atan_l = ln( 1 + 4 * hyl / ( x_sqr + sqr( 1 - hyl ) ) ) / 4.0; 03121 im_atan_l = Q_atan(abs_rez,hyl); // Blomquist 03122 times2pown(im_atan_l,-2); // Division by 4 03123 03124 if ( sm_zero(simz) ) 03125 // im_atan_u = -ln( 1 - 4 * hyu / ( x_sqr + sqr( 1 + hyu ) ) ) / 4.0; 03126 // im_atan_u = -lnp1(-4 * hyu / ( x_sqr + sqr( 1 + hyu ) ) ) / 4.0; // Blomquist 03127 im_atan_u = -Q_atan(abs_rez,-hyu); // Blomquist 03128 else 03129 // im_atan_u = ln( 1 + 4 * hyu / ( x_sqr + sqr( 1 - hyu ) ) ) / 4.0; 03130 im_atan_u = Q_atan(abs_rez,hyu); // Blomquist 03131 times2pown(im_atan_u,-2); // Division by 4 03132 03133 resyl = min( Inf( im_atan_l ), Inf( im_atan_u ) ); 03134 resyu = max( Sup( im_atan_l ), Sup( im_atan_u ) ); 03135 // 03136 // look for extremal values on vertical boundaries, 03137 // if vertical boundaries intersect extremal curves 03138 // 03139 lx_real abs_x_min = Inf( abs( rez ) ); 03140 lx_interval 03141 x_extr = lx_interval( abs_x_min ), 03142 // y_extr = sqrt( 1.0 + sqr( x_extr ) ); 03143 y_extr = sqrt1px2(x_extr); // Blomquist; 03144 03145 if( Inf( y_extr ) < simz && Sup( y_extr ) > iimz ) 03146 // extremal curve intersects left boundary of |x|+iy in quadrant I 03147 // update maximum 03148 // resyu = Sup( ln( 1 + 4 * y_extr / ( sqr( x_extr ) + sqr( 1 - y_extr ) ) ) ) / 4.0; 03149 // resyu = Sup( T_atan(abs_x_min)/4.0 ); // Blomquist 03150 { 03151 rez = T_atan(abs_x_min); 03152 times2pown(rez,-2); 03153 resyu = Sup(rez); 03154 } 03155 03156 if( -Sup(y_extr) < simz && -Inf(y_extr) > iimz ) 03157 // extremal curve intersects left boundary of |x|+iy in quadrant IV 03158 // update minimum 03159 // resyl = -Sup( ln( 1 + 4 * y_extr / ( sqr( x_extr ) + sqr( 1 - y_extr ) ) ) ) / 4.0; 03160 // resyl = -Sup( T_atan(abs_x_min)/4.0 ); // Blomquist 03161 { 03162 rez = T_atan(abs_x_min); 03163 times2pown(rez,-2); 03164 resyl = -Sup(rez); 03165 } 03166 03167 res = lx_cinterval( lx_interval(resxl,resxu),lx_interval(-resyu,-resyl) ); 03168 stagprec = stagsave; 03169 res = adjust(res); 03170 return res; 03171 } 03172 // 03173 //-- end acot ----------------------------------------------------------------- 03174 03175 //-- atanh -------------------------------------------------------------------- 03176 // 03177 lx_cinterval atanh( const lx_cinterval& z ) throw() 03178 // 03179 // atanh( z ) = - i * atan( i * z ) 03180 // 03181 { 03182 lx_cinterval res = atan( lx_cinterval( -Im(z), Re(z) ) ); 03183 return lx_cinterval( Im(res), -Re(res) ); 03184 } 03185 // 03186 //-- end atanh ---------------------------------------------------------------- 03187 03188 //-- acoth -------------------------------------------------------------------- 03189 // 03190 lx_cinterval acoth( const lx_cinterval& z ) throw() 03191 // 03192 // acoth( z ) = i * acot( i * z ) 03193 // 03194 { 03195 lx_cinterval res = acot( lx_cinterval( -Im(z), Re(z) ) ); 03196 return lx_cinterval( -Im(res), Re(res) ); 03197 } 03198 // 03199 //-- end acoth ---------------------------------------------------------------- 03200 03201 // ---- sqrt1px2 -------------------------------------------------------------- 03202 // 03203 lx_cinterval sqrt1px2(const lx_cinterval& z) throw() 03204 // sqrt(1+z^2); 03205 // Blomquist, 03.07.2008; 03206 { 03207 const lx_real c = lx_real(1600,l_real(1)); 03208 int stagsave = stagprec, 03209 stagmax = 30; // l_imath.cpp: sqrt() uses stagmax = 30; 03210 if (stagprec>stagmax) stagprec = stagmax; 03211 03212 lx_cinterval res; 03213 lx_interval absz(abs(z)); 03214 lx_real Inf_absz(Inf(absz)); 03215 03216 if (Inf_absz > c) 03217 { 03218 absz = 1 / lx_interval(Inf_absz); 03219 Inf_absz = Sup(absz); 03220 res = lx_cinterval(lx_interval(-Inf_absz,Inf_absz), 03221 lx_interval(-Inf_absz,Inf_absz)); 03222 // res is the correcture interval, i.e. 03223 // z + res or -z + res is the inclusionof sqrt{1+z^2} 03224 res = (Inf(Re(z))>=0)? z + res : -z + res; 03225 } 03226 else 03227 { 03228 res = lx_cinterval(lx_interval(0,l_interval(0)), 03229 lx_interval(0,l_interval(1))); // res = i 03230 if ( Sup(abs(z-res))<0.5 || Sup(abs(z+res))<0.5 ) 03231 { 03232 res = lx_cinterval(-Im(z),Re(z)); // Res = i*z; 03233 // (1 - i*z)*(1 + i*z) = 1+z^2; 03234 res = sqrt( (1-res)*(1+res) ); 03235 } 03236 else 03237 res = sqrt(1+sqr(z)); 03238 } 03239 if (sm_zero(Inf(Re(res)))) 03240 res = lx_cinterval(lx_interval(lx_real(0),Sup(Re(res))),Im(res)); 03241 stagprec = stagsave; 03242 res = adjust(res); 03243 return res; 03244 } 03245 // -- end sqrt1px2 ------------------------------------------------------------ 03246 03247 lx_cinterval sqrt1mx2(const lx_cinterval& z) throw() 03248 // sqrt(1-z^2); 03249 // Blomquist, 03.07.2008; 03250 { 03251 const lx_real c = lx_real(1600,l_real(1)); 03252 int stagsave = stagprec, 03253 stagmax = 30; // l_imath.cpp: sqrt() uses stagmax = 30; 03254 if (stagprec>stagmax) stagprec = stagmax; 03255 03256 lx_cinterval res,u; 03257 lx_interval absz(abs(z)); 03258 lx_real Inf_absz(Inf(absz)); 03259 03260 if (Inf_absz > c) 03261 { 03262 absz = 1 / lx_interval(Inf_absz); 03263 Inf_absz = Sup(absz); 03264 res = lx_cinterval(lx_interval(-Inf_absz,Inf_absz), 03265 lx_interval(-Inf_absz,Inf_absz)); // res = Delta 03266 u = lx_cinterval(-Im(z),Re(z)); // u = i*z; 03267 // res is the correcture interval, i.e. 03268 // i*z + res or -i*z + res is the inclusion of sqrt{1-z^2} 03269 res = (Inf(Im(z))>=0)? -u + res : u + res; 03270 } 03271 else 03272 { 03273 res = 1-z; u = 1+z; 03274 res = (Sup(abs(res))<0.5 || Sup(abs(u))<0.5)? sqrt(res*u) : sqrt(1-sqr(z)); 03275 } 03276 if (sm_zero(Inf(Re(res)))) 03277 res = lx_cinterval(lx_interval(lx_real(0),Sup(Re(res))),Im(res)); 03278 stagprec = stagsave; 03279 res = adjust(res); 03280 return res; 03281 } 03282 03283 lx_cinterval sqrtx2m1(const lx_cinterval& z) throw() 03284 // sqrt(z^2-1); 03285 // Blomquist, 03.07.2008; 03286 { 03287 const lx_real c = lx_real(1600,l_real(1)); 03288 int stagsave = stagprec, 03289 stagmax = 30; // l_imath.cpp: sqrt() uses stagmax = 30; 03290 if (stagprec>stagmax) stagprec = stagmax; 03291 03292 lx_cinterval res,u; 03293 lx_interval absz(abs(z)); 03294 lx_real Inf_absz(Inf(absz)); 03295 03296 if (Inf_absz > c) 03297 { 03298 absz = 1 / lx_interval(Inf_absz); 03299 Inf_absz = Sup(absz); 03300 res = lx_cinterval(lx_interval(-Inf_absz,Inf_absz), 03301 lx_interval(-Inf_absz,Inf_absz)); // res = Delta 03302 // res is the correcture interval, i.e. 03303 res = (Inf(Re(z))>=0)? z + res : -z + res; 03304 } 03305 else 03306 { 03307 res = z-1; u = z+1; 03308 res = (Sup(abs(res))<0.5 || Sup(abs(u))<0.5)? sqrt(res*u) : sqrt(sqr(z)-1); 03309 } 03310 03311 if (sm_zero(Inf(Re(res)))) 03312 res = lx_cinterval(lx_interval(lx_real(0),Sup(Re(res))),Im(res)); 03313 03314 stagprec = stagsave; 03315 res = adjust(res); 03316 return res; 03317 } 03318 03319 lx_cinterval sqrtp1m1(const lx_cinterval& z) throw() 03320 // sqrt(1+z)-1; 03321 // Blomquist, 08.07.2008; 03322 { 03323 const real c = 0.125; 03324 int stagsave = stagprec, 03325 stagmax = 30; // l_imath.cpp: sqrt() uses stagmax = 30; 03326 if (stagprec>stagmax) stagprec = stagmax; 03327 03328 lx_cinterval res; 03329 lx_interval absz(abs(z)); 03330 lx_real Sup_absz(Sup(absz)); 03331 03332 if (Sup_absz < c) 03333 res = z / (sqrt(z+1) + 1); 03334 else 03335 res = sqrt(z+1) - 1; 03336 03337 stagprec = stagsave; 03338 res = adjust(res); 03339 return res; 03340 } 03341 03342 lx_cinterval expm1(const lx_cinterval& z) throw() 03343 // exp(z) - 1; 03344 // Blomquist, 09.08.2008; 03345 { 03346 int stagsave = stagprec, 03347 stagmax = 30; // l_imath.cpp: sqrt() uses stagmax = 30; 03348 if (stagprec>stagmax) stagprec = stagmax; 03349 03350 const lx_interval cancl_test = lx_interval(0,l_interval(0.995,1.005)); 03351 lx_interval rez(Re(z)), imz(Im(z)); 03352 lx_interval exp_x, sin_y, cos_y, h, Xt; 03353 lx_cinterval res; 03354 03355 exp_x = exp(rez); 03356 sin_y = sin(imz); 03357 cos_y = cos(imz); 03358 03359 h = exp_x*cos_y; 03360 if (h < cancl_test && cos_y < cancl_test) 03361 { 03362 h = lnp1(-sqr(sin_y)); 03363 times2pown(h,-1); 03364 // h = 0.5 * ln(1-sqr( sin(y) )); 03365 h = expm1(rez+h); // Cancellation also possible here! 03366 } 03367 else h = h - 1; // Cancellation possible here (real part)! 03368 // h: Real part; 03369 imz = exp_x * sin_y; 03370 res = lx_cinterval(h,imz); 03371 03372 stagprec = stagsave; 03373 res = adjust(res); 03374 03375 return res; 03376 } 03377 03378 lx_cinterval lnp1(const lx_cinterval& z) throw() 03379 {// ln(1+z); 03380 // Blomquist, 11.08.2008; 03381 int stagsave = stagprec, 03382 stagmax = 30; 03383 if (stagprec > stagmax) stagprec = stagmax; 03384 const real c1 = 1e-128; 03385 lx_cinterval y; 03386 lx_interval abs_z(abs(z)); 03387 lx_real 03388 srez = Sup( Re(z) ), 03389 simz = Sup( Im(z) ), 03390 iimz = Inf( Im(z) ); 03391 03392 if (-1 <= z) // z contains -1 03393 cxscthrow(STD_FKT_OUT_OF_DEF( 03394 "lx_cinterval lnp1(const lx_cinterval& z); z contains -1")); 03395 if ( srez<-1 && iimz<0 && simz>=0 ) 03396 cxscthrow(STD_FKT_OUT_OF_DEF( 03397 "lx_cinterval lnp1(const lx_cinterval& z); z not allowed")); 03398 if (Sup(abs_z) < c1) 03399 { 03400 abs_z = Re(z); 03401 abs_z = lnp1( abs_z*(2+abs_z) + sqr(Im(z)) ); 03402 times2pown(abs_z,-1); 03403 y = lx_cinterval( abs_z, arg(1+z) ); 03404 } 03405 else 03406 y = Ln(1+z); 03407 stagprec = stagsave; 03408 y = adjust(y); 03409 03410 return y; 03411 } 03412 03413 //-- pow_all ------ 03414 // 03415 // Non-analytic power function for real l_interval exponent. 03416 // 03417 // If 0 \not\in z, then compute four rectangular intervals that comprehend 03418 // an annulus which contains all values zeta^phi, zeta in z, phi in p. 03419 // 03420 // If 0 in z and negative reals in p, then abort execution 03421 // (potential modification: return the entire complex plane). 03422 // 03423 std::list<lx_cinterval> pow_all( const lx_cinterval& z, 03424 const lx_interval& p ) throw() 03425 { 03426 lx_interval abs_z = abs(z); 03427 03428 if( 0.0 < Inf( abs_z ) ) 03429 { 03430 lx_interval abs_z_p = exp( p * ln( abs_z ) ); 03431 03432 // Inner and outer radii of the annulus are inf/sup( abs_z_n ) 03433 // Inscribed square has side length sqrt( 2 ) * rad_1 03434 lx_interval rad_1 = Sqrt2r_l_interval() * lx_interval(Inf( abs_z_p )); 03435 lx_interval rad_2 = lx_interval(Sup( abs_z_p )); 03436 03437 std::list<lx_cinterval> res; 03438 03439 // 4 intervals covering the annulus, counter-clockwise 03440 res.push_back(lx_cinterval(lx_interval( Inf( rad_1 ), Sup( rad_2 ) ), 03441 lx_interval( -Sup( rad_1 ), Sup( rad_2 ) ))); 03442 res.push_back(lx_cinterval(lx_interval( -Sup( rad_2 ), Sup( rad_1 ) ), 03443 lx_interval( Inf( rad_1 ), Sup( rad_2 ) ) )); 03444 res.push_back(lx_cinterval(lx_interval( -Sup( rad_2 ), -Inf( rad_1 ) ), 03445 lx_interval( -Sup( rad_2 ), Sup(rad_1) ) ) ); 03446 res.push_back(lx_cinterval(lx_interval( -Sup( rad_1 ), Sup( rad_2 ) ), 03447 lx_interval( -Sup( rad_2 ), -Inf(rad_1) ) ) ); 03448 03449 return res; 03450 } 03451 else 03452 // 0 in z 03453 { 03454 if( Inf( p ) > 0.0 ) 03455 // 03456 // All values zeta^phi, zeta in z, phi in p lie in a disk 03457 // 03458 { 03459 lx_interval abs_z_p = exp( p * ln( lx_interval( Sup( abs_z ) ) ) ); 03460 lx_real rad_p = Sup( abs_z_p ); 03461 03462 std::list<lx_cinterval> res; 03463 03464 res.push_back( lx_cinterval( lx_interval( -rad_p, rad_p ), 03465 lx_interval( -rad_p, rad_p ) ) ); 03466 return res; 03467 } 03468 else 03469 { 03470 // 03471 // The set zeta^phi, zeta in z, phi in p is unbounded 03472 // if inf( p ) < 0. 0^p is undefined for p <= 0. 03473 // 03474 cxscthrow(STD_FKT_OUT_OF_DEF( 03475 "pow_all(lx_cinterval, lx_interval); 0^p is undefined for p <= 0.")); 03476 std::list<lx_cinterval> res; 03477 return res; 03478 } 03479 } 03480 } 03481 // 03482 //-- end pow_all -------------------------------------------------------------- 03483 03484 } // namespace cxsc
1.7.5.1