C-XSC - A C++ Class Library for Extended Scientific Computing  2.5.4
imath.cpp
1 /*
2 ** CXSC is a C++ library for eXtended Scientific Computing (V 2.5.4)
3 **
4 ** Copyright (C) 1990-2000 Institut fuer Angewandte Mathematik,
5 ** Universitaet Karlsruhe, Germany
6 ** (C) 2000-2014 Wiss. Rechnen/Softwaretechnologie
7 ** Universitaet Wuppertal, Germany
8 **
9 ** This library is free software; you can redistribute it and/or
10 ** modify it under the terms of the GNU Library General Public
11 ** License as published by the Free Software Foundation; either
12 ** version 2 of the License, or (at your option) any later version.
13 **
14 ** This library is distributed in the hope that it will be useful,
15 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
16 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 ** Library General Public License for more details.
18 **
19 ** You should have received a copy of the GNU Library General Public
20 ** License along with this library; if not, write to the Free
21 ** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
22 */
23 
24 /* CVS $Id: imath.cpp,v 1.43 2014/01/30 17:23:45 cxsc Exp $ */
25 
26 #include "imath.hpp"
27 #include "rmath.hpp"
28 
29 // Auch fi_lib hat ein eigenes interval (wie RTS sein dotprecision)
30 //typedef struct fi_interval { double INF, SUP;} fi_interval;
31 #undef LINT_ARGS
32 #define CXSC_INCLUDE
33 #include <fi_lib.hpp>
34 
35 extern "C" {
36 #ifndef rfcth_included
37 #define rfcth_included
38 #include "r_fcth.h"
39 #endif
40 }
41 
42 namespace cxsc {
43 
44 
45 interval sqr (const interval &a) noexcept
46 {
47  interval res;
48  res= a*a;
49  if (Inf(res)<0) Inf(res)=0;
50  return res;
51 }
52 
53 interval sqrt (const interval &a, int n)
54 {
55  if ( ((n>0) && (Inf(a)>=0.0)) || ((n<0) && (Inf(a)>0.0)) )
56  return pow(a,interval(1.0,1.0)/n);
57  else {
58  cxscthrow(STD_FKT_OUT_OF_DEF("interval sqrt (const interval &a, int n)"));
59  return interval(-1.0); // dummy result
60  }
61 }
62 
63 interval sqrt1px2(const interval& x) noexcept
64 // Inclusion of sqrt(1+x^2); Blomquist, 13.12.02;
65 {
66  interval t = abs(x),y;
67  if (expo(Inf(t)) > 33)
68  {
69  y = t;
70  Sup(y) = succ(Sup(y));
71  } else if (expo(Sup(t)) > 33)
72  {
73  y = interval(Inf(t)); // interval(Inf(t)) is a point interval!
74  y = sqrt(1+y*y); // ---> y*y == sqr(y);
75  y = interval(Inf(y),succ(Sup(t)));
76  } else y = sqrt(1+sqr(t));
77  return y;
78 }
79 
80 interval sqrtx2y2(const interval& x, const interval& y) noexcept
81 // Inclusion of sqrt(x^2+y^2); Blomquist, 13.12.02;
82 {
83  interval a=abs(x), b=abs(y), r;
84  int exa=expo(Sup(a)), exb=expo(Sup(b)), ex;
85  if (exb > exa)
86  { // Permutation of a,b:
87  r = a; a = b; b = r;
88  ex = exa; exa = exb; exb = ex;
89  }
90  ex = 511 - exa;
91  times2pown(a,ex);
92  times2pown(b,ex);
93  r = sqrt(a*a + b*b);
94  times2pown(r,-ex);
95  return r;
96 } // sqrtx2y2
97 
98 //********************************************************************
99 // Constants for: interval sqrtp1m1(const interval& x) noexcept
100 // Blomquist 05.08.03
101 const real Delta_f = 2*minreal;
102 const real q_sqrtp1m1m = 9007199254740984.0 / 9007199254740992.0;
103 const real q_sqrtp1m1p = 4503599627370502.0 / 4503599627370496.0;
104 //********************************************************************
105 interval sqrtp1m1(const interval& x) noexcept
106 // interval(a,b) is an inclusion of sqrt(x+1)-1;
107 // Exported by imath.hpp; Blomquist, 05.08.03;
108 {
109  real a=0,b=0,ix=Inf(x),sx=Sup(x);
110  int ex_ix,ex_sx,sgn_ix,sgn_sx;
111 
112  // Calculating the lower bound:
113  ex_ix = expo(ix); sgn_ix = sign(ix);
114  if (ex_ix<=-1021) // ==> |ix| < 2^(-1021)
115  { if (sgn_ix) a = sqrtp1m1(ix) - Delta_f; }
116  else // |ix| >= 2^(-1021)
117  if (ex_ix<=-51)
118  {
119  times2pown(ix,-1); // exact division by 2
120  a = pred(ix);
121  }
122  else
123  if (sgn_ix>0) a = (ix>0.67) ?
124  Inf( sqrt(interval(ix)+1)-1 ) : sqrtp1m1(ix)*q_sqrtp1m1m;
125  else a = (ix<-0.25) ?
126  Inf( sqrt(interval(ix)+1)-1 ) : sqrtp1m1(ix)*q_sqrtp1m1p;
127  // Calculating the upper bound:
128  if (ix == sx) { ex_sx = ex_ix; sgn_sx = sgn_ix; }
129  else { ex_sx = expo(sx); sgn_sx = sign(sx); }
130  if (ex_sx<=-1021) // ==> |sx| < 2^(-1021)
131  { if (sgn_sx) b = sqrtp1m1(sx) + Delta_f; }
132  else // |sx| >= 2^(-1021)
133  if (ex_sx<=-47) { b = sx; times2pown(b,-1); }
134  else
135  if (sgn_sx>0) b = (sx>0.58) ?
136  Sup( sqrt(interval(sx)+1)-1 ) : sqrtp1m1(sx)*q_sqrtp1m1p;
137  else b = (sx<-0.32) ?
138  Sup( sqrt(interval(sx)+1)-1 ) : sqrtp1m1(sx)*q_sqrtp1m1m;
139  return interval(a,b);
140 } // sqrtp1m1()
141 
142 //********************************************************************
143 // Konstanten fr die Intervallfunktion sqrt(x^2-1):
144 real q_sqrtx2m1p(4503599627370501.0/4503599627370496.0); // (1+e(f))
145 real q_sqrtx2m1m(9007199254740986.0/9007199254740992.0); // (1-e(f))
146 //********************************************************************
147 
149 // sqrt(x^2-1); Blomquist, 13.04.04;
150 {
151  interval z = abs(x);
152  real r1,r2;
153  r1 = sqrtx2m1(Inf(z)) * q_sqrtx2m1m;
154  r2 = Sup(z);
155  if (expo(r2)<26)
156  r2 = sqrtx2m1(r2) * q_sqrtx2m1p;
157  // expo(r2) >= 26 --> r2 = Sup(z) ist optimale Oberschranke!
158  return interval(r1,r2);
159 } // sqrtx2m1 (Intervallargumente)
160 
161 
162 //********************************************************************
163 // Konstanten fr die Intervallfunktion sqrt(1-x2):
164 real q_sqrt1mx2p(4503599627370501.0/4503599627370496.0); // (1+e(f))
165 real q_sqrt1mx2m(9007199254740985.0/9007199254740992.0); // (1-e(f))
166 //********************************************************************
167 
169 // sqrt(1-x2); Blomquist, 13.04.04;
170 {
171  interval z = abs(x); // sqrt(1-t2) monoton fallend in t aus z;
172  real r1,r2,sz,iz;
173  sz = Sup(z); iz = Inf(z);
174  // Berechnung der Unterschranke r1:
175  r2 = sqrt1mx2(sz);
176  r1 = (sz==0)? 1 : r2 * q_sqrt1mx2m;
177  // Berechnung der Oberrschranke r2:
178  if (iz<4.81e-8) r2 = 1; // r2=1 ist immer korrekte Oberschranke!
179  else r2 = (sz==iz)? r2*q_sqrt1mx2p : sqrt1mx2(iz)*q_sqrt1mx2p;
180  return interval(r1,r2);
181 
182 } // sqrtx2m1 (Intervallargumente)
183 
184 
185 // Konstanten fuer die Intervallfunktion exp(-x^2):
186 const real q_exp_x2p = 4503599627370502.0 / 4503599627370496.0; // (1+e(f))
187 const real q_exp_x2m = 9007199254740984.0 / 9007199254740992.0; // (1-e(f))
188  // Oberschranke expmx2_UB fuer |x| > expmx2_x0:
189 const real expmx2_UB = 2.225073858507447856659E-308;
190 const real expmx2_x0 = 7491658466053896.0 / 281474976710656.0;
191 
193 // e^(-x^2); Blomquist, 05.07.04;
194 {
195  real y,r1,r2,Sz,Iz;
196  interval z = abs(x);
197  // Berechnung einer Unterschranke:
198  Sz = Sup(z); Iz = Inf(z);
199  y = expmx2(Sz);
200  r1 = y;
201  if (Sz<=expmx2_x0) r1 = r1 * q_exp_x2m; // Abrunden
202  if (Sz==0) r1 = y;
203  // Berechnung einer Oberschranke:
204  if (Iz>expmx2_x0) r2 = expmx2_UB;
205  else r2 = (Sz==Iz)? y * q_exp_x2p : expmx2(Iz) * q_exp_x2p;
206  if (r2>1) r2 = 1.0;
207  return interval(r1,r2);
208 
209 } // expmx2 (Intervallargumente)
210 
211 
213 // Interval function of exp(x)-1
214 {
215  real y,r1,r2,Sx,Ix;
216  Sx = Sup(x); Ix = Inf(x);
217  y = expm1(Ix);
218  // calculation of a lower bound r1:
219  r1 = (y>0)? y*q_exmm : y*q_exmp; // rounding downwards
220  if (r1<-1) r1 = -1;
221  // calculation of an upper bound r2:
222  if (Sx!=Ix) y = expm1(Sx);
223  r2 = (y>0)? y*q_exmp : y*q_exmm;
224 
225  return interval(r1,r2);
226 } // expm1(...)
227 
228 // Die folgenden Konstanten q_expx2_p, q_expx2_m beziehen sich
229 // auf eps(expx2) = 4.618958e-16; f(x) = e^{+x^2};
230 
231 const real q_expx2_p = 4503599627370500.0 / 4503599627370496.0; // (1+e(f))
232 const real q_expx2_m = 9007199254740985.0 / 9007199254740992.0; // (1-e(f))
233 
235 // e^(+x^2); Blomquist, 25.07.06;
236 {
237  real y,r1,r2,Sz,Iz;
238  interval z = abs(x);
239  Sz = Sup(z); Iz = Inf(z);
240  // Berechnung einer Unterschranke:
241  y = expx2(Iz);
242  r1 = y;
243  r1 = r1 * q_expx2_m; // Abrunden
244  if (r1<1.0) r1 = 1.0;
245  // Berechnung einer Oberschranke:
246  r2 = (Sz==Iz)? y * q_expx2_p : expx2(Sz) * q_expx2_p;
247  if (Sz==0) r2 = 1.0;
248 
249  return interval(r1,r2);
250 } // expx2 (Intervallfunktion)
251 
252 // Die folgenden Konstanten q_expx2m1_p, q_expx2m1_m beziehen sich
253 // auf eps(expx2m1) = 4.813220e-16; f(x) = e^{+x^2}-1;
254 
255 const real q_expx2m1_p = 4503599627370500.0 / 4503599627370496.0; // (1+e(f))
256 const real q_expx2m1_m = 9007199254740985.0 / 9007199254740992.0; // (1-e(f))
257 const real expx2m1_0 = comp(0.5,-510); // 2^(-511)
258 
259 void sqr2uv(const real&, real&, real&);
260 
261 real expx2m1_intv(const real& x)
262 // Zur Implementierung der Intervallfunktion;
263 // e^(+x^2)-1; rel. Fehlerschranke: eps = 4.813220E-16 = e(f) gilt
264 // fuer alle x, mit: 2^(-511) <= x <= x0 = 26.64174755704632....
265 // x0 = 7498985273150791.0 / 281474976710656.0;
266 // Fuer x > x0 --> Programmabbruch wegen Overflow;
267 // Fuer 0 <= x < 2^(-511) werden die Funktionswerte auf Null gesetzt.
268 // Ausfuehrlich getestet; Blomquist, 10.08.2006;
269 {
270  real t(x),u,v,y,res(0);
271  int ex;
272  if (t<0) t = -t; // t >= 0;
273 
274  if (t>=6.5) res = expx2(t);
275  else
276  {
277  ex = expo(t);
278  sqr2uv(x,u,v); // u := x*x und v aus S(2,53);
279  if (ex>=2) // g4(x)
280  {
281  y = exp(u);
282  res = 1 - v*y;
283  res = y - res;
284  }
285  else
286  if (ex>=-8) res = expm1(u) + v*exp(u); // g3(x)
287  else
288  if (ex>=-25) { // g2(x)
289  y = u*u;
290  times2pown(y,-1);
291  res = (1+u/3)*y + u;
292  }
293  else
294  if(ex>=-510) res = u; // g1(x)
295  }
296 
297  return res;
298 } // expx2m1_intv
299 
301 // e^(+x^2)-1; Blomquist, 10.08.06;
302 {
303  real y,r1,r2,Sz,Iz;
304  interval z = abs(x);
305  Sz = Sup(z); Iz = Inf(z);
306  // Berechnung einer Unterschranke:
307  y = expx2m1_intv(Iz);
308  r1 = y;
309  r1 = r1 * q_expx2m1_m; // Abrunden
310  // Berechnung einer Oberschranke:
311  if (Sz < expx2m1_0)
312  {
313  r2 = MinReal;
314  if (Sz==0) r2 = 0.0;
315  }
316  else r2 = (Sz==Iz)? y * q_expx2m1_p : expx2m1_intv(Sz) * q_expx2m1_p;
317 
318  return interval(r1,r2);
319 } // expx2m1 (Intervallfunktion)
320 
321 // ------ 1-eps and 1+eps for function lnp1, Blomquist 28,07,03; -----------
322 // ------------------------ eps = 2.5082e-16 -------------------------------
323 static real q_lnp1m = 9007199254740986.0 / 9007199254740992.0; // 1-eps
324 static real q_lnp1p = 4503599627370501.0 / 4503599627370496.0; // 1+eps
325 // ---------------------------------------------------------------------------
326 
327 interval lnp1(const interval& x) noexcept
328 // returns an inclusion of ln(1+t), with t in x; Blomquist 28.07.03;
329 {
330  real ix=Inf(x), sx=Sup(x),a,b; // ln(1+x) <= [a,b]
331  // Calculating the lower bound a:
332  int sgn_ix = sign(ix), ex_ix = expo(ix), sgn_sx,ex_sx;
333  if (!sgn_ix) a = 0; // optimal lower bound!
334  else if (sgn_ix>0)
335  a = (ex_ix<=-53) ? pred(ix) : lnp1(ix) * q_lnp1m;
336  else
337  a = (ex_ix<=-54) ? pred(ix) : lnp1(ix) * q_lnp1p;
338  // Calculating the upper bound b:
339  if (ix == sx) { sgn_sx = sgn_ix; ex_sx = ex_ix; }
340  else { sgn_sx = sign(sx); ex_sx = expo(sx); }
341  if (sgn_sx>=0)
342  b = (ex_sx<=-49) ? sx : lnp1(sx) * q_lnp1p;
343  else // sx < 0:
344  b = (ex_sx<=-50) ? sx : lnp1(sx) * q_lnp1m;
345  return interval(a,b); // ln(1+x) in [a,b]
346 } // lnp1
347 
354 interval erf (const interval &a) { return j_erf(a); }
361 interval erfc (const interval &a) { return j_erfc(a); }
362 
363 //interval pow (const interval &a, const interval &b)
364 //{
365 // if(Inf(a)>0)
366 // return j_exp(b*ln(a));
367 // else if(Inf(a)==0 && Inf(b)>=0)
368 // {
369 // if(Sup(a)>=1)
370 // return interval(0,pow(Sup(a),Sup(b)));
371 // else
372 // return interval(0,pow(Sup(a),Inf(b)));
373 // }
374 // else
375 // {
376 // cxscthrow(ERROR_INTERVAL_STD_FKT_OUT_OF_DEF("interval pow(const interval &,const interval &)"));
377 // return interval(0,0);
378 // }
379 //}
380 
381 inline a_intv _a_intv(const interval &x)
382 {
383  return *((const a_intv *)(&x));
384 }
385 inline interval _interval(const a_intv &x)
386 {
387  return *((const interval *)(&x));
388 }
389 
390 interval pow (const interval &a, const interval &b) noexcept
391  {
392  interval res;
393  if(Inf(a)==0 && Inf(b)>=0)
394  {
395  if(Sup(a)>=1)
396  res=interval(0,pow(Sup(a),Sup(b)));
397  else
398  res=interval(0,pow(Sup(a),Inf(b)));
399  }
400  else res = _interval(i_pow(_a_intv(a),_a_intv(b)));
401 
402  if (Inf(res) <= Sup(res))
403  return res;
404  else
405  return interval(Sup(res),Inf(res));
406  }
407 
408 //----------------------------------------------------------------------------
409 // Purpose: The local function 'Power()' is used to compute a lower or an
410 // upper bound for the power function with real argument and integer
411 // exponent, respectively.
412 // Parameters:
413 // In: 'x' : real argument
414 // 'n' : integer exponent
415 // 'RndMode': rounding mode,
416 // (-1 = downwardly directed, +1 = upwardly directed)
417 // Description:
418 // This function is used to speed up the interval power function defined
419 // below. The exponentiation is reduced to multiplications using the
420 // binary shift method. Depending on 'n', this function is up to 40 times
421 // as fast as the standard power function for real argument and real
422 // exponent. However, its accuracy is less than one ulp (unit in the last
423 // place of the mantissa) since about log2(n) multiplications are executed
424 // during computation. Since directed roundings are antisymmetric, one
425 // gets
426 //
427 // down(x^n) = -up((-x)^n) and up(x^n) = -down((-x)^n)
428 //
429 // for x < 0 and odd n, where 'down' and 'up' denote the downwardly and
430 // upwardly directed roundings, respectively.
431 //----------------------------------------------------------------------------
432 static real Power (const real & x, int n, int RndMode )
433 { // Signals change of the rounding mode
434  int ChangeRndMode; // for x < 0 and odd n
435  real p, z;
436 
437  ChangeRndMode = ( (x < 0.0) && (n % 2 == 1) );
438 
439  if (ChangeRndMode) {
440  z = -x;
441  RndMode = -RndMode;
442  }
443  else
444  z = x;
445 
446  p = 1.0;
447  switch (RndMode) { // Separate while-loops used
448  case -1 : while (n > 0) { // to gain speed at runtime
449  if (n % 2 == 1) p = muld(p,z); //--------------------------
450  n = n / 2;
451  if (n > 0) z = muld(z,z);
452  }
453  break;
454  case +1 : while (n > 0) {
455  if (n % 2 == 1) p = mulu(p,z);
456  n = n / 2;
457  if (n > 0) z = mulu(z,z);
458  }
459  break;
460  }
461 
462  if (ChangeRndMode)
463  return -p;
464  else
465  return p;
466 }
467 
468 interval power(const interval& a, int n)
469 // Calculating a^n;
470 // Examples: [-1,4]^2 = [0,16]; [-1,4]^3 = [-16,+64];
471 {
472  bool neg(n<0);
473  int N(n); //,k(-1),r;
474  interval res,h;
475  real Lower, Upper;
476  if (neg) N = -N;
477  if (N==0) res = 1;
478  else if (N==1) res = a;
479  else // N > 1:
480  if (Inf(a)>=MinReal) res = exp(N*ln(a));
481  else if (Sup(a)<=-MinReal) {
482  h = interval(-Sup(a),-Inf(a));
483  res = exp(N*ln(h));
484  if (N%2 == 1) res = -res;
485  }
486  else
487  {
488 // h = a;
489 // while(N>0)
490 // {
491 // k++;
492 // r = N % 2;
493 // if (k==0)
494 // if (r==1) res=a; else res=1;
495 // else {
496 // h = sqr(h);
497 // if (r==1) res *= h;
498 // }
499 // N = N / 2;
500 // }
501  if ( (0.0 < Inf(a)) || (N % 2 == 1) ) {
502  Lower = Power(Inf(a),N,-1);
503  Upper = Power(Sup(a),N,+1);
504  }
505  else if (0.0 > Sup(a)) {
506  Lower = Power(Sup(a),N,-1);
507  Upper = Power(Inf(a),N,+1);
508  }
509  else {
510  Lower = 0.0;
511  Upper = Power(AbsMax(a),N,+1);
512  }
513  res = interval(Lower,Upper);
514  }
515  if (neg) res = 1/res;
516  return res;
517 } // power
518 
519 
520 //----------------------------------------------------------------------------
521 // Purpose: This version of the function 'Power()' is used to compute an
522 // enclosure for the power function with interval argument and integer
523 // exponent.
524 // Parameters:
525 // In: 'x': interval argument
526 // 'n': integer exponent
527 // Description:
528 // In general, this implementation does not deliver a result of maximum
529 // accuracy, but it is about 30-40 times faster than the standard power
530 // function for interval arguments and interval exponents. The resulting
531 // interval has a width of approximately 2*log2(n) ulps. Since x^n is
532 // considered as a monomial, we define x^0 := 1. For negative exponents
533 // and 0 in 'x', the division at the end of the function will cause a
534 // runtime error (division by zero).
535 //----------------------------------------------------------------------------
536 interval Power (const interval & x, int n )
537 {
538  int m;
539  real Lower, Upper;
540 
541  if (n == 0) return(interval(1.0,1.0));
542 
543  if (n > 0) m = n; else m = -n;
544 
545  if ( (0.0 < Inf(x)) || (m % 2 == 1) ) {
546  Lower = Power(Inf(x),m,-1);
547  Upper = Power(Sup(x),m,+1);
548  }
549  else if (0.0 > Sup(x)) {
550  Lower = Power(Sup(x),m,-1);
551  Upper = Power(Inf(x),m,+1);
552  }
553  else {
554  Lower = 0.0;
555  Upper = Power(AbsMax(x),m,+1);
556  }
557 
558  if (n > 0)
559  return(interval(Lower,Upper));
560  else // Causes a runtime error
561  return(1.0/interval(Lower,Upper)); // if 0 in 'x'.
562 }
563 
564 //----------------------------------------------------------------------------
565 
566 interval Pi ( ) // Enclosure of constant pi
567  { return(4.0*atan(_interval(1.0,1.0))); } //-------------------------
568 
569 // Error bounds for the interval function ln_sqrtx2y2:
570 real ln_x2y2_abs(2.225076E-308); // Absolute error bond
571 real q_lnx2y2p(4503599627370502.0 / 4503599627370496.0); // 1+e(f)
572 real q_lnx2y2m(9007199254740984.0 / 9007199254740992.0); // 1-e(f)
573 // With the following b0
574 // real b0 = 6369051672525773.0 / 30191699398572330817932436647906151127335369763331523427009650401964993299137190816689013801421270140331747000246110759198164677039398341060491474011461568349195162615808.0;
575 real b0 = MakeHexReal(0,1022-510,0x0016A09E,0x667F3BCD);
576 // it holds:
577 // 1. b < bo ==> g(b) := (0.5*b)*b < MinReal with rounding downwards
578 // 2. b >= b0 ==> g(b) := (0.5*b)*b >= MinReal with arbitrary rounding
579 // modus by the two multiplications.
580 
581 interval ln_sqrtx2y2(const interval& x, const interval& y) noexcept
582 // ln( sqrt(x^2+y^2) ) == 0.5*ln(x^2+y^2); Blomquist, 22.11.03;
583 {
584  interval ax=abs(x), ay=abs(y);
585  real Ix=Inf(ax), Sx=Sup(ax), Iy=Inf(ay), Sy=Sup(ay),f,u1,u2;
586  // Calculating the lower bound u1:
587  f = ln_sqrtx2y2(Ix,Iy);
588  if ((Ix==1 && Iy<b0) || (Iy==1 && Ix<b0)) {
589  // f in the denormalized range!
590  u1 = f - ln_x2y2_abs; // directed rounding not necessary!
591  if (sign(u1)<0) u1 = 0;
592  } else u1 = (sign(f)<0) ? f*q_lnx2y2p : f*q_lnx2y2m;
593  // Calculating the upper bound u2:
594  if (Ix==Sx && Iy==Sy) // x and y are point-intervals
595  if ((Sx==1 && Sy<b0) || (Sy==1 && Sx<b0)) {
596  // f in the denormalized range!
597  u2 = (Sy==0 || Sx==0) ? f : f+ln_x2y2_abs;
598  } else u2 = (sign(f)<0) ? f*q_lnx2y2m : f*q_lnx2y2p;
599  else // x or y is no point-interval:
600  {
601  f = ln_sqrtx2y2(Sx,Sy);
602  if ((Sx==1 && Sy<b0) || (Sy==1 && Sx<b0))
603  // f in the denormalized range!
604  u2 = (sign(Sy)==0 || sign(Sx)==0) ? f : f+ln_x2y2_abs;
605  else u2 = (sign(f)<0) ? f*q_lnx2y2m : f*q_lnx2y2p;
606  }
607  return interval(u1,u2);
608 } // ln_sqrtx2y2
609 
610 
611 // Constants for the interval function acoshp1(x) = acosh(1+x):
612 // (1+e(f)):
613 static const real q_acoshp1p(4503599627370503.0/4503599627370496.0);
614 // (1-e(f))
615 static const real q_acoshp1m(9007199254740981.0/9007199254740992.0);
616 
618 // acoshp1; Blomquist, 28.03.2005;
619 {
620  real r1,r2,sx,ix;
621  sx = Sup(x); ix = Inf(x);
622  // Calculating of the lower bound r1:
623  r2 = acoshp1(ix);
624  r1 = r2 * q_acoshp1m;
625  // Calculating of the upper bound r2:
626  r2 = (sx==ix)? r2*q_acoshp1p : acoshp1(sx)*q_acoshp1p;
627  return interval(r1,r2);
628 
629 } // acoshp1 (interval arguments)
630 
631 // *****************************************************************************
632 // sin(Pi*x)/Pi *
633 // ********************************************************************************
634 
635 // Die folgenden Konstanten q_sinpix_p, q_sinpix_m beziehen sich
636 // auf eps(sinpix_pi) = 3.401895e-16; f(x) = sin(pi*x)/pi;
637 
638 const real q_sinpix_p = 4503599627370499.0 / 4503599627370496.0; // (1+e(f))
639 const real q_sinpix_m = 9007199254740986.0 / 9007199254740992.0; // (1-e(f))
640 
641 real rounded_up(const real& x)
642 {
643  real y;
644  y = (x>=0)? x*q_sinpix_p : x*q_sinpix_m;
645  return y;
646 }
647 
648 real rounded_down(const real& x)
649 {
650  real y;
651  y = (x>=0)? x*q_sinpix_m : x*q_sinpix_p;
652  return y;
653 }
654 
656 {
657  const real Pir = Sup(Pir_interval); // 1/pi upwards rounded
658  interval y;
659  int ma,mb;
660  real y1,y2,a(Inf(x)),b(Sup(x)),ya,yb;
661  bool bl, ya_klg_yb;
662 
663  ma = Round(a); mb = Round(b);
664  bl = (ma%2 != 0);
665  if (mb-ma>1)
666  {
667  y1 = -Pir; y2 = Pir;
668  }
669  else // 0 <= mb-ma <=1
670  if (mb==ma)
671  if (a==b) // x: Point interval
672  {
673  ya = sinpix_pi(a);
674  y1 = rounded_down(ya);
675  y2 = rounded_up(ya);
676  }
677  else // (mb==ma) and (a!=b)
678  {
679  ya = sinpix_pi(a);
680  yb = sinpix_pi(b);
681  if (!bl)
682  {
683  y1 = rounded_down(ya);
684  y2 = rounded_up(yb);
685  }
686  else
687  {
688  y1 = rounded_down(yb);
689  y2 = rounded_up(ya);
690  }
691  }
692  else // mb-ma=1;
693  {
694  ya = sinpix_pi(a);
695  yb = sinpix_pi(b);
696  ya_klg_yb = (ya <= yb);
697 
698  if (bl)
699  {
700  if (!ya_klg_yb)
701  yb = ya;
702  ya = -Pir;
703  }
704  else
705  {
706  if (!ya_klg_yb)
707  ya = yb;
708  yb = Pir;
709  }
710  y1 = rounded_down(ya);
711  if (y1<-Pir)
712  y1 = -Pir;
713  y2 = rounded_up(yb);
714  if (y2>Pir)
715  y2 = Pir;
716  }
717 
718  y = interval(y1,y2);
719  return y;
720 }
721 
722 // ********************************************************************************
723 // *************** Konstanten bez. Gamma(x) und 1/Gamma(x) ********************
724 // ********************************************************************************
725 
726 /* worst case relative error bound for 1/Gamma(x) */
727 /* eps(gammar) = 2.866906e-15; */
728 /* q_gammarm = 1 - eps(gammar) */
729 const real q_gammarm = 9007199254740964.0 / 9007199254740992.0;
730 /* q_gammarp = 1 + eps(gammar) */
731 const real q_gammarp = 4503599627370510.0 / 4503599627370496.0;
732 
733 // Inclusion of the extremes of 1/gamma(x):
734 interval pow2(const interval& x, int ex)
735 {
736  interval y(x);
737  times2pown(y,ex);
738  return y;
739 }
740 
741 const real Ne = 9007199254740992.0;
742 const real Ne1 = 1125899906842624.0;
743 const real Ne2 = 562949953421312.0;
744 const real Ne3 = 281474976710656.0;
745 const real Ne4 = 140737488355328.0;
746 const real Ne5 = 70368744177664.0;
747 const real Ne6 = 35184372088832.0;
748 
749 const interval gam_rxi[171] =
750 {
751  interval( 6582605572834349.0 / 4503599627370496.0,6582606400588712.0 /
752  4503599627370496.0 ),
753  interval( -4540376432147063.0 / 9007199254740992.0,-4540375772996112.0 /
754  9007199254740992.0 ),
755  interval( -7086407292338520.0 / 4503599627370496.0,-7086406981597106.0 /
756  4503599627370496.0 ),
757  interval( -5878820838740338.0 / 2251799813685248.0,-5878820690102701.0 /
758  2251799813685248.0 ),
759  interval( -8185952996852629.0 / 2251799813685248.0,-8185952850644519.0 /
760  2251799813685248.0 ),
761  interval( -5239079997162568.0 / Ne1,-5239079928185648.0 /
762  Ne1 ),
763  interval( -6380657697812205.0 / Ne1,-6380657632438250.0 /
764  Ne1 ),
765  interval( -7519230477777525.0 / Ne1,-7519230410301402.0 /
766  Ne1 ),
767  interval( -8655680190901081.0 / Ne1,-8655680125714323.0 /
768  Ne1 ),
769  interval( -4895280046470312.0 / Ne2,-4895280015191181.0 /
770  Ne2 ),
771  interval( -5462119069950045.0 / Ne2,-5462119039144308.0 /
772  Ne2 ),
773  interval( -6028485171921533.0 / Ne2,-6028485140574985.0 /
774  Ne2 ),
775  interval( -6594470676196825.0 / Ne2,-6594470646005838.0 /
776  Ne2 ),
777  interval( -7160144161412306.0 / Ne2,-7160144131821972.0 /
778  Ne2 ),
779  interval( -7725557826948019.0 / Ne2,-7725557796923813.0 /
780  Ne2 ),
781  interval( -8290752238810453.0 / Ne2,-8290752208368199.0 /
782  Ne2 ),
783  interval( -8855759486553113.0 / Ne2,-8855759457969009.0 /
784  Ne2 ),
785  interval( -4710302676530551.0 / Ne3,-4710302661730969.0 /
786  Ne3 ),
787  interval( -4992655414739558.0 / Ne3,-4992655400196254.0 /
788  Ne3 ),
789  interval( -5274946608174960.0 / Ne3,-5274946594248303.0 /
790  Ne3 ),
791  interval( -5557183461054268.0 / Ne3,-5557183446978818.0 /
792  Ne3 ),
793  interval( -5839372030862353.0 / Ne3,-5839372016359958.0 /
794  Ne3 ),
795  interval( -6121517453741464.0 / Ne3,-6121517439898965.0 /
796  Ne3 ),
797  interval( -6403624121061720.0 / Ne3,-6403624107245033.0 /
798  Ne3 ),
799  interval( -6685695811452843.0 / Ne3,-6685695797850064.0 /
800  Ne3 ),
801  interval( -6967735798799276.0 / Ne3,-6967735784842353.0 /
802  Ne3 ),
803  interval( -7249746935136834.0 / Ne3,-7249746921647546.0 /
804  Ne3 ),
805  interval( -7531731721183113.0 / Ne3,-7531731707547301.0 /
806  Ne3 ),
807  interval( -7813692357395941.0 / Ne3,-7813692343778691.0 /
808  Ne3 ),
809  interval( -8095630791428232.0 / Ne3,-8095630778040390.0 /
810  Ne3 ),
811  interval( -8377548753853165.0 / Ne3,-8377548740538224.0 /
812  Ne3 ),
813  interval( -8659447788741678.0 / Ne3,-8659447775645579.0 /
814  Ne3 ),
815  interval( -8941329278863280.0 / Ne3,-8941329265761967.0 /
816  Ne3 ),
817  interval( -4611597233397515.0 / Ne4,-4611597226897947.0 /
818  Ne4 ),
819  interval( -4752522236545681.0 / Ne4,-4752522230097182.0 /
820  Ne4 ),
821  interval( -4893440155674552.0 / Ne4,-4893440149271021.0 /
822  Ne4 ),
823  interval( -5034351450592848.0 / Ne4,-5034351444042898.0 /
824  Ne4 ),
825  interval( -5175256539394880.0 / Ne4,-5175256533033078.0 /
826  Ne4 ),
827  interval( -5316155803584886.0 / Ne4,-5316155797241740.0 /
828  Ne4 ),
829  interval( -5457049591969558.0 / Ne4,-5457049585698186.0 /
830  Ne4 ),
831  interval( -5597938224227692.0 / Ne4,-5597938218047102.0 /
832  Ne4 ),
833  interval( -5738821993949958.0 / Ne4,-5738821987610287.0 /
834  Ne4 ),
835  interval( -5879701171316002.0 / Ne4,-5879701164978161.0 /
836  Ne4 ),
837  interval( -6020576005634620.0 / Ne4,-6020575999374578.0 /
838  Ne4 ),
839  interval( -6161446727231484.0 / Ne4,-6161446721061003.0 /
840  Ne4 ),
841  interval( -6302313549465809.0 / Ne4,-6302313543147939.0 /
842  Ne4 ),
843  interval( -6443176669927728.0 / Ne4,-6443176663720294.0 /
844  Ne4 ),
845  interval( -6584036272403005.0 / Ne4,-6584036266213061.0 /
846  Ne4 ),
847  interval( -6724892527807298.0 / Ne4,-6724892521515390.0 /
848  Ne4 ),
849  interval( -6865745595463323.0 / Ne4,-6865745589218114.0 /
850  Ne4 ),
851  interval( -7006595624078029.0 / Ne4,-7006595617954995.0 /
852  Ne4 ),
853  interval( -7147442752315627.0 / Ne4,-7147442746313355.0 /
854  Ne4 ),
855  interval( -7288287110549528.0 / Ne4,-7288287104399241.0 /
856  Ne4 ),
857  interval( -7429128820262002.0 / Ne4,-7429128814290427.0 /
858  Ne4 ),
859  interval( -7569967996009183.0 / Ne4,-7569967989919521.0 /
860  Ne4 ),
861  interval( -7710804744971319.0 / Ne4,-7710804738911590.0 /
862  Ne4 ),
863  interval( -7851639168099204.0 / Ne4,-7851639162048862.0 /
864  Ne4 ),
865  interval( -7992471360380410.0 / Ne4,-7992471354478088.0 /
866  Ne4 ),
867  interval( -8133301411685542.0 / Ne4,-8133301405639704.0 /
868  Ne4 ),
869  interval( -8274129406251900.0 / Ne4,-8274129400330745.0 /
870  Ne4 ),
871  interval( -8414955423947592.0 / Ne4,-8414955418025784.0 /
872  Ne4 ),
873  interval( -8555779540237542.0 / Ne4,-8555779534343191.0 /
874  Ne4 ),
875  interval( -8696601826519818.0 / Ne4,-8696601820560983.0 /
876  Ne4 ),
877  interval( -8837422350374443.0 / Ne4,-8837422344547779.0 /
878  Ne4 ),
879  interval( -8978241175812537.0 / Ne4,-8978241170031881.0 /
880  Ne4 ),
881  interval( -4559529181955286.0 / Ne5,-4559529178973419.0 /
882  Ne5 ),
883  interval( -4629936985962592.0 / Ne5,-4629936983062621.0 /
884  Ne5 ),
885  interval( -4700344027557972.0 / Ne5,-4700344024658010.0 /
886  Ne5 ),
887  interval( -4770750332748164.0 / Ne5,-4770750329806175.0 /
888  Ne5 ),
889  interval( -4841155926312136.0 / Ne5,-4841155923514622.0 /
890  Ne5 ),
891  interval( -4911560831959402.0 / Ne5,-4911560829205074.0 /
892  Ne5 ),
893  interval( -4981965072279533.0 / Ne5,-4981965069314643.0 /
894  Ne5 ),
895  interval( -5052368668591817.0 / Ne5,-5052368665644419.0 /
896  Ne5 ),
897  interval( -5122771641448337.0 / Ne5,-5122771638497776.0 /
898  Ne5 ),
899  interval( -5193174010445884.0 / Ne5,-5193174007591888.0 /
900  Ne5 ),
901  interval( -5263575794318673.0 / Ne5,-5263575791488264.0 /
902  Ne5 ),
903  interval( -5333977010931459.0 / Ne5,-5333977008101358.0 /
904  Ne5 ),
905  interval( -5404377677505176.0 / Ne5,-5404377674566742.0 /
906  Ne5 ),
907  interval( -5474777810213221.0 / Ne5,-5474777807391522.0 /
908  Ne5 ),
909  interval( -5545177424917448.0 / Ne5,-5545177422065289.0 /
910  Ne5 ),
911  interval( -5615576536591204.0 / Ne5,-5615576533728944.0 /
912  Ne5 ),
913  interval( -5685975159660620.0 / Ne5,-5685975156829327.0 /
914  Ne5 ),
915  interval( -5756373307993982.0 / Ne5,-5756373305121310.0 /
916  Ne5 ),
917  interval( -5826770994840416.0 / Ne5,-5826770992034928.0 /
918  Ne5 ),
919  interval( -5897168232948637.0 / Ne5,-5897168230130823.0 /
920  Ne5 ),
921  interval( -5967565034702908.0 / Ne5,-5967565031698652.0 /
922  Ne5 ),
923  interval( -6037961411567024.0 / Ne5,-6037961408739972.0 /
924  Ne5 ),
925  interval( -6108357375160195.0 / Ne5,-6108357372312820.0 /
926  Ne5 ),
927  interval( -6178752936267084.0 / Ne5,-6178752933424618.0 /
928  Ne5 ),
929  interval( -6249148105390399.0 / Ne5,-6249148102629437.0 /
930  Ne5 ),
931  interval( -6319542892697349.0 / Ne5,-6319542889886519.0 /
932  Ne5 ),
933  interval( -6389937307844225.0 / Ne5,-6389937305053873.0 /
934  Ne5 ),
935  interval( -6460331360217357.0 / Ne5,-6460331357417915.0 /
936  Ne5 ),
937  interval( -6530725058889375.0 / Ne5,-6530725056107292.0 /
938  Ne5 ),
939  interval( -6601118412624317.0 / Ne5,-6601118409842826.0 /
940  Ne5 ),
941  interval( -6671511429748937.0 / Ne5,-6671511426971253.0 /
942  Ne5 ),
943  interval( -6741904118401788.0 / Ne5,-6741904115722630.0 /
944  Ne5 ),
945  interval( -6812296486576567.0 / Ne5,-6812296483762819.0 /
946  Ne5 ),
947  interval( -6882688541678557.0 / Ne5,-6882688538949795.0 /
948  Ne5 ),
949  interval( -6953080291125250.0 / Ne5,-6953080288362282.0 /
950  Ne5 ),
951  interval( -7023471742013192.0 / Ne5,-7023471739228496.0 /
952  Ne5 ),
953  interval( -7093862901097948.0 / Ne5,-7093862898353711.0 /
954  Ne5 ),
955  interval( -7164253775147266.0 / Ne5,-7164253772348037.0 /
956  Ne5 ),
957  interval( -7234644370421069.0 / Ne5,-7234644367656947.0 /
958  Ne5 ),
959  interval( -7305034693228153.0 / Ne5,-7305034690493625.0 /
960  Ne5 ),
961  interval( -7375424749560875.0 / Ne5,-7375424746816433.0 /
962  Ne5 ),
963  interval( -7445814545246703.0 / Ne5,-7445814542529259.0 /
964  Ne5 ),
965  interval( -7516204085881418.0 / Ne5,-7516204083182192.0 /
966  Ne5 ),
967  interval( -7586593377039902.0 / Ne5,-7586593374291972.0 /
968  Ne5 ),
969  interval( -7656982423878545.0 / Ne5,-7656982421139464.0 /
970  Ne5 ),
971  interval( -7727371231629585.0 / Ne5,-7727371228950751.0 /
972  Ne5 ),
973  interval( -7797759805243431.0 / Ne5,-7797759802601748.0 /
974  Ne5 ),
975  interval( -7868148149631789.0 / Ne5,-7868148146919441.0 /
976  Ne5 ),
977  interval( -7938536269389156.0 / Ne5,-7938536266702540.0 /
978  Ne5 ),
979  interval( -8008924169145439.0 / Ne5,-8008924166482991.0 /
980  Ne5 ),
981  interval( -8079311853288492.0 / Ne5,-8079311850638146.0 /
982  Ne5 ),
983  interval( -8149699326208468.0 / Ne5,-8149699323496057.0 /
984  Ne5 ),
985  interval( -8220086591946405.0 / Ne5,-8220086589180599.0 /
986  Ne5 ),
987  interval( -8290473654625668.0 / Ne5,-8290473651911160.0 /
988  Ne5 ),
989  interval( -8360860518207539.0 / Ne5,-8360860515490326.0 /
990  Ne5 ),
991  interval( -8431247186492146.0 / Ne5,-8431247183838485.0 /
992  Ne5 ),
993  interval( -8501633663256230.0 / Ne5,-8501633660616442.0 /
994  Ne5 ),
995  interval( -8572019952114999.0 / Ne5,-8572019949515073.0 /
996  Ne5 ),
997  interval( -8642406056614155.0 / Ne5,-8642406053936017.0 /
998  Ne5 ),
999  interval( -8712791980171444.0 / Ne5,-8712791977516241.0 /
1000  Ne5 ),
1001  interval( -8783177726114270.0 / Ne5,-8783177723474143.0 /
1002  Ne5 ),
1003  interval( -8853563297747803.0 / Ne5,-8853563295082565.0 /
1004  Ne5 ),
1005  interval( -8923948698211540.0 / Ne5,-8923948695603642.0 /
1006  Ne5 ),
1007  interval( -8994333930651793.0 / Ne5,-8994333928013544.0 /
1008  Ne5 ),
1009  interval( -4532359499005981.0 / Ne6,-4532359497669548.0 /
1010  Ne6 ),
1011  interval( -4567551951590474.0 / Ne6,-4567551950297444.0 /
1012  Ne6 ),
1013  interval( -4602744324609335.0 / Ne6,-4602744323253914.0 /
1014  Ne6 ),
1015  interval( -4637936619312858.0 / Ne6,-4637936618035482.0 /
1016  Ne6 ),
1017  interval( -4673128837201202.0 / Ne6,-4673128835890259.0 /
1018  Ne6 ),
1019  interval( -4708320979539414.0 / Ne6,-4708320978220509.0 /
1020  Ne6 ),
1021  interval( -4743513047586994.0 / Ne6,-4743513046289186.0 /
1022  Ne6 ),
1023  interval( -4778705042649671.0 / Ne6,-4778705041357978.0 /
1024  Ne6 ),
1025  interval( -4813896965939740.0 / Ne6,-4813896964638916.0 /
1026  Ne6 ),
1027  interval( -4849088818685547.0 / Ne6,-4849088817383583.0 /
1028  Ne6 ),
1029  interval( -4884280602023663.0 / Ne6,-4884280600723731.0 /
1030  Ne6 ),
1031  interval( -4919472317096448.0 / Ne6,-4919472315799988.0 /
1032  Ne6 ),
1033  interval( -4954663965077128.0 / Ne6,-4954663963749440.0 /
1034  Ne6 ),
1035  interval( -4989855546965370.0 / Ne6,-4989855545671197.0 /
1036  Ne6 ),
1037  interval( -5025047063935622.0 / Ne6,-5025047062637366.0 /
1038  Ne6 ),
1039  interval( -5060238516963350.0 / Ne6,-5060238515663821.0 /
1040  Ne6 ),
1041  interval( -5095429907060327.0 / Ne6,-5095429905761378.0 /
1042  Ne6 ),
1043  interval( -5130621235251224.0 / Ne6,-5130621233930590.0 /
1044  Ne6 ),
1045  interval( -5165812502482099.0 / Ne6,-5165812501185488.0 /
1046  Ne6 ),
1047  interval( -5201003709724055.0 / Ne6,-5201003708430144.0 /
1048  Ne6 ),
1049  interval( -5236194857910350.0 / Ne6,-5236194856593808.0 /
1050  Ne6 ),
1051  interval( -5271385947936424.0 / Ne6,-5271385946609048.0 /
1052  Ne6 ),
1053  interval( -5306576980681870.0 / Ne6,-5306576979376506.0 /
1054  Ne6 ),
1055  interval( -5341767957039505.0 / Ne6,-5341767955734187.0 /
1056  Ne6 ),
1057  interval( -5376958877810379.0 / Ne6,-5376958876528096.0 /
1058  Ne6 ),
1059  interval( -5412149743940434.0 / Ne6,-5412149742614602.0 /
1060  Ne6 ),
1061  interval( -5447340556059975.0 / Ne6,-5447340554811246.0 /
1062  Ne6 ),
1063  interval( -5482531315163654.0 / Ne6,-5482531313857339.0 /
1064  Ne6 ),
1065  interval( -5517722021905912.0 / Ne6,-5517722020635122.0 /
1066  Ne6 ),
1067  interval( -5552912677100886.0 / Ne6,-5552912675823002.0 /
1068  Ne6 ),
1069  interval( -5588103281480912.0 / Ne6,-5588103280171229.0 /
1070  Ne6 ),
1071  interval( -5623293835745456.0 / Ne6,-5623293834493516.0 /
1072  Ne6 ),
1073  interval( -5658484340702630.0 / Ne6,-5658484339399074.0 /
1074  Ne6 ),
1075  interval( -5693674796958315.0 / Ne6,-5693674795683212.0 /
1076  Ne6 ),
1077  interval( -5728865205222665.0 / Ne6,-5728865203973973.0 /
1078  Ne6 ),
1079  interval( -5764055566229013.0 / Ne6,-5764055564966520.0 /
1080  Ne6 ),
1081  interval( -5799245880597279.0 / Ne6,-5799245879305275.0 /
1082  Ne6 ),
1083  interval( -5834436148937784.0 / Ne6,-5834436147684294.0 /
1084  Ne6 ),
1085  interval( -5869626371953711.0 / Ne6,-5869626370690622.0 /
1086  Ne6 ),
1087  interval( -5904816550239413.0 / Ne6,-5904816548965896.0 /
1088  Ne6 ),
1089  interval( -5940006684383290.0 / Ne6,-5940006683093915.0 /
1090  Ne6 ),
1091  interval( -5975196775000579.0 / Ne6,-5975196773734016.0 /
1092  Ne6 ) };
1093 // Inclusion of the extremes of 1/gamma(x):
1094 const interval gam_ryi[171] = {
1095 pow2( interval( 5085347089749720.0 / Ne,5085347089749823.0 / Ne ) , 1 ) ,
1096 pow2( interval( -5082146609264467.0 / Ne,-5082146609264314.0 / Ne ) , -1 ) ,
1097 pow2( interval( 7824158147621733.0 / Ne,7824158147621966.0 / Ne ) , -1 ) ,
1098 pow2( interval( -5070842539852372.0 / Ne,-5070842539852221.0 / Ne ) , 1 ) ,
1099 pow2( interval( 4593118780547419.0 / Ne,4593118780547576.0 / Ne ) , 3 ) ,
1100 pow2( interval( -5333021955274733.0 / Ne,-5333021955274575.0 / Ne ) , 5 ) ,
1101 pow2( interval( 7546574203185105.0 / Ne,7546574203185319.0 / Ne ) , 7 ) ,
1102 pow2( interval( -6294628859031764.0 / Ne,-6294628859031469.0 / Ne ) , 10 ) ,
1103 pow2( interval( 6045310252810166.0 / Ne,6045310252811273.0 / Ne ) , 13 ) ,
1104 pow2( interval( -6568078652156336.0 / Ne,-6568078652156148.0 / Ne ) , 16 ) ,
1105 pow2( interval( 7963169065060572.0 / Ne,7963169065060801.0 / Ne ) , 19 ) ,
1106 pow2( interval( -5328217018030122.0 / Ne,-5328217018029960.0 / Ne ) , 23 ) ,
1107 pow2( interval( 7800142897041864.0 / Ne,7800142897042089.0 / Ne ) , 26 ) ,
1108 pow2( interval( -6199437664213474.0 / Ne,-6199437664213297.0 / Ne ) , 30 ) ,
1109 pow2( interval( 5316470282961123.0 / Ne,5316470282961284.0 / Ne ) , 34 ) ,
1110 pow2( interval( -4892929765135337.0 / Ne,-4892929765135165.0 / Ne ) , 38 ) ,
1111 pow2( interval( 4810107119289947.0 / Ne,4810107119290088.0 / Ne ) , 42 ) ,
1112 pow2( interval( -5030373421375086.0 / Ne,-5030373421374834.0 / Ne ) , 46 ) ,
1113 pow2( interval( 5576144001185310.0 / Ne,5576144001185479.0 / Ne ) , 50 ) ,
1114 pow2( interval( -6530685487420963.0 / Ne,-6530685487420774.0 / Ne ) , 54 ) ,
1115 pow2( interval( 8057940169576582.0 / Ne,8057940169576818.0 / Ne ) , 58 ) ,
1116 pow2( interval( -5223648494045513.0 / Ne,-5223648494045349.0 / Ne ) , 63 ) ,
1117 pow2( interval( 7099855957135674.0 / Ne,7099855957135885.0 / Ne ) , 67 ) ,
1118 pow2( interval( -5047359382236272.0 / Ne,-5047359382236084.0 / Ne ) , 72 ) ,
1119 pow2( interval( 7492585872478835.0 / Ne,7492585872479188.0 / Ne ) , 76 ) ,
1120 pow2( interval( -5795835662380422.0 / Ne,-5795835662380242.0 / Ne ) , 81 ) ,
1121 pow2( interval( 4664800910382651.0 / Ne,4664800910382790.0 / Ne ) , 86 ) ,
1122 pow2( interval( -7801058080117709.0 / Ne,-7801058080117472.0 / Ne ) , 90 ) ,
1123 pow2( interval( 6767162072327001.0 / Ne,6767162072327282.0 / Ne ) , 95 ) ,
1124 pow2( interval( -6082121514218736.0 / Ne,-6082121514218554.0 / Ne ) , 100 ) ,
1125 pow2( interval( 5656800000052189.0 / Ne,5656800000052359.0 / Ne ) , 105 ) ,
1126 pow2( interval( -5438268378952110.0 / Ne,-5438268378951951.0 / Ne ) , 110 ) ,
1127 pow2( interval( 5398375606367166.0 / Ne,5398375606367329.0 / Ne ) , 115 ) ,
1128 pow2( interval( -5527713447587841.0 / Ne,-5527713447587674.0 / Ne ) , 120 ) ,
1129 pow2( interval( 5833125895912623.0 / Ne,5833125895912799.0 / Ne ) , 125 ) ,
1130 pow2( interval( -6337936184674347.0 / Ne,-6337936184674153.0 / Ne ) , 130 ) ,
1131 pow2( interval( 7084743510515278.0 / Ne,7084743510515501.0 / Ne ) , 135 ) ,
1132 pow2( interval( -8141214882701327.0 / Ne,-8141214882701088.0 / Ne ) , 140 ) ,
1133 pow2( interval( 4804968547193877.0 / Ne,4804968547194018.0 / Ne ) , 146 ) ,
1134 pow2( interval( -5822137580509526.0 / Ne,-5822137580509355.0 / Ne ) , 151 ) ,
1135 pow2( interval( 7236772755227956.0 / Ne,7236772755228162.0 / Ne ) , 156 ) ,
1136 pow2( interval( -4610758665056508.0 / Ne,-4610758665056369.0 / Ne ) , 162 ) ,
1137 pow2( interval( 6019530845699084.0 / Ne,6019530845699266.0 / Ne ) , 167 ) ,
1138 pow2( interval( -8047036389398365.0 / Ne,-8047036389398123.0 / Ne ) , 172 ) ,
1139 pow2( interval( 5504580189086749.0 / Ne,5504580189086968.0 / Ne ) , 178 ) ,
1140 pow2( interval( -7703001513324420.0 / Ne,-7703001513324183.0 / Ne ) , 183 ) ,
1141 pow2( interval( 5510183009440391.0 / Ne,5510183009440581.0 / Ne ) , 189 ) ,
1142 pow2( interval( -8055535954952413.0 / Ne,-8055535954952173.0 / Ne ) , 194 ) ,
1143 pow2( interval( 6014315232803007.0 / Ne,6014315232803294.0 / Ne ) , 200 ) ,
1144 pow2( interval( -4584378555360492.0 / Ne,-4584378555360260.0 / Ne ) , 206 ) ,
1145 pow2( interval( 7132212380084113.0 / Ne,7132212380084326.0 / Ne ) , 211 ) ,
1146 pow2( interval( -5659549393054692.0 / Ne,-5659549393054526.0 / Ne ) , 217 ) ,
1147 pow2( interval( 4579461117155838.0 / Ne,4579461117155977.0 / Ne ) , 223 ) ,
1148 pow2( interval( -7554216840666713.0 / Ne,-7554216840666493.0 / Ne ) , 228 ) ,
1149 pow2( interval( 6348787715758027.0 / Ne,6348787715758222.0 / Ne ) , 234 ) ,
1150 pow2( interval( -5434979980476367.0 / Ne,-5434979980476204.0 / Ne ) , 240 ) ,
1151 pow2( interval( 4737681191908824.0 / Ne,4737681191908967.0 / Ne ) , 246 ) ,
1152 pow2( interval( -8407842664867513.0 / Ne,-8407842664867267.0 / Ne ) , 251 ) ,
1153 pow2( interval( 7592052521188700.0 / Ne,7592052521188935.0 / Ne ) , 257 ) ,
1154 pow2( interval( -6974119252551297.0 / Ne,-6974119252551090.0 / Ne ) , 263 ) ,
1155 pow2( interval( 6515520808385677.0 / Ne,6515520808385874.0 / Ne ) , 269 ) ,
1156 pow2( interval( -6188946869743481.0 / Ne,-6188946869743300.0 / Ne ) , 275 ) ,
1157 pow2( interval( 5975502808844840.0 / Ne,5975502808845020.0 / Ne ) , 281 ) ,
1158 pow2( interval( -5862842897072874.0 / Ne,-5862842897072704.0 / Ne ) , 287 ) ,
1159 pow2( interval( 5843967448508660.0 / Ne,5843967448508828.0 / Ne ) , 293 ) ,
1160 pow2( interval( -5916517001341501.0 / Ne,-5916517001341321.0 / Ne ) , 299 ) ,
1161 pow2( interval( 6082464626325325.0 / Ne,6082464626325503.0 / Ne ) , 305 ) ,
1162 pow2( interval( -6348157530347044.0 / Ne,-6348157530346858.0 / Ne ) , 311 ) ,
1163 pow2( interval( 6724699799057619.0 / Ne,6724699799057843.0 / Ne ) , 317 ) ,
1164 pow2( interval( -7228705737680202.0 / Ne,-7228705737679999.0 / Ne ) , 323 ) ,
1165 pow2( interval( 7883493269720206.0 / Ne,7883493269720561.0 / Ne ) , 329 ) ,
1166 pow2( interval( -8720834785364833.0 / Ne,-8720834785364561.0 / Ne ) , 335 ) ,
1167 pow2( interval( 4891722644546351.0 / Ne,4891722644546502.0 / Ne ) , 342 ) ,
1168 pow2( interval( -5564236710028970.0 / Ne,-5564236710028799.0 / Ne ) , 348 ) ,
1169 pow2( interval( 6416191129172903.0 / Ne,6416191129173091.0 / Ne ) , 354 ) ,
1170 pow2( interval( -7498890927628704.0 / Ne,-7498890927628487.0 / Ne ) , 360 ) ,
1171 pow2( interval( 8881515552460572.0 / Ne,8881515552460999.0 / Ne ) , 366 ) ,
1172 pow2( interval( -5328950915550370.0 / Ne,-5328950915550206.0 / Ne ) , 373 ) ,
1173 pow2( interval( 6478093314396794.0 / Ne,6478093314397089.0 / Ne ) , 379 ) ,
1174 pow2( interval( -7976303366065662.0 / Ne,-7976303366065426.0 / Ne ) , 385 ) ,
1175 pow2( interval( 4972846688449830.0 / Ne,4972846688450017.0 / Ne ) , 392 ) ,
1176 pow2( interval( -6278401907481090.0 / Ne,-6278401907480879.0 / Ne ) , 398 ) ,
1177 pow2( interval( 8024854758356088.0 / Ne,8024854758356345.0 / Ne ) , 404 ) ,
1178 pow2( interval( -5191277948909595.0 / Ne,-5191277948909444.0 / Ne ) , 411 ) ,
1179 pow2( interval( 6797621462551740.0 / Ne,6797621462551941.0 / Ne ) , 417 ) ,
1180 pow2( interval( -4503636668393666.0 / Ne,-4503636668393518.0 / Ne ) , 424 ) ,
1181 pow2( interval( 6037997262493341.0 / Ne,6037997262493523.0 / Ne ) , 430 ) ,
1182 pow2( interval( -8189485306115383.0 / Ne,-8189485306115130.0 / Ne ) , 436 ) ,
1183 pow2( interval( 5617805845426844.0 / Ne,5617805845427124.0 / Ne ) , 443 ) ,
1184 pow2( interval( -7795192616785187.0 / Ne,-7795192616784477.0 / Ne ) , 449 ) ,
1185 pow2( interval( 5469175405734180.0 / Ne,5469175405734422.0 / Ne ) , 456 ) ,
1186 pow2( interval( -7759929987383324.0 / Ne,-7759929987383086.0 / Ne ) , 462 ) ,
1187 pow2( interval( 5565727978288701.0 / Ne,5565727978288876.0 / Ne ) , 469 ) ,
1188 pow2( interval( -8070914994857895.0 / Ne,-8070914994857635.0 / Ne ) , 475 ) ,
1189 pow2( interval( 5914931467943193.0 / Ne,5914931467943373.0 / Ne ) , 482 ) ,
1190 pow2( interval( -8762204548045716.0 / Ne,-8762204548045455.0 / Ne ) , 488 ) ,
1191 pow2( interval( 6558513517606168.0 / Ne,6558513517606353.0 / Ne ) , 495 ) ,
1192 pow2( interval( -4960305627886271.0 / Ne,-4960305627886120.0 / Ne ) , 502 ) ,
1193 pow2( interval( 7580642983583672.0 / Ne,7580642983583897.0 / Ne ) , 508 ) ,
1194 pow2( interval( -5851844804194595.0 / Ne,-5851844804194367.0 / Ne ) , 515 ) ,
1195 pow2( interval( 4563038858728436.0 / Ne,4563038858728577.0 / Ne ) , 522 ) ,
1196 pow2( interval( -7187477492053316.0 / Ne,-7187477492052964.0 / Ne ) , 528 ) ,
1197 pow2( interval( 5716852908386950.0 / Ne,5716852908387214.0 / Ne ) , 535 ) ,
1198 pow2( interval( -4591808630269563.0 / Ne,-4591808630269411.0 / Ne ) , 542 ) ,
1199 pow2( interval( 7448102539955649.0 / Ne,7448102539955986.0 / Ne ) , 548 ) ,
1200 pow2( interval( -6098770429791387.0 / Ne,-6098770429791204.0 / Ne ) , 555 ) ,
1201 pow2( interval( 5041550443966798.0 / Ne,5041550443966946.0 / Ne ) , 562 ) ,
1202 pow2( interval( -8413996086583072.0 / Ne,-8413996086582821.0 / Ne ) , 568 ) ,
1203 pow2( interval( 7086939987269423.0 / Ne,7086939987269731.0 / Ne ) , 575 ) ,
1204 pow2( interval( -6024570065319942.0 / Ne,-6024570065319682.0 / Ne ) , 582 ) ,
1205 pow2( interval( 5168535487082451.0 / Ne,5168535487082609.0 / Ne ) , 589 ) ,
1206 pow2( interval( -8949051953781375.0 / Ne,-8949051953781115.0 / Ne ) , 595 ) ,
1207 pow2( interval( 7817344426895164.0 / Ne,7817344426895996.0 / Ne ) , 602 ) ,
1208 pow2( interval( -6889843867972878.0 / Ne,-6889843867972674.0 / Ne ) , 609 ) ,
1209 pow2( interval( 6126229646423302.0 / Ne,6126229646423484.0 / Ne ) , 616 ) ,
1210 pow2( interval( -5495122334906381.0 / Ne,-5495122334906222.0 / Ne ) , 623 ) ,
1211 pow2( interval( 4971972094727164.0 / Ne,4971972094727314.0 / Ne ) , 630 ) ,
1212 pow2( interval( -4537480959802395.0 / Ne,-4537480959802254.0 / Ne ) , 637 ) ,
1213 pow2( interval( 8352835047353300.0 / Ne,8352835047353555.0 / Ne ) , 643 ) ,
1214 pow2( interval( -7753443787904532.0 / Ne,-7753443787904298.0 / Ne ) , 650 ) ,
1215 pow2( interval( 7257653550749169.0 / Ne,7257653550749382.0 / Ne ) , 657 ) ,
1216 pow2( interval( -6850281165773769.0 / Ne,-6850281165773570.0 / Ne ) , 664 ) ,
1217 pow2( interval( 6519305845448896.0 / Ne,6519305845449168.0 / Ne ) , 671 ) ,
1218 pow2( interval( -6255266499085062.0 / Ne,-6255266499084872.0 / Ne ) , 678 ) ,
1219 pow2( interval( 6050802311308162.0 / Ne,6050802311308350.0 / Ne ) , 685 ) ,
1220 pow2( interval( -5900304762620398.0 / Ne,-5900304762620223.0 / Ne ) , 692 ) ,
1221 pow2( interval( 5799657649647993.0 / Ne,5799657649648165.0 / Ne ) , 699 ) ,
1222 pow2( interval( -5746047975553302.0 / Ne,-5746047975553134.0 / Ne ) , 706 ) ,
1223 pow2( interval( 5737835419331524.0 / Ne,5737835419331693.0 / Ne ) , 713 ) ,
1224 pow2( interval( -5774471890994117.0 / Ne,-5774471890993944.0 / Ne ) , 720 ) ,
1225 pow2( interval( 5856465763387432.0 / Ne,5856465763387600.0 / Ne ) , 727 ) ,
1226 pow2( interval( -5985387992102590.0 / Ne,-5985387992102406.0 / Ne ) , 734 ) ,
1227 pow2( interval( 6163919695584074.0 / Ne,6163919695584257.0 / Ne ) , 741 ) ,
1228 pow2( interval( -6395943042753787.0 / Ne,-6395943042753502.0 / Ne ) , 748 ) ,
1229 pow2( interval( 6686679647283150.0 / Ne,6686679647283350.0 / Ne ) , 755 ) ,
1230 pow2( interval( -7042883260256940.0 / Ne,-7042883260256730.0 / Ne ) , 762 ) ,
1231 pow2( interval( 7473096566380533.0 / Ne,7473096566380749.0 / Ne ) , 769 ) ,
1232 pow2( interval( -7987985534527481.0 / Ne,-7987985534527243.0 / Ne ) , 776 ) ,
1233 pow2( interval( 8600769311605383.0 / Ne,8600769311605633.0 / Ne ) , 783 ) ,
1234 pow2( interval( -4663884705694464.0 / Ne,-4663884705694325.0 / Ne ) , 791 ) ,
1235 pow2( interval( 5094554684614484.0 / Ne,5094554684614634.0 / Ne ) , 798 ) ,
1236 pow2( interval( -5604802840349871.0 / Ne,-5604802840349701.0 / Ne ) , 805 ) ,
1237 pow2( interval( 6209951739735886.0 / Ne,6209951739736072.0 / Ne ) , 812 ) ,
1238 pow2( interval( -6928963530888061.0 / Ne,-6928963530887851.0 / Ne ) , 819 ) ,
1239 pow2( interval( 7785368708274196.0 / Ne,7785368708274423.0 / Ne ) , 826 ) ,
1240 pow2( interval( -8808459126256481.0 / Ne,-8808459126256060.0 / Ne ) , 833 ) ,
1241 pow2( interval( 5017412797579486.0 / Ne,5017412797579638.0 / Ne ) , 841 ) ,
1242 pow2( interval( -5755173329981532.0 / Ne,-5755173329981361.0 / Ne ) , 848 ) ,
1243 pow2( interval( 6646385258439176.0 / Ne,6646385258439444.0 / Ne ) , 855 ) ,
1244 pow2( interval( -7727539896552529.0 / Ne,-7727539896552294.0 / Ne ) , 862 ) ,
1245 pow2( interval( 4522473425691912.0 / Ne,4522473425692052.0 / Ne ) , 870 ) ,
1246 pow2( interval( -5328812572761788.0 / Ne,-5328812572761623.0 / Ne ) , 877 ) ,
1247 pow2( interval( 6320558000502691.0 / Ne,6320558000502885.0 / Ne ) , 884 ) ,
1248 pow2( interval( -7546265781200776.0 / Ne,-7546265781200489.0 / Ne ) , 891 ) ,
1249 pow2( interval( 4534316912522546.0 / Ne,4534316912522688.0 / Ne ) , 899 ) ,
1250 pow2( interval( -5484491485989575.0 / Ne,-5484491485989407.0 / Ne ) , 906 ) ,
1251 pow2( interval( 6676632315202014.0 / Ne,6676632315202302.0 / Ne ) , 913 ) ,
1252 pow2( interval( -8180074398476253.0 / Ne,-8180074398476014.0 / Ne ) , 920 ) ,
1253 pow2( interval( 5042989707083422.0 / Ne,5042989707083666.0 / Ne ) , 928 ) ,
1254 pow2( interval( -6257379418480333.0 / Ne,-6257379418480019.0 / Ne ) , 935 ) ,
1255 pow2( interval( 7813097673618694.0 / Ne,7813097673619043.0 / Ne ) , 942 ) ,
1256 pow2( interval( -4908325621754370.0 / Ne,-4908325621754217.0 / Ne ) , 950 ) ,
1257 pow2( interval( 6205346227418363.0 / Ne,6205346227418597.0 / Ne ) , 957 ) ,
1258 pow2( interval( -7893590972392525.0 / Ne,-7893590972392227.0 / Ne ) , 964 ) ,
1259 pow2( interval( 5051411882876310.0 / Ne,5051411882876506.0 / Ne ) , 972 ) ,
1260 pow2( interval( -6504655602059905.0 / Ne,-6504655602059583.0 / Ne ) , 979 ) ,
1261 pow2( interval( 8426810051054742.0 / Ne,8426810051054986.0 / Ne ) , 986 ) ,
1262 pow2( interval( -5491407534973626.0 / Ne,-5491407534973452.0 / Ne ) , 994 ) ,
1263 pow2( interval( 7199960218142557.0 / Ne,7199960218142768.0 / Ne ) , 1001 ) ,
1264 pow2( interval( -4748178637044143.0 / Ne,-4748178637044000.0 / Ne ) , 1009 ) ,
1265 pow2( interval( 6299691458188962.0 / Ne,6299691458189149.0 / Ne ) , 1016 ) };
1266 
1267 // ****************************************************************************
1268 // ******************** Gamma(x), 1/Gamma(x) **********************************
1269 // ****************************************************************************
1270 
1271 inline int round_g(const real& x) noexcept
1272 // Only for the internal use ( interval gammar(x) )
1273 // For |x| < 2147483647.5 the assignment y = round_g(x) delivers:
1274 // y = round_g(-0.1); ---> y = 1;
1275 // y = round_g(-0.5); ---> y = 1;
1276 // y = round_g(-1.0); ---> y = 1;
1277 // y = round_g(-1.4); ---> y = 2;
1278 // y = round_g(-6.8); ---> y = 7;
1279 // y = round_g(+0.0); ---> y = 0;
1280 // y = round_g(+0.1); ---> y = 0;
1281 // y = round_g(+0.5); ---> y = 0;
1282 // y = round_g(+2.6); ---> y = 0;
1283 // Blomquist, 25.06.2009;
1284 {
1285  int n = ifloor(_double(x));
1286  n = (n>=0)? 0:-n;
1287 
1288  return n;
1289 }
1290 
1291 real gamr_even_Ma(const real& x1, const real& x2, int n1)
1292 {
1293  real y;
1294 
1295  if ( x2<Inf(gam_rxi[n1]) || Sup(gam_rxi[n1])<x1 )
1296  { // [x1,x2] & gam_rxi[n1] is the empty set:
1297  y = (x1<Inf(gam_rxi[n1]))? gammar(x2) : gammar(x1);
1298  y *= q_gammarp;
1299  }
1300  else // [x1,x2] & gam_rxi[n1] is not the empty set:
1301  y = Sup(gam_ryi[n1]);
1302 
1303  return y;
1304 }
1305 
1306 real gamr_even_Mi(const real& x1, const real& x2, int n1)
1307 {
1308  real y,y1;
1309 
1310  if ( x2<Inf(gam_rxi[n1]) || Sup(gam_rxi[n1])<x1 )
1311  { // [x1,x2] & gam_rxi[n1] is the empty set:
1312  std::cout << "Leere Menge:" << std::endl;
1313  y = (x1<Inf(gam_rxi[n1]))? gammar(x1) : gammar(x2);
1314  y *= q_gammarm;
1315  }
1316  else // [x1,x2] & gam_rxi[n1] is not the empty set:
1317  {
1318  y1 = gammar(x1)*q_gammarm;
1319  y = gammar(x2)*q_gammarm;
1320  if (y1<y) y=y1;
1321  }
1322 
1323  return y;
1324 }
1325 
1326 real gamr_odd_Mi(const real& x1, const real& x2, int n1)
1327 {
1328  real y;
1329 
1330  if ( x2<Inf(gam_rxi[n1]) || Sup(gam_rxi[n1])<x1 )
1331  { // [x1,x2] & gam_rxi[n1] is the empty set:
1332  y = (x1<Inf(gam_rxi[n1]))? gammar(x2) : gammar(x1);
1333  y *= q_gammarp;
1334  }
1335  else // [x1,x2] & gam_rxi[n1] is not the empty set:
1336  y = Inf(gam_ryi[n1]);
1337 
1338  return y;
1339 }
1340 
1341 real gamr_odd_Ma(const real& x1, const real& x2, int n1)
1342 {
1343  real y,y1;
1344 
1345  if ( x2<Inf(gam_rxi[n1]) || Sup(gam_rxi[n1])<x1 )
1346  { // [x1,x2] & gam_rxi[n1] is the empty set:
1347  std::cout << "Leere Menge:" << std::endl;
1348  y = (x1<Inf(gam_rxi[n1]))? gammar(x1) : gammar(x2);
1349  y *= q_gammarm;
1350  }
1351  else // [x1,x2] & gam_rxi[n1] is not the empty set:
1352  {
1353  y1 = gammar(x1)*q_gammarm;
1354  y = gammar(x2)*q_gammarm;
1355  if (y1>y) y=y1;
1356  }
1357 
1358  return y;
1359 }
1360 
1362 // Calculating inclusions of 1/Gamma(x);
1363 // Blomquist, 01.07.09
1364 {
1365  interval y;
1366  real x0, x1(Inf(x)), x2(Sup(x)), y0,y1(0),y2;
1367  int n1,n2;
1368 
1369  n1 = round_g(x1);
1370  n2 = round_g(x2);
1371  if (x1==x2) // x is point interval
1372  if (x1==-n1) y2 = y1; // y = [0,0];
1373  else
1374  {
1375  y1 = gammar(x1);
1376  y2 = y1;
1377  // Lower bound y1, Upper bound y2:
1378  if (y1<0)
1379  {
1380  y1 = y1*q_gammarp;
1381  y2 = y2*q_gammarm;
1382  }
1383  else
1384  {
1385  y1 = y1*q_gammarm;
1386  y2 = y2*q_gammarp;
1387  }
1388  }
1389  else // x2>x1:
1390  {
1391  if (n1%2==0) // n1 even, i.e. n1=0,2,4,6,...
1392  {
1393  if (n1==n2)
1394  {
1395  y1 = gamr_even_Mi(x1,x2,n1);
1396  y2 = gamr_even_Ma(x1,x2,n1);
1397  }
1398  else
1399  if (n2==n1-1)
1400  {
1401  y1 = gamr_odd_Mi(-n2,x2,n2);
1402  y2 = gamr_even_Ma(x1,-n2,n1);
1403  }
1404  else // n2 <= n1-2
1405  {
1406  y1 = Inf(gam_ryi[n1-1]); // Minimum OK
1407  y2 = gamr_even_Ma(x1,-n1+1,n1);
1408  x0 = x2;
1409  if (x2>n1-3 && x2<0)
1410  x0 = n1-3;
1411  y0 = gamr_even_Ma(-n1+2,x0,n1-2);
1412  if (y0>y2) y2 = y0;
1413 
1414  if (n1==4 && n2==0)
1415  {
1416  y0 = gamr_even_Ma(0,x2,0);
1417  if (y0>y2) y2=y0;
1418  }
1419  }
1420  } // n1 even;
1421  else // n1 odd:
1422  {
1423  if (n1==n2)
1424  {
1425  y1 = gamr_odd_Mi(x1,x2,n1);
1426  y2 = gamr_odd_Ma(x1,x2,n1);
1427  }
1428  else
1429  if (n2==n1-1)
1430  {
1431  y1 = gamr_odd_Mi(x1,-n2,n1);
1432  y2 = gamr_even_Ma(-n2,x2,n2);
1433  }
1434  else
1435  if (n2==n1-2)
1436  {
1437  y1 = gamr_odd_Mi(x1,n1-1,n1);
1438  y2 = gamr_odd_Mi(-n1+2,x2,n1-2);
1439  if (y2<y1) y1 = y2; // Minimum calculated
1440  y2 = Sup( gam_ryi[n1-1] );
1441  }
1442  else // 0 <= n2 <= n1-3;
1443  { // Calculating the minimum y1:
1444  y1 = gamr_odd_Mi(x1,n1-1,n1);
1445  y2 = Inf( gam_ryi[n1-2] );
1446  if (y2<y1) y1 = y2; // Minimum y1 calculated
1447  // Now calculating the maximum y2:
1448  if (n1==3)
1449  {
1450  y2 = Sup( gam_ryi[n1-1]);
1451  y0 = gamr_even_Ma(0,x2,0);
1452  if (y0>y2) y2=y0;
1453  }
1454  else // n1 = 5,7,9,....
1455  y2 = Sup( gam_ryi[n1-1] );
1456  }
1457 
1458  } // n1 odd
1459  }
1460 
1461  y = interval(y1,y2);
1462  return y;
1463 }
1464 
1466 // Calculating inclusions of 1/Gamma(x);
1467 // Blomquist, 01.07.09
1468 {
1469  return 1/gammar(x);
1470 }
1471 
1472 
1473 } // namespace cxsc
1474 
cxsc::AbsMax
real AbsMax(const interval &x)
Computes the greatest absolute value .
Definition: interval.cpp:303
cxsc::power
cinterval power(const cinterval &z, int n) noexcept
Calculates .
Definition: cimath.cpp:1941
cxsc::MinReal
const real MinReal
Smallest normalized representable floating-point number.
Definition: real.cpp:62
cxsc::gamma
interval gamma(const interval &x)
The Gamma function.
Definition: imath.cpp:1465
cxsc::lnp1
cinterval lnp1(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:867
cxsc::sqrt1px2
cinterval sqrt1px2(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:1071
cxsc::minreal
const real minreal
Smallest positive denormalized representable floating-point number.
Definition: real.cpp:63
cxsc::expmx2
interval expmx2(const interval &x)
Calculates .
Definition: imath.cpp:192
cxsc::expx2
interval expx2(const interval &x)
Calculates .
Definition: imath.cpp:234
cxsc::MakeHexReal
const real & MakeHexReal(int sign, unsigned int expo, a_btyp manthigh, a_btyp mantlow)
Produces an IEEE 64-bit floating-point number from given binary coded parts of an IEEE 64-bit floatin...
Definition: real.cpp:52
cxsc::interval
The Scalar Type interval.
Definition: interval.hpp:55
cxsc::sqrt
cinterval sqrt(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:1007
cxsc::ln
cinterval ln(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:851
cxsc::expx2m1
interval expx2m1(const interval &x)
Calculates .
Definition: imath.cpp:300
cxsc::abs
ivector abs(const cimatrix_subv &mv) noexcept
Returns the absolute value of the matrix.
Definition: cimatrix.inl:737
cxsc::Pir_interval
const interval Pir_interval
Enclosure-Interval for .
Definition: interval.cpp:366
cxsc::sqrtx2m1
cinterval sqrtx2m1(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:1109
cxsc::erfc
interval erfc(const interval &a)
The complementary Gauss error function .
Definition: imath.cpp:361
cxsc::acoshp1
interval acoshp1(const interval &x)
Calculates .
Definition: imath.cpp:617
cxsc::exp
cinterval exp(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:159
cxsc::Round
int Round(const real &x) noexcept
Rouding to the next integer; |x| < 2147483647.5.
Definition: rmath.cpp:536
cxsc::erf
interval erf(const interval &a)
The Gauss error function .
Definition: imath.cpp:354
cxsc::ln_sqrtx2y2
interval ln_sqrtx2y2(const interval &x, const interval &y) noexcept
Calculates .
Definition: imath.cpp:581
cxsc::pow
cinterval pow(const cinterval &z, const interval &p) noexcept
Calculates .
Definition: cimath.cpp:2074
cxsc::sinpix_pi
interval sinpix_pi(const interval &x)
Calculates ;.
Definition: imath.cpp:655
cxsc
The namespace cxsc, providing all functionality of the class library C-XSC.
Definition: cdot.cpp:29
cxsc::sqrtp1m1
cinterval sqrtp1m1(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:1054
cxsc::ifloor
int ifloor(const real &x) noexcept
Rounding to the greates integer smaller or equal x; -2147483649 < x <= 2147483647....
Definition: rmath.cpp:570
cxsc::expm1
cinterval expm1(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:177
cxsc::times2pown
void times2pown(cinterval &x, int n) noexcept
Fast multiplication of reference parameter [z] with .
Definition: cimath.cpp:2059
cxsc::atan
cinterval atan(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:2938
cxsc::gammar
interval gammar(const interval &x)
The inverse Gamma function: 1/Gamma(x)
Definition: imath.cpp:1361
cxsc::sqrt1mx2
cinterval sqrt1mx2(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:1140
cxsc::sqr
cinterval sqr(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:3342
cxsc::Pi
interval Pi()
Enclosure-Interval for .
Definition: imath.cpp:566
cxsc::sqrtx2y2
interval sqrtx2y2(const interval &x, const interval &y) noexcept
Calculates .
Definition: imath.cpp:80
cxsc::real
The Scalar Type real.
Definition: real.hpp:114