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