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
35extern "C" {
36#ifndef rfcth_included
37#define rfcth_included
38#include "r_fcth.h"
39#endif
40}
41
42namespace cxsc {
43
44
45interval 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
53interval 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
63interval 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
80interval 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
101const real Delta_f = 2*minreal;
102const real q_sqrtp1m1m = 9007199254740984.0 / 9007199254740992.0;
103const real q_sqrtp1m1p = 4503599627370502.0 / 4503599627370496.0;
104//********************************************************************
105interval 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):
144real q_sqrtx2m1p(4503599627370501.0/4503599627370496.0); // (1+e(f))
145real 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):
164real q_sqrt1mx2p(4503599627370501.0/4503599627370496.0); // (1+e(f))
165real 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):
186const real q_exp_x2p = 4503599627370502.0 / 4503599627370496.0; // (1+e(f))
187const real q_exp_x2m = 9007199254740984.0 / 9007199254740992.0; // (1-e(f))
188 // Oberschranke expmx2_UB fuer |x| > expmx2_x0:
189const real expmx2_UB = 2.225073858507447856659E-308;
190const 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
231const real q_expx2_p = 4503599627370500.0 / 4503599627370496.0; // (1+e(f))
232const 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
255const real q_expx2m1_p = 4503599627370500.0 / 4503599627370496.0; // (1+e(f))
256const real q_expx2m1_m = 9007199254740985.0 / 9007199254740992.0; // (1-e(f))
257const real expx2m1_0 = comp(0.5,-510); // 2^(-511)
258
259void sqr2uv(const real&, real&, real&);
260
261real 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 -------------------------------
323static real q_lnp1m = 9007199254740986.0 / 9007199254740992.0; // 1-eps
324static real q_lnp1p = 4503599627370501.0 / 4503599627370496.0; // 1+eps
325// ---------------------------------------------------------------------------
326
327interval 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
354interval erf (const interval &a) { return j_erf(a); }
361interval 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
381inline a_intv _a_intv(const interval &x)
382{
383 return *((const a_intv *)(&x));
384}
385inline interval _interval(const a_intv &x)
386{
387 return *((const interval *)(&x));
388}
389
390interval 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//----------------------------------------------------------------------------
432static 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
468interval 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//----------------------------------------------------------------------------
536interval 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
566interval 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:
570real ln_x2y2_abs(2.225076E-308); // Absolute error bond
571real q_lnx2y2p(4503599627370502.0 / 4503599627370496.0); // 1+e(f)
572real q_lnx2y2m(9007199254740984.0 / 9007199254740992.0); // 1-e(f)
573// With the following b0
574// real b0 = 6369051672525773.0 / 30191699398572330817932436647906151127335369763331523427009650401964993299137190816689013801421270140331747000246110759198164677039398341060491474011461568349195162615808.0;
575real 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
581interval 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)):
613static const real q_acoshp1p(4503599627370503.0/4503599627370496.0);
614// (1-e(f))
615static 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
638const real q_sinpix_p = 4503599627370499.0 / 4503599627370496.0; // (1+e(f))
639const real q_sinpix_m = 9007199254740986.0 / 9007199254740992.0; // (1-e(f))
640
641real 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
648real 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) */
729const real q_gammarm = 9007199254740964.0 / 9007199254740992.0;
730/* q_gammarp = 1 + eps(gammar) */
731const real q_gammarp = 4503599627370510.0 / 4503599627370496.0;
732
733// Inclusion of the extremes of 1/gamma(x):
734interval pow2(const interval& x, int ex)
735{
736 interval y(x);
737 times2pown(y,ex);
738 return y;
739}
740
741const real Ne = 9007199254740992.0;
742const real Ne1 = 1125899906842624.0;
743const real Ne2 = 562949953421312.0;
744const real Ne3 = 281474976710656.0;
745const real Ne4 = 140737488355328.0;
746const real Ne5 = 70368744177664.0;
747const real Ne6 = 35184372088832.0;
748
749const 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):
1094const interval gam_ryi[171] = {
1095pow2( interval( 5085347089749720.0 / Ne,5085347089749823.0 / Ne ) , 1 ) ,
1096pow2( interval( -5082146609264467.0 / Ne,-5082146609264314.0 / Ne ) , -1 ) ,
1097pow2( interval( 7824158147621733.0 / Ne,7824158147621966.0 / Ne ) , -1 ) ,
1098pow2( interval( -5070842539852372.0 / Ne,-5070842539852221.0 / Ne ) , 1 ) ,
1099pow2( interval( 4593118780547419.0 / Ne,4593118780547576.0 / Ne ) , 3 ) ,
1100pow2( interval( -5333021955274733.0 / Ne,-5333021955274575.0 / Ne ) , 5 ) ,
1101pow2( interval( 7546574203185105.0 / Ne,7546574203185319.0 / Ne ) , 7 ) ,
1102pow2( interval( -6294628859031764.0 / Ne,-6294628859031469.0 / Ne ) , 10 ) ,
1103pow2( interval( 6045310252810166.0 / Ne,6045310252811273.0 / Ne ) , 13 ) ,
1104pow2( interval( -6568078652156336.0 / Ne,-6568078652156148.0 / Ne ) , 16 ) ,
1105pow2( interval( 7963169065060572.0 / Ne,7963169065060801.0 / Ne ) , 19 ) ,
1106pow2( interval( -5328217018030122.0 / Ne,-5328217018029960.0 / Ne ) , 23 ) ,
1107pow2( interval( 7800142897041864.0 / Ne,7800142897042089.0 / Ne ) , 26 ) ,
1108pow2( interval( -6199437664213474.0 / Ne,-6199437664213297.0 / Ne ) , 30 ) ,
1109pow2( interval( 5316470282961123.0 / Ne,5316470282961284.0 / Ne ) , 34 ) ,
1110pow2( interval( -4892929765135337.0 / Ne,-4892929765135165.0 / Ne ) , 38 ) ,
1111pow2( interval( 4810107119289947.0 / Ne,4810107119290088.0 / Ne ) , 42 ) ,
1112pow2( interval( -5030373421375086.0 / Ne,-5030373421374834.0 / Ne ) , 46 ) ,
1113pow2( interval( 5576144001185310.0 / Ne,5576144001185479.0 / Ne ) , 50 ) ,
1114pow2( interval( -6530685487420963.0 / Ne,-6530685487420774.0 / Ne ) , 54 ) ,
1115pow2( interval( 8057940169576582.0 / Ne,8057940169576818.0 / Ne ) , 58 ) ,
1116pow2( interval( -5223648494045513.0 / Ne,-5223648494045349.0 / Ne ) , 63 ) ,
1117pow2( interval( 7099855957135674.0 / Ne,7099855957135885.0 / Ne ) , 67 ) ,
1118pow2( interval( -5047359382236272.0 / Ne,-5047359382236084.0 / Ne ) , 72 ) ,
1119pow2( interval( 7492585872478835.0 / Ne,7492585872479188.0 / Ne ) , 76 ) ,
1120pow2( interval( -5795835662380422.0 / Ne,-5795835662380242.0 / Ne ) , 81 ) ,
1121pow2( interval( 4664800910382651.0 / Ne,4664800910382790.0 / Ne ) , 86 ) ,
1122pow2( interval( -7801058080117709.0 / Ne,-7801058080117472.0 / Ne ) , 90 ) ,
1123pow2( interval( 6767162072327001.0 / Ne,6767162072327282.0 / Ne ) , 95 ) ,
1124pow2( interval( -6082121514218736.0 / Ne,-6082121514218554.0 / Ne ) , 100 ) ,
1125pow2( interval( 5656800000052189.0 / Ne,5656800000052359.0 / Ne ) , 105 ) ,
1126pow2( interval( -5438268378952110.0 / Ne,-5438268378951951.0 / Ne ) , 110 ) ,
1127pow2( interval( 5398375606367166.0 / Ne,5398375606367329.0 / Ne ) , 115 ) ,
1128pow2( interval( -5527713447587841.0 / Ne,-5527713447587674.0 / Ne ) , 120 ) ,
1129pow2( interval( 5833125895912623.0 / Ne,5833125895912799.0 / Ne ) , 125 ) ,
1130pow2( interval( -6337936184674347.0 / Ne,-6337936184674153.0 / Ne ) , 130 ) ,
1131pow2( interval( 7084743510515278.0 / Ne,7084743510515501.0 / Ne ) , 135 ) ,
1132pow2( interval( -8141214882701327.0 / Ne,-8141214882701088.0 / Ne ) , 140 ) ,
1133pow2( interval( 4804968547193877.0 / Ne,4804968547194018.0 / Ne ) , 146 ) ,
1134pow2( interval( -5822137580509526.0 / Ne,-5822137580509355.0 / Ne ) , 151 ) ,
1135pow2( interval( 7236772755227956.0 / Ne,7236772755228162.0 / Ne ) , 156 ) ,
1136pow2( interval( -4610758665056508.0 / Ne,-4610758665056369.0 / Ne ) , 162 ) ,
1137pow2( interval( 6019530845699084.0 / Ne,6019530845699266.0 / Ne ) , 167 ) ,
1138pow2( interval( -8047036389398365.0 / Ne,-8047036389398123.0 / Ne ) , 172 ) ,
1139pow2( interval( 5504580189086749.0 / Ne,5504580189086968.0 / Ne ) , 178 ) ,
1140pow2( interval( -7703001513324420.0 / Ne,-7703001513324183.0 / Ne ) , 183 ) ,
1141pow2( interval( 5510183009440391.0 / Ne,5510183009440581.0 / Ne ) , 189 ) ,
1142pow2( interval( -8055535954952413.0 / Ne,-8055535954952173.0 / Ne ) , 194 ) ,
1143pow2( interval( 6014315232803007.0 / Ne,6014315232803294.0 / Ne ) , 200 ) ,
1144pow2( interval( -4584378555360492.0 / Ne,-4584378555360260.0 / Ne ) , 206 ) ,
1145pow2( interval( 7132212380084113.0 / Ne,7132212380084326.0 / Ne ) , 211 ) ,
1146pow2( interval( -5659549393054692.0 / Ne,-5659549393054526.0 / Ne ) , 217 ) ,
1147pow2( interval( 4579461117155838.0 / Ne,4579461117155977.0 / Ne ) , 223 ) ,
1148pow2( interval( -7554216840666713.0 / Ne,-7554216840666493.0 / Ne ) , 228 ) ,
1149pow2( interval( 6348787715758027.0 / Ne,6348787715758222.0 / Ne ) , 234 ) ,
1150pow2( interval( -5434979980476367.0 / Ne,-5434979980476204.0 / Ne ) , 240 ) ,
1151pow2( interval( 4737681191908824.0 / Ne,4737681191908967.0 / Ne ) , 246 ) ,
1152pow2( interval( -8407842664867513.0 / Ne,-8407842664867267.0 / Ne ) , 251 ) ,
1153pow2( interval( 7592052521188700.0 / Ne,7592052521188935.0 / Ne ) , 257 ) ,
1154pow2( interval( -6974119252551297.0 / Ne,-6974119252551090.0 / Ne ) , 263 ) ,
1155pow2( interval( 6515520808385677.0 / Ne,6515520808385874.0 / Ne ) , 269 ) ,
1156pow2( interval( -6188946869743481.0 / Ne,-6188946869743300.0 / Ne ) , 275 ) ,
1157pow2( interval( 5975502808844840.0 / Ne,5975502808845020.0 / Ne ) , 281 ) ,
1158pow2( interval( -5862842897072874.0 / Ne,-5862842897072704.0 / Ne ) , 287 ) ,
1159pow2( interval( 5843967448508660.0 / Ne,5843967448508828.0 / Ne ) , 293 ) ,
1160pow2( interval( -5916517001341501.0 / Ne,-5916517001341321.0 / Ne ) , 299 ) ,
1161pow2( interval( 6082464626325325.0 / Ne,6082464626325503.0 / Ne ) , 305 ) ,
1162pow2( interval( -6348157530347044.0 / Ne,-6348157530346858.0 / Ne ) , 311 ) ,
1163pow2( interval( 6724699799057619.0 / Ne,6724699799057843.0 / Ne ) , 317 ) ,
1164pow2( interval( -7228705737680202.0 / Ne,-7228705737679999.0 / Ne ) , 323 ) ,
1165pow2( interval( 7883493269720206.0 / Ne,7883493269720561.0 / Ne ) , 329 ) ,
1166pow2( interval( -8720834785364833.0 / Ne,-8720834785364561.0 / Ne ) , 335 ) ,
1167pow2( interval( 4891722644546351.0 / Ne,4891722644546502.0 / Ne ) , 342 ) ,
1168pow2( interval( -5564236710028970.0 / Ne,-5564236710028799.0 / Ne ) , 348 ) ,
1169pow2( interval( 6416191129172903.0 / Ne,6416191129173091.0 / Ne ) , 354 ) ,
1170pow2( interval( -7498890927628704.0 / Ne,-7498890927628487.0 / Ne ) , 360 ) ,
1171pow2( interval( 8881515552460572.0 / Ne,8881515552460999.0 / Ne ) , 366 ) ,
1172pow2( interval( -5328950915550370.0 / Ne,-5328950915550206.0 / Ne ) , 373 ) ,
1173pow2( interval( 6478093314396794.0 / Ne,6478093314397089.0 / Ne ) , 379 ) ,
1174pow2( interval( -7976303366065662.0 / Ne,-7976303366065426.0 / Ne ) , 385 ) ,
1175pow2( interval( 4972846688449830.0 / Ne,4972846688450017.0 / Ne ) , 392 ) ,
1176pow2( interval( -6278401907481090.0 / Ne,-6278401907480879.0 / Ne ) , 398 ) ,
1177pow2( interval( 8024854758356088.0 / Ne,8024854758356345.0 / Ne ) , 404 ) ,
1178pow2( interval( -5191277948909595.0 / Ne,-5191277948909444.0 / Ne ) , 411 ) ,
1179pow2( interval( 6797621462551740.0 / Ne,6797621462551941.0 / Ne ) , 417 ) ,
1180pow2( interval( -4503636668393666.0 / Ne,-4503636668393518.0 / Ne ) , 424 ) ,
1181pow2( interval( 6037997262493341.0 / Ne,6037997262493523.0 / Ne ) , 430 ) ,
1182pow2( interval( -8189485306115383.0 / Ne,-8189485306115130.0 / Ne ) , 436 ) ,
1183pow2( interval( 5617805845426844.0 / Ne,5617805845427124.0 / Ne ) , 443 ) ,
1184pow2( interval( -7795192616785187.0 / Ne,-7795192616784477.0 / Ne ) , 449 ) ,
1185pow2( interval( 5469175405734180.0 / Ne,5469175405734422.0 / Ne ) , 456 ) ,
1186pow2( interval( -7759929987383324.0 / Ne,-7759929987383086.0 / Ne ) , 462 ) ,
1187pow2( interval( 5565727978288701.0 / Ne,5565727978288876.0 / Ne ) , 469 ) ,
1188pow2( interval( -8070914994857895.0 / Ne,-8070914994857635.0 / Ne ) , 475 ) ,
1189pow2( interval( 5914931467943193.0 / Ne,5914931467943373.0 / Ne ) , 482 ) ,
1190pow2( interval( -8762204548045716.0 / Ne,-8762204548045455.0 / Ne ) , 488 ) ,
1191pow2( interval( 6558513517606168.0 / Ne,6558513517606353.0 / Ne ) , 495 ) ,
1192pow2( interval( -4960305627886271.0 / Ne,-4960305627886120.0 / Ne ) , 502 ) ,
1193pow2( interval( 7580642983583672.0 / Ne,7580642983583897.0 / Ne ) , 508 ) ,
1194pow2( interval( -5851844804194595.0 / Ne,-5851844804194367.0 / Ne ) , 515 ) ,
1195pow2( interval( 4563038858728436.0 / Ne,4563038858728577.0 / Ne ) , 522 ) ,
1196pow2( interval( -7187477492053316.0 / Ne,-7187477492052964.0 / Ne ) , 528 ) ,
1197pow2( interval( 5716852908386950.0 / Ne,5716852908387214.0 / Ne ) , 535 ) ,
1198pow2( interval( -4591808630269563.0 / Ne,-4591808630269411.0 / Ne ) , 542 ) ,
1199pow2( interval( 7448102539955649.0 / Ne,7448102539955986.0 / Ne ) , 548 ) ,
1200pow2( interval( -6098770429791387.0 / Ne,-6098770429791204.0 / Ne ) , 555 ) ,
1201pow2( interval( 5041550443966798.0 / Ne,5041550443966946.0 / Ne ) , 562 ) ,
1202pow2( interval( -8413996086583072.0 / Ne,-8413996086582821.0 / Ne ) , 568 ) ,
1203pow2( interval( 7086939987269423.0 / Ne,7086939987269731.0 / Ne ) , 575 ) ,
1204pow2( interval( -6024570065319942.0 / Ne,-6024570065319682.0 / Ne ) , 582 ) ,
1205pow2( interval( 5168535487082451.0 / Ne,5168535487082609.0 / Ne ) , 589 ) ,
1206pow2( interval( -8949051953781375.0 / Ne,-8949051953781115.0 / Ne ) , 595 ) ,
1207pow2( interval( 7817344426895164.0 / Ne,7817344426895996.0 / Ne ) , 602 ) ,
1208pow2( interval( -6889843867972878.0 / Ne,-6889843867972674.0 / Ne ) , 609 ) ,
1209pow2( interval( 6126229646423302.0 / Ne,6126229646423484.0 / Ne ) , 616 ) ,
1210pow2( interval( -5495122334906381.0 / Ne,-5495122334906222.0 / Ne ) , 623 ) ,
1211pow2( interval( 4971972094727164.0 / Ne,4971972094727314.0 / Ne ) , 630 ) ,
1212pow2( interval( -4537480959802395.0 / Ne,-4537480959802254.0 / Ne ) , 637 ) ,
1213pow2( interval( 8352835047353300.0 / Ne,8352835047353555.0 / Ne ) , 643 ) ,
1214pow2( interval( -7753443787904532.0 / Ne,-7753443787904298.0 / Ne ) , 650 ) ,
1215pow2( interval( 7257653550749169.0 / Ne,7257653550749382.0 / Ne ) , 657 ) ,
1216pow2( interval( -6850281165773769.0 / Ne,-6850281165773570.0 / Ne ) , 664 ) ,
1217pow2( interval( 6519305845448896.0 / Ne,6519305845449168.0 / Ne ) , 671 ) ,
1218pow2( interval( -6255266499085062.0 / Ne,-6255266499084872.0 / Ne ) , 678 ) ,
1219pow2( interval( 6050802311308162.0 / Ne,6050802311308350.0 / Ne ) , 685 ) ,
1220pow2( interval( -5900304762620398.0 / Ne,-5900304762620223.0 / Ne ) , 692 ) ,
1221pow2( interval( 5799657649647993.0 / Ne,5799657649648165.0 / Ne ) , 699 ) ,
1222pow2( interval( -5746047975553302.0 / Ne,-5746047975553134.0 / Ne ) , 706 ) ,
1223pow2( interval( 5737835419331524.0 / Ne,5737835419331693.0 / Ne ) , 713 ) ,
1224pow2( interval( -5774471890994117.0 / Ne,-5774471890993944.0 / Ne ) , 720 ) ,
1225pow2( interval( 5856465763387432.0 / Ne,5856465763387600.0 / Ne ) , 727 ) ,
1226pow2( interval( -5985387992102590.0 / Ne,-5985387992102406.0 / Ne ) , 734 ) ,
1227pow2( interval( 6163919695584074.0 / Ne,6163919695584257.0 / Ne ) , 741 ) ,
1228pow2( interval( -6395943042753787.0 / Ne,-6395943042753502.0 / Ne ) , 748 ) ,
1229pow2( interval( 6686679647283150.0 / Ne,6686679647283350.0 / Ne ) , 755 ) ,
1230pow2( interval( -7042883260256940.0 / Ne,-7042883260256730.0 / Ne ) , 762 ) ,
1231pow2( interval( 7473096566380533.0 / Ne,7473096566380749.0 / Ne ) , 769 ) ,
1232pow2( interval( -7987985534527481.0 / Ne,-7987985534527243.0 / Ne ) , 776 ) ,
1233pow2( interval( 8600769311605383.0 / Ne,8600769311605633.0 / Ne ) , 783 ) ,
1234pow2( interval( -4663884705694464.0 / Ne,-4663884705694325.0 / Ne ) , 791 ) ,
1235pow2( interval( 5094554684614484.0 / Ne,5094554684614634.0 / Ne ) , 798 ) ,
1236pow2( interval( -5604802840349871.0 / Ne,-5604802840349701.0 / Ne ) , 805 ) ,
1237pow2( interval( 6209951739735886.0 / Ne,6209951739736072.0 / Ne ) , 812 ) ,
1238pow2( interval( -6928963530888061.0 / Ne,-6928963530887851.0 / Ne ) , 819 ) ,
1239pow2( interval( 7785368708274196.0 / Ne,7785368708274423.0 / Ne ) , 826 ) ,
1240pow2( interval( -8808459126256481.0 / Ne,-8808459126256060.0 / Ne ) , 833 ) ,
1241pow2( interval( 5017412797579486.0 / Ne,5017412797579638.0 / Ne ) , 841 ) ,
1242pow2( interval( -5755173329981532.0 / Ne,-5755173329981361.0 / Ne ) , 848 ) ,
1243pow2( interval( 6646385258439176.0 / Ne,6646385258439444.0 / Ne ) , 855 ) ,
1244pow2( interval( -7727539896552529.0 / Ne,-7727539896552294.0 / Ne ) , 862 ) ,
1245pow2( interval( 4522473425691912.0 / Ne,4522473425692052.0 / Ne ) , 870 ) ,
1246pow2( interval( -5328812572761788.0 / Ne,-5328812572761623.0 / Ne ) , 877 ) ,
1247pow2( interval( 6320558000502691.0 / Ne,6320558000502885.0 / Ne ) , 884 ) ,
1248pow2( interval( -7546265781200776.0 / Ne,-7546265781200489.0 / Ne ) , 891 ) ,
1249pow2( interval( 4534316912522546.0 / Ne,4534316912522688.0 / Ne ) , 899 ) ,
1250pow2( interval( -5484491485989575.0 / Ne,-5484491485989407.0 / Ne ) , 906 ) ,
1251pow2( interval( 6676632315202014.0 / Ne,6676632315202302.0 / Ne ) , 913 ) ,
1252pow2( interval( -8180074398476253.0 / Ne,-8180074398476014.0 / Ne ) , 920 ) ,
1253pow2( interval( 5042989707083422.0 / Ne,5042989707083666.0 / Ne ) , 928 ) ,
1254pow2( interval( -6257379418480333.0 / Ne,-6257379418480019.0 / Ne ) , 935 ) ,
1255pow2( interval( 7813097673618694.0 / Ne,7813097673619043.0 / Ne ) , 942 ) ,
1256pow2( interval( -4908325621754370.0 / Ne,-4908325621754217.0 / Ne ) , 950 ) ,
1257pow2( interval( 6205346227418363.0 / Ne,6205346227418597.0 / Ne ) , 957 ) ,
1258pow2( interval( -7893590972392525.0 / Ne,-7893590972392227.0 / Ne ) , 964 ) ,
1259pow2( interval( 5051411882876310.0 / Ne,5051411882876506.0 / Ne ) , 972 ) ,
1260pow2( interval( -6504655602059905.0 / Ne,-6504655602059583.0 / Ne ) , 979 ) ,
1261pow2( interval( 8426810051054742.0 / Ne,8426810051054986.0 / Ne ) , 986 ) ,
1262pow2( interval( -5491407534973626.0 / Ne,-5491407534973452.0 / Ne ) , 994 ) ,
1263pow2( interval( 7199960218142557.0 / Ne,7199960218142768.0 / Ne ) , 1001 ) ,
1264pow2( interval( -4748178637044143.0 / Ne,-4748178637044000.0 / Ne ) , 1009 ) ,
1265pow2( interval( 6299691458188962.0 / Ne,6299691458189149.0 / Ne ) , 1016 ) };
1266
1267// ****************************************************************************
1268// ******************** Gamma(x), 1/Gamma(x) **********************************
1269// ****************************************************************************
1270
1271inline 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
1291real 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
1306real 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
1326real 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
1341real 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
The Scalar Type interval.
Definition interval.hpp:55
The Scalar Type real.
Definition real.hpp:114
The namespace cxsc, providing all functionality of the class library C-XSC.
Definition cdot.cpp:29
cinterval sqrtp1m1(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:1054
cinterval sqrt1mx2(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:1140
const real MinReal
Smallest normalized representable floating-point number.
Definition real.cpp:62
interval erf(const interval &a)
The Gauss error function .
Definition imath.cpp:354
interval gamma(const interval &x)
The Gamma function.
Definition imath.cpp:1465
int ifloor(const real &x) noexcept
Rounding to the greates integer smaller or equal x; -2147483649 < x <= 2147483647....
Definition rmath.cpp:570
int Round(const real &x) noexcept
Rouding to the next integer; |x| < 2147483647.5.
Definition rmath.cpp:536
cinterval power(const cinterval &z, int n) noexcept
Calculates .
Definition cimath.cpp:1941
interval sinpix_pi(const interval &x)
Calculates ;.
Definition imath.cpp:655
interval expx2(const interval &x)
Calculates .
Definition imath.cpp:234
cinterval ln(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:851
const real minreal
Smallest positive denormalized representable floating-point number.
Definition real.cpp:63
cinterval pow(const cinterval &z, const interval &p) noexcept
Calculates .
Definition cimath.cpp:2074
interval expmx2(const interval &x)
Calculates .
Definition imath.cpp:192
interval acoshp1(const interval &x)
Calculates .
Definition imath.cpp:617
const interval Pir_interval
Enclosure-Interval for .
Definition interval.cpp:366
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
interval Pi()
Enclosure-Interval for .
Definition imath.cpp:566
cinterval sqrtx2m1(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:1109
cinterval sqrt1px2(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:1071
cinterval exp(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:159
interval gammar(const interval &x)
The inverse Gamma function: 1/Gamma(x)
Definition imath.cpp:1361
interval erfc(const interval &a)
The complementary Gauss error function .
Definition imath.cpp:361
interval expx2m1(const interval &x)
Calculates .
Definition imath.cpp:300
interval ln_sqrtx2y2(const interval &x, const interval &y) noexcept
Calculates .
Definition imath.cpp:581
real AbsMax(const interval &x)
Computes the greatest absolute value .
Definition interval.cpp:303
cinterval expm1(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:177
ivector abs(const cimatrix_subv &mv) noexcept
Returns the absolute value of the matrix.
Definition cimatrix.inl:737
cinterval sqrt(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:1007
void times2pown(cinterval &x, int n) noexcept
Fast multiplication of reference parameter [z] with .
Definition cimath.cpp:2059
cinterval sqr(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:3342
cinterval lnp1(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:867
cinterval atan(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:2938
interval sqrtx2y2(const interval &x, const interval &y) noexcept
Calculates .
Definition imath.cpp:80