C-XSC - A C++ Class Library for Extended Scientific Computing
2.5.4
|
00001 /* 00002 ** CXSC is a C++ library for eXtended Scientific Computing (V 2.5.4) 00003 ** 00004 ** Copyright (C) 1990-2000 Institut fuer Angewandte Mathematik, 00005 ** Universitaet Karlsruhe, Germany 00006 ** (C) 2000-2014 Wiss. Rechnen/Softwaretechnologie 00007 ** Universitaet Wuppertal, Germany 00008 ** 00009 ** This library is free software; you can redistribute it and/or 00010 ** modify it under the terms of the GNU Library General Public 00011 ** License as published by the Free Software Foundation; either 00012 ** version 2 of the License, or (at your option) any later version. 00013 ** 00014 ** This library is distributed in the hope that it will be useful, 00015 ** but WITHOUT ANY WARRANTY; without even the implied warranty of 00016 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00017 ** Library General Public License for more details. 00018 ** 00019 ** You should have received a copy of the GNU Library General Public 00020 ** License along with this library; if not, write to the Free 00021 ** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00022 */ 00023 00024 /* CVS $Id: l_cimath.cpp,v 1.21 2014/01/30 17:23:46 cxsc Exp $ */ 00025 00026 /* 00027 ** 00028 ** File: l_cimath.cpp, 2007/03/05. 00029 ** 00030 ** Copyright (C) Markus Neher, markus.neher@math.uni-karlsruhe.de 00031 ** Ingo Eble, ingoeble@web.de 00032 ** Frithjof Blomquist, Blomquist@math.uni-wuppertal.de 00033 */ 00034 00035 //Include header files 00036 #include <l_real.hpp> 00037 #include <l_cimath.hpp> // Complex functions in staggered format 00038 #include <l_imath.hpp> // "l_interval" standard functions 00039 00040 namespace cxsc{ 00041 00042 inline l_real min(const l_real& a, const l_real& b) 00043 { 00044 return (a<b)? a : b; 00045 } 00046 00047 inline l_real max(const l_real& a, const l_real& b) 00048 { 00049 return (a>b)? a : b; 00050 } 00051 00052 inline const l_interval& ONE_INTERVAL() 00053 { 00054 static const l_interval one = l_interval(1.0); 00055 return one; 00056 } 00057 00058 bool disjoint(const l_interval& x, const l_interval& y) 00059 { 00060 l_real ix( Inf(x) ),iy( Inf(y) ),sx( Sup(x) ), sy( Sup(y) ); 00061 l_real inf( ( ix > iy )? ix : iy ); 00062 l_real sup( ( sx < sy )? sx : sy ); 00063 00064 return ( inf > sup ); 00065 } 00066 00067 00068 /* ***************************************************************************/ 00069 /* ***************************************************************************/ 00070 /* *** Single-valued functions *** */ 00071 /* ***************************************************************************/ 00072 /* ***************************************************************************/ 00073 00074 00075 /* ***************************************************************************/ 00076 /* *** Power operator pow is not listed here, since it relies on the ****/ 00077 /* *** (multi-valued) logarithm ****/ 00078 /* ***************************************************************************/ 00079 00080 00081 /* ***************************************************************************/ 00082 /* *** The hyperbolic functions exp, sin, cos, sinh, cosh are separable: ****/ 00083 /* *** Their real and imaginary parts are products of real functions ****/ 00084 /* ***************************************************************************/ 00085 /* *** With Re(z)=x, Im(z)=y : ****/ 00086 /* *** ****/ 00087 /* *** exp : Re(exp(z)) = exp(x) * cos(y) ****/ 00088 /* *** Im(exp(z)) = exp(x) * sin(y) ****/ 00089 /* *** ****/ 00090 /* *** sin : Re(sin(z)) = sin(x) * cosh(y) ****/ 00091 /* *** Im(sin(x)) = cos(x) * sinh(y) ****/ 00092 /* *** ****/ 00093 /* *** cos : Re(cos(z)) = cos(x) * cosh(y) ****/ 00094 /* *** Im(sin(x)) = -sin(x) * sinh(y) ****/ 00095 /* *** ****/ 00096 /* *** sinh : Re(sinh(z)) = sinh(x) * cos(y) ****/ 00097 /* *** Im(sinh(z)) = cosh(x) * sin(y) ****/ 00098 /* *** ****/ 00099 /* *** cosh : Re(cosh(z)) = cosh(x) * cos(y) ****/ 00100 /* *** Im(cosh(z)) = sinh(x) * sin(y) ****/ 00101 /* *** ****/ 00102 /* ***************************************************************************/ 00103 00104 00105 l_cinterval exp(const l_cinterval& z) throw() 00106 { 00107 int stagsave = stagprec, 00108 stagmax = 19; 00109 l_interval lreal(Re(z)), limg(Im(z)); 00110 cinterval dz(z); 00111 l_cinterval y; 00112 bool grob = ( (Sup(Re(dz)) > succ(succ(Inf(Re(dz))))) || 00113 (Sup(Im(dz)) > succ(succ(Inf(Im(dz))))) ); 00114 00115 if (stagprec==1 || grob) y = exp(dz); // schnelle Auswertung! 00116 else 00117 { 00118 if (stagprec < stagmax) stagprec++; 00119 else stagprec = stagmax; 00120 l_interval 00121 A( exp(lreal) ), 00122 B( limg ); 00123 y = l_cinterval( A*cos( B ) , A*sin( B ) ); 00124 stagprec = stagsave; 00125 y = adjust(y); 00126 } 00127 return y; 00128 } 00129 00130 l_cinterval exp2(const l_cinterval& z) throw() 00131 { 00132 return exp(z*Ln2_l_interval()); 00133 } 00134 00135 l_cinterval exp10(const l_cinterval& z) throw() 00136 { 00137 return exp(z*Ln10_l_interval()); 00138 } 00139 00140 l_cinterval expm1(const l_cinterval& z) throw() 00141 // exp(z) - 1; 00142 // Calculates nearly optimal inclusions for not too wide intervals z. 00143 // Blomquist, 03.12.2008; 00144 { 00145 int stagsave = stagprec, 00146 stagmax = 30; // l_imath.cpp: sqrt() uses stagmax = 30; 00147 if (stagprec>stagmax) stagprec = stagmax; 00148 00149 const l_interval cancl_test = l_interval(0.995,1.005); 00150 l_interval rez(Re(z)), imz(Im(z)); 00151 l_interval exp_x, sin_y, cos_y, h, Xt; 00152 l_cinterval res; 00153 00154 exp_x = exp(rez); 00155 sin_y = sin(imz); 00156 cos_y = cos(imz); 00157 00158 h = exp_x*cos_y; 00159 if (h < cancl_test && cos_y < cancl_test) 00160 { 00161 h = lnp1(-sqr(sin_y)); 00162 times2pown(h,-1); 00163 // h = 0.5 * ln(1-sqr( sin(y) )); 00164 h = expm1(rez+h); // Cancellation also possible here! 00165 } 00166 else h = h - 1; // Cancellation possible here (real part)! 00167 // h: Real part; 00168 imz = exp_x * sin_y; 00169 res = l_cinterval(h,imz); 00170 00171 stagprec = stagsave; 00172 res = adjust(res); 00173 00174 return res; 00175 } 00176 00177 l_cinterval cos(const l_cinterval& z) throw() 00178 { 00179 int stagsave = stagprec, 00180 stagmax = 19; 00181 l_interval lreal(Re(z)),limg(Im(z)); 00182 cinterval dz(z); 00183 l_cinterval y; 00184 bool grob = ( (Sup(Re(dz)) > succ(succ(Inf(Re(dz))))) || 00185 (Sup(Im(dz)) > succ(succ(Inf(Im(dz))))) ); 00186 00187 if (stagprec==1 || grob) y = cos(dz); // schnelle Auswertung! 00188 else 00189 { 00190 if (stagprec < stagmax) stagprec++; 00191 else stagprec = stagmax; 00192 l_interval 00193 A( lreal ), 00194 B( limg ); 00195 y = l_cinterval( cos( A )*cosh( B ) , -sin( A )*sinh( B ) ); 00196 stagprec = stagsave; 00197 y = adjust(y); 00198 } 00199 return y; 00200 } 00201 00202 l_cinterval sin(const l_cinterval& z) throw() 00203 { 00204 int stagsave = stagprec, 00205 stagmax = 19; 00206 l_interval lreal(Re(z)),limg(Im(z)); 00207 cinterval dz(z); 00208 l_cinterval y; 00209 bool grob = ( (Sup(Re(dz)) > succ(succ(Inf(Re(dz))))) || 00210 (Sup(Im(dz)) > succ(succ(Inf(Im(dz))))) ); 00211 00212 if (stagprec==1 || grob) y = sin(dz); // schnelle Auswertung! 00213 else 00214 { 00215 if (stagprec < stagmax) stagprec++; 00216 else stagprec = stagmax; 00217 l_interval 00218 A( lreal ), 00219 B( limg ); 00220 y = l_cinterval( sin( A )*cosh( B ) , cos( A )*sinh( B ) ); 00221 stagprec = stagsave; 00222 y = adjust(y); 00223 } 00224 return y; 00225 } 00226 00227 l_cinterval cosh(const l_cinterval& z) throw() 00228 { 00229 int stagsave = stagprec, 00230 stagmax = 19; 00231 l_interval lreal(Re(z)),limg(Im(z)); 00232 cinterval dz(z); 00233 l_cinterval y; 00234 bool grob = ( (Sup(Re(dz)) > succ(succ(Inf(Re(dz))))) || 00235 (Sup(Im(dz)) > succ(succ(Inf(Im(dz))))) ); 00236 00237 if (stagprec==1 || grob) y = cosh(dz); // schnelle Auswertung! 00238 else 00239 { 00240 if (stagprec < stagmax) stagprec++; 00241 else stagprec = stagmax; 00242 l_interval 00243 A( lreal ), 00244 B( limg ); 00245 y = l_cinterval( cos( B )*cosh( A ) , sin( B )*sinh( A ) ); 00246 stagprec = stagsave; 00247 y = adjust(y); 00248 } 00249 return y; 00250 } 00251 00252 l_cinterval sinh(const l_cinterval& z) throw() 00253 { 00254 int stagsave = stagprec, 00255 stagmax = 19; 00256 l_interval lreal(Re(z)),limg(Im(z)); 00257 cinterval dz(z); 00258 l_cinterval y; 00259 bool grob = ( (Sup(Re(dz)) > succ(succ(Inf(Re(dz))))) || 00260 (Sup(Im(dz)) > succ(succ(Inf(Im(dz))))) ); 00261 00262 if (stagprec==1 || grob) y = sinh(dz); // schnelle Auswertung! 00263 else 00264 { 00265 if (stagprec < stagmax) stagprec++; 00266 else stagprec = stagmax; 00267 l_interval 00268 A( lreal ), 00269 B( limg ); 00270 y = l_cinterval( cos( B )*sinh( A ) , sin( B )*cosh( A ) ); 00271 stagprec = stagsave; 00272 y = adjust(y); 00273 } 00274 return y; 00275 } 00276 00277 l_cinterval sqr(const l_cinterval& z) throw() 00278 // Blomquist, 21.11.2006; 00279 { 00280 dotprecision akku; 00281 l_interval rez(Re(z)), reza(abs(rez)), 00282 imz(Im(z)), imza(abs(imz)); 00283 l_real 00284 irez = Inf(reza), 00285 srez = Sup(reza), 00286 iimz = Inf(imza), 00287 simz = Sup(imza); 00288 00289 akku = 0; 00290 accumulate(akku,irez,irez); 00291 accumulate(akku,-simz,simz); 00292 irez = cxsc::rnd_down(akku); 00293 00294 akku = 0; 00295 accumulate(akku,srez,srez); 00296 accumulate(akku,-iimz,iimz); 00297 srez = cxsc::rnd_up(akku); 00298 00299 rez = rez * imz; 00300 times2pown(rez,1); // fast multiplikation with 2; 00301 00302 return l_cinterval( l_interval(irez,srez), rez ); 00303 } 00304 00305 //----------------------------------------------------------------------------- 00306 // 00307 // Section 2: tan, cot, tanh, coth 00308 // 00309 // The implementation of cot, tanh, and coth is based on tan: 00310 // 00311 // cot( z ) = tan( pi/2 - z ) 00312 // tanh( z ) = transp( i * tan( transp( i * z ) ) 00313 // coth( z ) = i * cot( i * z ) = i * tan( pi/2 - i * z ) 00314 // 00315 //----------------------------------------------------------------------------- 00316 00317 //-- tan ------------------------------------------------------------ 040827 -- 00318 // 00319 // Complex tangent function 00320 // 00321 00322 void horizontal_check( //-------------------------------------------- 040726 -- 00323 const l_interval& hy, const l_interval& cosh_2y, l_real irez, l_real srez, 00324 const l_interval& hxl, 00325 const l_interval& hxu, l_real& resxl, l_real& resxu ) 00326 // 00327 // Subroutine of tangent function. 00328 // Check intersections with extremal curves of tan on a horizontal boundary. 00329 // This subroutine is only called if an intersection occurs. 00330 // In this case, sinh( 2 * hy ) <> 0.0 (poles are handled before). 00331 // 00332 // There may be 1 or 2 intersections. 00333 // If intersections lie next to a boundary of rez, then it is impossible to 00334 // decide if there are 1 or 2 intersections. 00335 // In this case, 2 intersections are assumed 00336 // (valid enclosure, at the expense of a potential slight overestimation). 00337 // 00338 { 00339 bool both = false, left = false, right = false; 00340 00341 if (srez - irez > Inf( Pid2_l_interval() )) 00342 // 2 intersections 00343 both = true; 00344 else 00345 { 00346 l_interval 00347 res_l = cos( 2 * hxl ) - cosh_2y, 00348 res_u = cos( 2 * hxu ) - cosh_2y; 00349 00350 if( Inf( res_l * res_u ) > 0.0 ) 00351 // 2 intersections 00352 both = true; 00353 else 00354 { 00355 if( Sup( res_l * res_u ) < 0.0 ) 00356 { 00357 // 1 intersection (3 intersections are PI() apart) 00358 // neither of the intervals res_l and res_u contains zero 00359 if( Inf( res_l ) > 0.0 ) 00360 // "left" intersection 00361 left = true; 00362 else 00363 // "right" intersection 00364 right = true; 00365 } 00366 else 00367 // 00368 // 1 (or both) intersections lie next to a boundary point 00369 // here, the task is to decide if 2 intersections occurs 00370 // if only 1 intersection takes place, then which one? 00371 // 00372 { 00373 l_interval 00374 sin_2xl = sin( 2 * hxl ), 00375 sin_2xu = sin( 2 * hxu ); 00376 00377 if( !disjoint( l_interval(0), res_l ) ) 00378 // intersection on the left boundary 00379 { 00380 if( Inf( sin_2xl ) >= 0.0 ) 00381 // "left" intersection 00382 { 00383 left = true; 00384 // remove the intersection by changing res_l! 00385 res_l = -l_interval(1); 00386 } 00387 else 00388 { 00389 if( Sup( sin_2xl ) <= 0.0 ) 00390 // "right" intersection 00391 { 00392 right = true; 00393 // remove the intersection by changing res_l! 00394 res_l = l_interval(1); 00395 } 00396 else 00397 // zero is interior point of sin_2xl 00398 // if the real sine function has optimal precision, 00399 // this case should never happen 00400 both = true; 00401 } 00402 } 00403 00404 if( !disjoint( l_interval(0), res_u ) ) 00405 // intersection on the right boundary 00406 { 00407 if( Inf( sin_2xu ) >= 0.0 ) 00408 // "left" intersection 00409 { 00410 left = true; 00411 // remove the intersection by changing res_u! 00412 res_u = l_interval(1); 00413 } 00414 else 00415 { 00416 if( Sup( sin_2xu ) <= 0.0 ) 00417 // "right" intersection 00418 { 00419 right = true; 00420 // remove the intersection by changing res_u! 00421 res_u = -l_interval(1); 00422 } 00423 else 00424 // zero is interior point of sin_2xu 00425 // if the real sine function has optimal precision, 00426 // this case should never happen 00427 both = true; 00428 } 00429 } 00430 // 00431 // finally, check if there is a second intersection 00432 // 00433 if( Inf( res_l * res_u ) < 0.0 ) 00434 both = true; 00435 } 00436 } 00437 } 00438 // 00439 // Calculate extremal values 00440 // 00441 l_interval re_tan = 1 / sinh( 2 * abs( hy ) ); 00442 00443 // "left" intersection, sin( 2x ) > 0 00444 if( left || both ) 00445 { 00446 resxl = min( resxl, Inf( re_tan ) ); 00447 resxu = max( resxu, Sup( re_tan ) ); 00448 } 00449 00450 // "right" intersection, sin( 2x ) < 0 00451 if( right || both ) 00452 { 00453 resxl = min( resxl, -Sup( re_tan ) ); 00454 resxu = max( resxu, -Inf( re_tan ) ); 00455 } 00456 } // end horizontal_check 00457 00458 00459 l_cinterval Tan( const l_cinterval& z ) throw() // -------------- 040827 -- 00460 { 00461 l_cinterval y; 00462 00463 l_interval 00464 rez = Re(z), // rez = z.re(), 00465 imz = Im(z); // imz = z.im(); 00466 00467 l_real 00468 irez = Inf(rez), 00469 srez = Sup(rez), 00470 iimz = Inf(imz), 00471 simz = Sup(imz); 00472 00473 l_interval 00474 hxl(irez), hxu(srez), hyl(iimz), hyu(simz); 00475 00476 l_real 00477 resxl, resxu, resyl, resyu; 00478 // 00479 // 1st: check for poles 00480 // 00481 if( ( !disjoint( l_interval(0), imz ) ) && 00482 ( !disjoint( l_interval(0), cos( rez ) ) ) ) 00483 cxscthrow (STD_FKT_OUT_OF_DEF("l_cinterval tan( const l_cinterval& z); Pole(s) in z")); 00484 // 00485 // 2nd: real part 00486 // 00487 // evaluate tan on vertical boundaries 00488 // 00489 l_interval 00490 cos_2rez = cos( 2 * rez ), 00491 sinh_imz_2 = sqr( sinh( imz ) ); 00492 00493 l_interval 00494 re_tan_l=sin( 2*hxl ) / ( 2*( sqr( cos( hxl ) ) + sinh_imz_2 ) ), 00495 re_tan_u=sin( 2*hxu ) / ( 2*( sqr( cos( hxu ) ) + sinh_imz_2 ) ); 00496 00497 resxl = min( Inf( re_tan_l ), Inf( re_tan_u ) ); 00498 resxu = max( Sup( re_tan_l ), Sup( re_tan_u ) ); 00499 00500 // 00501 // look for extremal values on horizontal boundaries 00502 // if one of the horizontal boundaries is the x-axes, 00503 // then the complex tangent is the real tangent on this 00504 // boundary, and due to monotonicity, its range 00505 // is already included in the ranges of the vertical boundaries 00506 // 00507 if( irez < srez ) 00508 { 00509 l_interval 00510 cosh_2yl = - 1 / cosh( 2 * hyl ), 00511 cosh_2yu = - 1 / cosh( 2 * hyu ); 00512 00513 if( !disjoint( cos_2rez, cosh_2yl ) && iimz != 0.0 ) 00514 //extremal curve intersects lower boundary 00515 horizontal_check(hyl,cosh_2yl,irez,srez,hxl,hxu,resxl,resxu); 00516 00517 if( !disjoint( cos_2rez, cosh_2yu ) && simz != 0.0 ) 00518 //extremal curve intersects upper boundary 00519 horizontal_check(hyu,cosh_2yu,irez,srez,hxl,hxu,resxl,resxu); 00520 } 00521 // 00522 // 3rd: imaginary part 00523 // 00524 // evaluate tan on horizontal boundaries 00525 // 00526 l_interval 00527 cos_rez_2 = sqr( cos( rez ) ); 00528 00529 l_interval 00530 im_tan_l = sinh( 2*hyl ) / (2*(cos_rez_2 + sqr( sinh( hyl ) ))), 00531 im_tan_u = sinh( 2*hyu ) / (2*(cos_rez_2 + sqr( sinh( hyu ) ))); 00532 00533 resyl = min( Inf( im_tan_l ), Inf( im_tan_u ) ); 00534 resyu = max( Sup( im_tan_l ), Sup( im_tan_u ) ); 00535 00536 // 00537 // look for extremal values on vertical boundaries 00538 // here, the situation is simpler than for the real part 00539 // if 2 intersections with extremal curves appear , 00540 // one lies in the lower half plane, the other in the upper half plane 00541 // 00542 l_interval 00543 cos_2xl = cos( 2 * hxl ), 00544 cos_2xu = cos( 2 * hxu ); 00545 l_interval im_tan; 00546 00547 if( iimz < 0.0 ) 00548 // intersection in lower half plane? 00549 { 00550 l_interval 00551 imz_h = l_interval( iimz, min( simz, l_real(0.0) ) ), 00552 cosh_2imz = - 1 / cosh( 2 * imz_h ); 00553 00554 if( ( !disjoint( cosh_2imz, cos_2xl ) ) ) 00555 //extremal curve intersects left boundary 00556 //in this case, sin( 2 * xl ) <> 0.0 (no poles here!) 00557 { 00558 im_tan = - 1 / abs( sin( 2 * hxl ) ); 00559 resyl = min( resyl, Inf( im_tan ) ); 00560 resyu = max( resyu, Sup( im_tan ) ); 00561 } 00562 if( ( !disjoint( cosh_2imz, cos_2xu ) ) ) 00563 //extremal curve intersects right boundary 00564 //in this case, sin( 2 * xu ) <> 0.0 (no poles here!) 00565 { 00566 im_tan = - 1 / abs( sin( 2 * hxu ) ); 00567 resyl = min( resyl, Inf( im_tan ) ); 00568 resyu = max( resyu, Sup( im_tan ) ); 00569 } 00570 } 00571 if( simz > 0.0 ) 00572 // intersection in upper half plane? 00573 { 00574 l_interval 00575 imz_h = l_interval( max( iimz, l_real(0.0) ), simz ), 00576 cosh_2imz = - 1 / cosh( 2 * imz_h ); 00577 00578 if( ( !disjoint( cosh_2imz, cos_2xl ) ) ) 00579 //extremal curve intersects left boundary 00580 //in this case, sin( 2 * xl ) <> 0.0 (no poles here!) 00581 { 00582 im_tan = + 1 / abs( sin( 2 * hxl ) ); 00583 resyl = min( resyl, Inf( im_tan ) ); 00584 resyu = max( resyu, Sup( im_tan ) ); 00585 } 00586 if( ( !disjoint( cosh_2imz, cos_2xu ) ) ) 00587 //extremal curve intersects right boundary 00588 //in this case, sin( 2 * xu ) <> 0.0 (no poles here!) 00589 { 00590 im_tan = + 1 / abs( sin( 2 * hxu ) ); 00591 resyl = min( resyl, Inf( im_tan ) ); 00592 resyu = max( resyu, Sup( im_tan ) ); 00593 } 00594 } 00595 00596 y = l_cinterval( l_interval(resxl,resxu ),l_interval(resyl,resyu ) ); 00597 00598 return y; 00599 00600 } // Tan 00601 00602 l_cinterval tan( const l_cinterval& z ) throw() 00603 { 00604 // tan(z) has the poles z_k = pi*(1+2k)/2; |k| in {0,1,2,3,...}. 00605 // z = z_k + eps = pi*(1+2k)/2 + eps; With |eps|<<1 k can be calculated 00606 // by: k = Re(z)/pi - 0.5; With this k the komplex value eps is given 00607 // by: eps = z - pi*(1+2k)/2; pi = 3.1415926... ; 00608 // It holds: 00609 // tan(z) = tan(z_k+eps) = tan[pi*(1+2k)/2 + eps] 00610 // = tan[pi/2 + pi*k + eps] = tan[pi/2 + eps] = -1 / tan(eps); 00611 // Definitions: u = Re(eps); u = abs(u); v = Im(eps); v = abs(v); 00612 // if (Sup(u)<S && Sup(v)<S) tan(z) = -1 / Tan(eps); 00613 // else tan(z) = Tan(z); S = 1e-15; 00614 // Thus, near the poles tan(z) is calculated in higher accuracy with 00615 // -1 / Tan(eps); 00616 // Blomquist, 28.09.2007; 00617 const real S = 1e-15; 00618 int stagsave = stagprec, 00619 stagmax = 19; 00620 if (stagprec < stagmax) stagprec++; 00621 else stagprec = stagmax; 00622 00623 l_cinterval y,eps; 00624 l_interval rez = Re(z); 00625 00626 interval re(rez),u,v; 00627 real x(mid(re)); 00628 double dbr = _double(x), pi(3.14159265358979323); 00629 int k,s; 00630 dbr = dbr/pi - 0.5; 00631 s = sign(dbr); 00632 k = (s>=0)? CUTINT(dbr+0.5) : CUTINT(dbr-0.5); 00633 if (k<-2147483647) 00634 cxscthrow (STD_FKT_OUT_OF_DEF( 00635 "l_cinterval tan(const l_cinterval& z); z out of range")); 00636 eps = z - Pid2_l_interval()*(1+2*k); 00637 00638 u = Re(eps); u = abs(u); 00639 v = Im(eps); v = abs(v); 00640 00641 if (Sup(u)<S && Sup(v)<S) 00642 y = -l_cinterval(1) / Tan(eps); 00643 else y = Tan(z); 00644 00645 stagprec = stagsave; 00646 y = adjust(y); 00647 00648 return y; 00649 } // tan() 00650 00651 l_cinterval cot( const l_cinterval& z ) throw() 00652 { 00653 // cot(z) has the poles z_k = k*pi; |k| in {0,1,2,3,...}. 00654 // z = z_k + eps = k*pi + eps; With |eps|<<1 k can be calculated 00655 // by: k = Re(z)/pi; With this k the komplex value eps is given 00656 // by: eps = z - k*pi = z - (pi/2)*2k; pi = 3.1415926... ; 00657 // It holds: 00658 // cot(z) = cot(z_k+eps) = cot(k*pi + eps) 00659 // = cot(eps) = 1 / tan(eps); 00660 // Definitions: u = Re(eps); u = abs(u); v = Im(eps); v = abs(v); 00661 // if (Sup(u)<S && Sup(v)<S) cot(z) = 1 / tan(eps); 00662 // else cot(z) = tan(pi/2 - z); S = 1e-15; 00663 // Thus, near the poles cot(z) is calculated in higher accuracy with 00664 // 1 / tan(eps); 00665 // Blomquist, 29.09.2007; 00666 const real S = 1e-15; 00667 int stagsave = stagprec, 00668 stagmax = 19; 00669 if (stagprec < stagmax) stagprec++; 00670 else stagprec = stagmax; 00671 00672 l_cinterval y,eps; 00673 00674 l_interval rez = Re(z); 00675 interval re(rez),u,v; 00676 real x(mid(re)); 00677 00678 double dbr = _double(x), pi(3.14159265358979323); 00679 int k,s; 00680 dbr = dbr/pi; 00681 s = sign(dbr); 00682 k = (s>=0)? CUTINT(dbr+0.5) : CUTINT(dbr-0.5); 00683 if (k<-2147483647) 00684 cxscthrow (STD_FKT_OUT_OF_DEF( 00685 "l_cinterval cot(const l_cinterval& z); z out of range")); 00686 eps = z - Pid2_l_interval()*2*k; 00687 u = Re(eps); u = abs(u); 00688 v = Im(eps); v = abs(v); 00689 if (Sup(u)<S && Sup(v)<S) 00690 y = cinterval(1) / Tan(eps); 00691 else y = Tan(Pid2_l_interval() - z); 00692 00693 stagprec = stagsave; 00694 y = adjust(y); 00695 00696 return y; 00697 } // cot 00698 00699 00700 00701 //-- tanh --------------------------------------------------------- 00702 // 00703 // tanh( z ) = transp( i * tan( transp( i * z ) ) 00704 // 00705 l_cinterval tanh( const l_cinterval& z ) throw() 00706 { 00707 int stagsave = stagprec, 00708 stagmax = 19; 00709 if (stagprec < stagmax) stagprec++; 00710 else stagprec = stagmax; 00711 00712 l_cinterval res = tan( l_cinterval( Im(z), Re(z) ) ),y; 00713 y = l_cinterval( Im(res), Re(res) ); 00714 00715 stagprec = stagsave; 00716 y = adjust(y); 00717 00718 return y; 00719 } 00720 // 00721 //-- end tanh ------------------------------------------------------- 00722 00723 //-- coth ----------------------------------------------------------- 00724 // 00725 // coth( z ) = i * cot( i * z ); 00726 // 00727 00728 l_cinterval coth( const l_cinterval& z ) throw() 00729 { // coth( z ) = i * cot( i * z ); 00730 int stagsave = stagprec, 00731 stagmax = 19; 00732 if (stagprec < stagmax) stagprec++; 00733 else stagprec = stagmax; 00734 00735 l_cinterval zh = l_cinterval( -Im(z), Re(z) ); // zh = i*z; 00736 l_cinterval res = cot(zh); 00737 zh = l_cinterval( -Im(res), Re(res) ); 00738 00739 stagprec = stagsave; 00740 zh = adjust(zh); 00741 00742 return zh; 00743 } 00744 // 00745 //-- end coth ----------------------------------------------------------------- 00746 00747 00748 //----------------------------------------------------------------------------- 00749 // 00750 // Part II: Multi-valued functions 00751 // 00752 //----------------------------------------------------------------------------- 00753 // 00754 // 00755 00756 l_interval Atan(const l_interval& y, const l_interval& x) throw() 00757 // Calculating an inclusion of atan(y/x) with x<>[0,0]; 00758 // This help function must only be used for POINT intervals y,x !! 00759 // This function avoids internal overflow by calculating y/x. 00760 // atan(t), t in R, is a point symmetrical function, so we only need 00761 // to consider t >= 0; 00762 // For sufficiently great t the values atan(t) will be included by 00763 // Pid2_l_interval(), declared in l_imath.hpp. Thus, for t > t0, we must 00764 // fulfill atan(t) > INF := Inf( Pid2_l_interval() ). For stagprec = 19 00765 // it holds: INF = pi/2 - d, d := 2.6817655261... * 10^-308; 00766 // ---> atan(t) > INF = pi/2 - d must be fulfilled for t > t0 = ? 00767 // atan(t) is a monotonic function, so t0 is determined by 00768 // atan(t0) = INF = pi/2 - d ---> t0 = tan(pi/2 - d) = cot(d); 00769 // With cot(d) = 1/d - d/3 - (d^3)/45 - ... < 1/d, 0 < d < pi, we get: 00770 // 00771 // t > 1/d = 3.7288... * 10^+307 --> atan(t) > INF := Inf(Pid2_l_interval()); 00772 // 00773 // In the following algorithm with ex_x,ex_y it holds: 00774 // Infy[1] = my * 2^ex_y; Infx[1] = mx * 2^ex_x; 0.5 <= my,mx < 1; 00775 // For ex_y - ex_x >= N0:=1023 in this algorithm atan(y/x) is included by 00776 // Pid2_l_interval(); It is easy to show, that for ex_y - ex_x = N0:=1023 00777 // the minimal value of y/x is given by u := 2^1022 = 4.49423...*10^+307, 00778 // which is greater than 1/d = 3.7288... * 10^+307, so the inclusion of 00779 // atan(y/x) with Pid2_l_interval() is correct! 00780 // On the other side, for ex_y - ex_x <= 1022 the quotient y/x is smaller 00781 // than 2^1023 < MaxReal = 0.9999...*2^1024, so that an overflow by 00782 // calculating y/x is not possible! 00783 // Blomquist, 04.12.2006; 00784 { 00785 const int c = 1022; 00786 l_interval res(0); 00787 l_real Infx(Inf(x)), 00788 Infy(Inf(y)); 00789 int ex_x(expo_gr(Infx)), 00790 ex_y(expo_gr(Infy)), 00791 signx(sign(Infx)), 00792 signy(sign(Infy)), 00793 signq; 00794 if (signy!=0) { 00795 signq = signx * signy; 00796 if (ex_y-ex_x > c) 00797 res = signq>0 ? Pid2_l_interval() : -Pid2_l_interval(); 00798 else res = atan(y/x); 00799 } 00800 00801 return res; 00802 } 00803 00804 l_interval Atan(const l_interval& y, const l_real& x) throw() 00805 // Calculating an inclusion of atan(y/x) with x<>0.0; 00806 // This help function must only be used for POINT intervals y !! 00807 // This function avoids internal overflow by calculating y/x. 00808 { 00809 l_interval xi(x); 00810 return Atan(y,xi); 00811 } 00812 00813 // 00814 // For the Multi-valued functions, two different procedures are 00815 // implemented. 00816 // 00817 // First, there is a subroutine for computing the principal value 00818 // of the respective function. The principal value is usually 00819 // analytic in a subset of the complex plane and undefined otherwise. 00820 // 00821 // Second, there are procedures for computing interval supersets 00822 // of all function values of the respective function, where feasible. 00823 // 00824 //----------------------------------------------------------------------------- 00825 // 00826 // Section 1: Argument functions 00827 // 00828 //----------------------------------------------------------------------------- 00829 00830 //-- Arg: analytic argument function -------------------------------- 040916 -- 00831 // 00832 // (i) Arg(Z) subset (-pi,pi). 00833 // (ii) Arg([0,0]) = 0. 00834 // (iii) Arg(Z) is undefined if Z contains negative real numbers. 00835 // (iv) Otherwise, Arg(Z) is the interval hull of { Arg(z) | z in Z, z<>0 }. 00836 // 00837 // atan is the inverse function of tan(t), t in (-pi/2,pi/2). 00838 // 00839 l_interval Arg( const l_cinterval& z ) throw() 00840 { 00841 l_real 00842 srez = Sup( Re(z) ), 00843 irez = Inf( Re(z) ), 00844 simz = Sup( Im(z) ), 00845 iimz = Inf( Im(z) ); 00846 00847 l_interval 00848 hxl(irez), hxu(srez), hyl(iimz), hyu(simz); 00849 00850 l_real resl, resu; 00851 00852 if( iimz > 0.0 ) 00853 // case I: Im(z) > 0 00854 { 00855 resl = ( srez > 0.0 ? Inf( Atan( hyl,hxu ) ) : 00856 ( srez < 0.0 ? Inf( Atan( hyu,hxu ) + Pi_l_interval() ) : 00857 Inf( Pid2_l_interval() ) ) ); 00858 resu = ( irez > 0.0 ? Sup( Atan( hyu,hxl ) ) : 00859 ( irez < 0.0 ? Sup( Atan( hyl,hxl ) + Pi_l_interval() ) : 00860 Sup( Pid2_l_interval() ) ) ); 00861 return l_interval( resl, resu ); 00862 } 00863 else 00864 { 00865 if( simz < 0.0 ) 00866 // case II: Im(z) < 0 00867 { 00868 resl = ( irez < 0.0 ? Inf( Atan( hyu,hxl ) - Pi_l_interval() ) : 00869 ( irez > 0.0 ? Inf( Atan( hyl,hxl ) ) : 00870 -Sup( Pid2_l_interval() ) ) ); 00871 resu = ( srez < 0.0 ? Sup( Atan( hyl,hxu ) - Pi_l_interval() ) : 00872 ( srez > 0.0 ? Sup( Atan( hyu,hxu ) ) : 00873 -Inf( Pid2_l_interval() ) ) ); 00874 return l_interval( resl, resu ); 00875 } 00876 else 00877 // 0 in Im(z) 00878 { 00879 if( irez > 0.0 ) 00880 // case III: Re(z) > 0 00881 // z contains positive real values 00882 { 00883 resl = iimz < 0.0 ? Inf( Atan( hyl,hxl ) ) : l_real(0.0); 00884 return l_interval( resl, Sup( Atan( hyu,hxl ) ) ); 00885 } 00886 else 00887 // z contains nonpositive real numbers 00888 { 00889 if( irez < 0.0 ) 00890 { 00891 // case IV: z contains negative real numbers 00892 cxscthrow (STD_FKT_OUT_OF_DEF("l_interval Arg( const l_cinterval& z ); z contains negative real numbers")); 00893 return l_interval(0.0); 00894 } 00895 else 00896 // case V: 0 in z, but z doesn't contain negative real numbers 00897 { 00898 if( srez > 0.0 ) 00899 // diam( Re(z) > 0.0 ) 00900 { 00901 resl = iimz < 0.0 ? -Sup(Pid2_l_interval()) : 00902 l_real(0.0); 00903 resu = simz > 0.0 ? Sup(Pid2_l_interval()) : 00904 l_real(0.0); 00905 return l_interval( resl, resu ); 00906 } 00907 else 00908 // Re(z) == 0.0 00909 { 00910 if( iimz == 0.0 && simz == 0.0 ) 00911 // Z == 0 00912 return l_interval(0.0); 00913 else 00914 { 00915 resl = ( iimz < 0.0 ? - Sup( Pid2_l_interval() ) : 00916 Inf( Pid2_l_interval() ) ); 00917 resu = ( simz > 0.0 ? Sup( Pid2_l_interval() ) : 00918 -Inf( Pid2_l_interval() ) ); 00919 return l_interval( resl, resu ); 00920 } 00921 } 00922 } 00923 } 00924 } 00925 } 00926 } 00927 // 00928 //-- end Arg ------------------------------------------------------------------ 00929 00930 //-- arg: non-analytic argument function ---------------------------- 040916 -- 00931 // 00932 // (i) arg(Z) is defined for all Z in IC. 00933 // (ii) arg(Z) subset [-pi,3*pi/2]. 00934 // (iii) arg(Z) == Arg(Z) if Arg(Z) is well-defined. 00935 // 00936 // atan is the inverse function of tan(t), t in (-pi/2,pi/2). 00937 // 00938 00939 l_interval arg( const l_cinterval& z ) throw() 00940 { 00941 l_real 00942 srez = Sup( Re(z) ), 00943 irez = Inf( Re(z) ), 00944 simz = Sup( Im(z) ), 00945 iimz = Inf( Im(z) ); 00946 00947 l_real resl, resu; 00948 00949 if( irez < 0.0 && iimz <= 0.0 && simz >= 0.0 ) 00950 // z contains negative real values 00951 { 00952 if( srez > 0.0 ) 00953 // 0 in z and 0 interior point of Re(z) 00954 { 00955 resl = ( iimz < 0.0 ? - Sup( Pi_l_interval() ) : l_real(0.0) ); 00956 resu = ( ( iimz < 0.0 && simz == 0.0 ) ? l_real(0.0) : 00957 Sup( Pi_l_interval() ) ); 00958 return l_interval( resl, resu ); 00959 } 00960 else 00961 { // srez <= 0.0 00962 if( iimz == simz ) 00963 // z is real interval containing no positive values 00964 return Pi_l_interval(); 00965 else 00966 // sup( Re(z) ) <= 0, diam( Im(z) ) > 0 00967 { 00968 if( srez == 0.0 ) 00969 { 00970 resl = ( simz > 0.0 ? Inf( Pid2_l_interval() ) : 00971 -Sup( Pi_l_interval() ) ); 00972 resu = ( iimz < 0.0 ? 00973 ( simz > 0.0 ? Sup( 3 * Pid2_l_interval() ) : 00974 -Inf( Pid2_l_interval() ) ) : 00975 Sup( Pi_l_interval() ) ); 00976 return l_interval( resl, resu ); 00977 } 00978 else 00979 // sup( Re(z) ) < 0, diam( Im(z) ) > 0 00980 { 00981 l_interval hyl(iimz), hyu(simz); 00982 resl = ( simz > 0.0 ? 00983 Inf( Atan( hyu,srez ) + Pi_l_interval() ) : 00984 -Sup( Pi_l_interval() ) ); 00985 resu = ( iimz < 0.0 ? ( simz > 0.0 ? 00986 Sup( Atan( hyl,srez ) + 00987 Pi_l_interval() ) : 00988 Sup( Atan( hyl,srez ) - 00989 Pi_l_interval() ) ) : 00990 Sup( Pi_l_interval() ) ); 00991 return l_interval( resl, resu ); 00992 } 00993 } 00994 } 00995 } 00996 else 00997 // Arg(z) is well-defined 00998 return Arg( z ); 00999 } 01000 // 01001 //-- end arg ------------------------------------------------------------------ 01002 01003 01004 // *************************************************************************** 01005 // *************************************************************************** 01006 // *** Multi-valued functions **** 01007 // *************************************************************************** 01008 // *************************************************************************** 01009 01010 01011 //-- arg_inclmon: non-analytic inclusion-monotone argument function - 040617 -- 01012 // 01013 // (i) arg_inclmon(Z) is defined for all Z in IC. 01014 // (ii) arg_inclmon(Z) = [-pi,pi] if Arg(Z) is not defined. 01015 // 01016 01017 l_interval arg_inclmon( const l_cinterval& z ) throw() 01018 { 01019 if( Inf( Re(z) ) < 0.0 && Inf( Im(z) ) <= 0.0 && Sup( Im(z) ) >= 0.0 ) 01020 return l_interval( -Sup( Pi_l_interval() ),Sup( Pi_l_interval() ) ); 01021 else return Arg(z); 01022 } 01023 // 01024 //-- end arg_inclmon ---------------------------------------------------------- 01025 01026 //----------------------------------------------------------------------------- 01027 // 01028 // Section 2: Logarithms 01029 // 01030 //----------------------------------------------------------------------------- 01031 01032 //-- Ln: analytic natural logarithm --------------------------------- 040625 -- 01033 // 01034 // Ln(z) is undefined if z contains zero; z must not touch the negative real 01035 // axis from below; 01036 // 01037 l_cinterval Ln( const l_cinterval& z ) throw() 01038 { // Blomquist; 01039 int stagsave = stagprec, 01040 stagmax = 19; 01041 if (stagprec < stagmax) stagprec++; 01042 else stagprec = stagmax; 01043 01044 l_cinterval y; 01045 l_real 01046 srez = Sup( Re(z) ), 01047 simz = Sup( Im(z) ), 01048 iimz = Inf( Im(z) ); 01049 l_interval a1( abs(Re(z)) ), 01050 a2( abs(Im(z)) ); 01051 if ( Inf(a1) == 0.0 && Inf(a2) == 0.0 ) 01052 // z contains 0 01053 cxscthrow(STD_FKT_OUT_OF_DEF("l_cinterval Ln( const l_cinterval& z ); z contains 0")); 01054 if ( srez<0 && iimz<0 && simz>=0 ) 01055 cxscthrow(STD_FKT_OUT_OF_DEF("l_cinterval Ln( const l_cinterval& z ); z not allowed")); 01056 y = l_cinterval( ln_sqrtx2y2(Re(z),Im(z)), arg(z) ); 01057 01058 stagprec = stagsave; 01059 y = adjust(y); 01060 01061 return y; 01062 } 01063 // 01064 //-- end Ln ------------------------------------------------------------------- 01065 01066 //-- ln: non-analytic natural logarithm ----------------------------- 040923 -- 01067 // 01068 // ln(z) is undefined if z contains zero. 01069 // 01070 l_cinterval ln( const l_cinterval& z ) throw() 01071 { 01072 int stagsave = stagprec, 01073 stagmax = 19; 01074 if (stagprec < stagmax) stagprec++; 01075 else stagprec = stagmax; 01076 01077 l_cinterval y; 01078 l_interval a1( abs(Re(z)) ), 01079 a2( abs(Im(z)) ); 01080 if ( Inf(a1) == 0.0 && Inf(a2) == 0.0 ) 01081 // z contains 0 01082 cxscthrow(STD_FKT_OUT_OF_DEF 01083 ("l_cinterval ln( const l_cinterval& z ); z contains 0")); 01084 y = l_cinterval( ln_sqrtx2y2(Re(z),Im(z)), arg(z) ); 01085 01086 stagprec = stagsave; 01087 y = adjust(y); 01088 01089 return y; 01090 } 01091 // 01092 //-- end ln ------------------------------------------------------------------- 01093 01094 l_cinterval lnp1(const l_cinterval& z) throw() 01095 { // ln(1+z); 01096 // Calculates nearly optimal inclusions for not too wide intervals z. 01097 // Blomquist, 03.12.2008; 01098 int stagsave = stagprec, 01099 stagmax = 30; 01100 if (stagprec > stagmax) 01101 stagprec = stagmax; 01102 01103 const real c1 = 1.0; 01104 l_cinterval y; 01105 l_interval abs_z(abs(z)); 01106 l_real 01107 srez = Sup( Re(z) ), 01108 simz = Sup( Im(z) ), 01109 iimz = Inf( Im(z) ); 01110 01111 if (-1 <= z) // z contains -1 01112 cxscthrow(STD_FKT_OUT_OF_DEF( 01113 "l_cinterval lnp1(const l_cinterval& z); z contains -1")); 01114 if ( srez<-1 && iimz<0 && simz>=0 ) 01115 cxscthrow(STD_FKT_OUT_OF_DEF( 01116 "l_cinterval lnp1(const l_cinterval& z); z not allowed")); 01117 01118 if (Sup(abs_z) < c1) 01119 { 01120 abs_z = Re(z); 01121 abs_z = lnp1( abs_z*(2+abs_z) + sqr(Im(z)) ); 01122 times2pown(abs_z,-1); 01123 y = l_cinterval( abs_z, arg(1+z) ); 01124 } 01125 else 01126 y = Ln(1+z); 01127 stagprec = stagsave; 01128 y = adjust(y); 01129 01130 return y; 01131 } 01132 01133 l_cinterval log2( const l_cinterval& z ) throw() 01134 { 01135 return Ln(z) / Ln2_l_interval(); 01136 } 01137 01138 l_cinterval log10( const l_cinterval& z ) throw() 01139 { 01140 return Ln(z) / Ln10_l_interval(); 01141 } 01142 01143 //----------------------------------------------------------------------------- 01144 // 01145 // Section 3: Root functions 01146 // 01147 //----------------------------------------------------------------------------- 01148 01149 l_interval Sqrt_zpx( const l_interval& x, const l_interval& y, int& d ) 01150 // Calculating: 01151 // res = sqrt(2^(-d)*|z|+2^(-d)*|x|) = 2^(-d/2) * sqrt(|z| + |x|) 01152 // without any internal overflow; 01153 // Notice: |z| = sqrt(x^2+y^2); sqrt(|z|+|x|) < Maxreal; 01154 // Blomquist, 05.03.2007; 01155 { 01156 const int c1 = 1021; 01157 l_real Infx(Inf(x)), Infy(Inf(y)); 01158 int ex_x(expo_gr(Infx)), ex_y(expo_gr(Infy)); 01159 l_interval xc(abs(x)),yc(y),res; 01160 bool yeq0(Infy==0); 01161 d = 0; 01162 01163 if ((ex_x>=c1) || (ex_y>=c1)) 01164 { // To avoid overflow, scaling is necessary: 01165 times2pown(xc,-2); 01166 if (yeq0) { 01167 times2pown(xc,1); 01168 res = sqrt(xc); 01169 } else { 01170 times2pown(yc,-2); 01171 res = sqrt( sqrtx2y2(xc,yc)+xc); 01172 } 01173 times2pown(res,1); // Backscaling 01174 } else // Normal calculation without scaling: 01175 if (yeq0) 01176 { 01177 times2pown(xc,1); 01178 res = sqrt(xc); 01179 } else 01180 { 01181 if (ex_y>ex_x) ex_x = ex_y; // ex_x = Max{ex_x,ex_y} 01182 if (ex_x%2 != 0) ex_x--; 01183 if (ex_x<-50) 01184 { 01185 d = ex_x; // d is divisible by 2; 01186 Times2pown(xc,-d); 01187 Times2pown(yc,-d); 01188 res = sqrt( sqrtx2y2(xc,yc)+xc); 01189 } else res = sqrt( sqrtx2y2(xc,yc)+xc); 01190 } 01191 return res; 01192 } // Sqrt_zpx() 01193 01194 //-- sqrt: analytic square root ------------------------------------- 040626 -- 01195 // 01196 l_interval Re_Sqrt_point( const l_interval& rez, const l_interval& imz ) 01197 // 01198 // Real part of analytic square root of A POINT INTERVAL ONLY. 01199 // Do not use this as a general function 01200 // - it's only a subroutine for sqrt and sqrt_all. 01201 // The calculation is void if (rez+I*imz) is not a complex number. 01202 // Blomquist, 05.03.2007; 01203 { 01204 int d; 01205 l_interval res; 01206 l_real 01207 irez = Inf( rez ), 01208 iimz = Inf( imz ); 01209 01210 if( iimz == 0.0 ) 01211 { 01212 if( irez >= 0.0 ) 01213 return sqrt( rez ); 01214 else return l_interval(0.0); 01215 } 01216 else 01217 { // iimz <> 0 01218 res = Sqrt_zpx(rez,imz,d); 01219 if (irez >= 0.0) 01220 { 01221 if (d!=0) times2pown(res,d/2); 01222 return Sqrt2r_l_interval() * res; 01223 } 01224 else 01225 { 01226 iimz = abs(iimz); 01227 if (d!=0) times2pown(iimz,-d/2); 01228 return Sqrt2r_l_interval() * iimz / res; 01229 } 01230 } 01231 } 01232 01233 l_interval Im_Sqrt_point( const l_interval& rez, const l_interval& imz ) 01234 // 01235 // Imaginary part of analytic square root of A POINT INTERVAL ONLY 01236 // Do not use this as a general function 01237 // - it's only a subroutine for sqrt and sqrt_all 01238 // The calculation is void if (rez+I*imz) is not a complex number 01239 // Blomquist, 05.03.2007; 01240 { 01241 int d; 01242 l_interval res; 01243 l_real 01244 irez = Inf( rez ), 01245 iimz = Inf( imz ); 01246 01247 if( iimz == 0.0 ) 01248 { 01249 if( irez >= 0.0 ) return l_interval(0.0); 01250 else return sqrt(-rez); 01251 } 01252 else 01253 { 01254 res = Sqrt_zpx(rez,imz,d); 01255 if( irez >= 0.0 ) 01256 { 01257 if (d!=0) times2pown(iimz,-d/2); 01258 return Sqrt2r_l_interval() * iimz / res; 01259 } 01260 else 01261 { 01262 if (d!=0) times2pown(res,d/2); 01263 res = Sqrt2r_l_interval() * res; 01264 if( iimz > 0.0 ) return res; 01265 else return -res; 01266 } 01267 } 01268 } 01269 01270 l_cinterval sqrt( const l_cinterval& z ) throw() // ---------------------- 01271 // 01272 // Analytic square root function with z in the principal branch. 01273 // The branch cut is the negative real axis. On the branch cut 01274 // the function values are defined by comming from the II quadrant. 01275 // Blomquist, 23.06.2005; 01276 // 01277 { 01278 int stagsave = stagprec, 01279 stagmax = 19; 01280 if (stagprec < stagmax) stagprec++; 01281 else stagprec = stagmax; 01282 01283 l_cinterval y; 01284 l_real 01285 irez = Inf( Re(z) ), 01286 srez = Sup( Re(z) ), 01287 iimz = Inf( Im(z) ), 01288 simz = Sup( Im(z) ); 01289 l_interval 01290 hxl( irez ), hxu( srez ), hyl( iimz ), hyu( simz ); 01291 l_real 01292 resxl, resxu, resyl, resyu; 01293 01294 if( irez < 0.0 && iimz < 0.0 && simz >= 0.0 ) 01295 // if z is not in the principal branch then the inclusion monotony 01296 // is violated! 01297 cxscthrow(STD_FKT_OUT_OF_DEF( 01298 "l_cinterval sqrt(const l_cinterval& z); z not in the principal branch.")); 01299 01300 if (iimz>=0.0) 01301 { 01302 resxl = Inf( Re_Sqrt_point( hxl,hyl ) ); 01303 resxu = Sup( Re_Sqrt_point( hxu,hyu ) ); 01304 01305 resyl = Inf( Im_Sqrt_point( hxu,hyl ) ); 01306 resyu = Sup( Im_Sqrt_point( hxl,hyu ) ); 01307 } else 01308 if (simz<=0.0) { 01309 resxl = Inf( Re_Sqrt_point( hxl,hyu ) ); 01310 resxu = Sup( Re_Sqrt_point( hxu,hyl ) ); 01311 01312 resyl = Inf( Im_Sqrt_point( hxl,hyl ) ); 01313 resyu = Sup( Im_Sqrt_point( hxu,hyu ) ); 01314 } else { 01315 resxl = Inf( sqrt( hxl ) ); 01316 resxu = ( -iimz > simz ? Sup( Re_Sqrt_point( hxu, hyl ) ) 01317 : Sup( Re_Sqrt_point( hxu, hyu ) ) ); 01318 resyl = Inf( Im_Sqrt_point( hxl,hyl ) ); 01319 resyu = Sup( Im_Sqrt_point( hxl,hyu ) ); 01320 } 01321 y = l_cinterval( l_interval( resxl,resxu ), l_interval( resyl,resyu ) ); 01322 01323 stagprec = stagsave; 01324 y = adjust(y); 01325 01326 return y; 01327 } 01328 01329 l_cinterval sqrtp1m1(const l_cinterval& z) throw() 01330 // sqrt(1+z)-1; 01331 // Calculates nearly optimal inclusions for not too wide intervals z. 01332 // Blomquist, 08.07.2008; 01333 { 01334 const real c = 0.125; 01335 int stagsave = stagprec, 01336 stagmax = 30; // l_imath.cpp: sqrt() uses stagmax = 30; 01337 if (stagprec>stagmax) stagprec = stagmax; 01338 01339 l_cinterval res; 01340 l_interval absz(abs(z)); 01341 l_real Sup_absz(Sup(absz)); 01342 01343 if (Sup_absz < c) 01344 res = z / (sqrt(z+1) + 1); 01345 else 01346 res = sqrt(z+1) - 1; 01347 01348 stagprec = stagsave; 01349 res = adjust(res); 01350 return res; 01351 } 01352 01353 l_cinterval sqrt1px2(const l_cinterval& z) throw() 01354 // sqrt(1+z^2); 01355 // Calculates nearly optimal inclusions for not too wide intervals z. 01356 // Blomquist, 03.07.2008; 01357 { 01358 const l_real c = l_real(1e152); 01359 int stagsave = stagprec, 01360 stagmax = 30; // l_imath.cpp: sqrt() uses stagmax = 30; 01361 if (stagprec>stagmax) stagprec = stagmax; 01362 01363 l_cinterval res; 01364 l_interval absz(abs(z)); 01365 l_real Inf_absz(Inf(absz)); 01366 01367 if (Inf_absz > c) 01368 { 01369 absz = 1 / l_interval(Inf_absz); 01370 Inf_absz = Sup(absz); 01371 res = l_cinterval( l_interval(-Inf_absz,Inf_absz), 01372 l_interval(-Inf_absz,Inf_absz) ); 01373 // res is the correcture interval, i.e. 01374 // z + res or -z + res is the inclusionof sqrt{1+z^2} 01375 res = (Inf(Re(z))>=0)? z + res : -z + res; 01376 } 01377 else 01378 { 01379 res = l_cinterval( l_interval(0), l_interval(1) ); // res = i 01380 if ( Sup(abs(z-res))<0.5 || Sup(abs(z+res))<0.5 ) 01381 { 01382 res = l_cinterval(-Im(z),Re(z)); // Res = i*z; 01383 // (1 - i*z)*(1 + i*z) = 1+z^2; 01384 res = sqrt( (1-res)*(1+res) ); 01385 } 01386 else 01387 res = sqrt(1+sqr(z)); 01388 } 01389 if (Inf(Re(res))<0) 01390 res = l_cinterval( l_interval(l_real(0),Sup(Re(res))) , Im(res) ); 01391 stagprec = stagsave; 01392 res = adjust(res); 01393 return res; 01394 } 01395 // -- end sqrt1px2 --------------------------------------------------------- 01396 01397 l_cinterval sqrtx2m1(const l_cinterval& z) throw() 01398 // sqrt(z^2-1); 01399 // Calculates nearly optimal inclusions for not too wide intervals z. 01400 // Blomquist, 04.12.2008; 01401 { 01402 const l_real c = l_real(1e152); 01403 int stagsave = stagprec, 01404 stagmax = 30; // l_imath.cpp: sqrt() uses stagmax = 30; 01405 if (stagprec>stagmax) stagprec = stagmax; 01406 01407 l_cinterval res,u; 01408 l_interval absz(abs(z)); 01409 l_real Inf_absz(Inf(absz)); 01410 01411 if (Inf_absz > c) 01412 { 01413 absz = 1 / l_interval(Inf_absz); 01414 Inf_absz = Sup(absz); 01415 res = l_cinterval(l_interval(-Inf_absz,Inf_absz), 01416 l_interval(-Inf_absz,Inf_absz)); // res = Delta 01417 // res is the correcture interval, i.e. 01418 res = (Inf(Re(z))>=0)? z + res : -z + res; 01419 } 01420 else 01421 { 01422 res = z-1; u = z+1; 01423 res = (Sup(abs(res))<0.5 || Sup(abs(u))<0.5)? sqrt(res*u) : sqrt(sqr(z)-1); 01424 } 01425 01426 if (Inf(Re(res))<0) 01427 res = l_cinterval(l_interval(l_real(0),Sup(Re(res))),Im(res)); 01428 01429 stagprec = stagsave; 01430 res = adjust(res); 01431 return res; 01432 } 01433 01434 l_cinterval sqrt1mx2(const l_cinterval& z) throw() 01435 // sqrt(1-z^2); 01436 // Calculates nearly optimal inclusions for not too wide intervals z. 01437 // Blomquist, 04.12.2008; 01438 { 01439 const l_real c = l_real(1e152); 01440 int stagsave = stagprec, 01441 stagmax = 30; // l_imath.cpp: sqrt() uses stagmax = 30; 01442 if (stagprec>stagmax) stagprec = stagmax; 01443 01444 l_cinterval res,u; 01445 l_interval absz(abs(z)); 01446 l_real Inf_absz(Inf(absz)); 01447 01448 if (Inf_absz > c) 01449 { 01450 absz = 1 / l_interval(Inf_absz); 01451 Inf_absz = Sup(absz); 01452 res = l_cinterval(l_interval(-Inf_absz,Inf_absz), 01453 l_interval(-Inf_absz,Inf_absz)); // res = Delta 01454 u = l_cinterval(-Im(z),Re(z)); // u = i*z; 01455 // res is the correcture interval, i.e. 01456 // i*z + res or -i*z + res is the inclusion of sqrt{1-z^2} 01457 res = (Inf(Im(z))>=0)? -u + res : u + res; 01458 } 01459 else 01460 { 01461 res = 1-z; u = 1+z; 01462 res = (Sup(abs(res))<0.5 || Sup(abs(u))<0.5)? sqrt(res*u) : sqrt(1-sqr(z)); 01463 } 01464 if (Inf(Re(res))<0) 01465 res = l_cinterval( l_interval(l_real(0),Sup(Re(res))) , Im(res) ); 01466 stagprec = stagsave; 01467 res = adjust(res); 01468 return res; 01469 } 01470 01471 01472 //-- sqrt_all ------------------------------------------------------- 040621 -- 01473 // 01474 // sqrt_all(z) computes a list of 2 intervals containing all square roots of z 01475 // 01476 std::list<l_cinterval> sqrt_all( const l_cinterval& z ) 01477 { 01478 l_real 01479 irez = Inf( Re(z) ), 01480 srez = Sup( Re(z) ), 01481 iimz = Inf( Im(z) ), 01482 simz = Sup( Im(z) ); 01483 l_interval 01484 hxl( irez ), hxu( srez ), hyl( iimz ), hyu( simz ); 01485 l_real 01486 resxl, resxu, resyl, resyu; 01487 l_cinterval w; 01488 01489 if( irez < 0.0 && iimz <= 0.0 && simz >= 0.0 ) 01490 // z contains negative real values 01491 { 01492 if( iimz == 0.0 ) 01493 // z in upper half plane 01494 // principal values can be used 01495 { 01496 // min( Re ( sqrt( z ) ) ) in lower left corner 01497 // max( Re ( sqrt( z ) ) ) in upper right corner 01498 resxl = Inf( Re_Sqrt_point( hxl, hyl ) ); 01499 resxu = Sup( Re_Sqrt_point( hxu, hyu ) ); 01500 // min( Im ( sqrt( z ) ) ) in lower right corner 01501 // max( Im ( sqrt( z ) ) ) in upper left corner 01502 resyl = Inf( Im_Sqrt_point( hxu, hyl ) ); 01503 resyu = Sup( Im_Sqrt_point( hxl, hyu ) ); 01504 } 01505 else 01506 { 01507 if( simz == 0.0 ) 01508 // z in lower half plane 01509 // principal values can be used in lower corners 01510 { 01511 // z in lower half plane 01512 // min( Re ( sqrt( z ) ) ) in upper left corner 01513 // max( Re ( sqrt( z ) ) ) in lower right corner 01514 resxl = 0; ; 01515 resxu = Sup( Re_Sqrt_point( hxu, hyl ) ); 01516 // min( Im ( sqrt( z ) ) ) in lower left corner 01517 // max( Im ( sqrt( z ) ) ) in upper right corner 01518 resyl = Inf( Im_Sqrt_point( hxl, hyl ) ); 01519 resyu = ( srez > 0.0 ? l_real(0.0) : -Inf( sqrt( -hxu ) ) ); 01520 } 01521 else 01522 // 0 is interior point of Im( z ) 01523 { 01524 if( srez <= 0.0 ) 01525 { 01526 // 0 is no interior point of Re (z ) 01527 // in quadrant III, Re( z ) = Im_Sqrt_point(|x|,y), 01528 // Im( z ) = Re_Sqrt_point(|x|,y) 01529 // min( Re ( sqrt( z ) ) ) in lower right corner = quadrant II/III 01530 // max( Re ( sqrt( z ) ) ) in upper right corner = quadrant II 01531 resxl = Inf( Im_Sqrt_point( -hxu, hyl ) ); 01532 resxu = Sup( Re_Sqrt_point( hxu, hyu ) ); 01533 // min( Im ( sqrt( z ) ) ) on real line 01534 // max( Im ( sqrt( z ) ) ) in lower or upper left corner 01535 resyl = Inf( sqrt( -hxu ) ); 01536 resyu = ( - iimz > simz ? Sup( Re_Sqrt_point( -hxl, hyl ) ) : Sup( Im_Sqrt_point( hxl, hyu ) ) ); 01537 } 01538 else 01539 // 0 is interior point of z 01540 // here, the principal values apply in all corner points 01541 { 01542 resxl = 0; 01543 // max( Re ( sqrt( z ) ) ) in lower or upper right corner 01544 resxu = ( - iimz > simz ? Sup( Re_Sqrt_point( hxu, hyl ) ) : Sup( Re_Sqrt_point( hxu, hyu ) ) ); 01545 // min( Im ( sqrt( z ) ) ) in lower left corner 01546 // max( Im ( sqrt( z ) ) ) in upper left corner 01547 resyl = Inf( Im_Sqrt_point( hxl, hyl ) ); 01548 resyu = Sup( Im_Sqrt_point( hxl, hyu ) ); 01549 } 01550 } 01551 } 01552 w = l_cinterval( l_interval(resxl,resxu), l_interval(resyl,resyu ) ); 01553 } 01554 else 01555 // sqrt( z ) is well-defined 01556 w = sqrt( z ); 01557 01558 std::list<l_cinterval> res; 01559 res.push_back( w ); 01560 res.push_back( -w ); 01561 01562 return res; 01563 } 01564 // 01565 //-- end sqrt_all ------------------------------------------------------------- 01566 01567 //-- sqrt(z,n): analytic n-th root ---------------------------------- 040624 -- 01568 // 01569 l_interval Re_Sqrt_point( const l_interval& rez, const l_interval& imz, 01570 int n ) // before: unsigned int n ---------- 040624 -- 01571 // 01572 // Real part of analytic n-th root. 01573 // Do not use this as a general function 01574 // - it's only a subroutine for sqrt(z,n) and sqrt_all(z,n). 01575 // The calculation is validated but largely overestimating 01576 // if (rez+I*imz) is not a complex number. 01577 // 01578 { 01579 l_interval abs_z_2 = sqr( rez ) + sqr( imz ); 01580 if( Sup( abs_z_2 ) == 0.0 ) 01581 // z == 0 01582 return l_interval(0); 01583 else 01584 return sqrt( abs_z_2, 2 * n ) * 01585 cos( Arg( l_cinterval( rez, imz ) ) / n ); 01586 } 01587 01588 l_interval Im_Sqrt_point( const l_interval& rez, const l_interval& imz, 01589 int n ) // before: unsigned int n --- 040624 -- 01590 // 01591 // Imaginary part of analytic n-th root. 01592 // Do not use this as a general function 01593 // - it's only a subroutine for sqrt(z,n) and sqrt_all(z,n). 01594 // The calculation is validated but largely overestimating 01595 // if (rez+I*imz) is not a complex number. 01596 // 01597 { 01598 l_interval abs_z_2 = sqr( rez ) + sqr( imz ); 01599 if( Sup( abs_z_2 ) == 0.0 ) 01600 // z == 0 01601 return l_interval(0); 01602 else 01603 return sqrt( abs_z_2, 2 * n ) * 01604 sin( Arg( l_cinterval( rez, imz ) ) / n ); 01605 } 01606 01607 l_cinterval sqrt( const l_cinterval& z, int n ) throw() // ----- 040915 -- 01608 // 01609 // Analytic n-th root function 01610 // sqrt(z,n) is undefined if z contains negative real numbers. 01611 // 01612 { 01613 if( n == 0 ) 01614 return l_cinterval(l_interval(1)); 01615 if( n == 1 ) 01616 return z; 01617 if( n == 2 ) 01618 return sqrt( z ); 01619 else 01620 { 01621 l_real 01622 irez = Inf( Re(z) ), 01623 srez = Sup( Re(z) ), 01624 iimz = Inf( Im(z) ), 01625 simz = Sup( Im(z) ); 01626 l_interval 01627 hxl( irez ), hxu( srez ), hyl( iimz ), hyu( simz ); 01628 l_real 01629 resxl, resxu, resyl, resyu; 01630 01631 if( irez < 0.0 && iimz <= 0.0 && simz >= 0.0 ) 01632 { 01633 // z contains negative real values 01634 cxscthrow(STD_FKT_OUT_OF_DEF("l_cinterval sqrt(const l_cinterval& z, int n ); z contains negative real values.")); 01635 return z; 01636 } 01637 else 01638 { 01639 if( simz < 0.0 ) 01640 { 01641 // z in lower half plane 01642 l_cinterval hres = sqrt( l_cinterval( Re(z), -Im(z) ), n ); 01643 return l_cinterval( Re(hres), -Im(hres) ); 01644 } 01645 else 01646 { 01647 if( iimz > 0.0 ) 01648 { 01649 // z in upper half plane 01650 l_interval tangle = tan( ( Pi_l_interval() * n ) 01651 / ( 2 * ( n-1 ) ) ); 01652 l_real tanglel = Inf( tangle ), 01653 tangleu = Sup( tangle ); 01654 // 01655 // min( Re( Root( z ) ) ) 01656 // 01657 if ( irez >= 0.0 || Sup( hyl / irez ) <= tanglel ) 01658 // lower boundary right of phi = n*Pi/(2n-2) 01659 // min( Re( Root( z ) ) ) in lower left corner 01660 resxl = Inf( Re_Sqrt_point( hxl, hyl, n ) ); 01661 else 01662 { 01663 if( srez < 0.0 && Inf( hyl / srez ) >= tangleu ) 01664 // lower boundary left of phi = n*Pi/(2n-2) 01665 // min( Re( Root( z ) ) ) in lower right corner 01666 resxl = Inf( Re_Sqrt_point( hxu, hyl, n ) ); 01667 else 01668 // lower boundary intersects phi = n*Pi/(2n-2) 01669 // min( Re( Root( z ) ) ) on intersection 01670 resxl = Inf( Re_Sqrt_point( iimz / tangle , 01671 hyl, n ) ); 01672 } 01673 // 01674 // max( Re( Root( z ) ) ) 01675 // 01676 if ( irez >= 0.0 || Sup( hyu / irez ) <= tanglel ) 01677 // upper boundary right of phi = n*Pi/(2n-2) 01678 // max( Re( Root( z ) ) ) in upper right corner 01679 resxu = Sup( Re_Sqrt_point( l_interval(srez), 01680 l_interval(simz), n ) ); 01681 else 01682 { 01683 if ( srez < 0.0 && Inf( hyu / srez ) >= tangleu ) 01684 // upper boundary left of phi = n*Pi/(2n-2) 01685 // max( Re( Root( z ) ) ) in upper left corner 01686 resxu = Sup( Re_Sqrt_point( hxl, hyu, n ) ); 01687 else 01688 // upper boundary intersects phi = n*Pi/(2n-2) 01689 // max( Re(Root( z )) ) on upper left or right corner 01690 resxu = max( Sup( Re_Sqrt_point( hxl, hyu, n ) ), 01691 Sup( Re_Sqrt_point( hxu, hyu, n ) ) ); 01692 } 01693 // 01694 // min( Im( Root( z ) ) ) 01695 // 01696 if ( srez >= 0.0 || Sup( hyl / srez ) <= tanglel ) 01697 // right boundary right of phi = n*Pi/(2n-2) 01698 // min( Im( Root( z ) ) ) in lower right corner 01699 resyl = Inf( Im_Sqrt_point( hxu, hyl, n ) ); 01700 else 01701 { 01702 if( Inf( hyu / srez ) >= tangleu ) 01703 // right boundary left of phi = n*Pi/(2n-2) 01704 // min( Im( Root( z ) ) ) in upper right corner 01705 resyl = Inf( Im_Sqrt_point( hxu, hyu, n ) ); 01706 else 01707 // right boundary intersects phi = n*Pi/(2n-2) 01708 // min( Im( Root( z ) ) ) on intersection 01709 resyl = Inf(Im_Sqrt_point( hxu, srez * tangle, n )); 01710 } 01711 // 01712 // max( Im( Root( z ) ) ) 01713 // 01714 if( irez >= 0.0 || Sup( hyl / irez ) <= tanglel ) 01715 // left boundary right of phi = n*Pi/(2n-2) 01716 // max( Im( Root( z ) ) ) in upper left corner 01717 resyu = Sup( Im_Sqrt_point( hxl, hyu, n ) ); 01718 else 01719 { 01720 if( Inf( hyu / irez ) >= tangleu ) 01721 // left boundary left of phi = n*Pi/(2n-2) 01722 // max( Im( Root( z ) ) ) in lower left corner 01723 resyu = Sup( Im_Sqrt_point( hxl, hyl, n ) ); 01724 else 01725 // left boundary intersects phi = n*Pi/(2n-2) 01726 // max( Im(Root( z )) ) on lower or upper left corner 01727 resyu = max( Sup( Im_Sqrt_point( hxl, hyl, n ) ), 01728 Sup( Im_Sqrt_point( hxl, hyu, n ) ) ); 01729 } 01730 } 01731 else 01732 { 01733 // z intersects positive real axis 01734 // min( Re( Root( z ) ) ) on real line 01735 // max( Re( Root( z ) ) ) in lower or upper right corner 01736 resxl = ( irez == 0.0 ? l_real(0.0) : 01737 Inf( sqrt( hxl, n ) ) ); 01738 resxu = ( - iimz > simz ? 01739 Sup( Re_Sqrt_point( hxu, hyl, n ) ) : 01740 Sup( Re_Sqrt_point( hxu, hyu, n ) ) ); 01741 // min( Im ( sqrt( z ) ) ) in lower left corner 01742 // max( Im ( sqrt( z ) ) ) in upper left corner 01743 resyl = Inf( Im_Sqrt_point( hxl, hyl, n ) ); 01744 resyu = Sup( Im_Sqrt_point( hxl, hyu, n ) ); 01745 } 01746 return l_cinterval( l_interval( resxl, resxu ), 01747 l_interval( resyl, resyu ) ); 01748 } 01749 } 01750 } 01751 } 01752 // 01753 //-- end sqrt ----------------------------------------------------------------- 01754 01755 //-- sqrt_all ------------------------------------------------------- 040628 -- 01756 // 01757 std::list<l_cinterval> sqrt_all( const l_cinterval& z, int n ) 01758 // 01759 // sqrt_all(z,n) computes a list of n intervals containing all n-th roots of z 01760 // 01761 // For n >=3, computing the optimal interval bounds is very expensive 01762 // and thus not considered cost-effective. 01763 // 01764 // Hence, the polar form is used to calculate sqrt_all(z,n) 01765 // (observed overestimation of Re() and Im() in test cases: 5-15%). 01766 // 01767 // z is enclosed into an interval in polar coordinates 01768 // (i.e. a sector of an annulus), from which the roots 01769 // (also polar intervals) are computed error-free 01770 // (save roundoff, which is enclosed). 01771 // 01772 // The the final inclusion of the roots into a rectangular complex interval 01773 // involves some additional overestimation. 01774 // 01775 { 01776 std::list<l_cinterval> res; 01777 01778 if( n == 0 ) 01779 { 01780 res.push_back( l_cinterval( l_interval(1), l_interval(0) ) ); 01781 return res; 01782 } 01783 else if( n == 1 ) 01784 { 01785 res.push_back(z); 01786 return res; 01787 } 01788 else if( n == 2 ) return sqrt_all( z ); 01789 else 01790 { 01791 l_interval 01792 arg_z = arg( z ), root_abs_z = sqrt( abs( z ), n ); 01793 01794 for(int k = 0; k < n; k++) 01795 { 01796 l_interval arg_k = ( arg_z + 2 * k * Pi_l_interval() ) / n; 01797 01798 res.push_back( l_cinterval( root_abs_z * cos( arg_k ), 01799 root_abs_z * sin( arg_k ) ) ); 01800 } 01801 return res; 01802 } 01803 } 01804 // 01805 //-- end sqrt_all ------------------------------------------------------------- 01806 01807 01808 //----------------------------------------------------------------------------- 01809 // 01810 // Section 4: Power functions 01811 // 01812 //----------------------------------------------------------------------------- 01813 01814 //-- power_fast ----------------------------------------------------- 040714 -- 01815 // 01816 // Fast, validated power function for integer powers, based on exp and ln. 01817 // Medium amount of overestimation. 01818 // 01819 l_cinterval power_fast( const l_cinterval& z, int n ) throw() 01820 { 01821 if( n == 0 ) 01822 return l_cinterval( l_interval(1) ); 01823 else if( n == 1 ) 01824 return z; 01825 else if( n == -1 ) 01826 return 1 / z; 01827 else if( n == 2 ) 01828 return sqr(z); 01829 else 01830 // n >= 3 or n <= -2 01831 // If z is a large interval, z^n is a distorted set, for which the 01832 // inclusion into a complex rectangle contains large overestimation. 01833 // For example, if n * Arg( z ) intersects the cartesian axes 01834 // more than once then 0 is contained in the rectangular enclosure 01835 // of z^n. 01836 // For n <= -2, also inversion of z or z^n is required; 01837 // both inversions lead to large overestimation of the resulting interval. 01838 // 01839 // The computation of an optimal rectangular enclosure is implemented 01840 // in power(z,n). In power_fast(z,n), z is enclosed into a sector S of 01841 // an annulus. S^n is again some sector S' of a (different) annulus. 01842 // S' is calculated exactly (apart from roundoff), and then enclosed 01843 // into a rectangle. There is a certain amount of overestimation 01844 // compared with the optimal rectangular enclosure of z^n, but the 01845 // calculations are much cheaper here. 01846 { 01847 l_interval abs_z = abs(z); 01848 01849 if( n < 0 && Inf( abs_z ) == 0.0 ) 01850 // z contains 0 01851 cxscthrow (STD_FKT_OUT_OF_DEF("l_cinterval power_fast(const l_cinterval& z, int n ); z contains 0.")); 01852 if( Sup( abs_z ) == 0.0 ) 01853 return l_cinterval( l_interval(0)); 01854 else 01855 { 01856 l_interval arg_z = arg( z ); 01857 l_interval abs_z_n = exp( n * ln( abs_z ) ); 01858 01859 return l_cinterval( abs_z_n * cos( n * arg_z ), 01860 abs_z_n * sin( n * arg_z ) ); 01861 } 01862 } 01863 } 01864 // 01865 //-- end power_fast ----------------------------------------------------------- 01866 01867 //-- power ---------------------------------------------------------- 040720 -- 01868 // 01869 // Power function for integer powers with optimal (save roundoff) accuracy. 01870 // 01871 // The extremal curves of Re(z^n) and Im(z^n) are straight lines defined 01872 // by the rays arg(z) = k pi/ ( 2n - 2 ), k = 0,...4n-5. 01873 // 01874 // Intersections of these rays with the boundary of z are called 01875 // "ray intersections" in the following. 01876 // 01877 l_cinterval power_point( const l_cinterval& z, int n ) //----------- 040715 -- 01878 // 01879 // z^n for A POINT INTERVAL z ONLY. 01880 // Do not use this as a general function 01881 // - it's only a subroutine for power. 01882 // The calculation may break down otherwise. 01883 // The case 0 in z for negative n is handled in power(z,n). 01884 // 01885 { 01886 if( Inf( Re(z) ) == 0.0 && Inf( Im(z) ) == 0.0 ) 01887 // if z is a valid point interval, it must be 0 01888 return l_cinterval( l_interval(0) ); 01889 else 01890 { 01891 l_interval ln_absz = ln_sqrtx2y2(Re(z),Im(z)); 01892 l_interval arg_z = arg(z); 01893 l_interval abs_z_n = exp( n * ln_absz ); 01894 01895 return l_cinterval( abs_z_n * cos( n * arg_z ), 01896 abs_z_n * sin( n * arg_z ) ); 01897 } 01898 } 01899 01900 void update_res( const l_cinterval& res, //------------------------- 040716 -- 01901 l_real& resxl, l_real& resxu, l_real& resyl, l_real& resyu ) 01902 // Subroutine of power(z,n). 01903 // Update extremal values of power function. 01904 { 01905 resxl = min( resxl, Inf( Re(res) ) ); 01906 resxu = max( resxu, Sup( Re(res) ) ); 01907 resyl = min( resyl, Inf( Im(res) ) ); 01908 resyu = max( resyu, Sup( Im(res) ) ); 01909 } 01910 01911 int trunc(const l_real& x) 01912 { 01913 l_real y(x); y += 0; 01914 bool neg(y[1]<0); 01915 if (neg) y = -y; // y >= 0; 01916 real r(y[1]); 01917 int res = int( _double(r) ); 01918 y = y - res; 01919 if (y[1]<0.0 && res>0) res -= 1; 01920 if (neg) res = -res; 01921 return res; 01922 } 01923 01924 void horizontal_check( //-------------------------------------------- 040720 -- 01925 const l_interval& hy, l_real hyy, const l_interval& arg_h, 01926 l_real irez, l_real srez, 01927 l_real& resxl, l_real& resxu, l_real& resyl, l_real& resyu, int n ) 01928 // Subroutine of power(z,n). 01929 // Check all relevant ray intersections on a horizontal boundary. 01930 { 01931 int r_int, il1, il2, iu1, iu2; 01932 l_cinterval res; 01933 // 01934 // Here, the most difficult task is 01935 // to determine the relevant ray intersections 01936 // Both arg_hl and n can be negative, which makes it very complicated 01937 // to decide the correct indexes for the "rightmost" intersections, etc. 01938 // 01939 // To simplify the distinction of cases, we introduce the variable 01940 // nofrays (number of rays) = abs( n-1 ) 01941 // 01942 // Note that we still have to pass n to power_point 01943 // 01944 iu1 = n-1; if (iu1<0) iu1 = -iu1; 01945 int nofrays = iu1; 01946 l_real arg_hlR = Inf( 2 * nofrays * arg_h / Pi_l_interval() ); 01947 if( arg_hlR[1] >= 0.0 ) 01948 r_int = trunc(arg_hlR); 01949 else 01950 r_int = trunc(arg_hlR - 1.0); 01951 il1 = r_int + 1; 01952 01953 l_real arg_huR = Sup( 2 * nofrays * arg_h / Pi_l_interval() ); 01954 if( arg_huR[1] >= 0.0 ) 01955 r_int = trunc(arg_huR); 01956 else 01957 r_int = trunc(arg_huR - 1.0); 01958 01959 iu1 = r_int; 01960 01961 if( iu1 >= il1 ) 01962 // at least one ray intersection 01963 // maybe more? 01964 { 01965 if( iu1 > il1 + 3 ) { 01966 // 01967 // we're in trouble: 01968 // there are more than 4 ray intersections 01969 // now we must decide which of these are relevant 01970 // depending on il1, iu1, n and on the quadrants, 01971 // 4 to 7 intersections must be considered 01972 // 01973 if( n > 0 ) 01974 // outmost intersections are relevant 01975 { 01976 if( irez >= 0.0 ) 01977 // z in quadrants I or IV: 01978 // only 4 rightmost intersections are relevant 01979 { 01980 if( hyy < 0.0 ) 01981 // quadrant IV 01982 il1 = iu1 - 3; 01983 else 01984 // quadrant I 01985 iu1 = il1 + 3; 01986 } 01987 else if( srez <= 0.0 ) 01988 // z in quadrants II or III: 01989 // only 4 leftmost intersections are relevant 01990 { 01991 if( hyy > 0.0 ) 01992 // quadrant II 01993 il1 = iu1 - 3; 01994 else 01995 // quadrant III 01996 iu1 = il1 + 3; 01997 } 01998 else 01999 // z intersects imaginary axes 02000 // we may have two lists of intersections! 02001 { 02002 iu2 = iu1; 02003 il2 = iu2 - 3; 02004 iu1 = il1 + 3; 02005 // remove intersection points that are contained in both lists 02006 if( il2 <= iu1 ) 02007 il2 = iu1 + 1; 02008 // 02009 // here, list 2 is processed 02010 // list 1 is processed below 02011 // 02012 for(int i = il2; i <= iu2; i++) 02013 { 02014 l_interval cotangle = cot( ( Pi_l_interval() * i ) / 02015 ( 2 * nofrays ) ); 02016 res = power_point( l_cinterval( hy * cotangle , hy ), n ); 02017 update_res( res, resxl, resxu, resyl, resyu ); 02018 } 02019 } 02020 } 02021 else 02022 // n < 0: 02023 // innermost intersections are relevant 02024 { 02025 if( irez >= 0.0 ) 02026 // z in quadrants I or IV: 02027 // only 4 leftmost intersections are relevant 02028 { 02029 if( hyy > 0.0 ) 02030 // quadrant I 02031 il1 = iu1 - 3; 02032 else 02033 // quadrant IV 02034 iu1 = il1 + 3; 02035 } 02036 else if( srez <= 0.0 ) 02037 // z in quadrants II or III: 02038 // only 4 rightmost intersections are relevant 02039 { 02040 if( hyy < 0.0 ) 02041 // quadrant III 02042 il1 = iu1 - 3; 02043 else 02044 // quadrant II 02045 iu1 = il1 + 3; 02046 } 02047 else 02048 // z intersects imaginary axes 02049 // we have one lists of 5 to 7 intersections 02050 { 02051 if( hyy > 0.0 ) 02052 { 02053 il2 = nofrays - 3; 02054 iu2 = nofrays + 3; 02055 } 02056 else 02057 { 02058 il2 = -nofrays - 3; 02059 iu2 = -nofrays + 3; 02060 } 02061 // list 1 contains all intersections 02062 // list 2 is a filter for all relevant intersections 02063 // store intersection of lists 1 and 2 in list 1 02064 il1 = ( il1 > il2 ? il1 : il2 ); 02065 iu1 = ( iu1 < iu2 ? iu1 : iu2 ); 02066 } 02067 } 02068 } 02069 // 02070 // list 1 has been left for processing 02071 // 02072 for(int i = il1; i <= iu1; i++) 02073 { 02074 l_interval cotangle = cot( ( Pi_l_interval() * i ) / 02075 ( 2 * nofrays ) ); 02076 res = power_point( l_cinterval( hy * cotangle , hy ), n ); 02077 update_res( res, resxl, resxu, resyl, resyu ); 02078 } 02079 } 02080 } 02081 02082 void vertical_check( //---------------------------------------------- 040720 -- 02083 const l_interval& hx, l_real hxx, const l_interval& arg_h, 02084 l_real iimz, l_real simz, 02085 l_real& resxl, l_real& resxu, l_real& resyl, l_real& resyu, int n ) 02086 // Subroutine of power(z,n). 02087 // Check all relevant ray intersections on a vertical boundary. 02088 { 02089 int r_int, il1, il2, iu1, iu2; 02090 l_cinterval res; 02091 // 02092 // Here, the most difficult task is 02093 // to determine the relevant ray intersections 02094 // Both arg_hl and n can be negative, which makes it very complicated 02095 // to decide the correct indexes for the topmost intersections, etc. 02096 // 02097 // To simplify the distinction of cases, we introduce the variable 02098 // nofrays (number of rays) = abs( n-1 ) 02099 // 02100 // Note that we still have to pass n to power_point 02101 // 02102 iu1 = n-1; if (iu1<0) iu1 = -iu1; 02103 int nofrays = iu1; 02104 02105 l_real arg_hlR = Inf( 2 * nofrays * arg_h / Pi_l_interval() ); 02106 if( arg_hlR[1] >= 0.0 ) 02107 r_int = trunc(arg_hlR); 02108 else 02109 r_int = trunc(arg_hlR - 1.0); 02110 il1 = r_int + 1; 02111 02112 l_real arg_huR = Sup( 2 * nofrays * arg_h / Pi_l_interval() ); 02113 if( arg_huR[1] >= 0.0 ) 02114 r_int = trunc(arg_huR); 02115 else 02116 r_int = trunc(arg_huR - 1.0); 02117 iu1 = r_int; 02118 02119 if( iu1 >= il1 ) 02120 // at least one ray intersection 02121 // maybe more? 02122 { 02123 if( iu1 > il1 + 3 ) { 02124 // 02125 // we're in trouble: 02126 // there are more than 4 ray intersections 02127 // now we must decide which of these are relevant 02128 // depending on il1, iu1, n and on the quadrants, 02129 // 4 to 7 intersections must be considered 02130 // 02131 if( n > 0 ) 02132 // outmost intersections are relevant 02133 { 02134 if( iimz >= 0.0 ) 02135 // z in quadrants I or II: 02136 // only 4 topmost intersections are relevant 02137 { 02138 if( hxx > 0.0 ) 02139 // quadrant I 02140 il1 = iu1 - 3; 02141 else 02142 // quadrant II 02143 iu1 = il1 + 3; 02144 } 02145 else if( simz <= 0.0 ) 02146 // z in quadrants III or IV: 02147 // only 4 lowest intersections are relevant 02148 { 02149 if( hxx < 0.0 ) 02150 // quadrant III 02151 il1 = iu1 - 3; 02152 else 02153 // quadrant IV 02154 iu1 = il1 + 3; 02155 } 02156 else 02157 // z intersects real axes 02158 // we may have two lists of intersections! 02159 { 02160 iu2 = iu1; 02161 il2 = iu2 - 3; 02162 iu1 = il1 + 3; 02163 // remove intersection points that are contained in both lists 02164 if( il2 <= iu1 ) 02165 il2 = iu1 + 1; 02166 // 02167 // here, list 2 is processed 02168 // list 1 is processed below 02169 // 02170 for(int i = il2; i <= iu2; i++) 02171 { 02172 l_interval tangle = tan( ( Pi_l_interval() * i ) / 02173 ( 2 * nofrays ) ); 02174 res = power_point( l_cinterval( hx, hx * tangle ), n ); 02175 update_res( res, resxl, resxu, resyl, resyu ); 02176 } 02177 } 02178 } 02179 else 02180 // n < 0: 02181 // innermost intersections are relevant 02182 { 02183 if( iimz >= 0.0 ) 02184 // z in quadrants I or II: 02185 // only 4 lowest intersections are relevant 02186 { 02187 if( hxx < 0.0 ) 02188 // quadrant II 02189 il1 = iu1 - 3; 02190 else 02191 // quadrant I 02192 iu1 = il1 + 3; 02193 } 02194 else if( simz <= 0.0 ) 02195 // z in quadrants III or IV: 02196 // only 4 topmost intersections are relevant 02197 { 02198 if( hxx > 0.0 ) 02199 // quadrant IV 02200 il1 = iu1 - 3; 02201 else 02202 // quadrant III 02203 iu1 = il1 + 3; 02204 } 02205 else 02206 // z intersects real axes 02207 // we have one lists of 5 to 7 intersections 02208 { 02209 if( hxx > 0.0 ) 02210 { 02211 il2 = -3; 02212 iu2 = 3; 02213 } 02214 else 02215 { 02216 il2 = 2 * nofrays - 3; 02217 iu2 = 2 * nofrays + 3; 02218 } 02219 // list 1 contains all intersections 02220 // list 2 is a filter for all relevant intersections 02221 // store intersection of lists 1 and 2 in list 1 02222 il1 = ( il1 > il2 ? il1 : il2 ); 02223 iu1 = ( iu1 < iu2 ? iu1 : iu2 ); 02224 } 02225 } 02226 } 02227 // 02228 // list 1 has been left for processing 02229 // 02230 for(int i = il1; i <= iu1; i++) 02231 { 02232 l_interval tangle = tan( ( Pi_l_interval() * i ) / 02233 ( 2 * nofrays ) ); 02234 res = power_point( l_cinterval( hx, hx * tangle ), n ); 02235 update_res( res, resxl, resxu, resyl, resyu ); 02236 } 02237 } 02238 } 02239 02240 l_cinterval power( const l_cinterval& z, int n ) throw() //---- 040720 -- 02241 // 02242 // Power function for integer powers with optimal (save roundoff) accuracy. 02243 // 02244 { 02245 if( n == 0 ) 02246 return l_cinterval( l_interval(1) ); 02247 else if( n == 1 ) 02248 return z; 02249 else if( n == -1 ) 02250 return 1 / z; 02251 else if( n == 2 ) 02252 return sqr(z); 02253 else 02254 // n >= 3 or n <= -2 02255 // 02256 // n may have a large absolute value and z may be a large interval. 02257 // In this case, we calculate function values at specific points 02258 // and look for min and max. 02259 // 02260 // An implementation with fewer function evaluations would be possible, 02261 // at the expense of an unreadable source code. 02262 { 02263 l_interval abs_z = abs(z); 02264 02265 if( n < 0 && Inf( abs_z ) == 0.0 ) 02266 // z contains 0 02267 cxscthrow (STD_FKT_OUT_OF_DEF("l_cinterval power(const l_cinterval& z, int n ); z contains 0.")); 02268 if( Sup( abs_z ) == 0.0 ) 02269 return l_cinterval( l_interval(0), l_interval(0) ); 02270 else 02271 { 02272 l_real 02273 irez = Inf(Re(z)), 02274 srez = Sup(Re(z)), 02275 iimz = Inf(Im(z)), 02276 simz = Sup(Im(z)); 02277 l_interval 02278 hxl( irez ), hxu( srez ), hyl( iimz ), hyu( simz ); 02279 l_real 02280 resxl, resxu, resyl, resyu; 02281 l_cinterval res; 02282 // 02283 // extremal values lie in corner points or on intersections of the 02284 // boundary of z with any ray arg( w ) = k*Pi/(2*(n-1)) 02285 // 02286 // first check the corners of z 02287 // 02288 res = power_point( l_cinterval( hxl, hyl ), n ); 02289 resxl = Inf( Re(res) ); 02290 resxu = Sup( Re(res) ); 02291 resyl = Inf( Im(res) ); 02292 resyu = Sup( Im(res) ); 02293 res = power_point( l_cinterval( hxu, hyl ), n ); 02294 update_res( res, resxl, resxu, resyl, resyu ); 02295 res = power_point( l_cinterval( hxl, hyu ), n ); 02296 update_res( res, resxl, resxu, resyl, resyu ); 02297 res = power_point( l_cinterval( hxu, hyu ), n ); 02298 update_res( res, resxl, resxu, resyl, resyu ); 02299 // 02300 // consider the origin, if it is a boundary point 02301 // (negative n have been taken care of above) 02302 // 02303 // if( irez * srez <= 0.0 and iimz * simz <= 0.0 02304 // and irez * srez * iimz * simz == 0.0 ) 02305 if ( 0<=Re(z) && 0<=Im(z) && 02306 (irez==0 || srez==0 || iimz==0 || simz==0 ) ) 02307 update_res( l_cinterval( l_interval(0), l_interval(0)), 02308 resxl, resxu, resyl, resyu ); 02309 // 02310 // now check for ray intersections 02311 // for each quadrant and each boundary, we must consider 02312 // the function values at at most 4 consecutive intersections 02313 // 02314 // Ia: lower boundary 02315 // 02316 if( iimz != 0.0 ) 02317 { 02318 l_interval arg_h = arg( l_cinterval( Re(z), hyl ) ); 02319 02320 // check if there are ray intersections 02321 horizontal_check( hyl, iimz, arg_h, irez, srez, 02322 resxl, resxu, resyl, resyu, n ); 02323 } 02324 // Ib: upper boundary 02325 if( simz != 0.0 ) 02326 { 02327 l_interval arg_h = arg( l_cinterval( Re(z), hyu ) ); 02328 02329 // check if there are ray intersections 02330 horizontal_check( hyu, simz, arg_h, irez, srez, 02331 resxl, resxu, resyl, resyu, n ); 02332 } 02333 // IIa: left boundary 02334 if( irez != 0.0 ) 02335 { 02336 l_interval arg_h = arg( l_cinterval( hxl, Im(z) ) ); 02337 02338 // check if there are ray intersections 02339 vertical_check( hxl, irez, arg_h, iimz, simz, 02340 resxl, resxu, resyl, resyu, n ); 02341 } 02342 // IIb: right boundary 02343 if( srez != 0.0 ) 02344 { 02345 l_interval arg_h = arg( l_cinterval( hxu, Im(z) ) ); 02346 02347 // check if there are ray intersections 02348 vertical_check( hxu, srez, arg_h, iimz, simz, 02349 resxl, resxu, resyl, resyu, n ); 02350 } 02351 return l_cinterval( l_interval( resxl, resxu ), 02352 l_interval( resyl, resyu ) ); 02353 } 02354 } 02355 } 02356 // 02357 //-- end power ---------------------------------------------------------------- 02358 02359 02360 //-- pow ------------------------------------------------------------ 040627 -- 02361 // 02362 // Analytic power function for real interval exponent, based on Ln. 02363 // 02364 02365 l_cinterval pow( const l_cinterval& z, const l_interval& p ) throw() 02366 // Neue Version von pow(...) von 040627 02367 { 02368 return exp( p*Ln(z) ); 02369 } 02370 02371 // 02372 //-- end pow ------------------------------------------------------------------ 02373 02374 //-- pow ------------------------------------------------------------ 040627 -- 02375 // 02376 // Analytic power function for complex interval exponent, based on Ln. 02377 // 02378 02379 l_cinterval pow( const l_cinterval& z, const l_cinterval& p ) throw() 02380 { 02381 return exp( p*Ln(z) ); 02382 } 02383 02384 // 02385 //-- end pow ------------------------------------------------------------------ 02386 02387 //-- pow_all -------------------------------------------------------- 041013 -- 02388 // 02389 // Non-analytic power function for real l_interval exponent. 02390 // 02391 // If 0 \not\in z, then compute four rectangular intervals that comprehend 02392 // an annulus which contains all values zeta^phi, zeta in z, phi in p. 02393 // 02394 // If 0 in z and negative reals in p, then abort execution 02395 // (potential modification: return the entire complex plane). 02396 // 02397 std::list<l_cinterval> pow_all( const l_cinterval& z, 02398 const l_interval& p ) throw() 02399 { 02400 l_interval abs_z = abs(z); 02401 02402 if( 0.0 < Inf( abs_z ) ) 02403 { 02404 l_interval abs_z_p = exp( p * ln( abs_z ) ); 02405 02406 // Inner and outer radii of the annulus are inf/sup( abs_z_n ) 02407 // Inscribed square has side length sqrt( 2 ) * rad_1 02408 l_interval rad_1 = Sqrt2r_l_interval() * Inf( abs_z_p ); 02409 l_interval rad_2 = l_interval(Sup( abs_z_p )); 02410 02411 std::list<l_cinterval> res; 02412 02413 // 4 intervals covering the annulus, counter-clockwise 02414 res.push_back(l_cinterval(l_interval( Inf( rad_1 ), Sup( rad_2 ) ), 02415 l_interval( -Sup( rad_1 ), Sup( rad_2 ) ))); 02416 res.push_back(l_cinterval(l_interval( -Sup( rad_2 ), Sup( rad_1 ) ), 02417 l_interval( Inf( rad_1 ), Sup( rad_2 ) ) )); 02418 res.push_back(l_cinterval(l_interval( -Sup( rad_2 ), -Inf( rad_1 ) ), 02419 l_interval( -Sup( rad_2 ), Sup(rad_1) ) ) ); 02420 res.push_back(l_cinterval(l_interval( -Sup( rad_1 ), Sup( rad_2 ) ), 02421 l_interval( -Sup( rad_2 ), -Inf(rad_1) ) ) ); 02422 02423 return res; 02424 } 02425 else 02426 // 0 in z 02427 { 02428 if( Inf( p ) > 0.0 ) 02429 // 02430 // All values zeta^phi, zeta in z, phi in p lie in a disk 02431 // 02432 { 02433 l_interval abs_z_p = exp( p * ln( l_interval( Sup( abs_z ) ) ) ); 02434 l_real rad_p = Sup( abs_z_p ); 02435 02436 std::list<l_cinterval> res; 02437 02438 res.push_back( l_cinterval( l_interval( -rad_p, rad_p ), 02439 l_interval( -rad_p, rad_p ) ) ); 02440 02441 return res; 02442 } 02443 else 02444 { 02445 // 02446 // The set zeta^phi, zeta in z, phi in p is unbounded 02447 // if inf( p ) < 0. 0^p is undefined for p <= 0. 02448 // 02449 cxscthrow(STD_FKT_OUT_OF_DEF("pow_all(l_cinterval, l_interval); 0^p is undefined for p <= 0.")); 02450 std::list<l_cinterval> res; 02451 return res; 02452 } 02453 } 02454 } 02455 // 02456 //-- end pow_all -------------------------------------------------------------- 02457 02458 02459 //----------------------------------------------------------------------------- 02460 // 02461 // Section 5: asin, acos, asinh, acosh 02462 // 02463 // The implementation of acos, asinh, and acosh is based on asin: 02464 // 02465 // acos( z ) = pi / 2 - asin( z ) 02466 // asinh( z ) = i * asin( - i * z ) 02467 // acosh( z ) = i * acos( z ) = +/- i * ( pi / 2 - asin( z ) ) 02468 // 02469 // Only principal values are computed. 02470 // 02471 //----------------------------------------------------------------------------- 02472 02473 //-- asin ----------------------------------------------------------- 040905 -- 02474 // 02475 // Analytic inverse sine function 02476 // 02477 l_interval f_aux_asin( const l_interval& x, const l_interval& y ) 02478 // 02479 // auxiliary function for evaluating the imaginary part of asin(z); 02480 // f_aux_asin(z) = ( |z+1| + |z-1| ) / 2; z = x + i*y; 02481 // = ( sqrt[(x+1)^2 + y^2] + sqrt[(x-1)^2 + y^2] )/2; 02482 // Blomquist, 23.01.07; 02483 { 02484 l_interval res; 02485 // interval sqr_y = sqr( y ); 02486 // interval res = ( sqrt( sqr( x + 1.0 ) + sqr_y ) + 02487 // sqrt( sqr( x - 1.0 ) + sqr_y ) ) / 2.0; 02488 res = abs(x); 02489 if (y != 0.0 || Inf(res) < 1.0) 02490 { 02491 res = sqrtx2y2(x+1.0,y) + sqrtx2y2(x-1.0,y); // Blomquist; 02492 times2pown(res,-1); 02493 } 02494 02495 if ( Sup(res)==Infinity ) // Blomquist: program stop, if overflow occurs. 02496 cxscthrow (STD_FKT_OUT_OF_DEF( 02497 "l_cinterval asin( const l_cinterval& z); z out of range")); 02498 02499 // Now we correct a possible overestimation of the lower bound 02500 // of res. 02501 // It holds: 02502 // (sqrt[(x+1)^2 + y^2] + sqrt[(x-1)^2 + y^2])/2 >= Max(1,|x|); 02503 l_real hlb = max( l_real(1.0), abs( Sup( x ) ) ); 02504 if( Inf( res ) < hlb ) // this is an invalid overestimation! 02505 res = l_interval( hlb, Sup(res) ); 02506 return res; 02507 } 02508 02509 l_interval ACOSH_f_aux( const l_interval& x, const l_interval& y ) 02510 // Hilfsfunktion zur Berechnung des Imagin�rteils von asin(z), 02511 // z = x + i*y; 02512 // Bezeichnungen: 02513 // f_aux_asin(x,y) = alpha(x,y); 02514 // alpha(x,y) := ( sqrt[(x+1)^2 + y^2] + sqrt[(x-1)^2 + y^2] )/2; 02515 // R�ckgabewert: 02516 // 1. res = acosh( f_aux_asin(x,y) ), falls |x|>2 oder |y|>2; 02517 // 2. res = ln[ alpha + sqrt(alpha^2 - 1) ]; d = alpha(x,y) -1 02518 // = ln[1 + sqrt(d)*(sqrt(d) + sqrt(2+d))], falls gilt: 02519 // |x|<=2 und |y|<=2; x,y sind nur Punktintervalle! 02520 // Blomquist, 03.02.2007; 02521 02522 { 02523 const int c = 1017, 02524 c1 = 1022; 02525 l_interval res,xa(abs(x)),ya(abs(y)),diff, diff1,S1,S2; 02526 l_real rx((Inf(xa))), ry(Inf(ya)); 02527 int d,p,ex,ex1,ex2,s(0); 02528 02529 if (rx>2.0 || ry>2.0) res = acosh( f_aux_asin(x,y) ); 02530 else 02531 { // Jetzt mit Verbesserungen:! 02532 if (rx == 1.0) 02533 { 02534 S1 = 0.5 + (ya/4) / (sqrt1px2(ya/2) + 1); 02535 S2 = sqrt(ya) * sqrt(S1); // S2 = sqrt(delta); 02536 res = lnp1(S2 * (S2 + sqrt( 2+sqr(S2) ))); 02537 } 02538 else 02539 if (rx<1.0) // 0 <= x < +1; 02540 { 02541 diff = 1 - xa; 02542 d = expo_gr(ry); 02543 ex = expo_gr(diff); 02544 if (d > ex) ex = d; 02545 ex = -ex; 02546 // ex = -Max( expo(ry[1]),expo(diff[1]) ) > 0; 02547 if (ex>LI_maxexpo1) 02548 { 02549 times2pown(diff,LI_maxexpo1); 02550 times2pown(diff,ex-LI_maxexpo1); 02551 times2pown(ya,LI_maxexpo1); 02552 times2pown(ya,ex-LI_maxexpo1); 02553 } 02554 else 02555 { 02556 times2pown(diff,ex); 02557 times2pown(ya,ex); 02558 } 02559 S2 = sqrtx2y2(diff,ya) + diff; // S2 '=' 1; 02560 if (ex>=LI_maxexpo1) 02561 { 02562 d = ex - c; // c = 1017; 02563 if (d%2 != 0) d += 1; 02564 } else d = 0; 02565 S2 = comp(0.5,ex-d+1) / S2; 02566 diff = 1+xa; 02567 S1 = comp(0.5,-d+1) / (sqrtx2y2(diff,y) + diff); 02568 S1 = S1 + S2; 02569 times2pown(S1,-1); // S1 = { ... } / 2; 02570 S1 = sqrt(S1); 02571 d = d/2; 02572 if (d != 0) times2pown(S1,d); 02573 S1 = abs(y) * S1; // S1 = sqrt(alpha - 1) 02574 res = lnp1(S1 * (S1 + sqrt( 2+sqr(S1) ))); 02575 } 02576 else // 1 < x <= 2; 02577 if (y == 0) 02578 { 02579 S1 = sqrt(xa-1); 02580 res = lnp1( S1*(S1 + sqrt(xa+1)) ); 02581 } 02582 else // 1 < x <= 2 and 0 < |y| <= 2 02583 { 02584 diff = xa - 1; 02585 ex1 = expo_gr(diff); 02586 ex2 = expo_gr(ry); 02587 ex = ex1; 02588 if (ex2>ex) ex = ex2; 02589 ex = -ex; // ex = -Max{expo(x-1),expo(y)} 02590 p = -2*ex2 + 1 + ex1; 02591 if (ex>p) p = ex; // p = Max{-2*ex2 + 1 + ex1,ex} 02592 if (p>c1) 02593 { 02594 d = p - c; // c = 1017; 02595 if (d%2 != 0) d += 1; 02596 } else d = 0; 02597 S1 = xa + 1; 02598 if (d>c1) 02599 { 02600 S2 = l_interval(MinReal); 02601 s = d - c1; // s > 0; 02602 times2pown(S1,s); 02603 times2pown(ya,s); 02604 } 02605 else S2 = l_interval(comp(0.5,-d+1)); 02606 S1 = S2 / (sqrtx2y2(S1,ya) + S1); // S1 = u * 2^-d; 02607 02608 diff = xa - 1; diff1 = diff; 02609 S2 = l_interval(1); 02610 Times2pown(S2,ex-d); 02611 ya = abs(y); 02612 if (ex>LI_maxexpo1) 02613 { 02614 times2pown(diff,LI_maxexpo1); 02615 times2pown(diff,ex-LI_maxexpo1); 02616 times2pown(ya,LI_maxexpo1); 02617 times2pown(ya,ex-LI_maxexpo1); 02618 } else 02619 { 02620 times2pown(diff,ex); 02621 times2pown(ya,ex); 02622 } 02623 S2 = S2 / (sqrtx2y2(diff,ya) + diff); 02624 // S2 = v * 2^-d; 02625 02626 s = -2*ex2-d+1; 02627 Times2pown(diff1,s); 02628 ya = abs(y); 02629 Times2pown(ya,-ex2); 02630 diff1 = diff1 / (ya*ya); 02631 // diff1 = w * 2^-d; 02632 02633 res = diff1 + S2; 02634 res = sqrt((res + S1)/2); 02635 s = expo_gr(res); 02636 if (d/2+s>c1) 02637 { // Overflow vermeiden durch: 02638 ex = c1 - s; 02639 times2pown(res,ex); 02640 ya = abs(y); 02641 times2pown(ya,d/2-ex); 02642 S1 = ya * res; 02643 } else 02644 { 02645 times2pown(res,d/2); 02646 S1 = abs(y) * res; 02647 } 02648 res = lnp1(S1 * (S1 + sqrt( 2 + sqr(S1) ))); 02649 } 02650 } 02651 return res; 02652 } // ACOSH_f_aux 02653 02654 l_interval Asin_beta( const l_interval& x, const l_interval& y ) 02655 // Calculating the improved real part of asin(z); Blomquist 22.01.2007; 02656 // Re(asin(z)) = asin[ 2x/(sqrt[(x+1)^2+y^2] + sqrt[(x-1)^2+y^2]) ]=asin[beta]; 02657 // Improvements for beta --> +1 and beta --> -1 are necessary because of 02658 // nearly vertical tangents of the real asin(t)-function for |t|-->+1; 02659 // This function should only be used for POINT intervals x,y! z = x + i*y; 02660 { 02661 const real c1 = 0.75; 02662 bool neg_x; 02663 l_real Infxa; 02664 int ex_d,ex_y,ex; 02665 l_interval res,beta,abs_beta,delta,w_delta,xa,Ne,Sqrt_Ne,Kla; 02666 Ne = sqrtx2y2(1+x,y) + sqrtx2y2(1-x,y); 02667 beta = x / ( Ne/2 ); 02668 if (Inf(beta)<-1) Inf(beta) = -1; 02669 if (Sup(beta)> 1) Sup(beta) = 1; 02670 abs_beta = abs(beta); 02671 if (Inf(abs_beta) < c1) res = asin(beta); // Normal calculation 02672 else { // Inf(abs_beta)>=c1; Calculation now with improvements: 02673 Ne = Ne + 2.0; // Ne = 2 + sqrt[(1+x)^2 + y^2] + sqrt[(1-x)^2 + y^2]; 02674 Sqrt_Ne = sqrt(Ne); 02675 xa = x; 02676 neg_x = Inf(x)<0; 02677 if (neg_x) xa = -xa; // Inf(xa) >0 : 02678 Infxa = Inf(xa); 02679 02680 if (Infxa > 1.0) 02681 { // x > +1; 02682 if (y == 0.0) 02683 { 02684 w_delta = 0.0; 02685 delta = 0.0; 02686 } 02687 else 02688 { 02689 beta = xa - 1; // beta > 0; 02690 Infxa = Inf(y); 02691 ex_y = expo_gr(Infxa); 02692 ex_d = expo_gr(beta); 02693 ex = ex_y-ex_d-50; 02694 if (ex > 0) 02695 { 02696 times2pown(beta,ex); 02697 res = l_interval(comp(0.5,-ex+1)); 02698 delta = abs(y)/beta; 02699 Kla = sqrtx2y2(res,delta) + res; 02700 w_delta = sqrt(2*beta)*delta/(Sqrt_Ne * sqrt(Kla)); 02701 } 02702 else 02703 { 02704 delta = abs(y)/beta; 02705 Kla = sqrt1px2(delta) + 1; 02706 w_delta = sqrt(2*beta)*delta / (Sqrt_Ne * sqrt(Kla)); 02707 } 02708 delta = sqr(w_delta); 02709 } 02710 } 02711 else 02712 if (Infxa == 1.0) 02713 { // x = 1; 02714 delta = 2*abs(y) / Ne; 02715 w_delta = sqrt(2*abs(y)) / Sqrt_Ne; 02716 } 02717 else 02718 { // 0.75 <= x < 1; 02719 if (y == 0.0) 02720 { 02721 delta = 1 - xa; // delta is a point interval! 02722 w_delta = sqrt(delta); 02723 } 02724 else // 0.75 <= x < 1 and y != 0: 02725 { 02726 beta = 1 - xa; // beta > 0; 02727 Infxa = Inf(y); 02728 ex_y = expo_gr(Infxa); 02729 ex_d = expo_gr(beta); 02730 ex = ex_y-ex_d-50; 02731 if (ex > 0) 02732 { 02733 times2pown(beta,ex); 02734 res = l_interval(comp(0.5,-ex+1)); 02735 Kla = sqrtx2y2(res,y/beta) + res; 02736 } 02737 else // Now without scaling 02738 Kla = sqrt1px2(y/beta) + 1; 02739 times2pown(beta,1); // beta * 2; 02740 delta = beta * Kla / Ne; 02741 w_delta = sqrt(beta) * sqrt(Kla) / Sqrt_Ne; 02742 } 02743 } 02744 res = Pid2_l_interval() - asin( w_delta*sqrt(2-delta) ); 02745 if (neg_x) res = -res; 02746 } 02747 return res; 02748 } // Asin_beta(...) 02749 02750 l_cinterval asin( const l_cinterval& z ) throw() //------------- 040730 -- 02751 { 02752 const real gr = 6.355804e307; // upper bound for abs(rez),abs(imz) 02753 l_interval 02754 rez = Re(z), 02755 imz = Im(z); 02756 02757 l_real irez = Inf(rez), 02758 srez = Sup(rez), 02759 iimz = Inf(imz), 02760 simz = Sup(imz); 02761 02762 l_interval hxl(irez), hxu(srez), hyl(iimz), hyu(simz); 02763 02764 l_real resxl, resxu, resyl, resyu; 02765 02766 bool bl = iimz< 0.0 && simz>0.0, 02767 raxis = iimz==0.0 && simz==0; 02768 02769 // 02770 // 1st: check for singularities 02771 // 02772 if( (irez<-1 && (bl || (iimz<0 && simz==0))) || 02773 (srez >1 && (bl || (iimz==0 && simz>0))) ) 02774 cxscthrow(STD_FKT_OUT_OF_DEF("l_cinterval asin( const l_cinterval& z ); z contains singularities.")); 02775 // 02776 // check for too large bounds of abs(rez) and abs(imz) to prevent 02777 // overflow by calculating f_aux_asin(...) 02778 // 02779 resxl = max(abs(irez),abs(srez)); 02780 resxu = max(abs(iimz),abs(simz)); 02781 if (resxl>gr || resxu>gr) 02782 cxscthrow(STD_FKT_OUT_OF_DEF("l_cinterval asin( const l_cinterval& z ); z with too large bounds.")); 02783 // 02784 // 2nd: real part 02785 // 02786 if( iimz < 0.0 && simz > 0.0 ) 02787 // z intersects [-1,1] 02788 { 02789 if( irez <= 0.0 ) 02790 resxl = Inf( asin( hxl ) ); 02791 else 02792 // resxl = Inf( asin( hxl / f_aux_asin( hxl, interval( max( - iimz, simz ) ) ) ) ); 02793 resxl = Inf( Asin_beta(hxl,l_interval( max(-iimz,simz) )) ); // Blomquist, 19.06.2005; 02794 if( srez < 0.0 ) 02795 // resxu = Sup( asin( hxu / f_aux_asin( hxu, interval( max( - iimz, simz ) ) ) ) ); 02796 resxu = Sup( Asin_beta(hxu,l_interval( max(-iimz,simz) )) ); // Blomquist, 19.06.2005; 02797 else 02798 resxu = Sup( asin( hxu ) ); 02799 } 02800 else 02801 { 02802 if( ( iimz >= 0.0 && irez >= 0.0 ) || ( simz <= 0.0 && irez <= 0.0 ) ) 02803 // left boundary in quadrants I or III 02804 // min( Re( z ) ) in upper left corner 02805 // resxl = Inf( asin( hxl / f_aux_asin( hxl, hyu ) ) ); 02806 resxl = Inf( Asin_beta(hxl,hyu) ); // Blomquist, 19.06.2005; 02807 else 02808 // left boundary in quadrants II or IV 02809 // min( Re( z ) ) in lower left corner 02810 // resxl = Inf( asin( hxl / f_aux_asin( hxl, hyl ) ) ); 02811 resxl = Inf( Asin_beta(hxl,hyl) ); // Blomquist, 19.06.2005; 02812 if( ( iimz >= 0.0 && srez >= 0.0 ) || ( simz <= 0.0 && srez <= 0.0 ) ) 02813 // right boundary in quadrants I or III 02814 // max( Re( z ) ) in lower right corner 02815 // resxu = Sup( asin( hxu / f_aux_asin( hxu, hyl ) ) ); 02816 resxu = Sup( Asin_beta(hxu,hyl) ); // Blomquist, 19.06.2005; 02817 else 02818 // right boundary in quadrants II or IV 02819 // max( Re( z ) ) in upper right corner 02820 // resxu = Sup( asin( hxu / f_aux_asin( hxu, hyu ) ) ); 02821 resxu = Sup( Asin_beta(hxu,hyu) ); // Blomquist, 19.06.2005; 02822 } 02823 // 02824 // 3rd: imaginary part 02825 // 02826 02827 if (raxis) { // Interval argument is now a subset of the real axis. 02828 // Blomquist, 16.06.2005; 02829 if (srez<0.0) resyl = Inf( ACOSH_f_aux( hxu, hyu )); 02830 else resyl = -Sup( ACOSH_f_aux( hxu, hyu )); 02831 if (irez>0.0) resyu = -Inf( ACOSH_f_aux( hxl, hyu )); 02832 else resyu = Sup( ACOSH_f_aux( hxl, hyu )); 02833 } else 02834 if( simz <= 0.0 ) 02835 // z in lower half plane 02836 // min( Im( z ) ) in point with max |z| 02837 // max( Im( z ) ) in point with min |z| 02838 { 02839 if( irez < -srez ) 02840 // most of z in quadrant III 02841 { 02842 resyl = - Sup( ACOSH_f_aux( hxl, hyl ) ); 02843 if( srez < 0.0 ) 02844 resyu = - Inf( ACOSH_f_aux( hxu, hyu ) ); 02845 else 02846 resyu = - Inf( ACOSH_f_aux( l_interval(0), hyu ) ); 02847 } 02848 else 02849 // most of z in quadrant IV 02850 { 02851 resyl = - Sup( ACOSH_f_aux( hxu, hyl ) ); 02852 if( irez > 0.0 ) 02853 resyu = - Inf( ACOSH_f_aux( hxl, hyu ) ); 02854 else 02855 resyu = - Inf( ACOSH_f_aux( l_interval(0), hyu ) ); 02856 } 02857 } 02858 else if( iimz >= 0.0 ) 02859 // z in upper half plane 02860 // min( Im( z ) ) in point with min |z| 02861 // max( Im( z ) ) in point with max |z| 02862 { 02863 if( irez < -srez ) // if( irez + srez < 0.0 ) 02864 // most of z in quadrant II 02865 { 02866 resyu = Sup( ACOSH_f_aux( hxl, hyu ) ); 02867 if( srez < 0.0 ) 02868 resyl = Inf( ACOSH_f_aux( hxu, hyl ) ); 02869 else 02870 resyl = Inf( ACOSH_f_aux( l_interval(0), hyl ) ); 02871 } 02872 else 02873 // most of z in quadrant I 02874 { 02875 resyu = Sup( ACOSH_f_aux( hxu, hyu ) ); 02876 if( irez > 0.0 ) 02877 resyl = Inf( ACOSH_f_aux( hxl, hyl ) ); 02878 else 02879 resyl = Inf( ACOSH_f_aux( l_interval(0), hyl ) ); 02880 } 02881 } 02882 else 02883 // z intersects imaginary axes 02884 // min( Im( z ) ) in point in lower half plane with max |z| 02885 // max( Im( z ) ) in point in upper half plane with max |z| 02886 { 02887 if( irez < -srez ) // if( irez + srez < 0.0 ) 02888 // most of z in quadrants II and IV 02889 { 02890 resyl = - Sup( ACOSH_f_aux( hxl, hyl ) ); 02891 resyu = Sup( ACOSH_f_aux( hxl, hyu ) ); 02892 } 02893 else 02894 { 02895 resyl = - Sup( ACOSH_f_aux( hxu, hyl ) ); 02896 resyu = Sup( ACOSH_f_aux( hxu, hyu ) ); 02897 } 02898 } 02899 02900 return l_cinterval( l_interval( resxl,resxu ), l_interval( resyl,resyu ) ); 02901 02902 } 02903 // 02904 //-- end asin ----------------------------------------------------------------- 02905 02906 l_interval Asin_arg( const l_interval& x, const l_interval& y ) 02907 // Asin_arg berechnet f�r Punktintervalle x,y mit 02908 // beta := 2*|x| / ( sqrt((x+1)^2+y^2) + sqrt((x-1)^2+y^2) ) 02909 // und delta := 1 - |beta| eine Einschlie�ung von: 02910 // arcsin[ sqrt(delta)*sqrt(2-delta) ]. 02911 // Das Punktintervall x darf auch negativ sein! 02912 // Blomquist, 06.02.2007; 02913 { 02914 l_interval res,Ne,Sqrt_Ne,xa,beta,delta,w_delta,Kla; 02915 l_real Infxa; 02916 int ex,ex_y,ex_d; 02917 bool neg_x; 02918 02919 Ne = 2 + sqrtx2y2(1+x,y) + sqrtx2y2(1-x,y); 02920 Sqrt_Ne = sqrt(Ne); 02921 xa = x; 02922 neg_x = Inf(x)<0; 02923 if (neg_x) xa = -xa; // Inf(xa) >0 : 02924 Infxa = Inf(xa); 02925 02926 if (Infxa > 1.0) { // x > +1; 02927 if (y == 0.0) { 02928 w_delta = 0.0; 02929 delta = 0.0; 02930 } else 02931 { 02932 beta = xa - 1; // beta > 0; 02933 Infxa = Inf(y); 02934 ex_y = expo_gr(Infxa); 02935 ex_d = expo_gr(beta); 02936 ex = ex_y-ex_d-50; 02937 if (ex > 0) { 02938 times2pown(beta,ex); 02939 res = l_interval(comp(0.5,-ex+1)); 02940 delta = abs(y)/beta; 02941 Kla = sqrtx2y2(res,delta) + res; 02942 w_delta = sqrt(2*beta)*delta/(Sqrt_Ne * sqrt(Kla)); 02943 } else 02944 { 02945 delta = abs(y)/beta; 02946 Kla = sqrt1px2(delta) + 1; 02947 w_delta = sqrt(2*beta)*delta / (Sqrt_Ne * sqrt(Kla)); 02948 } 02949 delta = sqr(w_delta); 02950 } 02951 } else 02952 if (Infxa == 1.0) { // x = 1; 02953 delta = 2*abs(y) / Ne; 02954 w_delta = sqrt(2*abs(y)) / Sqrt_Ne; 02955 } else { // 0.75 <= x < 1; 02956 if (y == 0.0){ 02957 beta = 1 - xa; // beta is a point interval! 02958 delta = 4*beta / Ne; 02959 w_delta = 2*sqrt(beta) / Sqrt_Ne; 02960 } else // 0.75 <= x < 1 and y != 0: 02961 { 02962 beta = 1 - xa; // beta > 0; 02963 Infxa = Inf(y); 02964 ex_y = expo_gr(Infxa); 02965 ex_d = expo_gr(beta); 02966 ex = ex_y-ex_d-50; 02967 if (ex > 0) { 02968 times2pown(beta,ex); 02969 res = l_interval(comp(0.5,-ex+1)); 02970 Kla = sqrtx2y2(res,y/beta) + res; 02971 } else // Now without scaling 02972 Kla = sqrt1px2(y/beta) + 1; 02973 times2pown(beta,1); // beta * 2; 02974 delta = beta * Kla / Ne; 02975 w_delta = sqrt(beta) * sqrt(Kla) / Sqrt_Ne; 02976 } 02977 } 02978 res = asin( w_delta*sqrt(2-delta) ); 02979 02980 return res; 02981 } // Asin_arg 02982 02983 l_interval Acos_beta( const l_interval& x, const l_interval& y ) 02984 // Calculating the improved real part of acos(z); Blomquist 05.06.2005; 02985 // Re(acos(z)) = acos[ 2x / (sqrt[(x+1)^2+y^2] + sqrt[(x-1)^2+y^2]) ] 02986 { 02987 const real c1 = 0.75; 02988 l_interval res(0),beta; 02989 beta = x / ( (sqrtx2y2(1+x,y) + sqrtx2y2(1-x,y))/2 ); 02990 if (Inf(beta)<-1) Inf(beta)=-1; 02991 if (Sup(beta)> 1) Sup(beta)= 1; 02992 02993 if (Sup(beta)<c1) 02994 if (Sup(beta)<-c1) { // Improvement for beta --> -1 02995 res = Pi_l_interval() - Asin_arg(x,y); 02996 } else res = acos(beta); // Normal calculation 02997 else // Sup(beta)>=c1 02998 res = Asin_arg(x,y); 02999 return res; 03000 } 03001 03002 03003 //-- acos ----------------------------------------------------------- 040730 -- 03004 // 03005 // l_cinterval acos( const l_cinterval& z ) 03006 // { 03007 // w := acos(z); 03008 // Re(w) in a new Version, 03009 // Im(w) = -Im(asin(z)); Blomquist, 14.06.2005; 03010 // } 03011 // 03012 //-- acos ----------------------------------------------------------------- 03013 03014 03015 //-- acos --------------------------------------------------- 03016 // 03017 l_cinterval acos( const l_cinterval& z ) throw() 03018 { 03019 const real gr = 6.355804e307; // upper bound for abs(rez),abs(imz) 03020 l_interval 03021 rez = Re(z), 03022 imz = Im(z); 03023 03024 l_real 03025 irez = Inf(rez), 03026 srez = Sup(rez), 03027 iimz = Inf(imz), 03028 simz = Sup(imz); 03029 03030 l_interval 03031 hxl(irez), hxu(srez), hyl(iimz), hyu(simz); 03032 03033 bool bl = iimz< 0.0 && simz>0.0, 03034 raxis = iimz==0.0 && simz==0; 03035 l_real 03036 resxl, resxu, resyl, resyu; 03037 // 03038 // 1st: check for singularities 03039 // 03040 if( (irez<-1 && (bl || (iimz<0 && simz==0))) || 03041 (srez >1 && (bl || (iimz==0 && simz>0))) ) 03042 cxscthrow(STD_FKT_OUT_OF_DEF("l_cinterval acos( const l_cinterval& z ); z contains singularities.")); 03043 // 03044 // check for too large bounds of abs(rez) and abs(imz) to prevent 03045 // overflow by calculating f_aux_asin(...) 03046 // 03047 resxl = max(abs(irez),abs(srez)); 03048 resxu = max(abs(iimz),abs(simz)); 03049 if (resxl>gr || resxu>gr) 03050 cxscthrow(STD_FKT_OUT_OF_DEF("l_cinterval acos( const l_cinterval& z ); z with too large bounds.")); 03051 // 03052 // 2nd: real part 03053 // 03054 // Blomquist, 05.02.2007; 03055 if( iimz < 0.0 && simz > 0.0 ) 03056 // z intersects [-1,1] on the x-axis 03057 { 03058 if( irez <= 0.0 ) resxu = Sup( acos( hxl ) ); 03059 else resxu = Sup( Acos_beta(hxl,l_interval( max(-iimz,simz) )) ); 03060 03061 if( srez < 0.0 ) 03062 resxl = Inf( Acos_beta(hxu,l_interval(max(-iimz,simz))) ); 03063 else resxl = Inf( acos( hxu ) ); 03064 } 03065 else 03066 { 03067 if (irez<0 && srez>0) 03068 // z intersects the posizive or negative y-axis 03069 if (iimz >= 0) { 03070 resxl = Inf( Acos_beta(hxu,hyl) ); 03071 resxu = Sup( Acos_beta(hxl,hyl) ); 03072 } else { 03073 resxl = Inf( Acos_beta(hxu,hyu) ); 03074 resxu = Sup( Acos_beta(hxl,hyu) ); 03075 } 03076 else 03077 { 03078 if( ( iimz >= 0.0 && irez >= 0.0 ) || ( simz <= 0.0 && irez < 0.0 ) ) 03079 // left boundary in quadrants I or III 03080 // min( Re( z ) ) in lower right corner 03081 resxl = Inf( Acos_beta(hxu,hyl) ); 03082 else 03083 // left boundary in quadrants II or IV 03084 // min( Re( z ) ) in upper right corner 03085 resxl = Inf( Acos_beta(hxu,hyu) ); 03086 03087 if( ( iimz >= 0.0 && srez > 0.0 ) || ( simz <= 0.0 && srez <= 0.0 ) ) 03088 // right boundary in quadrants I or III 03089 // max( Re( z ) ) in upper left corner 03090 resxu = Sup( Acos_beta(hxl,hyu) ); 03091 else 03092 // right boundary in quadrants II or IV 03093 // max( Re( z ) ) in lower left corner 03094 resxu = Sup( Acos_beta(hxl,hyl) ); 03095 } 03096 } 03097 // 03098 // 3rd: imaginary part 03099 // 03100 if (raxis) { // Interval argument is now a subset of the real axis. 03101 // Blomquist, 16.06.2005; 03102 if (srez<0.0) resyl = Inf( ACOSH_f_aux( hxu, hyu )); 03103 else resyl = -Sup( ACOSH_f_aux( hxu, hyu )); 03104 if (irez>0.0) resyu = -Inf( ACOSH_f_aux( hxl, hyu )); 03105 else resyu = Sup( ACOSH_f_aux( hxl, hyu )); 03106 } else 03107 if( simz <= 0.0 ) 03108 // z in lower half plane 03109 // min( Im( z ) ) in point with max |z| 03110 // max( Im( z ) ) in point with min |z| 03111 { 03112 if( irez + srez < 0.0 ) 03113 // most of z in quadrant III 03114 { 03115 resyl = -Sup( ACOSH_f_aux( hxl, hyl ) ); 03116 if( srez < 0.0 ) 03117 resyu = -Inf( ACOSH_f_aux( hxu, hyu ) ); 03118 else 03119 resyu = -Inf( ACOSH_f_aux( l_interval(0), hyu ) ); 03120 } 03121 else 03122 // most of z in quadrant IV 03123 { 03124 resyl = -Sup( ACOSH_f_aux( hxu, hyl ) ); 03125 if( irez > 0.0 ) 03126 resyu = -Inf( ACOSH_f_aux( hxl, hyu ) ); 03127 else 03128 resyu = -Inf( ACOSH_f_aux( l_interval(0), hyu ) ); 03129 } 03130 } 03131 else if( iimz >= 0.0 ) 03132 // z in upper half plane 03133 // min( Im( z ) ) in point with min |z| 03134 // max( Im( z ) ) in point with max |z| 03135 { 03136 if( irez < -srez ) // if( irez + srez < 0.0 ) 03137 // most of z in quadrant II 03138 { 03139 resyu = Sup( ACOSH_f_aux( hxl, hyu ) ); 03140 if( srez < 0.0 ) 03141 resyl = Inf( ACOSH_f_aux( hxu, hyl ) ); 03142 else 03143 resyl = Inf( ACOSH_f_aux( l_interval(0), hyl ) ); 03144 } 03145 else 03146 // most of z in quadrant I 03147 { 03148 resyu = Sup( ACOSH_f_aux( hxu, hyu ) ); 03149 if( irez > 0.0 ) 03150 resyl = Inf( ACOSH_f_aux( hxl, hyl ) ); 03151 else 03152 resyl = Inf( ACOSH_f_aux( l_interval(0), hyl ) ); 03153 } 03154 } 03155 else 03156 // z intersects imaginary axes 03157 // min( Im( z ) ) in point in lower half plane with max |z| 03158 // max( Im( z ) ) in point in upper half plane with max |z| 03159 { 03160 if( irez < -srez ) // if( irez + srez < 0.0 ) 03161 // most of z in quadrants II and IV 03162 { 03163 resyl = -Sup( ACOSH_f_aux( hxl, hyl ) ); 03164 resyu = Sup( ACOSH_f_aux( hxl, hyu ) ); 03165 } 03166 else 03167 { 03168 resyl = -Sup( ACOSH_f_aux( hxu, hyl ) ); 03169 resyu = Sup( ACOSH_f_aux( hxu, hyu ) ); 03170 } 03171 } 03172 03173 return l_cinterval( l_interval( resxl, resxu ),-l_interval( resyl, resyu ) ); 03174 03175 } 03176 // 03177 //-- end acos ----------------------------------------------------------------- 03178 03179 03180 //-- asinh ---------------------------------------------------------- 040730 -- 03181 // 03182 l_cinterval asinh( const l_cinterval& z ) throw() 03183 // 03184 // asinh( z ) = i * asin( -i * z ) 03185 // 03186 { 03187 l_cinterval res = asin( l_cinterval( Im(z), -Re(z) ) ); 03188 return l_cinterval( -Im(res), Re(res) ); 03189 } 03190 // 03191 //-- end asinh ---------------------------------------------------------------- 03192 03193 03194 //-- acosh ---------------------------------------------------------- 040908 -- 03195 // 03196 l_cinterval acosh( const l_cinterval& z ) throw() 03197 // 03198 // acosh( z ) = i * acos( z ) = +/- i * ( pi / 2 - asin( z ) ) 03199 // 03200 { 03201 l_interval 03202 rez = Re(z), 03203 imz = Im(z); 03204 03205 l_real 03206 irez = Inf(rez), 03207 srez = Sup(rez), 03208 iimz = Inf(imz), 03209 simz = Sup(imz); 03210 03211 l_interval 03212 hxl(irez), hxu(srez), hyl(iimz), hyu(simz); 03213 03214 l_real 03215 resxl, resxu, resyl, resyu; 03216 03217 // cinterval res; 03218 03219 // 03220 // 1st: check for singularities 03221 // 03222 if( ( iimz <= 0.0 && simz >= 0.0 ) && ( irez < 1.0 ) ) 03223 cxscthrow(STD_FKT_OUT_OF_DEF("l_cinterval acosh( const l_cinterval& z ); z contains singularities.")); 03224 // With this restriction the complex interval argument and the real axis must not have any common 03225 // point, if irez < +1; 03226 // So for example the negative real axis must not be touched from above if irez<1, although this 03227 // should be possible if the principal branch is considered! So the above restriction is too widely in 03228 // some cases; Blomquist, 21.06.2005; 03229 // 03230 // 2nd: z in upper half plane (or on the real axis) 03231 // acosh( z ) = + i * ( pi / 2 - asin( z ) ) 03232 // 03233 if( iimz > 0.0 ) 03234 { 03235 l_cinterval res = acos(z); 03236 return l_cinterval( -Im(res),Re(res) ); 03237 } 03238 // 03239 // 3rd: z in lower half plane 03240 // acosh( z ) = - i * ( pi / 2 - asin( z ) ) 03241 // 03242 if( simz < 0.0 ) 03243 { 03244 // cinterval res = HALFPI() - asin( z ); 03245 l_cinterval res = acos(z); // Blomquist, 14.06.2005 03246 return l_cinterval( Im(res), -Re(res) ); 03247 } 03248 // 03249 // z intersects [1,infinity) 03250 // 03251 // real part 03252 // minimum on the left on real axes, maximum in lower or upper right corner 03253 // 03254 resxl = Inf( acosh( hxl ) ); 03255 l_interval ytilde( max( -iimz, simz ) ); 03256 // resxu = Sup( acosh( f_aux_asin( hxu, ytilde ) ) ); 03257 resxu = Sup( ACOSH_f_aux(hxu,ytilde) ); // Blomquist, 14.06.2005; 03258 // 03259 // imaginary part 03260 // minimum in lower left corner, maximum in upper left corner 03261 // 03262 // resyl = -Sup( acos( hxl / f_aux_asin( hxl, hyl ) ) ); 03263 resyl = -Sup( Acos_beta(hxl,hyl) ); // Blomquist, 14.06.2005; 03264 // resyu = Sup( acos( hxl / f_aux_asin( hxl, hyu ) ) ); 03265 resyu = Sup( Acos_beta(hxl,hyu) ); 03266 03267 return l_cinterval(l_interval( resxl, resxu ),l_interval( resyl, resyu )); 03268 } 03269 // 03270 //-- end acosh ---------------------------------------------------------------- 03271 03272 03273 03274 //----------------------------------------------------------------------------- 03275 // 03276 // Section 6: atan, acot, atanh, acoth 03277 // 03278 // The implementation of acot, atanh, and acoth is based on atan: 03279 // 03280 // acot( z ) = atan( 1 / z ) 03281 // atanh( z ) = - i * atan( i * z ) 03282 // acoth( z ) = i * acot( i * z ) 03283 // 03284 // Only principal values are computed. 03285 // 03286 //----------------------------------------------------------------------------- 03287 03288 //-- atan ----------------------------------------------------------- 040912 -- 03289 // 03290 // Analytic inverse tangent function 03291 // 03292 03293 void re_atan( const l_interval& y, l_interval& x, l_interval& res ) 03294 // Calculating an inclusion res of 1 - y^2 - x^2; 03295 // If |1 - y^2 - x^2| is too small, then x = (2^d * x) and 03296 // res = 2^d * (1 - y^2 - x^2); 03297 // x is always a point interval. 03298 // Blomquist, 26.02.2007; 03299 { 03300 if ( x==0.0 ) res = 1.0; 03301 else // x <> 0: 03302 { 03303 interval y1,x1; 03304 int ex,ex_x,d; 03305 l_interval y_(y); 03306 l_real lr; 03307 y1 = y; x1 = x; 03308 y1 = 1 - sqr(y1) - sqr(x1); 03309 if ( expo(Sup(y1))<-2 ) 03310 { 03311 res = 1 - sqr(y) - sqr(x); 03312 lr = Sup(abs(res)); 03313 ex = expo_gr(lr); 03314 if (ex < -10000) res = 0; 03315 else 03316 { 03317 lr = Sup(x); 03318 ex_x = expo_gr(lr); 03319 d = 1 - ex; 03320 ex_x = 1022 - ex_x; 03321 if (ex_x<d) d = ex_x; // d: Minimum{1-ex,1022-ex_x} 03322 if (d>1022) d = 1022; // d: Minimum{1-ex,1022-ex_x,1022} 03323 if (d%2 != 0) d--; // d divisible by 2; 03324 ex = d/2; 03325 times2pown(x,ex); 03326 times2pown(y_,ex); 03327 res = comp(0.5,d+1) - sqr(x) - sqr(y_); 03328 times2pown(x,ex); 03329 } 03330 } 03331 else res = 1 - sqr(y) - sqr(x); 03332 } 03333 } // re_atan 03334 03335 void re_vert( const l_real& x, const l_interval& hx, 03336 const l_real& rew_inf, const l_real& rew_sup, 03337 l_real& resxl, l_real& resxu ) //---------------------- 040729 -- 03338 // 03339 // Subroutine of analytic inverse tangent function. 03340 // Evaluate real part on a vertical boundary. 03341 // 03342 { 03343 if( x == 0.0 ) 03344 // singularities have been handled before, hence Re( w ) > 0 03345 { 03346 resxl = 0.0; 03347 resxu = 0.0; 03348 } 03349 else 03350 { 03351 l_interval hx2(hx); 03352 times2pown(hx2,1); 03353 if( x > 0.0 ) 03354 // w in quadrants I and/or II 03355 // atan is the inverse function of tan(t), t in (-pi/2,pi/2). 03356 { 03357 resxl = rew_sup > 0.0 ? Inf( Atan( hx2,rew_sup )/2.0 ) 03358 : ( rew_sup < 0.0 ? Inf( (Atan( hx2,rew_sup ) + Pi_l_interval() )/2.0 ) 03359 : Inf( Pid4_l_interval() ) ); 03360 03361 resxu = rew_inf > 0.0 ? Sup( Atan( hx2,rew_inf )/2.0 ) 03362 : ( rew_inf < 0.0 ? Sup( (Atan( hx2,rew_inf ) + Pi_l_interval())/2.0 ) 03363 : Sup( Pid4_l_interval() ) ); 03364 } 03365 else 03366 // w in quadrants III and/or IV 03367 { 03368 resxl = rew_inf < 0.0 ? Inf( (Atan( hx2,rew_inf ) - Pi_l_interval())/2.0 ) 03369 : ( rew_inf > 0.0 ? Inf( Atan( hx2,rew_inf )/2.0 ) 03370 : -Sup( Pid2_l_interval()/2.0 ) ); 03371 resxu = rew_sup < 0.0 ? Sup( (Atan( hx2,rew_sup ) - Pi_l_interval())/2.0 ) 03372 : ( rew_sup > 0.0 ? Sup( Atan( hx2,rew_sup )/2.0 ) 03373 : -Inf( Pid4_l_interval() ) ); 03374 } 03375 } 03376 } // re_vert 03377 03378 l_interval Aux_1_atan(const l_real& x) 03379 // x > minreal; minreal = 2^(-1074); 03380 // Calculating: ln[ 1+2/(sqrt(1+x^2)-1) ], [x] = x, 03381 // [x] is a point interval ! 03382 // Tested in detail, Blomquist; 15.02.07; 03383 { 03384 const int exOv = +510; 03385 const int exUn = 1; 03386 03387 l_interval res, 03388 ix(x), // ix is point interval with x>0; 03389 r,t; 03390 int ex(expo_gr(x)); 03391 03392 if (ex>=exOv) 03393 { // preventing overflow 03394 r = 1/ix; 03395 res = lnp1( 2*r*( sqrt1px2(r) + r) ); 03396 } else 03397 if (ex<=exUn) 03398 { // x < 2^(1) 03399 res = sqr(ix); 03400 // Calculating r := 2*ln(ix) for x --> 0; 03401 if (ex<-1000) { 03402 Times2pown(ix,-ex); 03403 r = ln(ix) + ex*Ln2_l_interval(); 03404 } else r = ln(ix); 03405 times2pown(r,1); 03406 t = 1+sqrt(1+res); 03407 times2pown(t,1); // t = 2*{1+sqrt(1+x^2)} 03408 res = ln(t + res) - r; 03409 } else 03410 { // normal calculation 03411 res = sqrtp1m1( sqr(ix) ); // res = sqrt(1+x^2)-1 03412 stagprec++; 03413 res = 2/res; 03414 stagprec--; 03415 res = lnp1(res); // res = ln[1 + 2/(sqrt(1+x^2)-1) ] 03416 } 03417 03418 return res; 03419 } // Aux_1_atan 03420 03421 03422 l_interval Q_atan_UPSIGN(const l_interval& x, const l_interval& y) 03423 { 03424 // x: abs(Re(z)); x: a real interval x = [x1,x2], with 0<=x1<x2. 03425 // y: Inf(Im(z)); y is point interval. 03426 // Q_atan_UPSIGN returns an inclusion of ln[1 + 4y/(x^2+(1-y)^2)] 03427 // Tested in detail; Blomquist, 14.02.2007; 03428 const int n = 511; 03429 l_interval res(0.0),t,t1,t2; 03430 l_real lr; 03431 int ex_x,ex_y,ex,s; 03432 if (y==1.0) 03433 if(Sup(x)<1) 03434 { // Now: y = 1 and minreal <= x1 <= x2 < +1 : 03435 lr = Inf(x); 03436 // Calculation of t = ln(x) with scaling: 03437 // ex = expo(x); 03438 // ln(x) = ln(x*2^-ex) + ex * ln(2), 03439 // if Inf(x) is too small, i.e. expo(lr[1])<-600 03440 ex = expo_gr(lr); // ex < 0; 03441 if (ex<-600) 03442 { 03443 t1 = x; 03444 Times2pown(t1,-ex); // t1 = 2^(-ex) * x; 03445 t = ln(t1) + ex*Ln2_l_interval(); // t = ln(x) 03446 } else t = ln(x); // t: ln(x) 03447 times2pown(t,1); // t: 2*ln(x); 03448 res = ln(4+sqr(x)) - t; // Blomquist, 14.02.2007. 03449 } 03450 else res = lnp1( sqr(2/x) ); 03451 else 03452 { // Now we have: y>=0 and y<>1 03453 lr = Sup(x); 03454 ex_x = expo_gr(lr); 03455 lr = Sup(y); 03456 ex = expo_gr(lr); 03457 if (ex<-100000) return res; // y = 0 ---> res = 0.0; 03458 // Now we have: ( y>0 and y<>1 ) 03459 ex_y = ex; 03460 if (ex_x>ex) ex = ex_x; 03461 // ex = Maximum{expo(x),expo(y)} 03462 if (ex>n) 03463 { // Now scaling to prevent overflow by calculating 03464 // the denominator: x^2 + (1-y)^2 : 03465 s = n-ex-1; 03466 // s is chosen in such a way that x^2 + (1-y)^2 will be maximal 03467 // after scaling, but an overflow must not be realized! 03468 t = x; times2pown(t,s); // fast scaling with 2^(s) 03469 t1 = y; times2pown(t1,s); // fast scaling with 2^(s) 03470 t2 = sqr(t) + sqr(comp(0.5,s+1)-t1); // t2: denominator 03471 stagprec++; 03472 t = y / t2; // scaled quotient 03473 stagprec--; 03474 times2pown(t,2*s+2); // back-scaling with 2^(s+2); 03475 // '+2': factor 4 !! 03476 if (Inf(t)<0) SetInf(t,0.0); 03477 res = lnp1(t); 03478 } 03479 else 03480 { // Calculating x^2 + (1-y)^2 now without overflow!!! 03481 // Now we have: ( y>0 and y<>1 ) 03482 t = 1 - y; 03483 lr = Inf(t); 03484 ex = expo_gr(lr); 03485 if (ex_x>ex) ex = ex_x; // It holds: -1073 <= ex; 03486 ex--; 03487 // Now x^2 + (1-y)^2 is scaled, so that: 03488 // [x*2^(-ex)]^2 + [(1-y)*2^(-ex)]^2 '=' 1 03489 Times2pown(t,-ex); // t = (1-y)*2^(-ex) 03490 t1 = x; 03491 Times2pown(t1,-ex); // t = x * 2^(-ex) 03492 res = sqr(t1) + sqr(t); // res: scaled denominator: 03493 // res = (x*2^-ex)^2 + [(1-y)*2^-ex]^2; 03494 s = ex_y+2-2*ex; 03495 // For x --> 0 and y --> 1 it holds: 03496 // 2^s is the order of 2^(-2ex) * [x^2+(1+y)^2]. 03497 if (s < 1020) 03498 { // No overflow by calculating y * 2^(-2*ex+2) 03499 // is not possible: 03500 t = y; 03501 times2pown(t,-2*ex+2); // t = y * 2^(-2*ex+2) 03502 stagprec++; 03503 t = t / res; 03504 stagprec--; 03505 res = lnp1(t); 03506 } 03507 else // s >= 1020; 03508 { // For x --> 0 and y --> 1 it holds: 03509 // 2^s is the order of 2^(-2ex) * [x^2+(1+y)^2]. 03510 // Thus, for s >= 1020 an overflow by calculating 03511 // 2^(-2ex) * [x^2+(1+y)^2] is possible and an 03512 // appropriate scaling is necessary: 03513 ex_y = s-1020; 03514 if (ex_y%2 != 0) ex_y++; // ex_y = d, now divisible by 2. 03515 s = ex_y/2; // s = d/2; 03516 t1 = x; 03517 times2pown(t1,-ex-s); 03518 t2 = 1+y; 03519 times2pown(t2,-ex-s); 03520 t = sqr(t1) + sqr(t2); 03521 stagprec++; 03522 t = t / res; 03523 stagprec--; 03524 res = ex_y * Ln2_l_interval() + ln(t); 03525 } 03526 } 03527 } 03528 03529 return res; 03530 } // Q_atan_UPSIGN 03531 03532 03533 l_cinterval atan( const l_cinterval& z ) throw() //----- 040912 -- 03534 { 03535 l_interval 03536 rez = Re(z), 03537 imz = Im(z); 03538 03539 l_real 03540 irez = Inf(rez), 03541 srez = Sup(rez), 03542 iimz = Inf(imz), 03543 simz = Sup(imz); 03544 03545 const int n = 511; // For possible scaling 03546 03547 l_interval 03548 hxl(irez), hxu(srez), hyl(iimz), hyu(simz); // all theses intervals are point intervals! 03549 03550 l_real 03551 resxl, resxu, resyl, resyu; 03552 // 03553 // 1st: check for singularities 03554 // 03555 if( ( irez <= 0.0 && srez >= 0.0 ) && ( iimz <= -1.0 || simz >= 1.0 ) ) 03556 cxscthrow(STD_FKT_OUT_OF_DEF("l_cinterval atan( const l_cinterval& z ); z contains singularities.")); 03557 // 03558 // 2nd: real part 03559 // Re( atan( z ) ) = Arg( w ) / 2, where w = 1 - x^2 - y^2 + i * 2x ) 03560 // 03561 // evaluate atan on vertical boundaries 03562 // 03563 l_interval 03564 // y_sqr = sqr( imz ), 03565 // rew_l = (1 - y_sqr) - sqr( hxl ), // Blomquist; before: rew_l = 1 - sqr(hxl) - y_sqr, 03566 // rew_u = (1 - y_sqr) - sqr( hxu ); // Blomquist; before: rew_u = 1 - sqr(hxu) - y_sqr; 03567 rew_l, rew_u; 03568 03569 // ------------------------------ Blomquist --------------------------------- 03570 // ---------- Improvements for Im(z) = [1,1] or Im(z) = [-1,-1] 03571 bool sqrImz_1 = (iimz==simz) && (iimz==1.0 || iimz==-1.0); 03572 // Test for Im(z) = [1,1] or [-1,-1] 03573 03574 if (sqrImz_1) { 03575 rew_l = -abs(hxl); hxl = l_interval(sign(irez)); 03576 rew_u = -abs(hxu); hxu = l_interval(sign(srez)); 03577 } 03578 else { 03579 int ex,s; 03580 l_interval imz_, scf; 03581 int ex1 = expo_gr(iimz); 03582 int ex2 = expo_gr(simz); 03583 if (ex2>ex1) ex1 = ex2; 03584 03585 ex = expo_gr(irez); 03586 if(ex1>ex) ex = ex1; // ex = Maximum 03587 if (ex>n) { // Scaling necessary 03588 s = n - ex - 1; 03589 scf = l_interval(comp(0.5,s+1)); // scf: scaling factor 2^s 03590 times2pown(scf,s); // scf = 2^(2*s) 03591 times2pown(hxl,s); // hxl = hxl * 2^s 03592 imz_ = imz; 03593 times2pown(imz_,s); // imz_ = imz_ * 2^s 03594 rew_l = (scf - sqr(imz_)) - sqr(hxl); // here now without overflow!! 03595 times2pown(hxl,s); // hxl = hxl * 2^s 03596 } else // rew_l = (1 - sqr( imz )) - sqr( hxl ); 03597 re_atan(imz,hxl,rew_l); 03598 03599 ex = expo_gr(srez); // 03600 if(ex1>ex) ex = ex1; // Maximum 03601 if (ex>n) { // Scaling necessary 03602 s = n - ex - 1; 03603 scf = l_interval(comp(0.5,s+1)); // scf: scaling factor 2^s 03604 times2pown(scf,s); // scf = 2^(2*s) 03605 times2pown(hxu,s); // hxu = hxu * 2^s 03606 imz_ = imz; 03607 times2pown(imz_,s); // imz_ = imz_ * 2^s 03608 rew_u = (scf - sqr(imz_)) - sqr(hxu); // here now without overflow!! 03609 times2pown(hxu,s); // hxu = hxu * 2^s 03610 } else // rew_u = (1 - sqr( imz )) - sqr( hxu ); 03611 re_atan(imz,hxu,rew_u); 03612 } 03613 // ------------------------------ Blomquist; 22.02.05; -------------------- 03614 03615 // 03616 // left boundary 03617 // 03618 l_real rew_inf = Inf( rew_l ); 03619 l_real rew_sup = Sup( rew_l ); 03620 re_vert( irez, hxl, rew_inf, rew_sup, resxl, resxu ); 03621 03622 // 03623 // right boundary 03624 // 03625 rew_inf = Inf( rew_u ); 03626 rew_sup = Sup( rew_u ); 03627 l_real res_l, res_u; 03628 re_vert( srez, hxu, rew_inf, rew_sup, res_l, res_u ); 03629 03630 if (res_l<resxl) resxl = res_l; 03631 if (res_u>resxu) resxu = res_u; 03632 03633 // 03634 // look for extremal values on horizontal boundaries 03635 // since atan( x+iy ) = atan( x-iy ), 03636 // intersections can be considered in the upper half plane 03637 // 03638 l_real abs_y_min = Inf( abs( imz ) ); 03639 03640 if( abs_y_min > 1.0 ) 03641 { 03642 l_interval 03643 abs_hyl = l_interval( abs_y_min ), 03644 // abs_hxl = sqrt( sqr( abs_hyl ) - 1.0 ); 03645 abs_hxl = sqrtx2m1(abs_hyl); // Blomquist; 03646 03647 if( Sup( abs_hxl ) > irez && Inf( abs_hxl ) < srez ) 03648 // extremal curve intersects lower boundary of x+i|y| in quadrant I 03649 // intersection in Q I or Q IV: update minimum 03650 // resxl = inf( atan( abs_y_min / abs_hxl ) ) / 2.0; 03651 resxl = Inf( (Pi_l_interval() - atan( 1.0 / abs_hxl ))/2.0 ); 03652 else if( -Inf( abs_hxl ) > irez && -Sup( abs_hxl ) < srez ) 03653 // extremal curve intersects lower boundary of x+i|y| in quadrant II 03654 // intersection in Q II or Q III: update maximum 03655 resxu = Sup( (atan( 1.0 / abs_hxl ) - Pi_l_interval())/2.0 ); 03656 } 03657 03658 // 3rd: imaginary part 03659 // Im( atan( z ) ) = +/- Ln( 1 +/- 4y/( x^2 + (1 -/+ y)^2 ) ) / 4 03660 // 03661 // evaluate atan on horizontal boundaries 03662 l_interval 03663 abs_rez = abs(rez), 03664 im_atan_l, im_atan_u; 03665 03666 if( iimz < 0.0 ) 03667 // im_atan_l = -ln( 1 - 4 * hyl / ( x_sqr + sqr( 1 + hyl ) ) ) / 4.0; 03668 // im_atan_l = -lnp1(-4 * hyl / ( x_sqr + sqr( 1 + hyl ) )) / 4.0; // Blomquist 03669 im_atan_l = -Q_atan_UPSIGN(abs_rez,-hyl); // Blomquist 03670 else 03671 // im_atan_l = ln( 1 + 4 * hyl / ( x_sqr + sqr( 1 - hyl ) ) ) / 4.0; 03672 im_atan_l = Q_atan_UPSIGN(abs_rez,hyl); // Blomquist 03673 times2pown(im_atan_l,-2); // Division by 4 03674 03675 if( simz < 0.0 ) 03676 // im_atan_u = -ln( 1 - 4 * hyu / ( x_sqr + sqr( 1 + hyu ) ) ) / 4.0; 03677 // im_atan_u = -lnp1(-4 * hyu / ( x_sqr + sqr( 1 + hyu ) ) ) / 4.0; // Blomquist 03678 im_atan_u = -Q_atan_UPSIGN(abs_rez,-hyu); // Blomquist 03679 else 03680 // im_atan_u = ln( 1 + 4 * hyu / ( x_sqr + sqr( 1 - hyu ) ) ) / 4.0; 03681 im_atan_u = Q_atan_UPSIGN(abs_rez,hyu); // Blomquist 03682 times2pown(im_atan_u,-2); // Division by 4 03683 03684 resyl = min( Inf( im_atan_l ), Inf( im_atan_u ) ); 03685 resyu = max( Sup( im_atan_l ), Sup( im_atan_u ) ); 03686 // 03687 // look for extremal values on vertical boundaries, 03688 // if vertical boundaries intersect extremal curves 03689 // 03690 l_real abs_x_min = Inf( abs( rez ) ); 03691 l_interval 03692 x_extr = l_interval( abs_x_min ), 03693 // y_extr = sqrt( 1.0 + sqr( x_extr ) ); 03694 y_extr = sqrt1px2(x_extr); // Blomquist; 03695 03696 if( Inf( y_extr ) < simz && Sup( y_extr ) > iimz ) 03697 // extremal curve intersects left boundary of |x|+iy in quadrant I 03698 // update maximum 03699 // resyu = Sup( ln( 1 + 4 * y_extr / ( sqr( x_extr ) + sqr( 1 - y_extr ) ) ) ) / 4.0; 03700 // resyu = Sup( Aux_1_atan(abs_x_min)/4.0 ); // Blomquist 03701 { 03702 rez = Aux_1_atan(abs_x_min); 03703 times2pown(rez,-2); 03704 resyu = Sup(rez); 03705 } 03706 03707 if( -Sup( y_extr ) < simz && -Inf( y_extr ) > iimz ) 03708 // extremal curve intersects left boundary of |x|+iy in quadrant IV 03709 // update minimum 03710 // resyl = -Sup( ln( 1 + 4 * y_extr / ( sqr( x_extr ) + sqr( 1 - y_extr ) ) ) ) / 4.0; 03711 // resyl = -Sup( Aux_1_atan(abs_x_min)/4.0 ); // Blomquist 03712 { 03713 rez = Aux_1_atan(abs_x_min); 03714 times2pown(rez,-2); 03715 resyl = -Sup(rez); 03716 } 03717 03718 return l_cinterval( l_interval( resxl, resxu ), l_interval( resyl, resyu ) ); 03719 03720 } 03721 // 03722 //-- end atan ----------------------------------------------------------------- 03723 03724 03725 //-- acot ----------------------------------------------------------- 040912 -- 03726 // 03727 // Analytic inverse cotangent function 03728 // acot( z ) = atan( 1/z ) 03729 // The code of acot( z ) is almost identical to the code of atan( z ) 03730 // 03731 l_cinterval acot( const l_cinterval& z ) throw() 03732 { 03733 l_interval 03734 rez = Re(z), 03735 imz = Im(z); 03736 03737 l_real 03738 irez = Inf(rez), 03739 srez = Sup(rez), 03740 iimz = Inf(imz), 03741 simz = Sup(imz); 03742 03743 const int n = 511; // For possible scaling 03744 03745 l_interval 03746 hxl(irez), hxu(srez), hyl(iimz), hyu(simz); 03747 03748 l_real 03749 resxl, resxu, resyl, resyu; 03750 // 03751 // 1st: check for singularities 03752 // 03753 if( ( irez <= 0.0 && srez >= 0.0 ) && ( iimz <= 1.0 && simz >= -1.0 ) ) 03754 cxscthrow(STD_FKT_OUT_OF_DEF("l_cinterval acot( const l_cinterval& z ); z contains singularities.")); 03755 // 03756 // 2nd: real part 03757 // Re( atan( z ) ) = Arg( w ) / 2, where w = 1 - x^2 - y^2 + i * 2x ) 03758 // Re( atan( 1 / z ) ) = Arg( w ) / 2, where w = x^2 + y^2 - 1 + i * 2x ) 03759 // 03760 // evaluate acot on vertical boundaries 03761 // 03762 l_interval 03763 // y_sqr = sqr( imz ), 03764 // rew_l = (y_sqr - 1) + sqr(hxl), 03765 // rew_u = (y_sqr - 1) + sqr(hxu); 03766 // rew_l = (sqr( hxl )-1) + y_sqr, 03767 // rew_u = (sqr( hxu )-1) + y_sqr; 03768 rew_l, rew_u; 03769 // ------------------------------ Blomquist ------------------------------ 03770 // ---------- Improvements for Im(z) = [1,1] or Im(z) = [-1,-1] 03771 bool sqrImz_1 = (iimz==simz) && (iimz==1.0 || iimz==-1.0); 03772 // Test for Im(z) = [1,1] or [-1,-1] 03773 03774 if (sqrImz_1) { 03775 rew_l = abs(hxl); hxl = l_interval(sign(irez)); 03776 rew_u = abs(hxu); hxu = l_interval(sign(srez)); 03777 } 03778 else { 03779 int ex,s; 03780 // l_real scf; // Scaling factor 03781 l_interval imz_, scf; 03782 int ex1 = expo_gr(iimz); int ex2 = expo_gr(simz); 03783 if (ex2>ex1) ex1 = ex2; 03784 03785 ex = expo_gr(irez); 03786 if(ex1>ex) ex = ex1; // Maximum 03787 if (ex>n) { // Scaling necessary 03788 s = n - ex - 1; 03789 scf = l_interval(comp(0.5,s+1)); // scf: scaling factor 2^s 03790 times2pown(scf,s); // scf = 2^(2*s); 03791 times2pown(hxl,s); // hxl = hxl * 2^s 03792 imz_ = imz; 03793 times2pown(imz_,s); // imz_ = imz_ * 2^s 03794 rew_l = (sqr(imz_) - scf) + sqr(hxl); // here now without overflow!! 03795 times2pown(hxl,s); // hxl = hxl * 2^s 03796 } else // rew_l = (sqr( imz ) - 1.0) + sqr( hxl ); 03797 { 03798 re_atan(imz,hxl,rew_l); 03799 rew_l = -rew_l; 03800 } 03801 03802 ex = expo_gr(srez); 03803 if(ex1>ex) ex = ex1; // Maximum 03804 if (ex>n) { // Scaling necessary 03805 s = n - ex - 1; 03806 scf = l_interval(comp(0.5,s+1)); // scf: scaling factor 2^s 03807 times2pown(scf,s); // scf = 2^(2*s); 03808 times2pown(hxu,s); // hxu = hxu * 2^s 03809 imz_ = imz; 03810 times2pown(imz_,s); // imz_ = imz_ * 2^s 03811 rew_u = (sqr(imz_) - scf) + sqr(hxu); // here now without overflow!! 03812 times2pown(hxu,s); // hxu = hxu * 2^s 03813 } else // rew_u = (sqr( imz )-1.0) + sqr( hxu ); 03814 { 03815 re_atan(imz,hxu,rew_u); 03816 rew_u = -rew_u; 03817 } 03818 } 03819 // ------------------------------ Blomquist; 22.02.05; ------------------ 03820 03821 // 03822 // left boundary 03823 // 03824 l_real rew_inf = Inf( rew_l ); 03825 l_real rew_sup = Sup( rew_l ); 03826 re_vert( irez, hxl, rew_inf, rew_sup, resxl, resxu ); 03827 // 03828 // right boundary 03829 // 03830 rew_inf = Inf( rew_u ); 03831 rew_sup = Sup( rew_u ); 03832 l_real res_l, res_u; 03833 re_vert( srez, hxu, rew_inf, rew_sup, res_l, res_u ); 03834 03835 if (res_l<resxl) resxl = res_l; 03836 if (res_u>resxu) resxu = res_u; 03837 03838 // 03839 // look for extremal values on horizontal boundaries 03840 // since acot( x+iy ) = acot( x-iy ), 03841 // intersections can be considered in the upper half plane 03842 // 03843 l_real abs_y_min = Inf( abs( imz ) ); 03844 03845 if( abs_y_min > 1.0 ) 03846 { 03847 l_interval 03848 abs_hyl = l_interval( abs_y_min ), 03849 // abs_hxl = sqrt( sqr( abs_hyl ) - 1.0 ); 03850 abs_hxl = sqrtx2m1(abs_hyl); // Blomquist; 03851 if( Sup( abs_hxl ) > irez && Inf( abs_hxl ) < srez ) 03852 // extremal curve intersects lower boundary of x+i|y| in quadrant I 03853 // intersection in Q I or Q IV: update maximum 03854 resxu = Sup( atan( 1.0 / abs_hxl ) / 2.0 ); 03855 if( -Inf( abs_hxl ) > irez && -Sup( abs_hxl ) < srez ) 03856 // extremal curve intersects lower boundary of x+i|y| in quadrant II 03857 // intersection in Q II or Q III: update minimum 03858 resxl = -Sup( atan( 1.0 / abs_hxl ) / 2.0 ); 03859 } 03860 // 03861 // 3rd: imaginary part 03862 // Im( atan( z ) ) = +/- Ln( 1 +/- 4y/( x^2 + (1 -/+ y)^2 ) ) / 4 03863 // Im( acot ) = - Im ( atan ): We calculate Im( atan ) and return "-" 03864 // 03865 // evaluate atan on horizontal boundaries 03866 // 03867 l_interval 03868 // x_sqr = sqr( rez ), // overflow is avoided by calling Q_atan_UPSIGN(...) 03869 im_atan_l, im_atan_u, 03870 abs_rez = abs(rez); // Blomquist; 03871 if( iimz < 0.0 ) 03872 // im_atan_l = -ln( 1 - 4 * hyl / ( x_sqr + sqr( 1 + hyl ) ) ) / 4.0; 03873 im_atan_l = -Q_atan_UPSIGN(abs_rez,-hyl); // Blomquist 03874 else 03875 // im_atan_l = ln( 1 + 4 * hyl / ( x_sqr + sqr( 1 - hyl ) ) ) / 4.0; 03876 im_atan_l = Q_atan_UPSIGN(abs_rez,hyl); // Blomquist 03877 times2pown(im_atan_l,-2); // Division by 4; 03878 03879 if( simz < 0.0 ) 03880 // im_atan_u = -ln( 1 - 4 * hyu / ( x_sqr + sqr( 1 + hyu ) ) ) / 4.0; 03881 im_atan_u = -Q_atan_UPSIGN(abs_rez,-hyu); // Blomquist 03882 else 03883 // im_atan_u = ln( 1 + 4 * hyu / ( x_sqr + sqr( 1 - hyu ) ) ) / 4.0; 03884 im_atan_u = Q_atan_UPSIGN(abs_rez,hyu); // Blomquist 03885 times2pown(im_atan_u,-2); // Division by 4; 03886 03887 resyl = min( Inf( im_atan_l ), Inf( im_atan_u ) ); 03888 resyu = max( Sup( im_atan_l ), Sup( im_atan_u ) ); 03889 // 03890 // look for extremal values on vertical boundaries, 03891 // if vertical boundaries intersect extremal curves 03892 // 03893 l_real abs_x_min = Inf( abs( rez ) ); 03894 l_interval 03895 x_extr = l_interval( abs_x_min ), 03896 // y_extr = sqrt( 1.0 + sqr( x_extr ) ); 03897 y_extr = sqrt1px2(x_extr); // Blomquist 03898 if( Inf( y_extr ) < simz && Sup( y_extr ) > iimz ) 03899 // extremal curve intersects left boundary of |x|+iy in quadrant I 03900 // update maximum 03901 // resyu = Sup( ln( 1 + 4 * y_extr / ( sqr( x_extr ) + sqr( 1 - y_extr ) ) ) ) / 4.0; 03902 // resyu = Sup( Aux_1_atan(abs_x_min)/4.0 ); // Blomquist 03903 { 03904 rez = Aux_1_atan(abs_x_min); 03905 times2pown(rez,-2); 03906 resyu = Sup(rez); 03907 } 03908 03909 if( -Sup( y_extr ) < simz && -Inf( y_extr ) > iimz ) 03910 // extremal curve intersects left boundary of |x|+iy in quadrant IV 03911 // update minimum 03912 // resyl = -Sup( ln( 1 + 4 * y_extr / ( sqr( x_extr ) + sqr( 1 - y_extr ) ) ) ) / 4.0; 03913 // resyl = -Sup( Aux_1_atan(abs_x_min)/4.0 ); // Blomquist 03914 { 03915 rez = Aux_1_atan(abs_x_min); 03916 times2pown(rez,-2); 03917 resyl = -Sup(rez); 03918 } 03919 03920 return l_cinterval(l_interval( resxl, resxu ), l_interval( -resyu, -resyl ) ); 03921 03922 } 03923 // 03924 //-- end acot ----------------------------------------------------------------- 03925 03926 03927 //-- atanh ---------------------------------------------------------- 040912 -- 03928 // 03929 l_cinterval atanh( const l_cinterval& z ) throw() 03930 // 03931 // atanh( z ) = - i * atan( i * z ) 03932 // 03933 { 03934 l_cinterval res = atan( l_cinterval( -Im(z), Re(z) ) ); 03935 return l_cinterval( Im(res), -Re(res) ); 03936 } 03937 // 03938 //-- end atanh ---------------------------------------------------------------- 03939 03940 //-- acoth ---------------------------------------------------------- 040912 -- 03941 // 03942 l_cinterval acoth( const l_cinterval& z ) throw() 03943 // 03944 // acoth( z ) = i * acot( i * z ) 03945 // 03946 { 03947 l_cinterval res = acot( l_cinterval( -Im(z), Re(z) ) ); 03948 return l_cinterval( -Im(res), Re(res) ); 03949 } 03950 // 03951 //-- end acoth ---------------------------------------------------------------- 03952 03953 } // namespace cxsc 03954 03955 /* 03956 03957 End of File: l_cimath.cpp 03958 03959 */