C-XSC - A C++ Class Library for Extended Scientific Computing
2.5.4
|
00001 /* 00002 ** CXSC is a C++ library for eXtended Scientific Computing (V 2.5.4) 00003 ** 00004 ** Copyright (C) 1990-2000 Institut fuer Angewandte Mathematik, 00005 ** Universitaet Karlsruhe, Germany 00006 ** (C) 2000-2014 Wiss. Rechnen/Softwaretechnologie 00007 ** Universitaet Wuppertal, Germany 00008 ** 00009 ** This library is free software; you can redistribute it and/or 00010 ** modify it under the terms of the GNU Library General Public 00011 ** License as published by the Free Software Foundation; either 00012 ** version 2 of the License, or (at your option) any later version. 00013 ** 00014 ** This library is distributed in the hope that it will be useful, 00015 ** but WITHOUT ANY WARRANTY; without even the implied warranty of 00016 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00017 ** Library General Public License for more details. 00018 ** 00019 ** You should have received a copy of the GNU Library General Public 00020 ** License along with this library; if not, write to the Free 00021 ** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00022 */ 00023 00024 /* CVS $Id: 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