C-XSC - A C++ Class Library for Extended Scientific Computing  2.5.4
imath.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: imath.cpp,v 1.43 2014/01/30 17:23:45 cxsc Exp $ */
00025 
00026 #include "imath.hpp"
00027 #include "rmath.hpp"
00028 
00029 // Auch fi_lib hat ein eigenes interval (wie RTS sein dotprecision)
00030 //typedef struct fi_interval { double INF, SUP;} fi_interval;
00031 #undef LINT_ARGS
00032 #define CXSC_INCLUDE
00033 #include <fi_lib.hpp>
00034 
00035 extern "C" {
00036 #ifndef rfcth_included
00037 #define rfcth_included
00038 #include "r_fcth.h"
00039 #endif
00040 }
00041 
00042 namespace cxsc {
00043 
00044 
00045 interval sqr    (const interval &a) throw() 
00046 { 
00047   interval res;
00048   res= a*a;
00049   if (Inf(res)<0) Inf(res)=0;
00050   return res;
00051 }
00052 
00053 interval sqrt   (const interval &a, int n)  throw(STD_FKT_OUT_OF_DEF)
00054 { 
00055    if ( ((n>0) && (Inf(a)>=0.0)) || ((n<0) && (Inf(a)>0.0)) ) 
00056       return pow(a,interval(1.0,1.0)/n); 
00057    else { 
00058       cxscthrow(STD_FKT_OUT_OF_DEF("interval sqrt (const interval &a, int n)"));
00059       return interval(-1.0); // dummy result
00060    }
00061 }
00062 
00063 interval sqrt1px2(const interval& x) throw()
00064 // Inclusion of sqrt(1+x^2); Blomquist, 13.12.02;
00065 {
00066     interval t = abs(x),y;
00067     if (expo(Inf(t)) > 33)
00068     {
00069         y = t;
00070         Sup(y) = succ(Sup(y));
00071     } else if (expo(Sup(t)) > 33)
00072     {
00073         y = interval(Inf(t));  // interval(Inf(t)) is a point interval!
00074         y = sqrt(1+y*y);       // --->  y*y == sqr(y);
00075         y = interval(Inf(y),succ(Sup(t)));
00076      } else y = sqrt(1+sqr(t));
00077     return y;
00078 }
00079 
00080 interval sqrtx2y2(const interval& x, const interval& y) throw()
00081 // Inclusion of sqrt(x^2+y^2); Blomquist, 13.12.02;
00082 {
00083     interval a=abs(x), b=abs(y), r;
00084     int exa=expo(Sup(a)), exb=expo(Sup(b)), ex;
00085     if (exb > exa)
00086     {  // Permutation of a,b:
00087         r = a;  a = b;  b = r;
00088         ex = exa;  exa = exb;  exb = ex;
00089     }
00090     ex = 511 - exa;
00091     times2pown(a,ex);
00092     times2pown(b,ex);
00093     r = sqrt(a*a + b*b);
00094     times2pown(r,-ex);
00095     return r;
00096 } // sqrtx2y2
00097 
00098 //********************************************************************
00099 //  Constants for:   interval sqrtp1m1(const interval& x) throw()
00100 //  Blomquist 05.08.03
00101 const real Delta_f = 2*minreal;
00102 const real q_sqrtp1m1m = 9007199254740984.0 / 9007199254740992.0;
00103 const real q_sqrtp1m1p = 4503599627370502.0 / 4503599627370496.0;
00104 //********************************************************************
00105 interval sqrtp1m1(const interval& x) throw()
00106 // interval(a,b) is an inclusion of sqrt(x+1)-1;
00107 // Exported by imath.hpp;       Blomquist, 05.08.03;
00108 {
00109     real a=0,b=0,ix=Inf(x),sx=Sup(x);
00110     int ex_ix,ex_sx,sgn_ix,sgn_sx;
00111 
00112     // Calculating the lower bound:
00113     ex_ix = expo(ix);  sgn_ix = sign(ix);
00114     if (ex_ix<=-1021) //  ==> |ix| < 2^(-1021)
00115         { if (sgn_ix) a = sqrtp1m1(ix) - Delta_f; }
00116     else // |ix| >= 2^(-1021)
00117         if (ex_ix<=-51) 
00118         {
00119             times2pown(ix,-1); // exact division by 2
00120             a = pred(ix);
00121         }
00122         else 
00123             if (sgn_ix>0) a = (ix>0.67) ? 
00124                Inf( sqrt(interval(ix)+1)-1 ) : sqrtp1m1(ix)*q_sqrtp1m1m;
00125             else a = (ix<-0.25) ? 
00126                Inf( sqrt(interval(ix)+1)-1 ) : sqrtp1m1(ix)*q_sqrtp1m1p;
00127     // Calculating the upper bound:
00128     if (ix == sx) { ex_sx = ex_ix;  sgn_sx = sgn_ix; }
00129     else { ex_sx = expo(sx);  sgn_sx = sign(sx); }
00130     if (ex_sx<=-1021) //  ==> |sx| < 2^(-1021)
00131         { if (sgn_sx) b = sqrtp1m1(sx) + Delta_f; }
00132     else // |sx| >= 2^(-1021)
00133         if (ex_sx<=-47) { b = sx; times2pown(b,-1); }
00134         else 
00135             if (sgn_sx>0) b = (sx>0.58) ? 
00136                 Sup( sqrt(interval(sx)+1)-1 ) : sqrtp1m1(sx)*q_sqrtp1m1p;
00137             else b = (sx<-0.32) ? 
00138                  Sup( sqrt(interval(sx)+1)-1 ) : sqrtp1m1(sx)*q_sqrtp1m1m;
00139     return interval(a,b);
00140 } // sqrtp1m1()
00141 
00142 //********************************************************************
00143 // Konstanten fr die Intervallfunktion sqrt(x^2-1):
00144 real q_sqrtx2m1p(4503599627370501.0/4503599627370496.0);  // (1+e(f))
00145 real q_sqrtx2m1m(9007199254740986.0/9007199254740992.0);  // (1-e(f))
00146 //********************************************************************
00147 
00148 interval sqrtx2m1(const interval& x)
00149 // sqrt(x^2-1);  Blomquist, 13.04.04;
00150 {
00151     interval z = abs(x);
00152     real r1,r2;
00153     r1 = sqrtx2m1(Inf(z)) * q_sqrtx2m1m;
00154     r2 = Sup(z); 
00155     if (expo(r2)<26)  
00156         r2 = sqrtx2m1(r2) * q_sqrtx2m1p; 
00157     // expo(r2) >= 26 --> r2 = Sup(z) ist optimale Oberschranke!
00158     return interval(r1,r2);
00159 } // sqrtx2m1 (Intervallargumente)
00160 
00161 
00162 //********************************************************************
00163 // Konstanten fr die Intervallfunktion sqrt(1-x2):
00164 real q_sqrt1mx2p(4503599627370501.0/4503599627370496.0);  // (1+e(f))
00165 real q_sqrt1mx2m(9007199254740985.0/9007199254740992.0);  // (1-e(f))
00166 //********************************************************************
00167 
00168 interval sqrt1mx2(const interval& x)
00169 // sqrt(1-x2);  Blomquist, 13.04.04;
00170 {
00171     interval z = abs(x);  // sqrt(1-t2) monoton fallend in t aus z;
00172     real r1,r2,sz,iz;
00173     sz = Sup(z);  iz = Inf(z);
00174     // Berechnung der Unterschranke r1:
00175     r2 = sqrt1mx2(sz);
00176     r1 = (sz==0)? 1 : r2 * q_sqrt1mx2m;
00177     // Berechnung der Oberrschranke r2:
00178     if (iz<4.81e-8) r2 = 1; // r2=1 ist immer korrekte Oberschranke!
00179     else r2 = (sz==iz)? r2*q_sqrt1mx2p : sqrt1mx2(iz)*q_sqrt1mx2p;
00180     return interval(r1,r2);
00181 
00182 } // sqrtx2m1 (Intervallargumente)
00183 
00184 
00185 // Konstanten fuer die Intervallfunktion exp(-x^2):
00186 const real q_exp_x2p = 4503599627370502.0 / 4503599627370496.0;  // (1+e(f))
00187 const real q_exp_x2m = 9007199254740984.0 / 9007199254740992.0;  // (1-e(f))
00188   // Oberschranke expmx2_UB fuer |x| > expmx2_x0:  
00189 const real expmx2_UB = 2.225073858507447856659E-308;
00190 const real expmx2_x0 = 7491658466053896.0 / 281474976710656.0;
00191 
00192 interval expmx2(const interval& x)
00193 // e^(-x^2);  Blomquist, 05.07.04;
00194 {
00195     real y,r1,r2,Sz,Iz;
00196     interval z = abs(x); 
00197     // Berechnung einer Unterschranke:
00198     Sz = Sup(z);  Iz = Inf(z);
00199     y = expmx2(Sz);
00200     r1 = y;
00201     if (Sz<=expmx2_x0) r1 = r1 * q_exp_x2m; // Abrunden
00202     if (Sz==0) r1 = y;
00203     // Berechnung einer Oberschranke:
00204     if (Iz>expmx2_x0) r2 = expmx2_UB;
00205     else r2 = (Sz==Iz)? y * q_exp_x2p : expmx2(Iz) * q_exp_x2p; 
00206     if (r2>1) r2 = 1.0;
00207     return interval(r1,r2);
00208 
00209 } // expmx2 (Intervallargumente)
00210 
00211 
00212 interval expm1(const interval& x)
00213 // Interval function of exp(x)-1
00214 {
00215     real y,r1,r2,Sx,Ix;
00216     Sx = Sup(x);  Ix = Inf(x);
00217     y = expm1(Ix);
00218     // calculation of a lower bound r1:
00219     r1 = (y>0)? y*q_exmm : y*q_exmp;  // rounding downwards
00220     if (r1<-1) r1 = -1;
00221     // calculation of an upper bound r2:
00222     if (Sx!=Ix) y = expm1(Sx);
00223     r2 = (y>0)? y*q_exmp : y*q_exmm;
00224 
00225     return interval(r1,r2);
00226 } // expm1(...) 
00227 
00228 // Die folgenden Konstanten q_expx2_p, q_expx2_m beziehen sich 
00229 // auf eps(expx2) = 4.618958e-16;    f(x) = e^{+x^2};
00230 
00231 const real q_expx2_p = 4503599627370500.0 / 4503599627370496.0;  // (1+e(f))
00232 const real q_expx2_m = 9007199254740985.0 / 9007199254740992.0;  // (1-e(f))
00233 
00234 interval expx2(const interval& x)
00235 // e^(+x^2);  Blomquist, 25.07.06;
00236 {
00237     real y,r1,r2,Sz,Iz;
00238     interval z = abs(x);
00239     Sz = Sup(z);  Iz = Inf(z);
00240     // Berechnung einer Unterschranke:
00241     y = expx2(Iz);
00242     r1 = y;
00243     r1 = r1 * q_expx2_m; // Abrunden
00244     if (r1<1.0) r1 = 1.0;
00245     // Berechnung einer Oberschranke:
00246     r2 = (Sz==Iz)? y * q_expx2_p : expx2(Sz) * q_expx2_p;
00247     if (Sz==0) r2 = 1.0;
00248 
00249     return interval(r1,r2);
00250 } // expx2 (Intervallfunktion)
00251 
00252 // Die folgenden Konstanten q_expx2m1_p, q_expx2m1_m beziehen sich 
00253 // auf eps(expx2m1) = 4.813220e-16;    f(x) = e^{+x^2}-1;
00254 
00255 const real q_expx2m1_p = 4503599627370500.0 / 4503599627370496.0; // (1+e(f))
00256 const real q_expx2m1_m = 9007199254740985.0 / 9007199254740992.0; // (1-e(f))
00257 const real expx2m1_0   = comp(0.5,-510);  // 2^(-511)
00258 
00259 void sqr2uv(const real&, real&, real&);
00260 
00261 real expx2m1_intv(const real& x)
00262 // Zur Implementierung der Intervallfunktion;
00263 // e^(+x^2)-1; rel. Fehlerschranke: eps = 4.813220E-16 = e(f) gilt
00264 // fuer alle x, mit:  2^(-511) <= x <= x0 = 26.64174755704632....
00265 // x0 = 7498985273150791.0 / 281474976710656.0;
00266 // Fuer x > x0 --> Programmabbruch wegen Overflow;
00267 // Fuer 0 <= x < 2^(-511) werden die Funktionswerte auf Null gesetzt.
00268 // Ausfuehrlich getestet;  Blomquist, 10.08.2006;
00269 {
00270     real t(x),u,v,y,res(0);
00271     int ex;
00272     if (t<0) t = -t;  // t >= 0;
00273     
00274     if (t>=6.5) res = expx2(t);
00275     else 
00276     {
00277         ex = expo(t);
00278         sqr2uv(x,u,v);  // u := x*x und v aus S(2,53); 
00279         if (ex>=2) // g4(x)
00280         {
00281             y = exp(u); 
00282             res = 1 - v*y;
00283             res = y - res;
00284         }
00285         else 
00286             if (ex>=-8) res = expm1(u) + v*exp(u); // g3(x)
00287             else 
00288                 if (ex>=-25) { // g2(x)
00289                     y = u*u;
00290                     times2pown(y,-1);
00291                     res = (1+u/3)*y + u;
00292                 }
00293                 else 
00294                     if(ex>=-510) res = u;  // g1(x)
00295     }
00296 
00297     return res;
00298 } // expx2m1_intv
00299 
00300 interval expx2m1(const interval& x)
00301 // e^(+x^2)-1;  Blomquist, 10.08.06;
00302 {
00303     real y,r1,r2,Sz,Iz;
00304     interval z = abs(x);
00305     Sz = Sup(z);  Iz = Inf(z);
00306     // Berechnung einer Unterschranke:
00307     y = expx2m1_intv(Iz);
00308     r1 = y;
00309     r1 = r1 * q_expx2m1_m; // Abrunden
00310     // Berechnung einer Oberschranke:
00311     if (Sz < expx2m1_0) 
00312     { 
00313         r2 = MinReal;
00314         if (Sz==0) r2 = 0.0;
00315     }
00316     else r2 = (Sz==Iz)? y * q_expx2m1_p : expx2m1_intv(Sz) * q_expx2m1_p;
00317 
00318     return interval(r1,r2);
00319 } // expx2m1 (Intervallfunktion)
00320 
00321 // ------  1-eps and 1+eps for function lnp1, Blomquist 28,07,03;  -----------
00322 // ------------------------  eps = 2.5082e-16  -------------------------------
00323 static real q_lnp1m = 9007199254740986.0 / 9007199254740992.0; // 1-eps
00324 static real q_lnp1p = 4503599627370501.0 / 4503599627370496.0; // 1+eps
00325 // ---------------------------------------------------------------------------
00326 
00327 interval lnp1(const interval& x) throw()
00328 // returns an inclusion of ln(1+t), with t in x; Blomquist 28.07.03;
00329 { 
00330     real ix=Inf(x), sx=Sup(x),a,b;   // ln(1+x) <= [a,b]
00331     // Calculating the lower bound a:
00332     int sgn_ix = sign(ix), ex_ix = expo(ix), sgn_sx,ex_sx;
00333     if (!sgn_ix) a = 0;  // optimal lower bound!
00334     else if (sgn_ix>0)
00335              a = (ex_ix<=-53) ? pred(ix) : lnp1(ix) * q_lnp1m;
00336          else 
00337              a = (ex_ix<=-54) ? pred(ix) : lnp1(ix) * q_lnp1p;
00338     // Calculating the upper bound b:
00339     if (ix == sx) { sgn_sx = sgn_ix; ex_sx = ex_ix; }
00340     else { sgn_sx = sign(sx); ex_sx = expo(sx); }
00341     if (sgn_sx>=0)
00342         b = (ex_sx<=-49) ? sx : lnp1(sx) * q_lnp1p;
00343     else // sx < 0:
00344         b = (ex_sx<=-50) ? sx : lnp1(sx) * q_lnp1m;
00345     return interval(a,b); // ln(1+x) in [a,b]
00346 } // lnp1
00347 
00354 interval erf    (const interval &a)         { return j_erf(a);  }
00361 interval erfc   (const interval &a)         { return j_erfc(a); }
00362 
00363 //interval pow    (const interval &a, const interval &b) throw(ERROR_INTERVAL_STD_FKT_OUT_OF_DEF)
00364 //{
00365 //      if(Inf(a)>0)
00366 //              return j_exp(b*ln(a));
00367 //      else if(Inf(a)==0 && Inf(b)>=0)
00368 //      {
00369 //              if(Sup(a)>=1)
00370 //                      return interval(0,pow(Sup(a),Sup(b)));
00371 //              else
00372 //                      return interval(0,pow(Sup(a),Inf(b)));
00373 //      }
00374 //      else
00375 //      {
00376 //              cxscthrow(ERROR_INTERVAL_STD_FKT_OUT_OF_DEF("interval pow(const interval &,const interval &)"));
00377 //              return interval(0,0);
00378 //      }
00379 //}
00380 
00381 inline a_intv _a_intv(const interval &x)
00382 {       
00383        return *((const a_intv *)(&x));
00384 }
00385 inline interval _interval(const a_intv &x)
00386 {       
00387        return *((const interval *)(&x));
00388 }
00389 
00390 interval pow    (const interval &a, const interval &b) throw()
00391        { 
00392          interval res; 
00393          if(Inf(a)==0 && Inf(b)>=0)
00394          {
00395            if(Sup(a)>=1)
00396              res=interval(0,pow(Sup(a),Sup(b)));
00397            else
00398              res=interval(0,pow(Sup(a),Inf(b)));
00399          }
00400          else res = _interval(i_pow(_a_intv(a),_a_intv(b))); 
00401 
00402          if (Inf(res) <= Sup(res))
00403            return res;
00404          else
00405            return interval(Sup(res),Inf(res));
00406        }
00407        
00408 //----------------------------------------------------------------------------
00409 // Purpose: The local function 'Power()' is used to compute a lower or an
00410 //    upper bound for the power function with real argument and integer
00411 //    exponent, respectively.
00412 // Parameters:
00413 //    In: 'x'      : real argument
00414 //        'n'      : integer exponent
00415 //        'RndMode': rounding mode,
00416 //                   (-1 = downwardly directed, +1 = upwardly directed)
00417 // Description:
00418 //    This function is used to speed up the interval power function defined
00419 //    below. The exponentiation is reduced to multiplications using the
00420 //    binary shift method. Depending on 'n', this function is up to 40 times
00421 //    as fast as the standard power function for real argument and real
00422 //    exponent. However, its accuracy is less than one ulp (unit in the last
00423 //    place of the mantissa) since about log2(n) multiplications are executed
00424 //    during computation. Since directed roundings are antisymmetric, one
00425 //    gets
00426 //
00427 //       down(x^n) = -up((-x)^n)   and   up(x^n) = -down((-x)^n)
00428 //
00429 //    for x < 0 and odd n, where 'down' and 'up' denote the downwardly and
00430 //    upwardly directed roundings, respectively.
00431 //----------------------------------------------------------------------------
00432 static real Power (const real & x, int n, int RndMode )
00433 {                         // Signals change of the rounding mode
00434   int  ChangeRndMode;     // for x < 0 and odd n
00435   real p, z;
00436 
00437   ChangeRndMode = ( (x < 0.0) && (n % 2 == 1) );
00438 
00439   if (ChangeRndMode) {
00440     z = -x;
00441     RndMode = -RndMode;
00442   }
00443   else
00444     z = x;
00445 
00446   p = 1.0;
00447   switch (RndMode) {                             // Separate while-loops used
00448     case -1 : while (n > 0) {                    // to gain speed at runtime
00449                 if (n % 2 == 1) p = muld(p,z);   //--------------------------
00450                 n = n / 2;
00451                 if (n > 0) z = muld(z,z);
00452               }
00453               break;
00454     case +1 : while (n > 0) {
00455                 if (n % 2 == 1) p = mulu(p,z);
00456                 n = n / 2;
00457                 if (n > 0) z = mulu(z,z);
00458               }
00459               break;
00460   }
00461 
00462   if (ChangeRndMode)
00463     return -p;
00464   else
00465     return p;
00466 }
00467 
00468 interval power(const interval& a, int n)
00469 // Calculating a^n; 
00470 // Examples: [-1,4]^2 = [0,16];  [-1,4]^3 = [-16,+64];   
00471 {
00472     bool neg(n<0);
00473     int N(n);          //,k(-1),r; 
00474     interval res,h;
00475     real Lower, Upper;
00476     if (neg) N = -N;
00477     if (N==0) res = 1;
00478     else if (N==1) res = a;
00479     else // N > 1:
00480         if (Inf(a)>=MinReal) res = exp(N*ln(a));
00481         else if (Sup(a)<=-MinReal) {
00482             h = interval(-Sup(a),-Inf(a));
00483             res = exp(N*ln(h));
00484             if (N%2 == 1) res = -res;
00485         } 
00486         else 
00487         {
00488 //          h = a;
00489 //          while(N>0) 
00490 //          {
00491 //              k++;
00492 //              r = N % 2;
00493 //              if (k==0) 
00494 //                  if (r==1) res=a; else res=1;
00495 //              else {
00496 //                  h = sqr(h);
00497 //                  if (r==1) res *= h;
00498 //              }
00499 //              N = N / 2;
00500 //          }
00501           if ( (0.0 < Inf(a)) || (N % 2 == 1) ) {
00502             Lower = Power(Inf(a),N,-1);
00503             Upper = Power(Sup(a),N,+1);
00504           }
00505           else if (0.0 > Sup(a)) {
00506             Lower = Power(Sup(a),N,-1);
00507             Upper = Power(Inf(a),N,+1);
00508           }
00509           else {
00510           Lower = 0.0;
00511           Upper = Power(AbsMax(a),N,+1);
00512           }
00513           res = interval(Lower,Upper);
00514         }
00515     if (neg) res = 1/res;
00516     return res;
00517 } // power
00518 
00519 
00520 //----------------------------------------------------------------------------
00521 // Purpose: This version of the function 'Power()' is used to compute an
00522 //    enclosure for the power function with interval argument and integer
00523 //    exponent.
00524 // Parameters:
00525 //    In: 'x': interval argument
00526 //        'n': integer exponent
00527 // Description:
00528 //    In general, this implementation does not deliver a result of maximum
00529 //    accuracy, but it is about 30-40 times faster than the standard power
00530 //    function for interval arguments and interval exponents. The resulting
00531 //    interval has a width of approximately 2*log2(n) ulps. Since x^n is
00532 //    considered as a monomial, we define x^0 := 1. For negative exponents
00533 //    and 0 in 'x', the division at the end of the function will cause a
00534 //    runtime error (division by zero).
00535 //----------------------------------------------------------------------------
00536 interval Power (const interval & x, int n )
00537 {
00538   int  m;
00539   real Lower, Upper;
00540 
00541   if (n == 0) return(interval(1.0,1.0));
00542 
00543   if (n > 0)  m = n;  else  m = -n;
00544 
00545   if ( (0.0 < Inf(x)) || (m % 2 == 1) ) {
00546     Lower = Power(Inf(x),m,-1);
00547     Upper = Power(Sup(x),m,+1);
00548   }
00549   else if (0.0 > Sup(x)) {
00550     Lower = Power(Sup(x),m,-1);
00551     Upper = Power(Inf(x),m,+1);
00552   }
00553   else {
00554     Lower = 0.0;
00555     Upper = Power(AbsMax(x),m,+1);
00556   }
00557 
00558   if (n > 0)
00559     return(interval(Lower,Upper));
00560   else                                    // Causes a runtime error
00561     return(1.0/interval(Lower,Upper));    // if 0 in 'x'.
00562 }
00563 
00564 //----------------------------------------------------------------------------
00565 
00566 interval Pi ( )                                    // Enclosure of constant pi
00567   { return(4.0*atan(_interval(1.0,1.0))); }        //-------------------------
00568 
00569 // Error bounds for the interval function ln_sqrtx2y2:
00570 real ln_x2y2_abs(2.225076E-308); // Absolute error bond
00571 real q_lnx2y2p(4503599627370502.0 / 4503599627370496.0); // 1+e(f)
00572 real q_lnx2y2m(9007199254740984.0 / 9007199254740992.0); // 1-e(f)
00573 // With the following b0 
00574 // real b0 = 6369051672525773.0 / 30191699398572330817932436647906151127335369763331523427009650401964993299137190816689013801421270140331747000246110759198164677039398341060491474011461568349195162615808.0;
00575 real b0 = MakeHexReal(0,1022-510,0x0016A09E,0x667F3BCD);
00576 // it holds:
00577 // 1. b < bo  ==> g(b) := (0.5*b)*b < MinReal with rounding downwards
00578 // 2. b >= b0 ==> g(b) := (0.5*b)*b >= MinReal with arbitrary rounding
00579 //                                   modus by the two multiplications.
00580 
00581 interval ln_sqrtx2y2(const interval& x, const interval& y) throw()
00582 // ln( sqrt(x^2+y^2) ) == 0.5*ln(x^2+y^2);   Blomquist, 22.11.03;
00583 {
00584     interval ax=abs(x), ay=abs(y);
00585     real Ix=Inf(ax), Sx=Sup(ax), Iy=Inf(ay), Sy=Sup(ay),f,u1,u2;
00586     // Calculating the lower bound u1:
00587     f = ln_sqrtx2y2(Ix,Iy);
00588     if ((Ix==1 && Iy<b0) || (Iy==1 && Ix<b0)) {
00589     // f in the denormalized range!
00590         u1 = f - ln_x2y2_abs;   // directed rounding not necessary!
00591         if (sign(u1)<0) u1 = 0;
00592     } else  u1 = (sign(f)<0) ? f*q_lnx2y2p : f*q_lnx2y2m;
00593     // Calculating the upper bound u2:
00594     if (Ix==Sx && Iy==Sy) // x and y are point-intervals
00595         if ((Sx==1 && Sy<b0) || (Sy==1 && Sx<b0)) {
00596         // f in the denormalized range!
00597             u2 = (Sy==0 || Sx==0) ? f : f+ln_x2y2_abs; 
00598         } else  u2 = (sign(f)<0) ? f*q_lnx2y2m : f*q_lnx2y2p;
00599     else // x or y is no point-interval:
00600     {
00601         f = ln_sqrtx2y2(Sx,Sy);
00602         if ((Sx==1 && Sy<b0) || (Sy==1 && Sx<b0))
00603         // f in the denormalized range!
00604             u2 = (sign(Sy)==0 || sign(Sx)==0) ? f : f+ln_x2y2_abs; 
00605         else  u2 = (sign(f)<0) ? f*q_lnx2y2m : f*q_lnx2y2p;
00606     }
00607     return interval(u1,u2);
00608 } // ln_sqrtx2y2
00609 
00610 
00611 // Constants for the interval function acoshp1(x) = acosh(1+x):
00612 // (1+e(f)):
00613 static const real q_acoshp1p(4503599627370503.0/4503599627370496.0);  
00614 // (1-e(f))
00615 static const real q_acoshp1m(9007199254740981.0/9007199254740992.0);  
00616 
00617 interval acoshp1(const interval& x)
00618 // acoshp1;  Blomquist, 28.03.2005;
00619 {
00620     real r1,r2,sx,ix;
00621     sx = Sup(x);  ix = Inf(x);
00622     // Calculating of the lower bound r1:
00623     r2 = acoshp1(ix);
00624     r1 = r2 * q_acoshp1m;
00625     // Calculating of the upper bound r2:
00626     r2 = (sx==ix)? r2*q_acoshp1p : acoshp1(sx)*q_acoshp1p;
00627     return interval(r1,r2);
00628 
00629 } // acoshp1 (interval arguments)
00630 
00631 // *****************************************************************************
00632 //                               sin(Pi*x)/Pi                                  *
00633 // ********************************************************************************
00634 
00635 // Die folgenden Konstanten q_sinpix_p, q_sinpix_m beziehen sich 
00636 // auf eps(sinpix_pi) = 3.401895e-16;    f(x) = sin(pi*x)/pi;
00637 
00638 const real q_sinpix_p = 4503599627370499.0 / 4503599627370496.0;  // (1+e(f))
00639 const real q_sinpix_m = 9007199254740986.0 / 9007199254740992.0;  // (1-e(f))
00640 
00641 real rounded_up(const real& x)
00642 {
00643         real y;
00644         y = (x>=0)? x*q_sinpix_p : x*q_sinpix_m;
00645         return y;
00646 }
00647 
00648 real rounded_down(const real& x)
00649 {
00650         real y;
00651         y = (x>=0)? x*q_sinpix_m : x*q_sinpix_p;
00652         return y;
00653 }
00654 
00655 interval sinpix_pi(const interval& x)
00656 {
00657         const real Pir = Sup(Pir_interval); // 1/pi upwards rounded 
00658         interval y;
00659         int ma,mb;
00660         real y1,y2,a(Inf(x)),b(Sup(x)),ya,yb;
00661         bool bl, ya_klg_yb;
00662                 
00663         ma = Round(a);  mb = Round(b);
00664         bl = (ma%2 != 0);
00665         if (mb-ma>1)
00666         {
00667                 y1 = -Pir;  y2 = Pir;
00668         }
00669         else // 0 <= mb-ma <=1
00670                 if (mb==ma)
00671                         if (a==b)  // x: Point interval
00672                         {
00673                                 ya = sinpix_pi(a);
00674             y1 = rounded_down(ya);
00675                                 y2 = rounded_up(ya);
00676                         }
00677                         else // (mb==ma)  and  (a!=b)
00678                         {
00679                                 ya = sinpix_pi(a);
00680                                 yb = sinpix_pi(b);
00681                                 if (!bl)
00682                                 {
00683                                         y1 = rounded_down(ya);
00684                                         y2 = rounded_up(yb);
00685                                 }
00686                                 else
00687                                 {
00688                                         y1 = rounded_down(yb);
00689                                         y2 = rounded_up(ya);
00690                                 }
00691                         }
00692                 else // mb-ma=1;
00693                 {
00694                         ya = sinpix_pi(a);
00695                         yb = sinpix_pi(b);
00696                         ya_klg_yb = (ya <= yb);
00697                         
00698                         if (bl)
00699                         {
00700                                 if (!ya_klg_yb)
00701                                         yb = ya;
00702                                 ya = -Pir;
00703                         }
00704                         else
00705                         {
00706                                 if (!ya_klg_yb)
00707                                         ya = yb;
00708                                 yb = Pir;
00709                         }
00710                         y1 = rounded_down(ya);
00711                         if (y1<-Pir) 
00712                                 y1 = -Pir;
00713                         y2 = rounded_up(yb);
00714                         if (y2>Pir)
00715                                 y2 = Pir;
00716                 }
00717                 
00718         y = interval(y1,y2);
00719         return y;
00720 }
00721 
00722 // ********************************************************************************
00723 // ***************  Konstanten bez. Gamma(x) und 1/Gamma(x) ********************
00724 // ********************************************************************************
00725 
00726 /* worst case relative error bound for 1/Gamma(x)            */    
00727 /* eps(gammar)   = 2.866906e-15;                             */
00728 /*     q_gammarm  = 1 - eps(gammar)                          */
00729 const real   q_gammarm  = 9007199254740964.0 / 9007199254740992.0;
00730 /*     q_gammarp  = 1 + eps(gammar)                       */
00731 const   real   q_gammarp  = 4503599627370510.0 / 4503599627370496.0;
00732 
00733 // Inclusion of the extremes of 1/gamma(x):
00734 interval pow2(const interval& x, int ex)
00735 {
00736         interval y(x);
00737    times2pown(y,ex);
00738         return y;
00739 }
00740                                                                          
00741 const real Ne =  9007199254740992.0;
00742 const real Ne1 = 1125899906842624.0;
00743 const real Ne2 = 562949953421312.0;
00744 const real Ne3 = 281474976710656.0;  
00745 const real Ne4 = 140737488355328.0;  
00746 const real Ne5 = 70368744177664.0;   
00747 const real Ne6 = 35184372088832.0;
00748                     
00749 const interval gam_rxi[171] =
00750 {
00751         interval( 6582605572834349.0 / 4503599627370496.0,6582606400588712.0 / 
00752                        4503599627370496.0 ),
00753         interval( -4540376432147063.0 / 9007199254740992.0,-4540375772996112.0 /
00754                             9007199254740992.0 ),
00755         interval( -7086407292338520.0 / 4503599627370496.0,-7086406981597106.0 /
00756                             4503599627370496.0 ),
00757         interval( -5878820838740338.0 / 2251799813685248.0,-5878820690102701.0 /
00758                             2251799813685248.0 ),
00759         interval( -8185952996852629.0 / 2251799813685248.0,-8185952850644519.0 /
00760                             2251799813685248.0 ),
00761         interval( -5239079997162568.0 / Ne1,-5239079928185648.0 /
00762                             Ne1 ),
00763         interval( -6380657697812205.0 / Ne1,-6380657632438250.0 /
00764                             Ne1 ),
00765         interval( -7519230477777525.0 / Ne1,-7519230410301402.0 /
00766                             Ne1 ),
00767    interval( -8655680190901081.0 / Ne1,-8655680125714323.0 /
00768                             Ne1 ),
00769         interval( -4895280046470312.0 / Ne2,-4895280015191181.0 /
00770                             Ne2 ),
00771         interval( -5462119069950045.0 / Ne2,-5462119039144308.0 /
00772                             Ne2 ),
00773         interval( -6028485171921533.0 / Ne2,-6028485140574985.0 /
00774                             Ne2 ),
00775         interval( -6594470676196825.0 / Ne2,-6594470646005838.0 /
00776                             Ne2 ),
00777         interval( -7160144161412306.0 / Ne2,-7160144131821972.0 /
00778                             Ne2 ),
00779         interval( -7725557826948019.0 / Ne2,-7725557796923813.0 /
00780                             Ne2 ),
00781         interval( -8290752238810453.0 / Ne2,-8290752208368199.0 /
00782                             Ne2 ),
00783         interval( -8855759486553113.0 / Ne2,-8855759457969009.0 /
00784                             Ne2 ),
00785         interval( -4710302676530551.0 / Ne3,-4710302661730969.0 /
00786                             Ne3 ),
00787         interval( -4992655414739558.0 / Ne3,-4992655400196254.0 /
00788                             Ne3 ),
00789         interval( -5274946608174960.0 / Ne3,-5274946594248303.0 /
00790                             Ne3 ),
00791         interval( -5557183461054268.0 / Ne3,-5557183446978818.0 /
00792                             Ne3 ),
00793         interval( -5839372030862353.0 / Ne3,-5839372016359958.0 /
00794                             Ne3 ),
00795         interval( -6121517453741464.0 / Ne3,-6121517439898965.0 /
00796                             Ne3 ),
00797         interval( -6403624121061720.0 / Ne3,-6403624107245033.0 /
00798                             Ne3 ),
00799         interval( -6685695811452843.0 / Ne3,-6685695797850064.0 /
00800                             Ne3 ),
00801         interval( -6967735798799276.0 / Ne3,-6967735784842353.0 /
00802                             Ne3 ),
00803         interval( -7249746935136834.0 / Ne3,-7249746921647546.0 /
00804                             Ne3 ),
00805         interval( -7531731721183113.0 / Ne3,-7531731707547301.0 /
00806                             Ne3 ),
00807         interval( -7813692357395941.0 / Ne3,-7813692343778691.0 /
00808                             Ne3 ),
00809         interval( -8095630791428232.0 / Ne3,-8095630778040390.0 /
00810                             Ne3 ),
00811         interval( -8377548753853165.0 / Ne3,-8377548740538224.0 /
00812                             Ne3 ),
00813         interval( -8659447788741678.0 / Ne3,-8659447775645579.0 /
00814                             Ne3 ),
00815         interval( -8941329278863280.0 / Ne3,-8941329265761967.0 /
00816                             Ne3 ),
00817         interval( -4611597233397515.0 / Ne4,-4611597226897947.0 /
00818                             Ne4 ),
00819         interval( -4752522236545681.0 / Ne4,-4752522230097182.0 /
00820                             Ne4 ),
00821         interval( -4893440155674552.0 / Ne4,-4893440149271021.0 /
00822                             Ne4 ),
00823         interval( -5034351450592848.0 / Ne4,-5034351444042898.0 /
00824                             Ne4 ),
00825         interval( -5175256539394880.0 / Ne4,-5175256533033078.0 /
00826                             Ne4 ),
00827         interval( -5316155803584886.0 / Ne4,-5316155797241740.0 /
00828                             Ne4 ),
00829         interval( -5457049591969558.0 / Ne4,-5457049585698186.0 /
00830                             Ne4 ),
00831         interval( -5597938224227692.0 / Ne4,-5597938218047102.0 /
00832                             Ne4 ),
00833         interval( -5738821993949958.0 / Ne4,-5738821987610287.0 /
00834                             Ne4 ),
00835         interval( -5879701171316002.0 / Ne4,-5879701164978161.0 /
00836                             Ne4 ),
00837         interval( -6020576005634620.0 / Ne4,-6020575999374578.0 /
00838                             Ne4 ),
00839         interval( -6161446727231484.0 / Ne4,-6161446721061003.0 /
00840                             Ne4 ),
00841         interval( -6302313549465809.0 / Ne4,-6302313543147939.0 /
00842                             Ne4 ),
00843         interval( -6443176669927728.0 / Ne4,-6443176663720294.0 /
00844                             Ne4 ),
00845         interval( -6584036272403005.0 / Ne4,-6584036266213061.0 /
00846                             Ne4 ),
00847         interval( -6724892527807298.0 / Ne4,-6724892521515390.0 /
00848                             Ne4 ),
00849         interval( -6865745595463323.0 / Ne4,-6865745589218114.0 /
00850                             Ne4 ),
00851         interval( -7006595624078029.0 / Ne4,-7006595617954995.0 /
00852                             Ne4 ),
00853         interval( -7147442752315627.0 / Ne4,-7147442746313355.0 /
00854                             Ne4 ),
00855         interval( -7288287110549528.0 / Ne4,-7288287104399241.0 /
00856                             Ne4 ),
00857         interval( -7429128820262002.0 / Ne4,-7429128814290427.0 /
00858                             Ne4 ),
00859         interval( -7569967996009183.0 / Ne4,-7569967989919521.0 /
00860                             Ne4 ),
00861         interval( -7710804744971319.0 / Ne4,-7710804738911590.0 /
00862                             Ne4 ),
00863         interval( -7851639168099204.0 / Ne4,-7851639162048862.0 /
00864                             Ne4 ),
00865         interval( -7992471360380410.0 / Ne4,-7992471354478088.0 /
00866                             Ne4 ),
00867         interval( -8133301411685542.0 / Ne4,-8133301405639704.0 /
00868                             Ne4 ),
00869         interval( -8274129406251900.0 / Ne4,-8274129400330745.0 /
00870                             Ne4 ),
00871         interval( -8414955423947592.0 / Ne4,-8414955418025784.0 /
00872                             Ne4 ),
00873         interval( -8555779540237542.0 / Ne4,-8555779534343191.0 /
00874                             Ne4 ),
00875         interval( -8696601826519818.0 / Ne4,-8696601820560983.0 /
00876                             Ne4 ),
00877         interval( -8837422350374443.0 / Ne4,-8837422344547779.0 /
00878                             Ne4 ),
00879         interval( -8978241175812537.0 / Ne4,-8978241170031881.0 /
00880                             Ne4 ),
00881         interval( -4559529181955286.0 / Ne5,-4559529178973419.0 /
00882                             Ne5 ),
00883         interval( -4629936985962592.0 / Ne5,-4629936983062621.0 /
00884                             Ne5 ),
00885         interval( -4700344027557972.0 / Ne5,-4700344024658010.0 /
00886                             Ne5 ),
00887         interval( -4770750332748164.0 / Ne5,-4770750329806175.0 /
00888                             Ne5 ),
00889         interval( -4841155926312136.0 / Ne5,-4841155923514622.0 /
00890                             Ne5 ),
00891         interval( -4911560831959402.0 / Ne5,-4911560829205074.0 /
00892                             Ne5 ),
00893         interval( -4981965072279533.0 / Ne5,-4981965069314643.0 /
00894                             Ne5 ),
00895         interval( -5052368668591817.0 / Ne5,-5052368665644419.0 /
00896                             Ne5 ),
00897         interval( -5122771641448337.0 / Ne5,-5122771638497776.0 /
00898                             Ne5 ),
00899         interval( -5193174010445884.0 / Ne5,-5193174007591888.0 /
00900                             Ne5 ),
00901         interval( -5263575794318673.0 / Ne5,-5263575791488264.0 /
00902                             Ne5 ),
00903         interval( -5333977010931459.0 / Ne5,-5333977008101358.0 /
00904                             Ne5 ),
00905         interval( -5404377677505176.0 / Ne5,-5404377674566742.0 /
00906                             Ne5 ),
00907         interval( -5474777810213221.0 / Ne5,-5474777807391522.0 /
00908                             Ne5 ),
00909         interval( -5545177424917448.0 / Ne5,-5545177422065289.0 /
00910                             Ne5 ),
00911         interval( -5615576536591204.0 / Ne5,-5615576533728944.0 /
00912                             Ne5 ),
00913         interval( -5685975159660620.0 / Ne5,-5685975156829327.0 /
00914                             Ne5 ),
00915         interval( -5756373307993982.0 / Ne5,-5756373305121310.0 /
00916                             Ne5 ),
00917         interval( -5826770994840416.0 / Ne5,-5826770992034928.0 /
00918                             Ne5 ),
00919         interval( -5897168232948637.0 / Ne5,-5897168230130823.0 /
00920                             Ne5 ),
00921         interval( -5967565034702908.0 / Ne5,-5967565031698652.0 /
00922                             Ne5 ),
00923         interval( -6037961411567024.0 / Ne5,-6037961408739972.0 /
00924                             Ne5 ),
00925         interval( -6108357375160195.0 / Ne5,-6108357372312820.0 /
00926                             Ne5 ),
00927         interval( -6178752936267084.0 / Ne5,-6178752933424618.0 /
00928                             Ne5 ),
00929         interval( -6249148105390399.0 / Ne5,-6249148102629437.0 /
00930                             Ne5 ),
00931         interval( -6319542892697349.0 / Ne5,-6319542889886519.0 /
00932                             Ne5 ),
00933         interval( -6389937307844225.0 / Ne5,-6389937305053873.0 /
00934                             Ne5 ),
00935         interval( -6460331360217357.0 / Ne5,-6460331357417915.0 /
00936                             Ne5 ),
00937         interval( -6530725058889375.0 / Ne5,-6530725056107292.0 /
00938                             Ne5 ),
00939         interval( -6601118412624317.0 / Ne5,-6601118409842826.0 /
00940                             Ne5 ),
00941         interval( -6671511429748937.0 / Ne5,-6671511426971253.0 /
00942                             Ne5 ),
00943         interval( -6741904118401788.0 / Ne5,-6741904115722630.0 /
00944                             Ne5 ),
00945         interval( -6812296486576567.0 / Ne5,-6812296483762819.0 /
00946                             Ne5 ),
00947         interval( -6882688541678557.0 / Ne5,-6882688538949795.0 /
00948                             Ne5 ),
00949         interval( -6953080291125250.0 / Ne5,-6953080288362282.0 /
00950                             Ne5 ),
00951         interval( -7023471742013192.0 / Ne5,-7023471739228496.0 /
00952                             Ne5 ),
00953         interval( -7093862901097948.0 / Ne5,-7093862898353711.0 /
00954                             Ne5 ),
00955         interval( -7164253775147266.0 / Ne5,-7164253772348037.0 /
00956                             Ne5 ),
00957         interval( -7234644370421069.0 / Ne5,-7234644367656947.0 /
00958                             Ne5 ),
00959         interval( -7305034693228153.0 / Ne5,-7305034690493625.0 /
00960                             Ne5 ),
00961         interval( -7375424749560875.0 / Ne5,-7375424746816433.0 /
00962                             Ne5 ),
00963         interval( -7445814545246703.0 / Ne5,-7445814542529259.0 /
00964                             Ne5 ),
00965         interval( -7516204085881418.0 / Ne5,-7516204083182192.0 /
00966                             Ne5 ),
00967         interval( -7586593377039902.0 / Ne5,-7586593374291972.0 /
00968                             Ne5 ),
00969         interval( -7656982423878545.0 / Ne5,-7656982421139464.0 /
00970                             Ne5 ),
00971         interval( -7727371231629585.0 / Ne5,-7727371228950751.0 /
00972                             Ne5 ),
00973         interval( -7797759805243431.0 / Ne5,-7797759802601748.0 /
00974                             Ne5 ),
00975         interval( -7868148149631789.0 / Ne5,-7868148146919441.0 /
00976                             Ne5 ),
00977         interval( -7938536269389156.0 / Ne5,-7938536266702540.0 /
00978                             Ne5 ),
00979         interval( -8008924169145439.0 / Ne5,-8008924166482991.0 /
00980                             Ne5 ),
00981         interval( -8079311853288492.0 / Ne5,-8079311850638146.0 /
00982                             Ne5 ),
00983         interval( -8149699326208468.0 / Ne5,-8149699323496057.0 /
00984                             Ne5 ),
00985         interval( -8220086591946405.0 / Ne5,-8220086589180599.0 /
00986                             Ne5 ),
00987         interval( -8290473654625668.0 / Ne5,-8290473651911160.0 /
00988                             Ne5 ),
00989         interval( -8360860518207539.0 / Ne5,-8360860515490326.0 /
00990                             Ne5 ),
00991         interval( -8431247186492146.0 / Ne5,-8431247183838485.0 /
00992                             Ne5 ),
00993         interval( -8501633663256230.0 / Ne5,-8501633660616442.0 /
00994                             Ne5 ),
00995         interval( -8572019952114999.0 / Ne5,-8572019949515073.0 /
00996                             Ne5 ),
00997         interval( -8642406056614155.0 / Ne5,-8642406053936017.0 /
00998                             Ne5 ),
00999         interval( -8712791980171444.0 / Ne5,-8712791977516241.0 /
01000                             Ne5 ),
01001         interval( -8783177726114270.0 / Ne5,-8783177723474143.0 /
01002                             Ne5 ),
01003         interval( -8853563297747803.0 / Ne5,-8853563295082565.0 /
01004                             Ne5 ),
01005         interval( -8923948698211540.0 / Ne5,-8923948695603642.0 /
01006                             Ne5 ),
01007         interval( -8994333930651793.0 / Ne5,-8994333928013544.0 /
01008                             Ne5 ),
01009         interval( -4532359499005981.0 / Ne6,-4532359497669548.0 /
01010                             Ne6 ),
01011         interval( -4567551951590474.0 / Ne6,-4567551950297444.0 /
01012                             Ne6 ),
01013         interval( -4602744324609335.0 / Ne6,-4602744323253914.0 /
01014                             Ne6 ),
01015         interval( -4637936619312858.0 / Ne6,-4637936618035482.0 /
01016                             Ne6 ),
01017         interval( -4673128837201202.0 / Ne6,-4673128835890259.0 /
01018                             Ne6 ),
01019         interval( -4708320979539414.0 / Ne6,-4708320978220509.0 /
01020                             Ne6 ),
01021         interval( -4743513047586994.0 / Ne6,-4743513046289186.0 /
01022                             Ne6 ),
01023         interval( -4778705042649671.0 / Ne6,-4778705041357978.0 /
01024                             Ne6 ),
01025         interval( -4813896965939740.0 / Ne6,-4813896964638916.0 /
01026                             Ne6 ),
01027         interval( -4849088818685547.0 / Ne6,-4849088817383583.0 /
01028                             Ne6 ),
01029         interval( -4884280602023663.0 / Ne6,-4884280600723731.0 /
01030                             Ne6 ),
01031         interval( -4919472317096448.0 / Ne6,-4919472315799988.0 /
01032                             Ne6 ),
01033         interval( -4954663965077128.0 / Ne6,-4954663963749440.0 /
01034                             Ne6 ),
01035         interval( -4989855546965370.0 / Ne6,-4989855545671197.0 /
01036                             Ne6 ),
01037         interval( -5025047063935622.0 / Ne6,-5025047062637366.0 /
01038                             Ne6 ),
01039         interval( -5060238516963350.0 / Ne6,-5060238515663821.0 /
01040                             Ne6 ),
01041         interval( -5095429907060327.0 / Ne6,-5095429905761378.0 /
01042                             Ne6 ),
01043         interval( -5130621235251224.0 / Ne6,-5130621233930590.0 /
01044                             Ne6 ),
01045         interval( -5165812502482099.0 / Ne6,-5165812501185488.0 /
01046                             Ne6 ),
01047         interval( -5201003709724055.0 / Ne6,-5201003708430144.0 /
01048                             Ne6 ),
01049         interval( -5236194857910350.0 / Ne6,-5236194856593808.0 /
01050                             Ne6 ),
01051         interval( -5271385947936424.0 / Ne6,-5271385946609048.0 /
01052                             Ne6 ),
01053         interval( -5306576980681870.0 / Ne6,-5306576979376506.0 /
01054                             Ne6 ),
01055         interval( -5341767957039505.0 / Ne6,-5341767955734187.0 /
01056                             Ne6 ),
01057         interval( -5376958877810379.0 / Ne6,-5376958876528096.0 /
01058                             Ne6 ),
01059         interval( -5412149743940434.0 / Ne6,-5412149742614602.0 /
01060                             Ne6 ),
01061         interval( -5447340556059975.0 / Ne6,-5447340554811246.0 /
01062                             Ne6 ),
01063         interval( -5482531315163654.0 / Ne6,-5482531313857339.0 /
01064                             Ne6 ),
01065         interval( -5517722021905912.0 / Ne6,-5517722020635122.0 /
01066                             Ne6 ),
01067         interval( -5552912677100886.0 / Ne6,-5552912675823002.0 /
01068                             Ne6 ),
01069         interval( -5588103281480912.0 / Ne6,-5588103280171229.0 /
01070                             Ne6 ),
01071         interval( -5623293835745456.0 / Ne6,-5623293834493516.0 /
01072                             Ne6 ),
01073         interval( -5658484340702630.0 / Ne6,-5658484339399074.0 /
01074                             Ne6 ),
01075         interval( -5693674796958315.0 / Ne6,-5693674795683212.0 /
01076                             Ne6 ),
01077         interval( -5728865205222665.0 / Ne6,-5728865203973973.0 /
01078                             Ne6 ),
01079         interval( -5764055566229013.0 / Ne6,-5764055564966520.0 /
01080                             Ne6 ),
01081         interval( -5799245880597279.0 / Ne6,-5799245879305275.0 /
01082                             Ne6 ),
01083         interval( -5834436148937784.0 / Ne6,-5834436147684294.0 /
01084                             Ne6 ),
01085         interval( -5869626371953711.0 / Ne6,-5869626370690622.0 /
01086                             Ne6 ),
01087         interval( -5904816550239413.0 / Ne6,-5904816548965896.0 /
01088                             Ne6 ),
01089         interval( -5940006684383290.0 / Ne6,-5940006683093915.0 /
01090                             Ne6 ),
01091         interval( -5975196775000579.0 / Ne6,-5975196773734016.0 /
01092                             Ne6 ) };
01093 // Inclusion of the extremes of 1/gamma(x): 
01094 const interval gam_ryi[171] = { 
01095 pow2( interval( 5085347089749720.0 / Ne,5085347089749823.0 / Ne ) , 1 ) ,
01096 pow2( interval( -5082146609264467.0 / Ne,-5082146609264314.0 / Ne ) , -1 ) ,
01097 pow2( interval( 7824158147621733.0 / Ne,7824158147621966.0 / Ne ) , -1 ) ,
01098 pow2( interval( -5070842539852372.0 / Ne,-5070842539852221.0 / Ne ) , 1 ) ,
01099 pow2( interval( 4593118780547419.0 / Ne,4593118780547576.0 / Ne ) , 3 ) ,
01100 pow2( interval( -5333021955274733.0 / Ne,-5333021955274575.0 / Ne ) , 5 ) ,
01101 pow2( interval( 7546574203185105.0 / Ne,7546574203185319.0 / Ne ) , 7 ) ,
01102 pow2( interval( -6294628859031764.0 / Ne,-6294628859031469.0 / Ne ) , 10 ) ,
01103 pow2( interval( 6045310252810166.0 / Ne,6045310252811273.0 / Ne ) , 13 ) ,
01104 pow2( interval( -6568078652156336.0 / Ne,-6568078652156148.0 / Ne ) , 16 ) ,
01105 pow2( interval( 7963169065060572.0 / Ne,7963169065060801.0 / Ne ) , 19 ) ,
01106 pow2( interval( -5328217018030122.0 / Ne,-5328217018029960.0 / Ne ) , 23 ) ,
01107 pow2( interval( 7800142897041864.0 / Ne,7800142897042089.0 / Ne ) , 26 ) ,
01108 pow2( interval( -6199437664213474.0 / Ne,-6199437664213297.0 / Ne ) , 30 ) ,
01109 pow2( interval( 5316470282961123.0 / Ne,5316470282961284.0 / Ne ) , 34 ) ,
01110 pow2( interval( -4892929765135337.0 / Ne,-4892929765135165.0 / Ne ) , 38 ) ,
01111 pow2( interval( 4810107119289947.0 / Ne,4810107119290088.0 / Ne ) , 42 ) ,
01112 pow2( interval( -5030373421375086.0 / Ne,-5030373421374834.0 / Ne ) , 46 ) ,
01113 pow2( interval( 5576144001185310.0 / Ne,5576144001185479.0 / Ne ) , 50 ) ,
01114 pow2( interval( -6530685487420963.0 / Ne,-6530685487420774.0 / Ne ) , 54 ) ,
01115 pow2( interval( 8057940169576582.0 / Ne,8057940169576818.0 / Ne ) , 58 ) ,
01116 pow2( interval( -5223648494045513.0 / Ne,-5223648494045349.0 / Ne ) , 63 ) ,
01117 pow2( interval( 7099855957135674.0 / Ne,7099855957135885.0 / Ne ) , 67 ) ,
01118 pow2( interval( -5047359382236272.0 / Ne,-5047359382236084.0 / Ne ) , 72 ) ,
01119 pow2( interval( 7492585872478835.0 / Ne,7492585872479188.0 / Ne ) , 76 ) ,
01120 pow2( interval( -5795835662380422.0 / Ne,-5795835662380242.0 / Ne ) , 81 ) ,
01121 pow2( interval( 4664800910382651.0 / Ne,4664800910382790.0 / Ne ) , 86 ) ,
01122 pow2( interval( -7801058080117709.0 / Ne,-7801058080117472.0 / Ne ) , 90 ) ,
01123 pow2( interval( 6767162072327001.0 / Ne,6767162072327282.0 / Ne ) , 95 ) ,
01124 pow2( interval( -6082121514218736.0 / Ne,-6082121514218554.0 / Ne ) , 100 ) ,
01125 pow2( interval( 5656800000052189.0 / Ne,5656800000052359.0 / Ne ) , 105 ) ,
01126 pow2( interval( -5438268378952110.0 / Ne,-5438268378951951.0 / Ne ) , 110 ) ,
01127 pow2( interval( 5398375606367166.0 / Ne,5398375606367329.0 / Ne ) , 115 ) ,
01128 pow2( interval( -5527713447587841.0 / Ne,-5527713447587674.0 / Ne ) , 120 ) ,
01129 pow2( interval( 5833125895912623.0 / Ne,5833125895912799.0 / Ne ) , 125 ) ,
01130 pow2( interval( -6337936184674347.0 / Ne,-6337936184674153.0 / Ne ) , 130 ) ,
01131 pow2( interval( 7084743510515278.0 / Ne,7084743510515501.0 / Ne ) , 135 ) ,
01132 pow2( interval( -8141214882701327.0 / Ne,-8141214882701088.0 / Ne ) , 140 ) ,
01133 pow2( interval( 4804968547193877.0 / Ne,4804968547194018.0 / Ne ) , 146 ) ,
01134 pow2( interval( -5822137580509526.0 / Ne,-5822137580509355.0 / Ne ) , 151 ) ,
01135 pow2( interval( 7236772755227956.0 / Ne,7236772755228162.0 / Ne ) , 156 ) ,
01136 pow2( interval( -4610758665056508.0 / Ne,-4610758665056369.0 / Ne ) , 162 ) ,
01137 pow2( interval( 6019530845699084.0 / Ne,6019530845699266.0 / Ne ) , 167 ) ,
01138 pow2( interval( -8047036389398365.0 / Ne,-8047036389398123.0 / Ne ) , 172 ) ,
01139 pow2( interval( 5504580189086749.0 / Ne,5504580189086968.0 / Ne ) , 178 ) ,
01140 pow2( interval( -7703001513324420.0 / Ne,-7703001513324183.0 / Ne ) , 183 ) ,
01141 pow2( interval( 5510183009440391.0 / Ne,5510183009440581.0 / Ne ) , 189 ) ,
01142 pow2( interval( -8055535954952413.0 / Ne,-8055535954952173.0 / Ne ) , 194 ) ,
01143 pow2( interval( 6014315232803007.0 / Ne,6014315232803294.0 / Ne ) , 200 ) ,
01144 pow2( interval( -4584378555360492.0 / Ne,-4584378555360260.0 / Ne ) , 206 ) ,
01145 pow2( interval( 7132212380084113.0 / Ne,7132212380084326.0 / Ne ) , 211 ) ,
01146 pow2( interval( -5659549393054692.0 / Ne,-5659549393054526.0 / Ne ) , 217 ) ,
01147 pow2( interval( 4579461117155838.0 / Ne,4579461117155977.0 / Ne ) , 223 ) ,
01148 pow2( interval( -7554216840666713.0 / Ne,-7554216840666493.0 / Ne ) , 228 ) ,
01149 pow2( interval( 6348787715758027.0 / Ne,6348787715758222.0 / Ne ) , 234 ) ,
01150 pow2( interval( -5434979980476367.0 / Ne,-5434979980476204.0 / Ne ) , 240 ) ,
01151 pow2( interval( 4737681191908824.0 / Ne,4737681191908967.0 / Ne ) , 246 ) ,
01152 pow2( interval( -8407842664867513.0 / Ne,-8407842664867267.0 / Ne ) , 251 ) ,
01153 pow2( interval( 7592052521188700.0 / Ne,7592052521188935.0 / Ne ) , 257 ) ,
01154 pow2( interval( -6974119252551297.0 / Ne,-6974119252551090.0 / Ne ) , 263 ) ,
01155 pow2( interval( 6515520808385677.0 / Ne,6515520808385874.0 / Ne ) , 269 ) ,
01156 pow2( interval( -6188946869743481.0 / Ne,-6188946869743300.0 / Ne ) , 275 ) ,
01157 pow2( interval( 5975502808844840.0 / Ne,5975502808845020.0 / Ne ) , 281 ) ,
01158 pow2( interval( -5862842897072874.0 / Ne,-5862842897072704.0 / Ne ) , 287 ) ,
01159 pow2( interval( 5843967448508660.0 / Ne,5843967448508828.0 / Ne ) , 293 ) ,
01160 pow2( interval( -5916517001341501.0 / Ne,-5916517001341321.0 / Ne ) , 299 ) ,
01161 pow2( interval( 6082464626325325.0 / Ne,6082464626325503.0 / Ne ) , 305 ) ,
01162 pow2( interval( -6348157530347044.0 / Ne,-6348157530346858.0 / Ne ) , 311 ) ,
01163 pow2( interval( 6724699799057619.0 / Ne,6724699799057843.0 / Ne ) , 317 ) ,
01164 pow2( interval( -7228705737680202.0 / Ne,-7228705737679999.0 / Ne ) , 323 ) ,
01165 pow2( interval( 7883493269720206.0 / Ne,7883493269720561.0 / Ne ) , 329 ) ,
01166 pow2( interval( -8720834785364833.0 / Ne,-8720834785364561.0 / Ne ) , 335 ) ,
01167 pow2( interval( 4891722644546351.0 / Ne,4891722644546502.0 / Ne ) , 342 ) ,
01168 pow2( interval( -5564236710028970.0 / Ne,-5564236710028799.0 / Ne ) , 348 ) ,
01169 pow2( interval( 6416191129172903.0 / Ne,6416191129173091.0 / Ne ) , 354 ) ,
01170 pow2( interval( -7498890927628704.0 / Ne,-7498890927628487.0 / Ne ) , 360 ) ,
01171 pow2( interval( 8881515552460572.0 / Ne,8881515552460999.0 / Ne ) , 366 ) ,
01172 pow2( interval( -5328950915550370.0 / Ne,-5328950915550206.0 / Ne ) , 373 ) ,
01173 pow2( interval( 6478093314396794.0 / Ne,6478093314397089.0 / Ne ) , 379 ) ,
01174 pow2( interval( -7976303366065662.0 / Ne,-7976303366065426.0 / Ne ) , 385 ) ,
01175 pow2( interval( 4972846688449830.0 / Ne,4972846688450017.0 / Ne ) , 392 ) ,
01176 pow2( interval( -6278401907481090.0 / Ne,-6278401907480879.0 / Ne ) , 398 ) ,
01177 pow2( interval( 8024854758356088.0 / Ne,8024854758356345.0 / Ne ) , 404 ) ,
01178 pow2( interval( -5191277948909595.0 / Ne,-5191277948909444.0 / Ne ) , 411 ) ,
01179 pow2( interval( 6797621462551740.0 / Ne,6797621462551941.0 / Ne ) , 417 ) ,
01180 pow2( interval( -4503636668393666.0 / Ne,-4503636668393518.0 / Ne ) , 424 ) ,
01181 pow2( interval( 6037997262493341.0 / Ne,6037997262493523.0 / Ne ) , 430 ) ,
01182 pow2( interval( -8189485306115383.0 / Ne,-8189485306115130.0 / Ne ) , 436 ) ,
01183 pow2( interval( 5617805845426844.0 / Ne,5617805845427124.0 / Ne ) , 443 ) ,
01184 pow2( interval( -7795192616785187.0 / Ne,-7795192616784477.0 / Ne ) , 449 ) ,
01185 pow2( interval( 5469175405734180.0 / Ne,5469175405734422.0 / Ne ) , 456 ) ,
01186 pow2( interval( -7759929987383324.0 / Ne,-7759929987383086.0 / Ne ) , 462 ) ,
01187 pow2( interval( 5565727978288701.0 / Ne,5565727978288876.0 / Ne ) , 469 ) ,
01188 pow2( interval( -8070914994857895.0 / Ne,-8070914994857635.0 / Ne ) , 475 ) ,
01189 pow2( interval( 5914931467943193.0 / Ne,5914931467943373.0 / Ne ) , 482 ) ,
01190 pow2( interval( -8762204548045716.0 / Ne,-8762204548045455.0 / Ne ) , 488 ) ,
01191 pow2( interval( 6558513517606168.0 / Ne,6558513517606353.0 / Ne ) , 495 ) ,
01192 pow2( interval( -4960305627886271.0 / Ne,-4960305627886120.0 / Ne ) , 502 ) ,
01193 pow2( interval( 7580642983583672.0 / Ne,7580642983583897.0 / Ne ) , 508 ) ,
01194 pow2( interval( -5851844804194595.0 / Ne,-5851844804194367.0 / Ne ) , 515 ) ,
01195 pow2( interval( 4563038858728436.0 / Ne,4563038858728577.0 / Ne ) , 522 ) ,
01196 pow2( interval( -7187477492053316.0 / Ne,-7187477492052964.0 / Ne ) , 528 ) ,
01197 pow2( interval( 5716852908386950.0 / Ne,5716852908387214.0 / Ne ) , 535 ) ,
01198 pow2( interval( -4591808630269563.0 / Ne,-4591808630269411.0 / Ne ) , 542 ) ,
01199 pow2( interval( 7448102539955649.0 / Ne,7448102539955986.0 / Ne ) , 548 ) ,
01200 pow2( interval( -6098770429791387.0 / Ne,-6098770429791204.0 / Ne ) , 555 ) ,
01201 pow2( interval( 5041550443966798.0 / Ne,5041550443966946.0 / Ne ) , 562 ) ,
01202 pow2( interval( -8413996086583072.0 / Ne,-8413996086582821.0 / Ne ) , 568 ) ,
01203 pow2( interval( 7086939987269423.0 / Ne,7086939987269731.0 / Ne ) , 575 ) ,
01204 pow2( interval( -6024570065319942.0 / Ne,-6024570065319682.0 / Ne ) , 582 ) ,
01205 pow2( interval( 5168535487082451.0 / Ne,5168535487082609.0 / Ne ) , 589 ) ,
01206 pow2( interval( -8949051953781375.0 / Ne,-8949051953781115.0 / Ne ) , 595 ) ,
01207 pow2( interval( 7817344426895164.0 / Ne,7817344426895996.0 / Ne ) , 602 ) ,
01208 pow2( interval( -6889843867972878.0 / Ne,-6889843867972674.0 / Ne ) , 609 ) ,
01209 pow2( interval( 6126229646423302.0 / Ne,6126229646423484.0 / Ne ) , 616 ) ,
01210 pow2( interval( -5495122334906381.0 / Ne,-5495122334906222.0 / Ne ) , 623 ) ,
01211 pow2( interval( 4971972094727164.0 / Ne,4971972094727314.0 / Ne ) , 630 ) ,
01212 pow2( interval( -4537480959802395.0 / Ne,-4537480959802254.0 / Ne ) , 637 ) ,
01213 pow2( interval( 8352835047353300.0 / Ne,8352835047353555.0 / Ne ) , 643 ) ,
01214 pow2( interval( -7753443787904532.0 / Ne,-7753443787904298.0 / Ne ) , 650 ) ,
01215 pow2( interval( 7257653550749169.0 / Ne,7257653550749382.0 / Ne ) , 657 ) ,
01216 pow2( interval( -6850281165773769.0 / Ne,-6850281165773570.0 / Ne ) , 664 ) ,
01217 pow2( interval( 6519305845448896.0 / Ne,6519305845449168.0 / Ne ) , 671 ) ,
01218 pow2( interval( -6255266499085062.0 / Ne,-6255266499084872.0 / Ne ) , 678 ) ,
01219 pow2( interval( 6050802311308162.0 / Ne,6050802311308350.0 / Ne ) , 685 ) ,
01220 pow2( interval( -5900304762620398.0 / Ne,-5900304762620223.0 / Ne ) , 692 ) ,
01221 pow2( interval( 5799657649647993.0 / Ne,5799657649648165.0 / Ne ) , 699 ) ,
01222 pow2( interval( -5746047975553302.0 / Ne,-5746047975553134.0 / Ne ) , 706 ) ,
01223 pow2( interval( 5737835419331524.0 / Ne,5737835419331693.0 / Ne ) , 713 ) ,
01224 pow2( interval( -5774471890994117.0 / Ne,-5774471890993944.0 / Ne ) , 720 ) ,
01225 pow2( interval( 5856465763387432.0 / Ne,5856465763387600.0 / Ne ) , 727 ) ,
01226 pow2( interval( -5985387992102590.0 / Ne,-5985387992102406.0 / Ne ) , 734 ) ,
01227 pow2( interval( 6163919695584074.0 / Ne,6163919695584257.0 / Ne ) , 741 ) ,
01228 pow2( interval( -6395943042753787.0 / Ne,-6395943042753502.0 / Ne ) , 748 ) ,
01229 pow2( interval( 6686679647283150.0 / Ne,6686679647283350.0 / Ne ) , 755 ) ,
01230 pow2( interval( -7042883260256940.0 / Ne,-7042883260256730.0 / Ne ) , 762 ) ,
01231 pow2( interval( 7473096566380533.0 / Ne,7473096566380749.0 / Ne ) , 769 ) ,
01232 pow2( interval( -7987985534527481.0 / Ne,-7987985534527243.0 / Ne ) , 776 ) ,
01233 pow2( interval( 8600769311605383.0 / Ne,8600769311605633.0 / Ne ) , 783 ) ,
01234 pow2( interval( -4663884705694464.0 / Ne,-4663884705694325.0 / Ne ) , 791 ) ,
01235 pow2( interval( 5094554684614484.0 / Ne,5094554684614634.0 / Ne ) , 798 ) ,
01236 pow2( interval( -5604802840349871.0 / Ne,-5604802840349701.0 / Ne ) , 805 ) ,
01237 pow2( interval( 6209951739735886.0 / Ne,6209951739736072.0 / Ne ) , 812 ) ,
01238 pow2( interval( -6928963530888061.0 / Ne,-6928963530887851.0 / Ne ) , 819 ) ,
01239 pow2( interval( 7785368708274196.0 / Ne,7785368708274423.0 / Ne ) , 826 ) ,
01240 pow2( interval( -8808459126256481.0 / Ne,-8808459126256060.0 / Ne ) , 833 ) ,
01241 pow2( interval( 5017412797579486.0 / Ne,5017412797579638.0 / Ne ) , 841 ) ,
01242 pow2( interval( -5755173329981532.0 / Ne,-5755173329981361.0 / Ne ) , 848 ) ,
01243 pow2( interval( 6646385258439176.0 / Ne,6646385258439444.0 / Ne ) , 855 ) ,
01244 pow2( interval( -7727539896552529.0 / Ne,-7727539896552294.0 / Ne ) , 862 ) ,
01245 pow2( interval( 4522473425691912.0 / Ne,4522473425692052.0 / Ne ) , 870 ) ,
01246 pow2( interval( -5328812572761788.0 / Ne,-5328812572761623.0 / Ne ) , 877 ) ,
01247 pow2( interval( 6320558000502691.0 / Ne,6320558000502885.0 / Ne ) , 884 ) ,
01248 pow2( interval( -7546265781200776.0 / Ne,-7546265781200489.0 / Ne ) , 891 ) ,
01249 pow2( interval( 4534316912522546.0 / Ne,4534316912522688.0 / Ne ) , 899 ) ,
01250 pow2( interval( -5484491485989575.0 / Ne,-5484491485989407.0 / Ne ) , 906 ) ,
01251 pow2( interval( 6676632315202014.0 / Ne,6676632315202302.0 / Ne ) , 913 ) ,
01252 pow2( interval( -8180074398476253.0 / Ne,-8180074398476014.0 / Ne ) , 920 ) ,
01253 pow2( interval( 5042989707083422.0 / Ne,5042989707083666.0 / Ne ) , 928 ) ,
01254 pow2( interval( -6257379418480333.0 / Ne,-6257379418480019.0 / Ne ) , 935 ) ,
01255 pow2( interval( 7813097673618694.0 / Ne,7813097673619043.0 / Ne ) , 942 ) ,
01256 pow2( interval( -4908325621754370.0 / Ne,-4908325621754217.0 / Ne ) , 950 ) ,
01257 pow2( interval( 6205346227418363.0 / Ne,6205346227418597.0 / Ne ) , 957 ) ,
01258 pow2( interval( -7893590972392525.0 / Ne,-7893590972392227.0 / Ne ) , 964 ) ,
01259 pow2( interval( 5051411882876310.0 / Ne,5051411882876506.0 / Ne ) , 972 ) ,
01260 pow2( interval( -6504655602059905.0 / Ne,-6504655602059583.0 / Ne ) , 979 ) ,
01261 pow2( interval( 8426810051054742.0 / Ne,8426810051054986.0 / Ne ) , 986 ) ,
01262 pow2( interval( -5491407534973626.0 / Ne,-5491407534973452.0 / Ne ) , 994 ) ,
01263 pow2( interval( 7199960218142557.0 / Ne,7199960218142768.0 / Ne ) , 1001 ) ,
01264 pow2( interval( -4748178637044143.0 / Ne,-4748178637044000.0 / Ne ) , 1009 ) ,
01265 pow2( interval( 6299691458188962.0 / Ne,6299691458189149.0 / Ne ) , 1016 ) };
01266 
01267 // ****************************************************************************
01268 // ******************** Gamma(x), 1/Gamma(x) **********************************
01269 // ****************************************************************************
01270 
01271 inline int round_g(const real& x) throw() 
01272 // Only for the internal use ( interval gammar(x) )
01273 // For |x| < 2147483647.5 the assignment  y = round_g(x)  delivers:
01274 // y = round_g(-0.1);            --->  y = 1;
01275 // y = round_g(-0.5);            --->  y = 1;
01276 // y = round_g(-1.0);            --->  y = 1;
01277 // y = round_g(-1.4);            --->  y = 2;           
01278 // y = round_g(-6.8);            --->  y = 7;
01279 // y = round_g(+0.0);            --->  y = 0;
01280 // y = round_g(+0.1);            --->  y = 0;
01281 // y = round_g(+0.5);            --->  y = 0;
01282 // y = round_g(+2.6);            --->  y = 0;
01283 // Blomquist, 25.06.2009;
01284 {
01285         int n = ifloor(_double(x));
01286         n = (n>=0)? 0:-n;
01287         
01288         return n;
01289 }
01290 
01291 real gamr_even_Ma(const real& x1, const real& x2, int n1)
01292 {
01293         real y;
01294         
01295         if ( x2<Inf(gam_rxi[n1]) || Sup(gam_rxi[n1])<x1 ) 
01296         {  // [x1,x2] & gam_rxi[n1] is the empty set:
01297                 y = (x1<Inf(gam_rxi[n1]))? gammar(x2) : gammar(x1);
01298                 y *= q_gammarp;
01299         }
01300         else // [x1,x2] & gam_rxi[n1] is not the empty set:
01301                 y = Sup(gam_ryi[n1]); 
01302         
01303         return y;
01304 }
01305 
01306 real gamr_even_Mi(const real& x1, const real& x2, int n1)
01307 {
01308         real y,y1;
01309         
01310         if ( x2<Inf(gam_rxi[n1]) || Sup(gam_rxi[n1])<x1 ) 
01311         {  // [x1,x2] & gam_rxi[n1] is the empty set:
01312                 std::cout << "Leere Menge:" << std::endl;
01313                 y = (x1<Inf(gam_rxi[n1]))? gammar(x1) : gammar(x2);
01314                 y *= q_gammarm;
01315         }
01316         else // [x1,x2] & gam_rxi[n1] is not the empty set:
01317         {
01318                 y1 = gammar(x1)*q_gammarm;
01319                 y  = gammar(x2)*q_gammarm;
01320                 if (y1<y) y=y1;
01321         }
01322         
01323         return y;
01324 }
01325 
01326 real gamr_odd_Mi(const real& x1, const real& x2, int n1)
01327 {
01328         real y;
01329         
01330         if ( x2<Inf(gam_rxi[n1]) || Sup(gam_rxi[n1])<x1 ) 
01331         {  // [x1,x2] & gam_rxi[n1] is the empty set:
01332                 y = (x1<Inf(gam_rxi[n1]))? gammar(x2) : gammar(x1);
01333                 y *= q_gammarp;
01334         }
01335         else // [x1,x2] & gam_rxi[n1] is not the empty set:
01336                 y = Inf(gam_ryi[n1]); 
01337         
01338         return y;
01339 }
01340 
01341 real gamr_odd_Ma(const real& x1, const real& x2, int n1)
01342 {
01343         real y,y1;
01344         
01345         if ( x2<Inf(gam_rxi[n1]) || Sup(gam_rxi[n1])<x1 ) 
01346         {  // [x1,x2] & gam_rxi[n1] is the empty set:
01347                 std::cout << "Leere Menge:" << std::endl;
01348                 y = (x1<Inf(gam_rxi[n1]))? gammar(x1) : gammar(x2);
01349                 y *= q_gammarm;
01350         }
01351         else // [x1,x2] & gam_rxi[n1] is not the empty set:
01352         {
01353                 y1 = gammar(x1)*q_gammarm;
01354                 y  = gammar(x2)*q_gammarm;
01355                 if (y1>y) y=y1;
01356         }
01357         
01358         return y;
01359 }
01360 
01361 interval gammar(const interval& x)
01362 // Calculating inclusions of 1/Gamma(x);
01363 // Blomquist, 01.07.09
01364 {
01365         interval y;
01366         real x0, x1(Inf(x)), x2(Sup(x)), y0,y1(0),y2;
01367         int n1,n2;
01368         
01369         n1 = round_g(x1); 
01370         n2 = round_g(x2);
01371         if (x1==x2) // x is point interval
01372                 if (x1==-n1) y2 = y1; // y = [0,0];
01373         else 
01374         {
01375                 y1 = gammar(x1);
01376                 y2 = y1;
01377                         // Lower bound y1, Upper bound y2:
01378                 if (y1<0) 
01379                 {
01380                         y1 = y1*q_gammarp;
01381                         y2 = y2*q_gammarm;
01382                 }
01383                 else 
01384                 {
01385                         y1 = y1*q_gammarm;
01386                         y2 = y2*q_gammarp;
01387                 }
01388         }
01389         else // x2>x1:
01390         {
01391                 if (n1%2==0) // n1 even, i.e. n1=0,2,4,6,...
01392                 {
01393                         if (n1==n2)
01394                         {
01395                                 y1 = gamr_even_Mi(x1,x2,n1);
01396                                 y2 = gamr_even_Ma(x1,x2,n1);
01397                         }
01398                         else
01399                                 if (n2==n1-1)
01400                         {
01401                                 y1 = gamr_odd_Mi(-n2,x2,n2);
01402                                 y2 = gamr_even_Ma(x1,-n2,n1);
01403                         }
01404                         else // n2 <= n1-2
01405                         {       
01406                                 y1 = Inf(gam_ryi[n1-1]);  // Minimum OK
01407                                 y2 = gamr_even_Ma(x1,-n1+1,n1);
01408                                 x0 = x2;
01409                                 if (x2>n1-3 && x2<0) 
01410                                         x0 = n1-3;
01411                                 y0 = gamr_even_Ma(-n1+2,x0,n1-2);
01412                                 if (y0>y2) y2 = y0;
01413                                         
01414                                 if (n1==4 && n2==0)
01415                                 {
01416                                         y0 = gamr_even_Ma(0,x2,0);
01417                                         if (y0>y2) y2=y0;
01418                                 }
01419                         }
01420                 } // n1 even;
01421                 else // n1 odd:
01422                 {
01423                         if (n1==n2)
01424                         {
01425                                 y1 = gamr_odd_Mi(x1,x2,n1);
01426                                 y2 = gamr_odd_Ma(x1,x2,n1);
01427                         }
01428                         else
01429                                 if (n2==n1-1)
01430                         {
01431                                 y1 = gamr_odd_Mi(x1,-n2,n1);
01432                                 y2 = gamr_even_Ma(-n2,x2,n2);
01433                         }
01434                         else
01435                                 if (n2==n1-2)
01436                         {
01437                                 y1 = gamr_odd_Mi(x1,n1-1,n1);
01438                                 y2 = gamr_odd_Mi(-n1+2,x2,n1-2);
01439                                 if (y2<y1) y1 = y2; // Minimum calculated
01440                                 y2 = Sup( gam_ryi[n1-1] );
01441                         }
01442                         else // 0 <= n2 <= n1-3;
01443                         {  // Calculating the minimum y1:
01444                                 y1 = gamr_odd_Mi(x1,n1-1,n1);
01445                                 y2 = Inf( gam_ryi[n1-2] );
01446                                 if (y2<y1) y1 = y2; // Minimum y1 calculated
01447                                                 // Now calculating the maximum y2:
01448                                 if (n1==3)
01449                                 {
01450                                         y2 = Sup( gam_ryi[n1-1]);
01451                                         y0 = gamr_even_Ma(0,x2,0);
01452                                         if (y0>y2) y2=y0;
01453                                 }
01454                                 else // n1 = 5,7,9,....
01455                                         y2 = Sup( gam_ryi[n1-1] );
01456                         }
01457                                                 
01458                 } // n1 odd
01459         }
01460         
01461         y = interval(y1,y2);    
01462         return y;
01463 }
01464 
01465 interval gamma(const interval& x)
01466 // Calculating inclusions of 1/Gamma(x);
01467 // Blomquist, 01.07.09
01468 {
01469         return 1/gammar(x);
01470 }
01471 
01472 
01473 } // namespace cxsc
01474