C-XSC - A C++ Class Library for Extended Scientific Computing 2.5.4
rmath.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: rmath.cpp,v 1.30 2014/01/30 17:23:48 cxsc Exp $ */
25
26#include "rmath.hpp"
27
28#include "rtsrmath.h"
29
30namespace cxsc {
31// Constants and values for the function ln_sqrtx2y2:
32real ln2_1067(6505485212531678.0/8796093022208.0); // 1067*ln(2)
33// Exponents of the interval bounds:
34int B_lnx2y2_1[22] = {21, 31, 51, 101, 151, 201, 251, 301, 351, 401,
35 451, 501, 551, 601, 651, 701, 751, 801, 851,
36 901, 951, 1025};
37int B_lnx2y2_2[22] = {-1021,-949,-899,-849,-799,-749,-699,-649,-599,
38 -549,-499,-449,-399,-349,-299,-249,-199,-149,
39 -99,-49,-29,-19};
40// Optimal values N for N*ln(2):
41int B_lnx2y2_N1[21] ={20, 40, 61, 122, 160, 229, 259, 320, 366, 427,
42 488, 518, 549, 610, 671, 732, 763, 825, 885,
43 945, 976};
44real B_lnx2y2_c1[21] =
45{
46 7804143460206699.0 / 562949953421312.0, // N*ln(2) with N = 20;
47 7804143460206699.0 / 281474976710656.0, // N*ln(2) with N = 40;
48 5950659388407608.0 / 140737488355328.0, // N*ln(2) with N = 61;
49 5950659388407608.0 / 70368744177664.0, // N*ln(2) with N = 122;
50 7804143460206699.0 / 70368744177664.0, // N*ln(2) with N = 160;
51 5584840163710419.0 / 35184372088832.0, // N*ln(2) with N = 229;
52 6316478613104797.0 / 35184372088832.0, // N*ln(2) with N = 259;
53 7804143460206699.0 / 35184372088832.0, // N*ln(2) with N = 320;
54 8925989082611412.0 / 35184372088832.0, // N*ln(2) with N = 366;
55 5206826964856657.0 / 17592186044416.0, // N*ln(2) with N = 427;
56 5950659388407608.0 / 17592186044416.0, // N*ln(2) with N = 488;
57 6316478613104797.0 / 17592186044416.0, // N*ln(2) with N = 518;
58 6694491811958559.0 / 17592186044416.0, // N*ln(2) with N = 549;
59 7438324235509510.0 / 17592186044416.0, // N*ln(2) with N = 610;
60 8182156659060461.0 / 17592186044416.0, // N*ln(2) with N = 671;
61 8925989082611412.0 / 17592186044416.0, // N*ln(2) with N = 732;
62 4652001140732587.0 / 8796093022208.0, // N*ln(2) with N = 763;
63 5030014339586349.0 / 8796093022208.0, // N*ln(2) with N = 825;
64 5395833564283538.0 / 8796093022208.0, // N*ln(2) with N = 885;
65 5761652788980727.0 / 8796093022208.0, // N*ln(2) with N = 945;
66 5950659388407608.0 / 8796093022208.0 // N*ln(2) with N = 976;
67};
68
69int uint_trail(const unsigned int& n)
70// Liefert die Anzahl aufeinanderfolgender Null-Bits am Ende von n;
71// n=0 --> 32; n=1 --> 0; n=2 --> 1; n=4 --> 2; n=1024 --> 10;
72{
73 int m1,p=0;
74 unsigned int m=n;
75 if (m==0) p = 32;
76 else { // m != 0, d.h. die folgende while-Schleife bricht ab!
77 do
78 {
79 m1 = m & 1; // m1=0(1), wenn letztes Bit von m gleich 0(1) ist.
80 if (m1==0) p++;
81 m = m >> 1; // Bit-Muster um 1 nach rechts schieben.
82 } while(m1==0);
83 }
84 return p;
85} // uint_trail
86
87void sqr2uv(const real& x, real& u, real& v)
88// Liefert u,v für: x2 = u + v; EXAKTE Darstellung, falls kein overflow
89// auftritt und v im normalisierten Bereich liegt. u > |v|
90// Vorsicht: Funktioniert zunächst nur auf INTEL-Systemen!!!
91{
92 real a,b,t,y1,y2;
93 a = Cut26(x);
94 b = x-a; // x = a+b;
95 u = x*x;
96 t = u - a*a; // exakte Auswertung!
97 y2 = a*b;
98 times2pown(y2,1); // y2 = 2*a*b, exakt!
99 t -= y2;
100// Jetzt fehlt noch: t2 - b*b, aber b*b wird nicht immer korrekt berechnet,
101// daher nochmalige Aufspaltung von b in y1+y2!!
102 y1 = Cut25(b);
103 y2 = b - y1; // b = y1+y2, exakt;
104 t -= y1*y1;
105 if (sign(y2)!=0)
106 {
107 a = y1*y2;
108 times2pown(a,1); // a = 2*y1*y2, exakt!
109 t -= a;
110 t -= y2*y2;
111 }
112 v = -t;
113} // sqr2uv
114
115
116//------------------------------------------------------------------
117
118// real-Konstante expmx2_x0 fuer die Funktion e^(-x2);
119const real expmx2_x0 = 7491658466053896.0 / 281474976710656.0;
120// Fuer x > expmx2_x0 werden die Funktionswerte auf Null gesetzt.
121// Die relative Fehlerschranke e(f) := 4.618919E-16 gilt fuer
122// alle |x| <= expmx2_x0 = 26.61571750925....
123
124real expmx2(const real& x) noexcept
125// e^(-x^2); rel. Fehlerschranke: eps = 4.618919E-16 = e(f) gilt
126// fuer alle |x| <= expmx2_x0 = 26.61571750925....
127// Fuer |x| > expmx2_x0 --> expmx2(x) = 0;
128// Blomquist, 05.07.04;
129{
130 real t=x,u,v,res=0;
131 int ex;
132 if (t<0) t = -t; // t >= 0;
133 ex = expo(t);
134 if (ex<=-26) res = 1; // t < 2^(-26)
135 else if (ex<=-6) // t < 2^(-6)
136 {
137 u = t*t; v = u;
138 times2pown(v,-1); // v: 0.5*x2
139 res = 1-u*( 1-v*(1-u/3) );
140 } else if (t<=expmx2_x0) {
141 sqr2uv(x,u,v); // u:= x*x,v aus S(2,53); x2 = u+v (exakt!)
142 res = exp(-u);
143 if (v!=0) {
144 times2pown(res,500); // Die Skalierung verhindert, dass
145 v *= res; // v*exp(-u) in den denormalisierten Bereich faellt
146 res -= v;
147 times2pown(res,-500); // Rueckskalierung
148 }
149 }
150 return res;
151} // expmx2
152
153real expx2(const real& x)
154// e^(+x^2); rel. Fehlerschranke: eps = 4.618958E-16 = e(f) gilt
155// fuer alle |x| <= x0 = 26.64174755704632....
156// x0 = 7498985273150791.0 / 281474976710656.0;
157// Fuer |x| > x0 --> Programmabbruch
158// Ausfuehrlich getestet; Blomquist, 26.07.06;
159{
160 real t=x,u,v,res;
161 int ex;
162 if (t<0) t = -t; // t >= 0;
163 ex = expo(t);
164 if (ex<=-26) res = 1; // t < 2^(-26)
165 else
166 if (ex<=-6) // t < 2^(-6)
167 {
168 u = t*t; v = u;
169 times2pown(v,-1); // v: 0.5*x^2
170 res = 1+u*( 1+v*(1+u/3) );
171 }
172 else
173 {
174 sqr2uv(x,u,v); // u := x*x und v aus S(2,53);
175 // x^2 = u+v ist exakt!
176 res = exp(u);
177 v *= res; // v*exp(+u)
178 res += v;
179 }
180 return res;
181} // expx2
182
184// e^(+x^2)-1; rel. Fehlerschranke: eps = 4.813220E-16 = e(f) gilt
185// fuer alle x, mit: 2^(-511) <= x <= x0 = 26.64174755704632....
186// x0 = 7498985273150791.0 / 281474976710656.0;
187// Fuer x > x0 --> Programmabbruch wegen Overflow;
188// Fuer 0 < x < 2^(-511) --> Programmabbruch, denorm. Bereich!!
189// Ausfuehrlich getestet; Blomquist, 10.08.2006;
190{
191 real t(x),u,v,y,res(0);
192 int ex;
193 if (t<0) t = -t; // t >= 0;
194
195 if (t>=6.5) res = expx2(t);
196 else
197 {
198 ex = expo(t);
199 sqr2uv(x,u,v); // u := x*x und v aus S(2,53);
200 if (ex>=2) // g4(x)
201 {
202 y = exp(u);
203 res = 1 - v*y;
204 res = y - res;
205 }
206 else
207 if (ex>=-8) res = expm1(u) + v*exp(u); // g3(x)
208 else
209 if (ex>=-25) { // g2(x)
210 y = u*u;
211 times2pown(y,-1);
212 res = (1+u/3)*y + u;
213 }
214 else
215 if(ex>=-510) res = u; // g1(x)
216 else if (ex>=-1073)
217 {
218 std::cerr << "expx2m1: denormalized range!"
219 << std::endl; exit(1);
220 }
221 }
222
223 return res;
224} // expx2m1
225
226
227
228//------------------------------------------------------------------
229
230real sqrt1px2(const real& x) noexcept
231// sqrt(1+x^2); Blomquist 13.12.02;
232{
233 if (expo(x) > 33) return abs(x);
234 else return sqrt(1+x*x);
235}
236
237// Folgende Konstante sqrtp1m1_s wird gebraucht in
238// real sqrtp1m1(const real& x) noexcept; Blomquist, 05.08.03
239const real sqrtp1m1_s = 9007199254740984.0 / 1125899906842624.0;
240
241real sqrtp1m1(const real& x) noexcept
242// sqrtp1m1(x) = sqrt(x+1)-1;
243// Blomquist, 05.08.03;
244{
245 real y = x;
246 int ex = expo(x);
247 if (ex<=-50) times2pown(y,-1); // |x|<2^(-50); fast division by 2
248 else if (ex>=105) y = sqrt(x); // x >= 2^(+104) = 2.02824...e+31
249 else if (ex>=53) y = sqrt(x)-1; // x >= 2^(+52) = 4.50359...e+15
250 else if (x>-0.5234375 && x<=sqrtp1m1_s ) y = x / (sqrt(x+1) + 1);
251 else y = sqrt(x+1)-1;
252 return y;
253}
254
255
256//------------------------------------------------------------------
257
258real sqrtx2y2(const real& x, const real& y) noexcept
259// calculating sqrt(x^2 + y^2) in high accuracy. Blomquist 01.12.02
260{
261 real a,b,r;
262 dotprecision dot;
263 int exa,exb,ex;
264 a = x; b = y;
265 exa = expo(a); exb = expo(b);
266 if (exb > exa)
267 { // Permutation of a,b:
268 r = a; a = b; b = r;
269 ex = exa; exa = exb; exb = ex;
270 } // |a| >= |b|
271 // a = |a| >= |b| > 0:
272 ex = 511 - exa; // scaling a with 2^ex --> expo(a) == 511, --> a*a and
273 times2pown(a,ex); // b*b without overflow. An underflow of b*b will not
274 times2pown(b,ex); // affect the accuracy of the result!
275 dot = 0;
276 accumulate(dot,a,a);
277 accumulate(dot,b,b);
278 r = rnd(dot);
279 r = sqrt(r); // sqrt(...) declared in rmath.hpp
280 times2pown(r,-ex); // back-scaling
281 return r;
282} // sqrtx2y2
283
284//------------------------------------------------------------------
285
287// sqrt(x^2-1); rel. Fehlerschranke: eps = 2.221305E-16 = e(f)
288// Blomquist, 13.04.04;
289{ const real c1 = 1.000732421875,
290 c2 = 44000.0,
291 c3 = 1024.0; // c1,c2,c3 werden exakt gespeichert!
292 real res,ep,ep2,s1,s2,x1,x2,arg=x;
293 if (sign(arg)<0) arg = -arg; // arg = |x| >= 0
294 if (arg <= c1) { // x = 1+ep; x^2-1 = 2*ep + ep^2;
295 ep = x - 1; // Differenz rundungsfehlerfrei!
296 ep2 = ep*ep; // ep2 i.a. fehlerbehaftet!
297 times2pown(ep,1); // ep = 2*ep;
298 res = sqrt(ep+ep2); // res=y0: Startwert;
299 // x - y0^2 = (2*eps - s1^2) + [eps^2 - s2*(y0 + s1)]
300 s1 = Cut26(res); s2 = res - s1; // Startwert y0 = s1 + s2;
301 arg = ep - s1*s1; // arg = 2*eps - s1^2;
302 arg += (ep2 - s2*(res+s1)); // arg = x - y0^2
303 if (sign(arg)>0) {
304 arg = arg / res;
305 times2pown(arg,-1);
306 res += arg; // 1. Newton-Schritt beendet; eps = 2.221261E-16
307 }
308 } else if (arg<c2) {
309 // x-y0^2 = [(x1^2-1)-s1^2] + [x2*(x+x1)-s2*(y0+s1)]
310 x1 = Cut26(arg); x2 = arg - x1; // arg = x = x1 + x2;
311 ep2 = x2*(arg+x1); // ep2 ist fehlerbehaftet
312 x2 = x1*x1; ep = x2-1;
313 res = sqrt(ep+ep2); // res ist Startwert für Newton-Verfahren
314 s1 = Cut26(res); s2 = res - s1; // Startwert = s1 + s2;
315 ep2 = ep2 - s2 * (res+s1); // ep2 = [x2*(x+x1)-s2*(y0+s1)]
316 if (arg<c3) ep -= s1*s1; // ep = (x1^2-1) - s1^2;
317 else {
318 x2 -= s1*s1; // x2 = x1^2-s1^2
319 ep = x2 - 1; } // ep = (x1^2-s1^2) - 1;
320 ep += ep2; // ep = x - y0^2;
321 ep /= res;
322 times2pown(ep,-1);
323 res = res + ep; // 1. Newton-Schritt in hoher Genauigkeit
324 // beendet; eps = 2.221305E-16
325 } else { // arg = |x| >= 44000;
326 res = -1/arg;
327 times2pown(res,-1); // Multiplikation mit 0.5;
328 res += arg; // res = x - 1/(2*x); eps = 2.221114E-16
329 }
330 return res;
331} // sqrtx2m1 (Punktargumente)
332
333//------------------------------------------------------------------
334
336// sqrt(1-x2); rel. Fehlerschranke: eps = 3.700747E-16 = e(f)
337// Blomquist, 19.06.04;
338{
339 real t=x,res;
340 int ex;
341 if (sign(t)<0) t = -t; // Argument t >=0;
342 if (t>1) cxscthrow(STD_FKT_OUT_OF_DEF("real sqrt1mx2(const real&)"));
343 // For argument t now it holds: 0 <= t <=1;
344 ex = expo(t);
345 if (ex<=-26) res = 1; // t < 2^(-26) --> res = 1
346 else if (ex<=-15) { // t < 2^(-15) --> res = 1-x2/2
347 res = t*t;
348 times2pown(res,-1);
349 res = 1 - res;
350 } else {
351 if (ex>=0)
352 { // ex>=0 --> t>=0.5;
353 res = 1-t; // res: delta = 1-t;
354 t = res * res;
355 times2pown(res,1); // res: 2*delta
356 res = res - t; // res: 1-x2 = 2*delta - delta2
357 } else res = 1-t*t; // res: Maschinenwert von 1-x2
358 res = sqrt(res); // res: Nullte Naeherung
359 }
360 return res;
361} // sqrt1mx2
362
363
364//------------------------------------------------------------------
365
366int Interval_Nr(int* v, const int& n, const int& ex)
367// n>0 subintervals: |...\|...\|...\ ..... |...|
368// subinterval Nr.: 0 1 2 ..... n-1
369{
370 int i=0,j=n,k; // n>0: Number of subintervals
371 do {
372 k = (i+j)/2;
373 if (ex < v[k]) j = k-1;
374 else i = k+1;
375 } while(i<=j);
376 return j; // x with ex=expo(x) lies in the subinterval number j
377}
378
379//------------------------------------------------------------------
380
381real ln_sqrtx2y2(const real& x, const real& y)
382
383// ln( sqrt(x^2+y^2) ) == 0.5*ln(x^2+y^2); Blomquist, 21.11.03;
384// Relative error bound: 5.160563E-016;
385// Absolute error bound: 2.225075E-308; if x=1 and 0<=y<=b0;
386{
387 int j,N;
388 real a,b,r,r1;
389 dotprecision dot;
390
391 a = sign(x)<0 ? -x : x; // a = |x| >= 0;
392 b = sign(y)<0 ? -y : y; // b = |y| >= 0;
393 int exa=expo(a), exb=expo(b), ex;
394 if (b > a)
395 {
396 r = a; a = b; b = r;
397 ex = exa; exa = exb; exb = ex;
398 }
399 // It holds now: 0 <= b <= a
400 if (sign(a)==0)
401 cxscthrow(STD_FKT_OUT_OF_DEF
402 ("real ln_sqrtx2y2(const real&, const real&)"));
403 if (exa>20) // to avoid overflow by calculating a^2 + b^2
404 { // a>=2^(20):
405 j = Interval_Nr(B_lnx2y2_1,21,exa); // j: No. of subinterval
406 N = B_lnx2y2_N1[j]; // N: Optimal int value
407 if (exb-exa > -25)
408 { // For (exb-exa>-25) we use the complete term:
409 // N*ln(2) + [ln(2^(-N)*a)+0.5*ln(1+(b/a)^2)]
410 b = b/a; // a > 0
411 b = lnp1(b*b);
412 times2pown(b,-1); // exact division by 2
413 times2pown(a,-N);
414 r = b + ln(a); // [ ... ] calculated!
415 r += B_lnx2y2_c1[j];
416 }
417 else { // For (exb-exa<=-25) only two summands!:
418 times2pown(a,-N);
419 r = ln(a) + B_lnx2y2_c1[j];
420 }
421 }
422 else // exa<=20 or a<2^(20):
423 { // Now calculation of a^2+b^2 without overflow:
424 if (exa<=-20) // to avoid underflow by calculating a^2+b^2
425 if (exa<=-1022) // a in the denormalized range
426 {
427 r = b/a;
428 r = lnp1(r*r); times2pown(r,-1); // r: 0.5*ln(1+..)
429 times2pown(a,1067);
430 r += ln(a); // [ .... ] ready
431 r -= ln2_1067; // rel. error = 2.459639e-16;
432 }
433 else // MinReal=2^(-1022) <= a < 2^(-20)
434 { // Calculating the number j of the subinterval:
435 j = 20 - Interval_Nr(B_lnx2y2_2,21,exa);
436 r = a; times2pown(r,B_lnx2y2_N1[j]);
437 r = ln(r); // r: ln(2^N*a);
438 if (exb-exa > -25) { // calculating the complete term
439 b = b/a;
440 a = lnp1(b*b);
441 times2pown(a,-1);
442 r += a; // [ ... ] ready now
443 }
444 // We now have: exb-exa<=-25, ==> b/a <= 2^(-24);
445 r -= B_lnx2y2_c1[j]; // 0.5*ln(1+(b/a)^2) neglected!
446 // relative error = 4.524090e-16 in both cases;
447 }
448 else // calculation of a^2+b^2 without overflow or underflow:
449 { // exa>-20 respective a>=2^(-20):
450 dot = 0;
451 accumulate(dot,a,a);
452 accumulate(dot,b,b); // dot = a^2+b^2, exact!
453 real s = rnd(dot); // s = a^2 + b^2, rounded!
454 if (s>=0.25 && s<=1.75)
455 if (s>=0.828125 && s<=1.171875)
456 { // Series:
457 if (a==1 && exb<=-28)
458 {
459 r = b; times2pown(r,-1);
460 r *= b;
461 }
462 else {
463 dot -= 1;
464 r = rnd(dot); // r = a^2+b^2-1 rounded!
465 r = lnp1(r);
466 times2pown(r,-1);
467 }
468 }
469 else { // Reading dot = a^2+b^2 twice:
470 r = rnd(dot);
471 dot -= r;
472 r1 = rnd(dot); // a^2+b^2 = r+r1, rounded!
473 r1 = lnp1(r1/r);
474 r = ln(r) + r1;
475 times2pown(r,-1); // exact division by 2
476 }
477 else { // calculating straight from: 0.5*ln(x^2+y^2)
478 r = ln(s);
479 times2pown(r,-1);
480 }
481 }
482 }
483 return r;
484} // ln_sqrtx2y2
485
486typedef union { double f; char intern[8]; } help_real;
487
488real Cut24(const real& x){
489 // y = Cut24(x) liefert ein y, das mit den ersten 24 Mantissenbits
490 // von x übereinstimmt, das hidden bit ist dabei mitgezählt!
491 // Die restlichen 53-24=29 Mantissenbits werden auf Null gesetzt.
492 help_real y;
493 y.f = _double(x);
494#if INTEL
495 y.intern[3] &= 224;
496 y.intern[0] = y.intern[1] = y.intern[2] = 0;
497#else
498 y.intern[4] &= 224;
499 y.intern[7] = y.intern[6] = y.intern[5] = 0;
500#endif
501 return real(y.f);
502}
503
504real Cut25(const real& x){
505 // y = Cut25(x) liefert ein y, das mit den ersten 25 Mantissenbits
506 // von x übereinstimmt, das hidden bit ist dabei mitgezählt!
507 // Die restlichen 53-25=28 Mantissenbits werden auf Null gesetzt.
508 help_real y;
509 y.f = _double(x);
510#if INTEL
511 y.intern[3] &= 240;
512 y.intern[0] = y.intern[1] = y.intern[2] = 0;
513#else
514 y.intern[4] &= 240;
515 y.intern[7] = y.intern[6] = y.intern[5] = 0;
516#endif
517 return real(y.f);
518}
519
520real Cut26(const real& x){
521 // y = Cut26(x) liefert ein y, das mit den ersten 26 Mantissenbits
522 // von x übereinstimmt, das hidden bit ist dabei mitgezählt!
523 // Die restlichen 53-26=27 Mantissenbits werden auf Null gesetzt.
524 help_real y;
525 y.f = _double(x);
526#if INTEL
527 y.intern[3] &= 248;
528 y.intern[0] = y.intern[1] = y.intern[2] = 0;
529#else
530 y.intern[4] &= 248;
531 y.intern[7] = y.intern[6] = y.intern[5] = 0;
532#endif
533 return real(y.f);
534}
535
536int Round(const real& x) noexcept
537// y = Round(x) delivers the rounded value y of x.
538// For |x| < 2147483647.5 the assignment y = Round(x) delivers
539// the next integer value y.
540// Examples:
541// y = Round(-0.1); ---> y = 0;
542// y = Round(-123.5); ---> y = -124.0;
543// y = Round(+123.5); ---> y = +124.0;
544// y = Round(+123.499); ---> y = +123.0;
545// y = Round(+2147483647.499); ---> y = +2147483647.0;
546// y = Round(+2147483647.5); ---> y = -2147483648.0;
547// y = Round(-2147483647.499); ---> y = -2147483647.0;
548// y = Round(-2147483647.5); ---> y = -2147483648.0;
549// y = Round(-2147483648.501); ---> y = -2147483648.0; wrong rounding!
550// Blomquist, 29.10.2008;
551{
552 double dbl;
553
554 dbl = _double(x);
555 return dbl<0 ? int(dbl-0.5) : int(dbl+0.5);
556}
557
558int ceil(const real& x) noexcept
559{
560 real y(x);
561 bool neg(y<0);
562 if (neg) y = -y; // y >= 0;
563 int res = int( _double(y) );
564 y = y - res;
565 if (!neg && y>0) res += 1;
566 if (neg) res = -res;
567 return res;
568}
569
570int ifloor(const real& x) noexcept
571{
572 real y(x);
573 bool neg(y<0);
574 if (neg) y = -y; // y >= 0;
575 int res = int( _double(y) );
576 y = y - res;
577 if (neg)
578 {
579 res = -res;
580 if (y>0) res -= 1;
581 }
582 return res;
583}
584
585static real q_acoshp1[5] = // Polynomial coefficients of Q_4(x)
586 // roundet to nearest. acosh(1+x) = sqrt(2*x)*Q_4(x)
587{ 1.0 / 1.0, // q_acoshp1[0]
588 -6004799503160661.0 / 72057594037927936.0,
589 +5404319552844595.0 / 288230376151711744.0,
590 -6433713753386423.0 / 1152921504606846976.0,
591 +8756999275442631.0 / 4611686018427387904.0 // q_acoshp1[4]
592};
593
594static const real c_ln2_B = 6243314768165359.0 / 9007199254740992.0;
595// c_ln2_B < ln(2) is the nearest machine number for ln(2) with an
596// absolute error < 2.3190469E-17;
597
598real acoshp1(const real& x) noexcept
599// acoshp1(x) = acosh(1+x); rel. error: eps = 7.792706E-16 = e(f)
600// Ausfuehrlich getestet; Blomquist, 27.03.05;
601{
602 real res;
603 int ex(expo(x));
604 if (x<0)
605 cxscthrow(STD_FKT_OUT_OF_DEF("real acoshp1(const real&)"));
606 // For argument x now it holds: 0 <= x <= MaxReal;
607 if (ex<=-50) res = sqrt(2*x); // 0<=x<2^(-50): acoshp1(x)=sqrt(2x)
608 else if (ex<=-9) // 2^(-50)<=x<2^{-9}: acoshp1(x)=sqrt(2x)*Q_4(x)
609 res = sqrt(2*x)*((((q_acoshp1[4]*x+q_acoshp1[3])*x+q_acoshp1[2])
610 *x+q_acoshp1[1])*x + q_acoshp1[0]);
611 else if (ex<=0) res = lnp1(x+sqrt(2*x+x*x)); // range A_3
612 else if (ex<=50) res = lnp1(x*(1+sqrt(1+2/x))); // range A_4
613 else if (ex<=1022) res = ln(2*x); // range A_5
614 else res = ln(x) + c_ln2_B; // range A_6
615
616 return res;
617} // acoshp1
618
619
620// *****************************************************************************
621// sin(Pi*x)/Pi *
622// ********************************************************************************
623// Sinpix_pi: Teilpunkte des Intervalls [0,0.5]:
624// ********************************************************************************
625real a_sinpix_pi[8] = {0,
626 7.450580596923828125e-9, // 2^(-27)
627 0.05078125,
628 0.1015625,
629 0.15234375,
630 0.2578125,
631 0.3828125,
632 0.5 };
633// ********************************************************************************
634
635// ********************************************************************************
636// Sinpix_pi: a_k,b_k des Kettenbruchs K2 in A1 = [2^-27,0.05078125]:
637// ********************************************************************************
638 const real q_sin1_a[3] = {0.0,
639 -7408124450506663.0 / 4503599627370496.0,
640 +4595816915236340.0 / 36028797018963968.0 };
641
642 const real q_sin1_b[3] = {+1.0,
643 +8889749340277122.0 / 18014398509481984.0,
644 -6103596185201565.0 / 36028797018963968.0};
645// ********************************************************************************
646
647// ********************************************************************************
648// Sinpix_pi: a_k,b_k des Kettenbruchs K4 in A2 = [0.05078125,0.1015625]:
649// ********************************************************************************
650 const real q_sin2_a[5] = {0.0,
651 -4829370543434630.0 / 18014398509481984.0,
652 +5229154385037560.0 / 140737488355328.0,
653 -6000119407941526.0 / 576460752303423488.0,
654 +4592842702840413.0 / 1125899906842624.0 };
655
656 const real q_sin2_b[5] = {-6359670668682560.0 / 576460752303423488.0,
657 -6771299865719208.0 / 1125899906842624.0,
658 +6861350950372367.0 / 1125899906842624.0,
659 -4669728680615660.0 / 2251799813685248.0,
660 +4607746187351458.0 / 2251799813685248.0 };
661// ********************************************************************************
662
663// ********************************************************************************
664// Sinpix_pi: a_k,b_k des Kettenbruchs K4 in A3 = [0.1015625,0.15234375]:
665// ********************************************************************************
666 const real q_sin3_a[5] = {0.0,
667 -7954137925457403.0 / 18014398509481984.0,
668 +7534721222201852.0 / 562949953421312.0,
669 -8481413201661564.0 / 288230376151711744.0,
670 +6484759126571114.0 / 4503599627370496.0};
671 const real q_sin3_b[5] = {-8780869688499296.0 / 288230376151711744.0,
672 -7929689932517363.0 / 2251799813685248.0,
673 +8223190551048856.0 / 2251799813685248.0,
674 -5789057670884168.0 / 4503599627370496.0,
675 +5586166285064551.0 / 4503599627370496.0 };
676// ********************************************************************************
677
678// ********************************************************************************
679// Sinpix_pi: a_k,b_k des Kettenbruchs K5 in A4 = [0.15234375,0.2578125]:
680// ********************************************************************************
681 const real q_sin4_a[6] = {+0.0,
682 -5469389903424399.0 / 9007199254740992.0,
683 +7704524171392577.0 / 1125899906842624.0,
684 -8513632445057187.0 / 144115188075855872.0,
685 +6447776473700967.0 / 9007199254740992.0,
686 -7362728796775514.0 / 288230376151711744.0 };
687
688 const real q_sin4_b[6] = {-8529339040415872.0 / 144115188075855872.0,
689 -5452397600582528.0 / 2251799813685248.0,
690 +5848879160432416.0 / 2251799813685248.0,
691 -8609222354677976.0 / 9007199254740992.0,
692 +8031029418839465.0 / 9007199254740992.0,
693 -4904826237544497.0 / 9007199254740992.0 };
694// ********************************************************************************
695
696// ********************************************************************************
697// Sinpix_pi: a_k,b_k des Kettenbruchs K5 in A5 = [0.2578125,0.3828125]:
698// ********************************************************************************
699 const real q_sin5_a[6] = {0.0,
700 4*(5276398025110506.0 / 9007199254740992.0),
701 +7169490078506431.0 / 1125899906842624.0,
702 -7877917296929161.0 / 36028797018963968.0,
703 +5470058796541645.0 / 9007199254740992.0,
704 -7233501679212276.0 / 144115188075855872.0 };
705
706 const real q_sin5_b[6] = {+4598161403289824.0 / 144115188075855872.0,
707 +4893639363223553.0 / 2251799813685248.0,
708 -5525704702280927.0 / 2251799813685248.0,
709 +4608444394217766.0 / 4503599627370496.0,
710 -8018266172128683.0 / 9007199254740992.0,
711 +5822429161405618.0 / 9007199254740992.0 };
712// ********************************************************************************
713
714// ********************************************************************************
715// Sinpix_pi: a_k,b_k des Kettenbruchs K5 in A6 = [0.3828125,0.5]:
716// ********************************************************************************
717 const real q_sin6_a[6] = {0.0,
718 +7461973733147728.0 / 36028797018963968.0,
719 +7979710108639590.0 / 140737488355328.0,
720 -6289496525317218.0 / 288230376151711744.0,
721 +6964694082726429.0 / 1125899906842624.0,
722 -5522978779702360.0 / 1152921504606846976.0 };
723
724 const real q_sin6_b[6] = {+5609829449006976.0 / 18014398509481984.0,
725 +8354019229473476.0 / 1125899906842624.0,
726 -8475200820952730.0 / 1125899906842624.0,
727 +5829377160898302.0 / 2251799813685248.0,
728 -5714036688191727.0 / 2251799813685248.0,
729 +7900583259891052.0 / 4503599627370496.0 };
730// ********************************************************************************
731
732// interne Funktion, wird fuer Gammafunktion benoetigt
733int int_no(real *a,
734 const int n, // n: Anzahl der Elemente des Feldes a
735 const real& x) // x: Eine real-Zahl
736// Ein Intervall [A,B] wird durch ein Feld a mit n Elementen in
737// n-1 Teilintervalle unterteilt. Fuer ein x aus [A,B] wird eine
738// Intervall-Nummer zurueckgegeben, wobei das erste Teilintervall
739// die Nr. 0 und das letzte Teilintervall die Nr. n-2 erhaelt.
740// Fuer x<A ist Nr. = -1, und fuer x>=B gilt Nr. = n-1;
741{
742 int i,j,k;
743 i = 0;
744 j = n-1;
745 do
746 {
747 k = (i+j)/2;
748 if (x<a[k]) j = k-1;
749 else i = k+1;
750 }
751 while (i <= j);
752
753 return j;
754}
755// *****************************************************************************
756// Sinpix_pi: Approximation in A1 = [2^-27,0.05078125]:
757// *****************************************************************************
758
759real sinpi_A1(const real& x)
760{
761 real y,v;
762
763 v = 1/(x*x);
764 // Approximation: sin(pi*x)/pi = x*K2(v);
765 y = q_sin1_a[2] / ( v + q_sin1_b[2]);
766 y = q_sin1_a[1] / ( (v + q_sin1_b[1]) + y);
767 y *= x;
768 y += x;
769
770 return y;
771}
772
773// ********************************************************************************
774// Sinpix_pi: Approximation in A2 = [0.05078125,0.1015625]:
775// ********************************************************************************
776real sinpi_A2(const real& x)
777{
778 real y,v;
779
780 if (x == 0.08203125)
781 y = q_sin2_b[0];
782 else
783 {
784 v = 1/(x-0.08203125);
785 // Approximation: sin(pi*x)/pi = x*K4(v);
786 y = q_sin2_a[4] / ( v + q_sin2_b[4]);
787 y = q_sin2_a[3] / ( (v + q_sin2_b[3]) + y);
788 y = q_sin2_a[2] / ( (v + q_sin2_b[2]) + y);
789 y = q_sin2_a[1] / ( (v + q_sin2_b[1]) + y) + q_sin2_b[0];
790 }
791 y *= x;
792 y += x;
793
794 return y;
795}
796
797// ********************************************************************************
798// Sinpix_pi: Approximation in A3 = [0.1015625,0.15234375]:
799// ********************************************************************************
800real sinpi_A3(const real& x)
801{
802 real y,v;
803
804 if (x == 0.13671875)
805 y = q_sin3_b[0];
806 else
807 {
808 v = 1/(x-0.13671875);
809 // Approximation: sin(pi*x)/pi = x*K4(v);
810 y = q_sin3_a[4] / ( v + q_sin3_b[4]);
811 y = q_sin3_a[3] / ( (v + q_sin3_b[3]) + y);
812 y = q_sin3_a[2] / ( (v + q_sin3_b[2]) + y);
813 y = q_sin3_a[1] / ( (v + q_sin3_b[1]) + y) + q_sin3_b[0];
814 }
815 y *= x;
816 y += x;
817
818 return y;
819}
820
821// ********************************************************************************
822// Sinpix_pi: Approximation in A4 = [0.15234375,0.2578125]:
823// ********************************************************************************
824real sinpi_A4(const real& x)
825{
826 real y,v;
827
828 if (x == 0.19140625)
829 y = q_sin4_b[0];
830 else
831 {
832 v = 1/(x-0.19140625);
833 // Approximation: sin(pi*x)/pi = x*K5(v);
834 y = q_sin4_a[5] / ( v + q_sin4_b[5]);
835 y = q_sin4_a[4] / ( (v + q_sin4_b[4]) + y);
836 y = q_sin4_a[3] / ( (v + q_sin4_b[3]) + y);
837 y = q_sin4_a[2] / ( (v + q_sin4_b[2]) + y);
838 y = q_sin4_a[1] / ( (v + q_sin4_b[1]) + y) + q_sin4_b[0];
839 }
840 y *= x;
841 y += x;
842
843 return y;
844}
845
846// ********************************************************************************
847// Sinpix_pi: Approximation in A5 = [0.2578125,0.3828125]:
848// ********************************************************************************
849real sinpi_A5(const real& x)
850{
851 real y,v;
852
853 if (x == 0.30078125)
854 y = q_sin5_b[0];
855 else
856 {
857 v = 1/(x-0.30078125);
858 // Approximation: sin(pi*x)/pi = K5(v);
859 y = q_sin5_a[5] / ( v + q_sin5_b[5]);
860 y = q_sin5_a[4] / ( (v + q_sin5_b[4]) + y);
861 y = q_sin5_a[3] / ( (v + q_sin5_b[3]) + y);
862 y = q_sin5_a[2] / ( (v + q_sin5_b[2]) + y);
863 y = q_sin5_a[1] / ( (v + q_sin5_b[1]) + y) + q_sin5_b[0];
864 }
865 y +=1.0;
866 times2pown(y,-2); // Division durch 4
867
868 return y;
869}
870
871// ********************************************************************************
872// Sinpix_pi: Approximation in A6 = [0.3828125,0.5]:
873// ********************************************************************************
874real sinpi_A6(const real& x)
875{
876 real y,v;
877
878 if (x == 0.43359375)
879 y = q_sin6_b[0];
880 else
881 {
882 v = 1/(x-0.43359375);
883 // Approximation: sin(pi*x)/pi = K5(v);
884 y = q_sin6_a[5] / ( v + q_sin6_b[5]);
885 y = q_sin6_a[4] / ( (v + q_sin6_b[4]) + y);
886 y = q_sin6_a[3] / ( (v + q_sin6_b[3]) + y);
887 y = q_sin6_a[2] / ( (v + q_sin6_b[2]) + y);
888 y = q_sin6_a[1] / ( (v + q_sin6_b[1]) + y) + q_sin6_b[0];
889 }
890
891 return y;
892}
893
895{ // Sin(Pi*x)/Pi; rel. Fehlerschranke: 4.870167e-16;
896 real res(0),xr;
897 int m,nr;
898 bool neg;
899
900 m = Round(x);
901 if (m == -2147483648)
902 cxscthrow(STD_FKT_OUT_OF_DEF("real sinpix_pi(const real&)"));
903 xr = x - m; // xr: reduced argument, exactly calculated!
904 // |xr| <= 0.5;
905 neg = xr<0;
906 if (neg) xr = -xr;
907 // 0 <= xr <= 0.5;
908 nr = int_no(a_sinpix_pi,8,xr);
909
910 switch(nr)
911 {
912 case 0: res = xr; break;
913 case 1: res = sinpi_A1(xr); break;
914 case 2: res = sinpi_A2(xr); break;
915 case 3: res = sinpi_A3(xr); break;
916 case 4: res = sinpi_A4(xr); break;
917 case 5: res = sinpi_A5(xr); break;
918 case 6: res = sinpi_A6(xr); break;
919 case 7: res = 5734161139222659.0 / 18014398509481984.0;
920 }
921
922 if (neg) res = -res;
923 if (m%2 != 0) res = -res;
924
925 return res;
926}
927
928// ********************************************************************************
929// *************** Konstanten bez. Gamma(x) und 1/Gamma(x) ********************
930// ********************************************************************************
931
932// Interval bounds of the neighbour intervals with respect to the
933// intervals S_k for Gamma(x):
934real gam_f85[19] =
935{-0.5, 8.5, 16.5, 24.5, 35.5, 46.5, 57.5, 68.5, 79.5, 90.5,
936 101.5, 112.5, 122.5, 132.5, 142.5, 150.5, 157.5, 164.5, 171.5 };
937
938 // 1/Gamma(x): S0=[1.5,2.5], x0 = 2.0;
939 const real q_gams0_a[8] =
940 {0.0,
941 -7616205496030159.0 / 18014398509481984.0,
942 6808968414757234.0 / 9007199254740992.0,
943 -7179656013820698.0 / 144115188075855872.0,
944 7646522333236288.0 / 144115188075855872.0,
945 -7301751514589296.0 / 144115188075855872.0,
946 6062274112599236.0 / 288230376151711744.0,
947 -7580146497393670.0 / 576460752303423488.0 };
948
949 const real q_gams0_b[8] =
950 {1.0,
951 -4965940208012024.0 / 9007199254740992.0,
952 8627036189024519.0 / 9007199254740992.0,
953 -7819065539096245.0 / 72057594037927936.0,
954 7433219784613049.0 / 36028797018963968.0,
955 7389378817674059.0 / 72057594037927936.0,
956 -6060163955665826.0 / 144115188075855872.0,
957 5769165265609039.0 / 18014398509481984.0 };
958
959 // Gamma(x): S1=[10.5,11.5], x0 = 11.125;
960 const real q_gams1_a[7] =
961 {0.0,
962 5262165366811426.0 / 72057594037927936.0,
963 5574236285326272.0 / 9007199254740992.0,
964 -8823729115230752.0 / 9223372036854775808.0,
965 4639548867796070.0 / 72057594037927936.0,
966 -4784760610310013.0 / 4611686018427387904.0,
967 5900435828568799.0 / 288230376151711744.0 };
968
969 const real q_gams1_b[7] =
970 {7108556377252296.0 / 36028797018963968.0,
971 -7219014979973550.0 / 9007199254740992.0,
972 7224885284258947.0 / 9007199254740992.0,
973 -8569030245421939.0 / 36028797018963968.0,
974 4838653997792226.0 / 18014398509481984.0,
975 -4567015707194151.0 / 36028797018963968.0,
976 5705105513288442.0 / 36028797018963968.0 };
977
978 // Gamma(x): S2=[18.5,19.5], x0 = 18.96875;
979 const real q_gams2_a[6] =
980 {0.0,
981 7322606177885925.0 / 72057594037927936.0,
982 5856515731454725.0 / 144115188075855872.0,
983 -5420981868254561.0 / 1152921504606846976.0,
984 6209135328051357.0 / 1152921504606846976.0,
985 -8261990134869242.0 / 2305843009213693952.0 };
986
987 const real q_gams2_b[6] =
988 {-5267325780498998.0 / 18014398509481984.0,
989 -4688643295042637.0 / 18014398509481984.0,
990 6865435581764388.0 / 36028797018963968.0,
991 -5147410677045000.0 / 144115188075855872.0,
992 5878944809110092.0 / 144115188075855872.0,
993 5315270486582075.0 / 576460752303423488.0 };
994
995 // Gamma(x): S3=[29.5,30.5], x0 = 29.9375;
996 const real q_gams3_a[6] =
997 {0.0,
998 4608894695736103.0 / 9007199254740992.0,
999 4622056591672246.0 / 144115188075855872.0,
1000 7936321632956658.0 / 1152921504606846976.0,
1001 7533691379187725.0 / 2305843009213693952.0,
1002 8317076627520632.0 / 4611686018427387904.0 };
1003
1004 const real q_gams3_b[6] =
1005 {-5793090370845400.0 / 36028797018963968.0,
1006 -5993724737769479.0 / 18014398509481984.0,
1007 -7355055148590989.0 / 288230376151711744.0,
1008 -6811362569175008.0 / 288230376151711744.0,
1009 -5643865380181690.0 / 288230376151711744.0,
1010 -7257737463164532.0 / 288230376151711744.0 };
1011
1012 // Gamma(x): S4=[40.5,41.5], x0 = 41.140625;
1013 const real q_gams4_a[6] =
1014 {0.0,
1015 -7504135817814391.0 / 18014398509481984.0,
1016 4990516833351222.0 / 576460752303423488.0,
1017 7409944094502278.0 / 4611686018427387904.0,
1018 8239437277348520.0 / 2305843009213693952.0,
1019 -5361398870955224.0 / 4611686018427387904.0 };
1020
1021 const real q_gams4_b[6] =
1022 {7405548010914168.0 / 18014398509481984.0,
1023 6819421126036941.0 / 36028797018963968.0,
1024 8474998846370093.0 / 288230376151711744.0,
1025 7733944242652532.0 / 144115188075855872.0,
1026 -6547033063052812.0 / 144115188075855872.0,
1027 -6636555235967309.0 / 576460752303423488.0 };
1028
1029 // Gamma(x): S5=[51.5,52.5], x0 = 52.015625;
1030 const real q_gams5_a[6] =
1031 {0.0,
1032 -7282435522169731.0 / 144115188075855872.0,
1033 7812862759300002.0 / 288230376151711744.0,
1034 -7591789799457117.0 / 9223372036854775808.0,
1035 7290376444377792.0 / 2305843009213693952.0,
1036 -7051827298983325.0 / 9223372036854775808.0 };
1037
1038 const real q_gams5_b[6] =
1039 {-4692552597358020.0 / 36028797018963968.0,
1040 7065248316566075.0 / 36028797018963968.0,
1041 -5667667511562837.0 / 36028797018963968.0,
1042 7741409329322298.0 / 144115188075855872.0,
1043 -6371025742549657.0 / 144115188075855872.0,
1044 8045016330906311.0 / 288230376151711744.0 };
1045
1046 // Gamma(x): S6=[62.5,63.5], x0 = 63.015625;
1047 const real q_gams6_a[6] =
1048 {0.0,
1049 6719861857699617.0 / 36028797018963968.0,
1050 6146108740657636.0 / 1152921504606846976.0,
1051 -8996491371873855.0 / 4611686018427387904.0,
1052 5082150623010284.0 / 2305843009213693952.0,
1053 -4650554434211955.0 / 18446744073709551616.0 };
1054
1055 const real q_gams6_b[6] =
1056 {6795471792542416.0 / 18014398509481984.0,
1057 -4567367130077106.0 / 36028797018963968.0,
1058 8403023343856889.0 / 288230376151711744.0,
1059 5083390457138636.0 / 144115188075855872.0,
1060 -4614092542204999.0 / 144115188075855872.0,
1061 7104477795873427.0 / 72057594037927936.0 };
1062
1063 // Gamma(x): S7=[73.5,74.5], x0 = 74.16015625;
1064 const real q_gams7_a[6] =
1065 {0.0 / 1.0,
1066 5372276754422009.0 / 18014398509481984.0,
1067 4663470139182266.0 / 576460752303423488.0,
1068 8173857739398457.0 / 4611686018427387904.0,
1069 4932343261664244.0 / 4611686018427387904.0,
1070 -5822870305854791.0 / 295147905179352825856.0 };
1071
1072 const real q_gams7_b[6] =
1073 {-4806355488125952.0 / 1152921504606846976.0,
1074 -6211400907258787.0 / 36028797018963968.0,
1075 -5429997192499743.0 / 288230376151711744.0,
1076 -5653057298064211.0 / 288230376151711744.0,
1077 -4746908442350821.0 / 1152921504606846976.0,
1078 8670733620601602.0 / 18014398509481984.0 };
1079
1080 // Gamma(x): S8=[84.5,85.5], x0 = 85.1015625;
1081 const real q_gams8_a[5] =
1082 { 0.0 / 1.0,
1083 -8754830070807807.0 / 72057594037927936.0,
1084 7931975738629914.0 / 2305843009213693952.0,
1085 6414800170661986.0 / 73786976294838206464.0,
1086 4558538137025832.0 / 18014398509481984.0};
1087
1088 const real q_gams8_b[5] =
1089 {-4924952929015874.0 / 18014398509481984.0,
1090 8571263173618757.0 / 72057594037927936.0,
1091 7912939311540494.0 / 576460752303423488.0,
1092 8975599623306617.0 / 18014398509481984.0,
1093 -4569152133381960.0 / 9007199254740992.0 };
1094
1095 // Gamma(x): S9=[95.5,96.5], x0 = 95.984375;
1096 const real q_gams9_a[5] =
1097 {0.0 / 1.0,
1098 -6140046099729716.0 / 144115188075855872.0,
1099 7278997066304035.0 / 576460752303423488.0,
1100 -4759432412103962.0 / 9223372036854775808.0,
1101 6912711986126027.0 / 4611686018427387904.0 };
1102
1103 const real q_gams9_b[5] =
1104 {-5611185402650416.0 / 72057594037927936.0,
1105 4915638659703209.0 / 36028797018963968.0,
1106 -7692469721238990.0 / 72057594037927936.0,
1107 5045754589592299.0 / 144115188075855872.0,
1108 -8290188163752846.0 / 288230376151711744.0 };
1109
1110 // Gamma(x): S10=[106.5,107.5], x0 = 107.078125;
1111 const real q_gams10_a[5] =
1112 {0.0 / 1.0,
1113 4717576521494070.0 / 72057594037927936.0,
1114 6906626643101502.0 / 1152921504606846976.0,
1115 -8147816895338065.0 / 9223372036854775808.0,
1116 7960901442281121.0 / 9223372036854775808.0 };
1117
1118 const real q_gams10_b[5] =
1119 {7952058244493952.0 / 288230376151711744.0,
1120 -7601355771801470.0 / 72057594037927936.0,
1121 4923655207606594.0 / 72057594037927936.0,
1122 -7238499638430702.0 / 576460752303423488.0,
1123 5852764550899526.0 / 576460752303423488.0 };
1124
1125 // Gamma(x): S11=[117.5,118.5], x0 = 117.8671875;
1126 const real q_gams11_a[5] =
1127 {0.0 / 1.0,
1128 5000157748740957.0 / 36028797018963968.0,
1129 6733781271808263.0 / 2305843009213693952.0,
1130 4912448110743899.0 / 18446744073709551616.0,
1131 6917678469221762.0 / 576460752303423488.0 };
1132
1133 const real q_gams11_b[5] =
1134 {-4805173503400964.0 / 36028797018963968.0,
1135 -7686561799456867.0 / 72057594037927936.0,
1136 -6649023763933472.0 / 576460752303423488.0,
1137 -7324214555543306.0 / 72057594037927936.0,
1138 8284786467509721.0 / 72057594037927936.0 };
1139
1140 // Gamma(x): Stuetzinterval S12=[126.5,127.5], x0 = 126.7421875;
1141 const real q_gams12_a[5] =
1142 {0.0 / 1.0,
1143 5226845396890227.0 / 36028797018963968.0,
1144 5602232597459763.0 / 1152921504606846976.0,
1145 4920783552645014.0 / 4611686018427387904.0,
1146 5570041179715558.0 / 9223372036854775808.0 };
1147
1148 const real q_gams12_b[5] =
1149 {-6799658092196962.0 / 18014398509481984.0,
1150 -4810318281953418.0 / 36028797018963968.0,
1151 -8340064868499919.0 / 576460752303423488.0,
1152 -8375346434191481.0 / 576460752303423488.0,
1153 -6231299816839781.0 / 1152921504606846976.0 };
1154
1155 // Gamma(x): S13=[137.5,138.5], x0 = 138.0390625;
1156 const real q_gams13_a[5] =
1157 {0.0 / 1.0,
1158 5077602289066213.0 / 18014398509481984.0,
1159 4971389438544026.0 / 576460752303423488.0,
1160 8326326060413143.0 / 4611686018427387904.0,
1161 7582214946963028.0 / 9223372036854775808.0 };
1162
1163 const real q_gams13_b[5] =
1164 {-8336661118182000.0 / 72057594037927936.0,
1165 -6152827173580884.0 / 36028797018963968.0,
1166 -6306415240260295.0 / 576460752303423488.0,
1167 -5949724688158565.0 / 576460752303423488.0,
1168 -5625482009313801.0 / 576460752303423488.0 };
1169
1170 // Gamma(x): S14=[146.5,147.5], x0 = 146.94921875;
1171 const real q_gams14_a[6] =
1172 {0.0 / 1.0,
1173 8624518348171320.0 / 18014398509481984.0,
1174 7049908374954842.0 / 576460752303423488.0,
1175 5769548496333642.0 / 2305843009213693952.0,
1176 5110826386599631.0 / 4611686018427387904.0,
1177 5902099580937297.0 / 9223372036854775808.0 };
1178
1179 const real q_gams14_b[6] =
1180 {4591842870321392.0 / 18014398509481984.0,
1181 -7195103997153274.0 / 36028797018963968.0,
1182 -5062133254875404.0 / 576460752303423488.0,
1183 -4899776496053167.0 / 576460752303423488.0,
1184 -4703735605739871.0 / 576460752303423488.0,
1185 -8717622510435264.0 / 1152921504606846976.0 };
1186
1187 // Gamma(x): Stuetzinterval S15=[153.5,154.5], x0 = 153.90234375;
1188 const real q_gams15_a[6] =
1189 {0.0 / 1.0,
1190 5047093325633351.0 / 9007199254740992.0,
1191 8838568428828895.0 / 576460752303423488.0,
1192 7169875770591066.0 / 2305843009213693952.0,
1193 6277116142443993.0 / 4611686018427387904.0,
1194 7160380168510556.0 / 9223372036854775808.0 };
1195
1196 const real q_gams15_b[6] =
1197 {5575895624898616.0 / 18014398509481984.0,
1198 -7982725134478240.0 / 36028797018963968.0,
1199 -8683131154525203.0 / 1152921504606846976.0,
1200 -8504376590080638.0 / 1152921504606846976.0,
1201 -8272396567408614.0 / 1152921504606846976.0,
1202 -7370684458128734.0 / 1152921504606846976.0 };
1203
1204 // Gamma(x): S16=[160.5,161.5], x0 = 161.08984375;
1205 const real q_gams16_a[6] =
1206 {0.0 / 1.0,
1207 8928052213955455.0 / 18014398509481984.0,
1208 5405751921217612.0 / 288230376151711744.0,
1209 8725861134087760.0 / 2305843009213693952.0,
1210 7583718370137025.0 / 4611686018427387904.0,
1211 8577302812898546.0 / 9223372036854775808.0 };
1212
1213 const real q_gams16_b[6] =
1214 {6669445708872192.0 / 144115188075855872.0,
1215 -8769965967955460.0 / 36028797018963968.0,
1216 -7524042392635917.0 / 1152921504606846976.0,
1217 -7422547685779587.0 / 1152921504606846976.0,
1218 -7284571995868819.0 / 1152921504606846976.0,
1219 -7797018373646746.0 / 1152921504606846976.0 };
1220
1221 // Gamma(x): S17=[167.5,168.5], x0 = 168.0;
1222 const real q_gams17_a[6] =
1223 {0.0 / 1.0,
1224 4642907317603488.0 / 9007199254740992.0,
1225 6403634366308572.0 / 288230376151711744.0,
1226 5153510478021874.0 / 1152921504606846976.0,
1227 8919203526515161.0 / 4611686018427387904.0,
1228 5015558226948384.0 / 4611686018427387904.0 };
1229
1230 const real q_gams17_b[6] =
1231 {-6229620043098112.0 / 9223372036854775808.0,
1232 -4750296278401558.0 / 18014398509481984.0,
1233 -6639785613322219.0 / 1152921504606846976.0,
1234 -6577910457226093.0 / 1152921504606846976.0,
1235 -6491061443859352.0 / 1152921504606846976.0,
1236 -6372399542597081.0 / 1152921504606846976.0 };
1237
1238// ****************************************************************************
1239// ******************** Gamma(x), 1/Gamma(x) **********************************
1240// ****************************************************************************
1241
1242 real gam_S0(const real& x)
1243// Calculating approximations for 1/Gamma(x) in S0 = [1.5,2.5];
1244// Rel. error bound: 3.930138e-16;
1245// Blomquist, 19.05.2009;
1246{
1247 real y(0),v;
1248
1249 // Continued fraction: K_7(v), v = 1/(x-x0), x0=2;
1250 if (x==2)
1251 y = q_gams0_b[0];
1252 else
1253 {
1254 v = 1/(x-2);
1255
1256 y = q_gams0_a[7] / ( v + q_gams0_b[7]);
1257 y = q_gams0_a[6] / ( (v + q_gams0_b[6]) + y );
1258 y = q_gams0_a[5] / ( (v + q_gams0_b[5]) + y );
1259 y = q_gams0_a[4] / ( (v + q_gams0_b[4]) + y );
1260 y = q_gams0_a[3] / ( (v + q_gams0_b[3]) + y );
1261 y = q_gams0_a[2] / ( (v + q_gams0_b[2]) + y );
1262 y = q_gams0_a[1] / ( (v + q_gams0_b[1]) + y ) + q_gams0_b[0];
1263 }
1264 return y;
1265}
1266
1267 real gam_S0_n0(const real& x)
1268// Calculating approximations for K_7(v), v=1/x, x in [-0.5,+0.5];
1269// Rel. error bound: 3.930138e-16;
1270// Blomquist, 19.05.2009;
1271{
1272 real y(0),v;
1273
1274 // Continued fraction: K_7(v), v = 1/x, x0=0;
1275 if (x==0)
1276 y = q_gams0_b[0];
1277 else
1278 {
1279 v = 1/x;
1280
1281 y = q_gams0_a[7] / ( v + q_gams0_b[7]);
1282 y = q_gams0_a[6] / ( (v + q_gams0_b[6]) + y );
1283 y = q_gams0_a[5] / ( (v + q_gams0_b[5]) + y );
1284 y = q_gams0_a[4] / ( (v + q_gams0_b[4]) + y );
1285 y = q_gams0_a[3] / ( (v + q_gams0_b[3]) + y );
1286 y = q_gams0_a[2] / ( (v + q_gams0_b[2]) + y );
1287 y = q_gams0_a[1] / ( (v + q_gams0_b[1]) + y ) + q_gams0_b[0];
1288 }
1289
1290 return y;
1291}
1292
1293real gam_S0_n1(const real& x)
1294// Calculating approximations for K_7(v), v=1/(x-1), x in [0.5,1.5];
1295// Rel. error bound: 3.930138e-16;
1296// Blomquist, 21.05.2009;
1297{
1298 real y(0),v;
1299
1300 // Continued fraction: K_7(v), v = 1/(x-1), x0=1;
1301 if (x==1)
1302 y = q_gams0_b[0];
1303 else
1304 {
1305 v = 1/(x-1);
1306
1307 y = q_gams0_a[7] / ( v + q_gams0_b[7]);
1308 y = q_gams0_a[6] / ( (v + q_gams0_b[6]) + y );
1309 y = q_gams0_a[5] / ( (v + q_gams0_b[5]) + y );
1310 y = q_gams0_a[4] / ( (v + q_gams0_b[4]) + y );
1311 y = q_gams0_a[3] / ( (v + q_gams0_b[3]) + y );
1312 y = q_gams0_a[2] / ( (v + q_gams0_b[2]) + y );
1313 y = q_gams0_a[1] / ( (v + q_gams0_b[1]) + y ) + q_gams0_b[0];
1314 }
1315
1316 return y;
1317}
1318
1319real gammar_S0(const real& x)
1320// Calculating approximations for 1/Gamma(x), x in (-0.5,+8.5];
1321// eps = 1.725284e-15;
1322// Blomquist, 21.05.2009;
1323{
1324 real y,Ne;
1325 int n,p;
1326
1327 n = Round(x);
1328
1329 switch(n)
1330 {
1331 case 0: if (expo(x)<=-52) y = x;
1332 else y = x*(x+1)*gam_S0_n0(x); break; // x in (-0.5,+0.5);
1333 case 1: y = x*gam_S0_n1(x); break; // x in [0.5,1.5);
1334 case 2: y = gam_S0(x); break; // x in [1.5,2.5);
1335
1336 default: // n=3,4,...,8; x in [2.5,8.5];
1337 p = n-2; Ne = x-1;
1338 for (int k=2; k<=p; k++) Ne *= (x-k);
1339 y = gam_S0(x-p) / Ne;
1340 }
1341
1342 return y;
1343}
1344
1345real gamma_S0(const real& x)
1346// Calculating approximations for Gamma(x) in S0 = (-0.5,+8.5];
1347// eps = 1.947330E-15;
1348// Blomquist, 21.05.2009;
1349{
1350 real y = 1/gammar_S0(x);
1351 return y;
1352}
1353
1354real gam_S1(const real& x)
1355// Calculating approximations for Gamma(x) in S1 = [10.5,11.5];
1356// Rel. error bound: 7.696082e-16;
1357// Blomquist, 22.05.2009;
1358{
1359 real y(0),v;
1360
1361 // Continued fraction: K_6(v), v = 1/(x-x0), x0=11.125;
1362 if (x==11.125)
1363 y = q_gams1_b[0];
1364 else
1365 {
1366 v = 1/(x-11.125);
1367
1368 y = q_gams1_a[6] / ( v + q_gams1_b[6]);
1369 y = q_gams1_a[5] / ( (v + q_gams1_b[5]) + y );
1370 y = q_gams1_a[4] / ( (v + q_gams1_b[4]) + y );
1371 y = q_gams1_a[3] / ( (v + q_gams1_b[3]) + y );
1372 y = q_gams1_a[2] / ( (v + q_gams1_b[2]) + y );
1373 y = q_gams1_a[1] / ( (v + q_gams1_b[1]) + y ) + q_gams1_b[0];
1374 }
1375 y += 1;
1376 y *= q_ex10(x);
1377 times2pown(y,-15);
1378
1379 return y;
1380}
1381
1382real gamma_S1(const real& x)
1383// Calculating approximations for Gamma(x), x in [8.5,16.5];
1384// eps = 1.879833e-15;
1385// Blomquist, 24.05.2009;
1386{
1387 real y,Pr;
1388 int n,p;
1389
1390 n = Round(x);
1391 p = n-11;
1392
1393 if (n>11) // neighbour intervals to the right
1394 {
1395 Pr = x-1;
1396 for (int k=2; k<=p; k++)
1397 Pr *= x-k;
1398 y = Pr*gam_S1(x-p);
1399 }
1400 else // neighbour intervals to the left and S1 itself
1401 {
1402 p = -p; Pr = x;
1403 for (int k=1; k<=p-1; k++)
1404 Pr *= x+k;
1405 y = (p==0)? gam_S1(x) : gam_S1(x+p)/Pr;
1406 }
1407
1408 return y;
1409}
1410
1411real gam_S2(const real& x)
1412// Calculating approximations for Gamma(x) in S2 = [18.5,19.5];
1413// Rel. error bound: 7.411454e-16;
1414// Blomquist, 25.05.2009;
1415{
1416 real y(0),v;
1417
1418 // Continued fraction: K_5(v), v = 1/(x-x0), x0=18.96875;
1419 if (x==18.96875)
1420 y = q_gams2_b[0];
1421 else
1422 {
1423 v = 1/(x-18.96875);
1424
1425 y = q_gams2_a[5] / ( v + q_gams2_b[5]);
1426 y = q_gams2_a[4] / ( (v + q_gams2_b[4]) + y );
1427 y = q_gams2_a[3] / ( (v + q_gams2_b[3]) + y );
1428 y = q_gams2_a[2] / ( (v + q_gams2_b[2]) + y );
1429 y = q_gams2_a[1] / ( (v + q_gams2_b[1]) + y ) + q_gams2_b[0];
1430 }
1431 y += 1;
1432 y *= q_exp2(4*x);
1433 times2pown(y,-23);
1434
1435 return y;
1436}
1437
1438real gamma_S2(const real& x)
1439// Calculating approximations for Gamma(x), x in [16.5,24.5];
1440// eps = 1.851370e-15;
1441// Blomquist, 25.05.2009;
1442{
1443 real y,Pr;
1444 int n,p;
1445
1446 n = Round(x);
1447 p = n-19;
1448
1449 if (n>19) // neighbour intervals to the right
1450 {
1451 Pr = x-1;
1452 for (int k=2; k<=p; k++)
1453 Pr *= x-k;
1454 y = Pr*gam_S2(x-p);
1455 }
1456 else // neighbour intervals to the left and S2 itself
1457 {
1458 p = -p; Pr = x;
1459 for (int k=1; k<=p-1; k++)
1460 Pr *= x+k;
1461 y = (p==0)? gam_S2(x) : gam_S2(x+p)/Pr;
1462 }
1463
1464 return y;
1465}
1466
1467real gam_S3(const real& x)
1468// Calculating approximations for Gamma(x) in S3 = [29.5,30.5];
1469// Rel. error bound: 9.724029e-16;
1470// Blomquist, 26.05.2009;
1471{
1472 real y(0),v;
1473
1474 // Continued fraction: K_5(v), v = 1/(x-x0), x0=29.9375;
1475 if (x==29.9375)
1476 y = q_gams3_b[0];
1477 else
1478 {
1479 v = 1/(x-29.9375);
1480
1481 y = q_gams3_a[5] / ( v + q_gams3_b[5]);
1482 y = q_gams3_a[4] / ( (v + q_gams3_b[4]) + y );
1483 y = q_gams3_a[3] / ( (v + q_gams3_b[3]) + y );
1484 y = q_gams3_a[2] / ( (v + q_gams3_b[2]) + y );
1485 y = q_gams3_a[1] / ( (v + q_gams3_b[1]) + y ) + q_gams3_b[0];
1486 }
1487 y += 1;
1488 y *= q_exp2(4*x);
1489 times2pown(y,-17);
1490
1491 return y;
1492}
1493
1494real gamma_S3(const real& x)
1495// Calculating approximations for Gamma(x), x in [24.5,35.5];
1496// eps = 2.082627e-15;
1497// Blomquist, 26.05.2009;
1498{
1499 real y,Pr;
1500 int n,p;
1501
1502 n = Round(x);
1503 p = n-30;
1504
1505 if (n>30) // neighbour intervals to the right
1506 {
1507 Pr = x-1;
1508 for (int k=2; k<=p; k++)
1509 Pr *= x-k;
1510 y = Pr*gam_S3(x-p);
1511 }
1512 else // neighbour intervals to the left and S3 itself
1513 {
1514 p = -p; Pr = x;
1515 for (int k=1; k<=p-1; k++)
1516 Pr *= x+k;
1517 y = (p==0)? gam_S3(x) : gam_S3(x+p)/Pr;
1518 }
1519
1520 return y;
1521}
1522
1523real gam_S4(const real& x)
1524// Calculating approximations for Gamma(x) in S4 = [40.5,41.5];
1525// Rel. error bound: 8.170628e-16;
1526// Blomquist, 26.05.2009;
1527{
1528 real y(0),v;
1529
1530 // Continued fraction: K_5(v), v = 1/(x-x0), x0=41.140625;
1531 if (x==41.140625)
1532 y = q_gams4_b[0];
1533 else
1534 {
1535 v = 1/(x-41.140625);
1536
1537 y = q_gams4_a[5] / ( v + q_gams4_b[5]);
1538 y = q_gams4_a[4] / ( (v + q_gams4_b[4]) + y );
1539 y = q_gams4_a[3] / ( (v + q_gams4_b[3]) + y );
1540 y = q_gams4_a[2] / ( (v + q_gams4_b[2]) + y );
1541 y = q_gams4_a[1] / ( (v + q_gams4_b[1]) + y ) + q_gams4_b[0];
1542 }
1543 y += 1;
1544 y *= exp(4*x);
1545 times2pown(y,-78);
1546
1547 return y;
1548}
1549
1550real gamma_S4(const real& x)
1551// Calculating approximations for Gamma(x), x in [35.5,46.5];
1552// eps = 1.927286e-15;
1553// Blomquist, 28.05.2009;
1554{
1555 real y,Pr;
1556 int n,p;
1557
1558 n = Round(x);
1559 p = n-41;
1560
1561 if (n>41) // neighbour intervals to the right
1562 {
1563 Pr = x-1;
1564 for (int k=2; k<=p; k++)
1565 Pr *= x-k;
1566 y = Pr*gam_S4(x-p);
1567 }
1568 else // neighbour intervals to the left and S4 itself
1569 {
1570 p = -p; Pr = x;
1571 for (int k=1; k<=p-1; k++)
1572 Pr *= x+k;
1573 y = (p==0)? gam_S4(x) : gam_S4(x+p)/Pr;
1574 }
1575
1576 return y;
1577}
1578
1579real gam_S5(const real& x)
1580// Calculating approximations for Gamma(x) in S5 = [51.5,52.5];
1581// Rel. error bound: 6.351369e-16;
1582// Blomquist, 28.05.2009;
1583{
1584 real y(0),v;
1585
1586 // Continued fraction: K_5(v), v = 1/(x-x0), x0=52.015625;
1587 if (x==52.015625)
1588 y = q_gams5_b[0];
1589 else
1590 {
1591 v = 1/(x-52.015625);
1592
1593 y = q_gams5_a[5] / ( v + q_gams5_b[5]);
1594 y = q_gams5_a[4] / ( (v + q_gams5_b[4]) + y );
1595 y = q_gams5_a[3] / ( (v + q_gams5_b[3]) + y );
1596 y = q_gams5_a[2] / ( (v + q_gams5_b[2]) + y );
1597 y = q_gams5_a[1] / ( (v + q_gams5_b[1]) + y ) + q_gams5_b[0];
1598 }
1599 y += 1;
1600 y *= exp(4*x);
1601 times2pown(y,-80);
1602
1603 return y;
1604}
1605
1606real gamma_S5(const real& x)
1607// Calculating approximations for Gamma(x), x in [46.5,57.5];
1608// eps = 1.745360e-15;
1609// Blomquist, 29.05.2009;
1610{
1611 real y,Pr;
1612 int n,p;
1613
1614 n = Round(x);
1615 p = n-52;
1616
1617 if (n>52) // neighbour intervals to the right
1618 {
1619 Pr = x-1;
1620 for (int k=2; k<=p; k++)
1621 Pr *= x-k;
1622 y = Pr*gam_S5(x-p);
1623 }
1624 else // neighbour intervals to the left and S5 itself
1625 {
1626 p = -p; Pr = x;
1627 for (int k=1; k<=p-1; k++)
1628 Pr *= x+k;
1629 y = (p==0)? gam_S5(x) : gam_S5(x+p)/Pr;
1630 }
1631
1632 return y;
1633}
1634
1635real gam_S6(const real& x)
1636// Calculating approximations for Gamma(x) in S6 = [62.5,63.5];
1637// Rel. error bound: 8.310104e-16;
1638// Blomquist, 29.05.2009;
1639{
1640 real y(0),v;
1641
1642 // Continued fraction: K_5(v), v = 1/(x-x0), x0=63.015625;
1643 if (x==63.015625)
1644 y = q_gams6_b[0];
1645 else
1646 {
1647 v = 1/(x-63.015625);
1648
1649 y = q_gams6_a[5] / ( v + q_gams6_b[5]);
1650 y = q_gams6_a[4] / ( (v + q_gams6_b[4]) + y );
1651 y = q_gams6_a[3] / ( (v + q_gams6_b[3]) + y );
1652 y = q_gams6_a[2] / ( (v + q_gams6_b[2]) + y );
1653 y = q_gams6_a[1] / ( (v + q_gams6_b[1]) + y ) + q_gams6_b[0];
1654 }
1655 y += 1;
1656 y *= exp(4*x);
1657 times2pown(y,-80);
1658
1659 return y;
1660}
1661
1662real gamma_S6(const real& x)
1663// Calculating approximations for Gamma(x), x in [57.5,68.5];
1664// eps = 1.941234e-15;
1665// Blomquist, 06.06.2009;
1666{
1667 real y,Pr;
1668 int n,p;
1669
1670 n = Round(x);
1671 p = n-63;
1672
1673 if (n>63) // neighbour intervals to the right
1674 {
1675 Pr = x-1;
1676 for (int k=2; k<=p; k++)
1677 Pr *= x-k;
1678 y = Pr*gam_S6(x-p);
1679 }
1680 else // neighbour intervals to the left and S6 itself
1681 {
1682 p = -p; Pr = x;
1683 for (int k=1; k<=p-1; k++)
1684 Pr *= x+k;
1685 y = (p==0)? gam_S6(x) : gam_S6(x+p)/Pr;
1686 }
1687
1688 return y;
1689}
1690
1691
1692real gam_S7(const real& x)
1693// Calculating approximations for Gamma(x) in S7 = [73.5,74.5];
1694// Rel. error bound: 7.999142e-16;
1695// Blomquist, 05.06.2009;
1696{
1697 real y(0),v;
1698
1699 // Continued fraction: K_5(v), v = 1/(x-x0), x0=74.16015625;
1700 if (x==74.16015625)
1701 y = q_gams7_b[0];
1702 else
1703 {
1704 v = 1/(x-74.16015625);
1705
1706 y = q_gams7_a[5] / ( v + q_gams7_b[5]);
1707 y = q_gams7_a[4] / ( (v + q_gams7_b[4]) + y );
1708 y = q_gams7_a[3] / ( (v + q_gams7_b[3]) + y );
1709 y = q_gams7_a[2] / ( (v + q_gams7_b[2]) + y );
1710 y = q_gams7_a[1] / ( (v + q_gams7_b[1]) + y ) + q_gams7_b[0];
1711 }
1712 y += 1;
1713 y *= exp(4*x);
1714 times2pown(y,-76);
1715
1716 return y;
1717}
1718
1719real gamma_S7(const real& x)
1720// Calculating approximations for Gamma(x), x in [68.5,79.5];
1721// eps = 1.910138e-15;
1722// Blomquist, 08.06.2009;
1723{
1724 real y,Pr;
1725 int n,p;
1726
1727 n = Round(x);
1728 p = n-74;
1729
1730 if (n>74) // neighbour intervals to the right
1731 {
1732 Pr = x-1;
1733 for (int k=2; k<=p; k++)
1734 Pr *= x-k;
1735 y = Pr*gam_S7(x-p);
1736 }
1737 else // neighbour intervals to the left and S7 itself
1738 {
1739 p = -p; Pr = x;
1740 for (int k=1; k<=p-1; k++)
1741 Pr *= x+k;
1742 y = (p==0)? gam_S7(x) : gam_S7(x+p)/Pr;
1743 }
1744
1745 return y;
1746}
1747
1748real gam_S8(const real& x)
1749// Calculating approximations for Gamma(x) in S8 = [84.5,85.5];
1750// Rel. error bound: 7.192334e-16;
1751// Blomquist, 08.06.2009;
1752{
1753 real y(0),v;
1754
1755 // Continued fraction: K_4(v), v = 1/(x-x0), x0=85.1015625;
1756 if (x==85.1015625)
1757 y = q_gams8_b[0];
1758 else
1759 {
1760 v = 1/(x-85.1015625);
1761
1762 y = q_gams8_a[4] / ( v + q_gams8_b[4]);
1763 y = q_gams8_a[3] / ( (v + q_gams8_b[3]) + y );
1764 y = q_gams8_a[2] / ( (v + q_gams8_b[2]) + y );
1765 y = q_gams8_a[1] / ( (v + q_gams8_b[1]) + y ) + q_gams8_b[0];
1766 }
1767 y += 1;
1768 y *= q_ex10(2*x);
1769 times2pown(y,-144);
1770
1771 return y;
1772}
1773
1774real gamma_S8(const real& x)
1775// Calculating approximations for Gamma(x), x in [79.5,90.5];
1776// eps = 1.829457e-15;
1777// Blomquist, 09.06.2009;
1778{
1779 real y,Pr;
1780 int n,p;
1781
1782 n = Round(x);
1783 p = n-85;
1784
1785 if (n>85) // neighbour intervals to the right
1786 {
1787 Pr = x-1;
1788 for (int k=2; k<=p; k++)
1789 Pr *= x-k;
1790 y = Pr*gam_S8(x-p);
1791 }
1792 else // neighbour intervals to the left and S8 itself
1793 {
1794 p = -p; Pr = x;
1795 for (int k=1; k<=p-1; k++)
1796 Pr *= x+k;
1797 y = (p==0)? gam_S8(x) : gam_S8(x+p)/Pr;
1798 }
1799
1800 return y;
1801}
1802
1803
1804real gam_S9(const real& x)
1805// Calculating approximations for Gamma(x) in S9 = [95.5,96.5];
1806// Rel. error bound: 6.622145e-16;
1807// Blomquist, 09.06.2009;
1808{
1809 real y(0),v;
1810
1811 // Continued fraction: K_4(v), v = 1/(x-x0), x0=95.984375;
1812 if (x==95.984375)
1813 y = q_gams9_b[0];
1814 else
1815 {
1816 v = 1/(x-95.984375);
1817
1818 y = q_gams9_a[4] / ( v + q_gams9_b[4]);
1819 y = q_gams9_a[3] / ( (v + q_gams9_b[3]) + y );
1820 y = q_gams9_a[2] / ( (v + q_gams9_b[2]) + y );
1821 y = q_gams9_a[1] / ( (v + q_gams9_b[1]) + y ) + q_gams9_b[0];
1822 }
1823 y += 1;
1824 y *= q_ex10(2*x);
1825 times2pown(y,-146);
1826
1827 return y;
1828}
1829
1830real gamma_S9(const real& x)
1831// Calculating approximations for Gamma(x), x in [90.5,101.5];
1832// eps = 1.772438e-15;
1833// Blomquist, 10.06.2009;
1834{
1835 real y,Pr;
1836 int n,p;
1837
1838 n = Round(x);
1839 p = n-96;
1840
1841 if (n>96) // neighbour intervals to the right
1842 {
1843 Pr = x-1;
1844 for (int k=2; k<=p; k++)
1845 Pr *= x-k;
1846 y = Pr*gam_S9(x-p);
1847 }
1848 else // neighbour intervals to the left and S9 itself
1849 {
1850 p = -p; Pr = x;
1851 for (int k=1; k<=p-1; k++)
1852 Pr *= x+k;
1853 y = (p==0)? gam_S9(x) : gam_S9(x+p)/Pr;
1854 }
1855
1856 return y;
1857}
1858
1859
1860real gam_S10(const real& x)
1861// Calculating approximations for Gamma(x) in S10 = [106.5,107.5];
1862// Rel. error bound: 8.254965e-16;
1863// Blomquist, 10.06.2009;
1864{
1865 real y(0),v;
1866
1867 // Continued fraction: K_4(v), v = 1/(x-x0), x0=107.078125;
1868 if (x==107.078125)
1869 y = q_gams10_b[0];
1870 else
1871 {
1872 v = 1/(x-107.078125);
1873
1874 y = q_gams10_a[4] / ( v + q_gams10_b[4]);
1875 y = q_gams10_a[3] / ( (v + q_gams10_b[3]) + y );
1876 y = q_gams10_a[2] / ( (v + q_gams10_b[2]) + y );
1877 y = q_gams10_a[1] / ( (v + q_gams10_b[1]) + y ) + q_gams10_b[0];
1878 }
1879 y += 1;
1880 y *= q_ex10(2*x);
1881 times2pown(y,-146);
1882
1883 return y;
1884}
1885
1886real gamma_S10(const real& x)
1887// Calculating approximations for Gamma(x), x in [101.5,112.5];
1888// eps = 1.935720e-15;
1889// Blomquist, 12.06.2009;
1890{
1891 real y,Pr;
1892 int n,p;
1893
1894 n = Round(x);
1895 p = n-107;
1896
1897 if (n>107) // neighbour intervals to the right
1898 {
1899 Pr = x-1;
1900 for (int k=2; k<=p; k++)
1901 Pr *= x-k;
1902 y = Pr*gam_S10(x-p);
1903 }
1904 else // neighbour intervals to the left and S10 itself
1905 {
1906 p = -p; Pr = x;
1907 for (int k=1; k<=p-1; k++)
1908 Pr *= x+k;
1909 y = (p==0)? gam_S10(x) : gam_S10(x+p)/Pr;
1910 }
1911
1912 return y;
1913}
1914
1915real gam_S11(const real& x)
1916// Calculating approximations for Gamma(x) in S11 = [117.5,118.5];
1917// Rel. error bound: 6.766734e-16;
1918// Blomquist, 11.06.2009;
1919{
1920 real y(0),v;
1921
1922 // Continued fraction: K_4(v), v = 1/(x-x0), x0=117.8671875;
1923 if (x==117.8671875)
1924 y = q_gams11_b[0];
1925 else
1926 {
1927 v = 1/(x-117.8671875);
1928
1929 y = q_gams11_a[4] / ( v + q_gams11_b[4]);
1930 y = q_gams11_a[3] / ( (v + q_gams11_b[3]) + y );
1931 y = q_gams11_a[2] / ( (v + q_gams11_b[2]) + y );
1932 y = q_gams11_a[1] / ( (v + q_gams11_b[1]) + y ) + q_gams11_b[0];
1933 }
1934 y += 1;
1935 y *= q_ex10(2*x);
1936 times2pown(y,-144);
1937
1938 return y;
1939}
1940
1941real gamma_S11(const real& x)
1942// Calculating approximations for Gamma(x) in [112.5,122.5];
1943// eps = 1.786897e-15;
1944// Blomquist, 14.06.2009;
1945{
1946 real y,Pr;
1947 int n,p;
1948
1949 n = Round(x);
1950 p = n-118;
1951
1952 if (n>118) // neighbour intervals to the right
1953 {
1954 Pr = x-1;
1955 for (int k=2; k<=p; k++)
1956 Pr *= x-k;
1957 y = Pr*gam_S11(x-p);
1958 }
1959 else // left neighbour intervals and S11 itself
1960 {
1961 p = -p; Pr = x;
1962 for (int k=1; k<=p-1; k++)
1963 Pr *= x+k;
1964 y = (p==0)? gam_S11(x) : gam_S11(x+p)/Pr;
1965 }
1966
1967 return y;
1968}
1969
1970real gam_S12(const real& x)
1971// Calculating approximations for Gamma(x) in S12 = [126.5,127.5];
1972// Rel. error bound: 8.175777e-16;
1973// Blomquist, 13.06.2009;
1974{
1975 real y(0),v;
1976
1977 // Continued fraction: K_4(v), v = 1/(x-x0), x0=126.7421875;
1978 if (x==126.7421875)
1979 y = q_gams12_b[0];
1980 else
1981 {
1982 v = 1/(x-126.7421875);
1983
1984 y = q_gams12_a[4] / ( v + q_gams12_b[4]);
1985 y = q_gams12_a[3] / ( (v + q_gams12_b[3]) + y );
1986 y = q_gams12_a[2] / ( (v + q_gams12_b[2]) + y );
1987 y = q_gams12_a[1] / ( (v + q_gams12_b[1]) + y ) + q_gams12_b[0];
1988 }
1989 y += 1;
1990 y *= q_ex10(2*x);
1991 times2pown(y,-141);
1992
1993 return y;
1994}
1995
1996real gamma_S12(const real& x)
1997// Calculating approximations for Gamma(x) in [122.5,132.5];
1998// eps = 1.927801e-15;
1999// Blomquist, 15.06.2009;
2000{
2001 real y,Pr;
2002 int n,p;
2003
2004 n = Round(x);
2005 p = n-127;
2006 if (n>127) // neighbour intervals to the right
2007 {
2008 Pr = x-1;
2009 for (int k=2; k<=p; k++)
2010 Pr *= x-k;
2011 y = Pr*gam_S12(x-p);
2012 }
2013 else // left neighbour intervals and S12 itself
2014 {
2015 p = -p; Pr = x;
2016 for (int k=1; k<=p-1; k++)
2017 Pr *= x+k;
2018 y = (p==0)? gam_S12(x) : gam_S12(x+p)/Pr;
2019 }
2020
2021 return y;
2022}
2023
2024real gam_S13(const real& x)
2025// Calculating approximations for Gamma(x) in S13 = [137.5,138.5];
2026// Rel. error bound: 8.563188e-16;
2027// Blomquist, 14.06.2009;
2028{
2029 real y(0),v;
2030
2031 // Continued fraction: K_4(v), v = 1/(x-x0), x0=138.0390625;
2032 if (x==138.0390625)
2033 y = q_gams13_b[0];
2034 else
2035 {
2036 v = 1/(x-138.0390625);
2037
2038 y = q_gams13_a[4] / ( v + q_gams13_b[4]);
2039 y = q_gams13_a[3] / ( (v + q_gams13_b[3]) + y );
2040 y = q_gams13_a[2] / ( (v + q_gams13_b[2]) + y );
2041 y = q_gams13_a[1] / ( (v + q_gams13_b[1]) + y ) + q_gams13_b[0];
2042 }
2043 y += 1;
2044 y *= q_ex10(2*x);
2045 times2pown(y,-137);
2046
2047 return y;
2048}
2049
2050real gamma_S13(const real& x)
2051// Calculating approximations for Gamma(x) in [132.5,142.5];
2052// eps = 1.966542e-15;
2053// Blomquist, 15.06.2009;
2054{
2055 real y,Pr;
2056 int n,p;
2057
2058 n = Round(x);
2059 p = n-138;
2060 if (n>138) // neighbour intervals to the right
2061 {
2062 Pr = x-1;
2063 for (int k=2; k<=p; k++)
2064 Pr *= x-k;
2065 y = Pr*gam_S13(x-p);
2066 }
2067 else // left neighbour intervals and S13 itself
2068 {
2069 p = -p; Pr = x;
2070 for (int k=1; k<=p-1; k++)
2071 Pr *= x+k;
2072 y = (p==0)? gam_S13(x) : gam_S13(x+p)/Pr;
2073 }
2074
2075 return y;
2076}
2077
2078real gam_S14(const real& x)
2079// Calculating approximations for Gamma(x) in S14 = [146.5,147.5];
2080// Rel. error bound: 8.570174e-16;
2081// Blomquist, 15.06.2009;
2082{
2083 real y(0),v;
2084
2085 // Continued fraction: K_5(v), v = 1/(x-x0), x0=146.94921875;
2086 if (x==146.94921875)
2087 y = q_gams14_b[0];
2088 else
2089 {
2090 v = 1/(x-146.94921875);
2091
2092 y = q_gams14_a[5] / ( v + q_gams14_b[5]);
2093 y = q_gams14_a[4] / ( (v + q_gams14_b[4]) + y );
2094 y = q_gams14_a[3] / ( (v + q_gams14_b[3]) + y );
2095 y = q_gams14_a[2] / ( (v + q_gams14_b[2]) + y );
2096 y = q_gams14_a[1] / ( (v + q_gams14_b[1]) + y ) + q_gams14_b[0];
2097 }
2098 y += 1;
2099 y *= q_ex10(2*x);
2100 times2pown(y,-133);
2101
2102 return y;
2103}
2104
2105real gamma_S14(const real& x)
2106// Calculating approximations for Gamma(x) in [142.5,150.5];
2107// eps = 1.745196e-15;
2108// Blomquist, 16.06.2009;
2109{
2110 real y,Pr;
2111 int n,p;
2112
2113 n = Round(x);
2114 p = n-147;
2115 if (n>147) // neighbour intervals to the right:
2116 {
2117 Pr = x-1;
2118 for (int k=2; k<=p; k++)
2119 Pr *= x-k;
2120 y = Pr*gam_S14(x-p);
2121 }
2122 else // left neighbour intervals and S14 itself:
2123 {
2124 p = -p; Pr = x;
2125 for (int k=1; k<=p-1; k++)
2126 Pr *= x+k;
2127 y = (p==0)? gam_S14(x) : gam_S14(x+p)/Pr;
2128 }
2129
2130 return y;
2131}
2132
2133real gam_S15(const real& x)
2134// Calculating approximations for Gamma(x) in S15 = [153.5,154.5];
2135// Rel. error bound: 1.279047e-15;
2136// Blomquist, 16.06.2009;
2137{
2138 real y(0),v;
2139
2140 // Continued fraction: K_5(v), v = 1/(x-x0), x0=153.90234375;
2141 if (x==153.90234375)
2142 y = q_gams15_b[0];
2143 else
2144 {
2145 v = 1/(x-153.90234375);
2146
2147 y = q_gams15_a[5] / ( v + q_gams15_b[5]);
2148 y = q_gams15_a[4] / ( (v + q_gams15_b[4]) + y );
2149 y = q_gams15_a[3] / ( (v + q_gams15_b[3]) + y );
2150 y = q_gams15_a[2] / ( (v + q_gams15_b[2]) + y );
2151 y = q_gams15_a[1] / ( (v + q_gams15_b[1]) + y ) + q_gams15_b[0];
2152 }
2153 y += 1;
2154 v = q_ex10(x);
2155 times2pown(y,-129);
2156 y *= v;
2157 y *= v;
2158
2159 return y;
2160}
2161
2162real gamma_S15(const real& x)
2163// Calculating approximations for Gamma(x) in [150.5,157.5];
2164// eps = 1.945181e-15;
2165// Blomquist, 19.06.2009;
2166{
2167 real y,Pr;
2168 int n,p;
2169
2170 n = Round(x);
2171 p = n-154;
2172 if (n>154) // neighbour intervals to the right
2173 {
2174 Pr = x-1;
2175 for (int k=2; k<=p; k++)
2176 Pr *= x-k;
2177 y = Pr*gam_S15(x-p);
2178 }
2179 else // left neighbour intervals and S15 itself
2180 {
2181 p = -p; Pr = x;
2182 for (int k=1; k<=p-1; k++)
2183 Pr *= x+k;
2184 y = (p==0)? gam_S15(x) : gam_S15(x+p)/Pr;
2185 }
2186
2187 return y;
2188}
2189
2190real gam_S16(const real& x)
2191// Calculating approximations for Gamma(x) in S16 = [160.5,161.5];
2192// Rel. error bound: 1.381209e-15;
2193// Blomquist, 19.06.2009;
2194{
2195 real y(0),v;
2196
2197 // Continued fraction: K_5(v), v = 1/(x-x0), x0=161.08984375;
2198 if (x==161.08984375)
2199 y = q_gams16_b[0];
2200 else
2201 {
2202 v = 1/(x-161.08984375);
2203
2204 y = q_gams16_a[5] / ( v + q_gams16_b[5]);
2205 y = q_gams16_a[4] / ( (v + q_gams16_b[4]) + y );
2206 y = q_gams16_a[3] / ( (v + q_gams16_b[3]) + y );
2207 y = q_gams16_a[2] / ( (v + q_gams16_b[2]) + y );
2208 y = q_gams16_a[1] / ( (v + q_gams16_b[1]) + y ) + q_gams16_b[0];
2209 }
2210 y += 1;
2211 v = q_ex10(x);
2212 times2pown(v,-62);
2213 y *= sqr(v);
2214 return y;
2215}
2216
2217real gamma_S16(const real& x)
2218// Calculating approximations for Gamma(x) in [157.5,164.5];
2219// eps = 2.047343e-15;
2220// Blomquist, 20.06.2009;
2221{
2222 real y,Pr;
2223 int n,p;
2224
2225 n = Round(x);
2226 p = n-161;
2227 if (n>161) // neighbour intervals to the right
2228 {
2229 Pr = x-1;
2230 for (int k=2; k<=p; k++)
2231 Pr *= x-k;
2232 y = Pr*gam_S16(x-p);
2233 }
2234 else // left neighbour intervals and S16 itself
2235 {
2236 p = -p; Pr = x;
2237 for (int k=1; k<=p-1; k++)
2238 Pr *= x+k;
2239 y = (p==0)? gam_S16(x) : gam_S16(x+p)/Pr;
2240 }
2241
2242 return y;
2243}
2244
2245real gam_S17(const real& x)
2246// Calculating approximations for Gamma(x) in: S17 = [167.5,168.5];
2247// Rel. error bound: 1.390239e-15;
2248// Blomquist, 20.06.2009;
2249{
2250 real y(0),v;
2251
2252 // Continued fraction: K_5(v), v = 1/(x-x0), x0=168.0;
2253 if (x==168.0)
2254 y = q_gams17_b[0];
2255 else
2256 {
2257 v = 1/(x-168.0);
2258
2259 y = q_gams17_a[5] / ( v + q_gams17_b[5]);
2260 y = q_gams17_a[4] / ( (v + q_gams17_b[4]) + y );
2261 y = q_gams17_a[3] / ( (v + q_gams17_b[3]) + y );
2262 y = q_gams17_a[2] / ( (v + q_gams17_b[2]) + y );
2263 y = q_gams17_a[1] / ( (v + q_gams17_b[1]) + y ) + q_gams17_b[0];
2264 }
2265 y += 1;
2266 v = q_ex10(x);
2267 times2pown(y,-119);
2268 y *= v;
2269 y *= v;
2270 return y;
2271}
2272
2273real gamma_S17(const real& x)
2274// Calculating approximations for Gamma(x), x in [164.5,171.5];
2275// eps = 2.056373e-15;
2276// Blomquist, 20.06.2009;
2277{
2278 real y,Pr;
2279 int n,p;
2280
2281 n = Round(x);
2282 p = n-168;
2283 if (n>168) // neighbour intervals to the right
2284 {
2285 Pr = x-1;
2286 for (int k=2; k<=p; k++)
2287 Pr *= x-k;
2288 y = Pr*gam_S17(x-p);
2289 }
2290 else // left neighbour intervals and S17 itself
2291 {
2292 p = -p; Pr = x;
2293 for (int k=1; k<=p-1; k++)
2294 Pr *= x+k;
2295 y = (p==0)? gam_S17(x) : gam_S17(x+p)/Pr;
2296 }
2297
2298 return y;
2299}
2300
2301real gamma_05(const real& x)
2302// Calculating Gamma(x), x in (-0.5,+171.5]
2303// eps = 2.082627e-15;
2304// Blomquist, 21.06.2009;
2305{
2306 real y(0);
2307 int Nr;
2308
2309 Nr = int_no(gam_f85,19,x);
2310
2311 switch(Nr)
2312 {
2313 case 0: y = gamma_S0(x); break; // x in (-0.5,8.5);
2314 case 1: y = gamma_S1(x); break; // x in [8.5,16.5);
2315 case 2: y = gamma_S2(x); break; // x in [16.5,24.5);
2316 case 3: y = gamma_S3(x); break; // x in [24.5,35.5);
2317 case 4: y = gamma_S4(x); break; // x in [35.5,46.5);
2318 case 5: y = gamma_S5(x); break; // x in [46.5,57.5);
2319 case 6: y = gamma_S6(x); break; // x in [57.5,68.5);
2320 case 7: y = gamma_S7(x); break; // x in [68.5,79.5);
2321 case 8: y = gamma_S8(x); break; // x in [79.5,90.5);
2322 case 9: y = gamma_S9(x); break; // x in [90.5,101.5);
2323 case 10: y = gamma_S10(x); break; // x in [101.5,112.5);
2324 case 11: y = gamma_S11(x); break; // x in [112.5,122.5);
2325 case 12: y = gamma_S12(x); break; // x in [122.5,132.5);
2326 case 13: y = gamma_S13(x); break; // x in [132.5,142.5);
2327 case 14: y = gamma_S14(x); break; // x in [142.5,150.5);
2328 case 15: y = gamma_S15(x); break; // x in [150.5,157.5);
2329 case 16: y = gamma_S16(x); break; // x in [157.5,164.5);
2330 case 17: y = gamma_S17(x); break; // x in [164.5,171.5);
2331
2332 default: // n=18; x in [171.5, ...];
2333 y = gamma_S17(x);
2334 }
2335
2336 return y;
2337}
2338
2340// Calculating 1/Gamma(x) for x in [-170.0,+171.0];
2341// eps = 2.866906e-15;
2342// Blomquist, 21.06.2009;
2343{
2344 real y;
2345
2346 if (x<-170.0 || x>171.0)
2347 cxscthrow(STD_FKT_OUT_OF_DEF("real gammar(const real& x)"));
2348
2349 if (x<=-0.5)
2350 y = -sinpix_pi(x)*x*gamma_05(-x);
2351 else
2352 if (x<=8.5)
2353 y = gammar_S0(x);
2354 else y = 1/gamma_05(x);
2355
2356 return y;
2357}
2358
2359real gamma(const real& x)
2360// Calculating Gamma(x) for x in [-170.0,+171.5]
2361// eps = 3.088951e-15;
2362// Blomquist, 21.06.2009;
2363{
2364 real y;
2365
2366 if (x>171.5 || x<-170.0)
2367 cxscthrow(STD_FKT_OUT_OF_DEF("real gamma(const real& x)"));
2368
2369 if (x<=-0.5)
2370 y = -1/( sinpix_pi(x)*x*gamma_05(-x) );
2371 else
2372 y = gamma_05(x);
2373
2374 return y;
2375}
2376
2377
2378extern "C" {
2379 void r_lfsr(void) {;} // Siehe real.hpp in real_ari...?!?!
2380}
2381
2382} // namespace cxsc
2383
The Data Type dotprecision.
Definition dot.hpp:112
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
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
interval sinpix_pi(const interval &x)
Calculates ;.
Definition imath.cpp:655
interval expx2(const interval &x)
Calculates .
Definition imath.cpp:234
real Cut25(const real &x)
Returns a real value, which corresponds with the first 25 mantissa bits of x.
Definition rmath.cpp:504
cinterval ln(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:851
int ceil(const real &x) noexcept
Rounding to the smallest integer greater or equal x; -2147483649 < x <= 2147483647....
Definition rmath.cpp:558
interval expmx2(const interval &x)
Calculates .
Definition imath.cpp:192
real Cut24(const real &x)
Returns a real value, which corresponds with the first 24 mantissa bits of x.
Definition rmath.cpp:488
interval acoshp1(const interval &x)
Calculates .
Definition imath.cpp:617
interval arg(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:741
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 expx2m1(const interval &x)
Calculates .
Definition imath.cpp:300
interval ln_sqrtx2y2(const interval &x, const interval &y) noexcept
Calculates .
Definition imath.cpp:581
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
real Cut26(const real &x)
Returns a real value, which corresponds with the first 26 mantissa bits of x.
Definition rmath.cpp:520
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
interval sqrtx2y2(const interval &x, const interval &y) noexcept
Calculates .
Definition imath.cpp:80