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 
30 namespace cxsc {
31 // Constants and values for the function ln_sqrtx2y2:
32 real ln2_1067(6505485212531678.0/8796093022208.0); // 1067*ln(2)
33 // Exponents of the interval bounds:
34 int 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};
37 int 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):
41 int 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};
44 real 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 
69 int 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 
87 void 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);
119 const 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 
124 real 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 
153 real 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 
183 real expx2m1(const real& x)
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 
230 real 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
239 const real sqrtp1m1_s = 9007199254740984.0 / 1125899906842624.0;
240 
241 real 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 
258 real 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 
286 real sqrtx2m1(const real& x)
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 
335 real sqrt1mx2(const real& x)
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 
366 int 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 
381 real 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 
486 typedef union { double f; char intern[8]; } help_real;
487 
488 real 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 
504 real 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 
520 real 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 
536 int 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 
558 int 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 
570 int 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 
585 static 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 
594 static 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 
598 real 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 // ********************************************************************************
625 real 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
733 int 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 
759 real 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 // ********************************************************************************
776 real 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 // ********************************************************************************
800 real 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 // ********************************************************************************
824 real 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 // ********************************************************************************
849 real 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 // ********************************************************************************
874 real 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 
894 real sinpix_pi(const real& x)
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):
934 real 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 
1293 real 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 
1319 real 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 
1345 real 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 
1354 real 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 
1382 real 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 
1411 real 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 
1438 real 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 
1467 real 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 
1494 real 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 
1523 real 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 
1550 real 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 
1579 real 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 
1606 real 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 
1635 real 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 
1662 real 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 
1692 real 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 
1719 real 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 
1748 real 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 
1774 real 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 
1804 real 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 
1830 real 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 
1860 real 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 
1886 real 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 
1915 real 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 
1941 real 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 
1970 real 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 
1996 real 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 
2024 real 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 
2050 real 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 
2078 real 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 
2105 real 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 
2133 real 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 
2162 real 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 
2190 real 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 
2217 real 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 
2245 real 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 
2273 real 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 
2301 real 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 
2339 real gammar(const real& x)
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 
2359 real 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 
2378 extern "C" {
2379  void r_lfsr(void) {;} // Siehe real.hpp in real_ari...?!?!
2380 }
2381 
2382 } // namespace cxsc
2383 
cxsc::gamma
interval gamma(const interval &x)
The Gamma function.
Definition: imath.cpp:1465
cxsc::lnp1
cinterval lnp1(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:867
cxsc::sqrt1px2
cinterval sqrt1px2(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:1071
cxsc::expmx2
interval expmx2(const interval &x)
Calculates .
Definition: imath.cpp:192
cxsc::expx2
interval expx2(const interval &x)
Calculates .
Definition: imath.cpp:234
cxsc::arg
interval arg(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:741
cxsc::Cut26
real Cut26(const real &x)
Returns a real value, which corresponds with the first 26 mantissa bits of x.
Definition: rmath.cpp:520
cxsc::sqrt
cinterval sqrt(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:1007
cxsc::ln
cinterval ln(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:851
cxsc::expx2m1
interval expx2m1(const interval &x)
Calculates .
Definition: imath.cpp:300
cxsc::abs
ivector abs(const cimatrix_subv &mv) noexcept
Returns the absolute value of the matrix.
Definition: cimatrix.inl:737
cxsc::sqrtx2m1
cinterval sqrtx2m1(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:1109
cxsc::Cut25
real Cut25(const real &x)
Returns a real value, which corresponds with the first 25 mantissa bits of x.
Definition: rmath.cpp:504
cxsc::acoshp1
interval acoshp1(const interval &x)
Calculates .
Definition: imath.cpp:617
cxsc::exp
cinterval exp(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:159
cxsc::Cut24
real Cut24(const real &x)
Returns a real value, which corresponds with the first 24 mantissa bits of x.
Definition: rmath.cpp:488
cxsc::Round
int Round(const real &x) noexcept
Rouding to the next integer; |x| < 2147483647.5.
Definition: rmath.cpp:536
cxsc::dotprecision
The Data Type dotprecision.
Definition: dot.hpp:112
cxsc::ceil
int ceil(const real &x) noexcept
Rounding to the smallest integer greater or equal x; -2147483649 < x <= 2147483647....
Definition: rmath.cpp:558
cxsc::ln_sqrtx2y2
interval ln_sqrtx2y2(const interval &x, const interval &y) noexcept
Calculates .
Definition: imath.cpp:581
cxsc::sinpix_pi
interval sinpix_pi(const interval &x)
Calculates ;.
Definition: imath.cpp:655
cxsc
The namespace cxsc, providing all functionality of the class library C-XSC.
Definition: cdot.cpp:29
cxsc::sqrtp1m1
cinterval sqrtp1m1(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:1054
cxsc::ifloor
int ifloor(const real &x) noexcept
Rounding to the greates integer smaller or equal x; -2147483649 < x <= 2147483647....
Definition: rmath.cpp:570
cxsc::expm1
cinterval expm1(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:177
cxsc::times2pown
void times2pown(cinterval &x, int n) noexcept
Fast multiplication of reference parameter [z] with .
Definition: cimath.cpp:2059
cxsc::gammar
interval gammar(const interval &x)
The inverse Gamma function: 1/Gamma(x)
Definition: imath.cpp:1361
cxsc::sqrt1mx2
cinterval sqrt1mx2(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:1140
cxsc::sqr
cinterval sqr(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:3342
cxsc::sqrtx2y2
interval sqrtx2y2(const interval &x, const interval &y) noexcept
Calculates .
Definition: imath.cpp:80
cxsc::real
The Scalar Type real.
Definition: real.hpp:114