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