C-XSC - A C++ Class Library for Extended Scientific Computing  2.5.4
rmath.cpp
00001 /*
00002 **  CXSC is a C++ library for eXtended Scientific Computing (V 2.5.4)
00003 **
00004 **  Copyright (C) 1990-2000 Institut fuer Angewandte Mathematik,
00005 **                          Universitaet Karlsruhe, Germany
00006 **            (C) 2000-2014 Wiss. Rechnen/Softwaretechnologie
00007 **                          Universitaet Wuppertal, Germany   
00008 **
00009 **  This library is free software; you can redistribute it and/or
00010 **  modify it under the terms of the GNU Library General Public
00011 **  License as published by the Free Software Foundation; either
00012 **  version 2 of the License, or (at your option) any later version.
00013 **
00014 **  This library is distributed in the hope that it will be useful,
00015 **  but WITHOUT ANY WARRANTY; without even the implied warranty of
00016 **  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00017 **  Library General Public License for more details.
00018 **
00019 **  You should have received a copy of the GNU Library General Public
00020 **  License along with this library; if not, write to the Free
00021 **  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00022 */
00023 
00024 /* CVS $Id: rmath.cpp,v 1.30 2014/01/30 17:23:48 cxsc Exp $ */
00025 
00026 #include "rmath.hpp"
00027 
00028 #include "rtsrmath.h"
00029 
00030 namespace cxsc {
00031 // Constants and values for the function ln_sqrtx2y2:
00032 real ln2_1067(6505485212531678.0/8796093022208.0); // 1067*ln(2)
00033 // Exponents of the interval bounds:
00034 int B_lnx2y2_1[22] = {21, 31, 51, 101, 151, 201, 251, 301, 351, 401,
00035                       451, 501, 551, 601, 651, 701, 751, 801, 851, 
00036                       901, 951, 1025}; 
00037 int B_lnx2y2_2[22] = {-1021,-949,-899,-849,-799,-749,-699,-649,-599,
00038                       -549,-499,-449,-399,-349,-299,-249,-199,-149,
00039                       -99,-49,-29,-19}; 
00040 // Optimal values N for N*ln(2):
00041 int B_lnx2y2_N1[21] ={20, 40, 61, 122, 160, 229, 259, 320, 366, 427,
00042                       488, 518, 549, 610, 671, 732, 763, 825, 885, 
00043                       945, 976};
00044 real B_lnx2y2_c1[21] = 
00045 {
00046     7804143460206699.0 / 562949953421312.0, // N*ln(2) with N = 20;
00047     7804143460206699.0 / 281474976710656.0, // N*ln(2) with N = 40;
00048     5950659388407608.0 / 140737488355328.0, // N*ln(2) with N = 61;
00049     5950659388407608.0 / 70368744177664.0,  // N*ln(2) with N = 122;
00050     7804143460206699.0 / 70368744177664.0,  // N*ln(2) with N = 160;
00051     5584840163710419.0 / 35184372088832.0,  // N*ln(2) with N = 229;
00052     6316478613104797.0 / 35184372088832.0,  // N*ln(2) with N = 259;
00053     7804143460206699.0 / 35184372088832.0,  // N*ln(2) with N = 320;
00054     8925989082611412.0 / 35184372088832.0,  // N*ln(2) with N = 366;
00055     5206826964856657.0 / 17592186044416.0,  // N*ln(2) with N = 427;
00056     5950659388407608.0 / 17592186044416.0,  // N*ln(2) with N = 488;
00057     6316478613104797.0 / 17592186044416.0,  // N*ln(2) with N = 518;
00058     6694491811958559.0 / 17592186044416.0,  // N*ln(2) with N = 549;
00059     7438324235509510.0 / 17592186044416.0,  // N*ln(2) with N = 610;
00060     8182156659060461.0 / 17592186044416.0,  // N*ln(2) with N = 671;
00061     8925989082611412.0 / 17592186044416.0,  // N*ln(2) with N = 732;
00062     4652001140732587.0 / 8796093022208.0,   // N*ln(2) with N = 763;
00063     5030014339586349.0 / 8796093022208.0,   // N*ln(2) with N = 825;
00064     5395833564283538.0 / 8796093022208.0,   // N*ln(2) with N = 885;
00065     5761652788980727.0 / 8796093022208.0,   // N*ln(2) with N = 945;
00066     5950659388407608.0 / 8796093022208.0    // N*ln(2) with N = 976;
00067 };
00068 
00069 int uint_trail(const unsigned int& n)
00070 //  Liefert die Anzahl aufeinanderfolgender Null-Bits am Ende von n;
00071 //  n=0 --> 32;   n=1 --> 0;   n=2 --> 1;   n=4 --> 2;   n=1024 --> 10;
00072 {
00073     int m1,p=0;
00074     unsigned int m=n;
00075     if (m==0) p = 32;
00076     else { // m != 0, d.h. die folgende while-Schleife bricht ab!
00077        do
00078        {
00079            m1 = m & 1;  // m1=0(1), wenn letztes Bit von m gleich 0(1) ist.
00080            if (m1==0) p++;
00081            m = m >> 1;  // Bit-Muster um 1 nach rechts schieben.
00082        } while(m1==0);
00083     }
00084     return p;
00085 } // uint_trail
00086 
00087 void sqr2uv(const real& x, real& u, real& v)
00088 // Liefert u,v für: x2 = u + v; EXAKTE Darstellung, falls kein overflow 
00089 // auftritt und v im normalisierten Bereich liegt. u > |v|
00090 // Vorsicht: Funktioniert zunächst nur auf INTEL-Systemen!!!
00091 {
00092     real a,b,t,y1,y2;
00093     a = Cut26(x);
00094     b = x-a;  // x = a+b;
00095     u = x*x;
00096     t = u - a*a; // exakte Auswertung!
00097     y2 = a*b;
00098     times2pown(y2,1); // y2 = 2*a*b, exakt!
00099     t -= y2;
00100 // Jetzt fehlt noch:  t2 - b*b, aber b*b wird nicht immer korrekt berechnet,
00101 // daher nochmalige Aufspaltung von b in y1+y2!!
00102     y1 = Cut25(b);
00103     y2 = b - y1;   // b = y1+y2, exakt;
00104     t -= y1*y1;
00105     if (sign(y2)!=0)
00106     {
00107         a = y1*y2;
00108         times2pown(a,1); // a = 2*y1*y2, exakt!
00109         t -= a;
00110         t -= y2*y2;
00111     }
00112     v = -t;
00113 } // sqr2uv
00114 
00115 
00116 //------------------------------------------------------------------
00117 
00118 // real-Konstante  expmx2_x0  fuer die Funktion e^(-x2);
00119 const real expmx2_x0 = 7491658466053896.0 / 281474976710656.0;
00120 // Fuer x > expmx2_x0 werden die Funktionswerte auf Null gesetzt.
00121 // Die relative Fehlerschranke e(f) := 4.618919E-16 gilt fuer
00122 // alle |x| <= expmx2_x0 = 26.61571750925.... 
00123 
00124 real expmx2(const real& x) throw()
00125 // e^(-x^2);  rel. Fehlerschranke:  eps = 4.618919E-16 = e(f) gilt
00126 // fuer alle |x| <= expmx2_x0 = 26.61571750925....
00127 // Fuer |x| > expmx2_x0 --> expmx2(x) = 0;
00128 // Blomquist, 05.07.04;
00129 { 
00130     real t=x,u,v,res=0;
00131     int ex;
00132     if (t<0) t = -t;  // t >= 0;
00133     ex = expo(t);
00134     if (ex<=-26) res = 1;  // t < 2^(-26)
00135     else if (ex<=-6)       // t < 2^(-6)
00136     {
00137         u = t*t;  v = u; 
00138         times2pown(v,-1);  // v: 0.5*x2
00139         res = 1-u*( 1-v*(1-u/3) );
00140     } else if (t<=expmx2_x0) {
00141         sqr2uv(x,u,v);  // u:= x*x,v aus S(2,53); x2 = u+v (exakt!)
00142         res = exp(-u); 
00143         if (v!=0) { 
00144             times2pown(res,500);  // Die Skalierung verhindert, dass 
00145             v *= res; // v*exp(-u) in den denormalisierten Bereich faellt
00146             res -= v;
00147             times2pown(res,-500); // Rueckskalierung
00148         }
00149     }
00150     return res;
00151 } // expmx2
00152 
00153 real expx2(const real& x)
00154 // e^(+x^2);  rel. Fehlerschranke:  eps = 4.618958E-16 = e(f) gilt
00155 // fuer alle |x| <= x0 = 26.64174755704632....
00156 // x0 = 7498985273150791.0 / 281474976710656.0;
00157 // Fuer |x| > x0 --> Programmabbruch
00158 // Ausfuehrlich getestet;  Blomquist, 26.07.06;
00159 {
00160     real t=x,u,v,res;
00161     int ex;
00162     if (t<0) t = -t;  // t >= 0;
00163     ex = expo(t);
00164     if (ex<=-26) res = 1;  // t < 2^(-26)
00165     else 
00166         if (ex<=-6)       // t < 2^(-6)
00167         {
00168             u = t*t;  v = u;
00169             times2pown(v,-1);  // v: 0.5*x^2
00170             res = 1+u*( 1+v*(1+u/3) );
00171         } 
00172         else 
00173         {
00174             sqr2uv(x,u,v);  // u := x*x und v aus S(2,53); 
00175                             // x^2 = u+v ist exakt!
00176             res = exp(u);
00177             v *= res;       // v*exp(+u)
00178             res += v;
00179         }
00180     return res;
00181 } // expx2
00182 
00183 real expx2m1(const real& x)
00184 // e^(+x^2)-1; rel. Fehlerschranke: eps = 4.813220E-16 = e(f) gilt
00185 // fuer alle x, mit:  2^(-511) <= x <= x0 = 26.64174755704632....
00186 // x0 = 7498985273150791.0 / 281474976710656.0;
00187 // Fuer x > x0 --> Programmabbruch wegen Overflow;
00188 // Fuer 0 < x < 2^(-511) --> Programmabbruch, denorm. Bereich!!
00189 // Ausfuehrlich getestet;  Blomquist, 10.08.2006;
00190 {
00191     real t(x),u,v,y,res(0);
00192     int ex;
00193     if (t<0) t = -t;  // t >= 0;
00194 
00195     if (t>=6.5) res = expx2(t);
00196     else 
00197     {
00198         ex = expo(t);
00199         sqr2uv(x,u,v);  // u := x*x und v aus S(2,53); 
00200         if (ex>=2) // g4(x)
00201         {
00202             y = exp(u); 
00203             res = 1 - v*y;
00204             res = y - res;
00205         }
00206         else 
00207             if (ex>=-8) res = expm1(u) + v*exp(u); // g3(x)
00208             else 
00209                 if (ex>=-25) { // g2(x)
00210                     y = u*u;
00211                     times2pown(y,-1);
00212                     res = (1+u/3)*y + u;
00213                 }
00214                 else 
00215                     if(ex>=-510) res = u;  // g1(x)
00216                     else if (ex>=-1073) 
00217                     {
00218                         std::cerr << "expx2m1: denormalized range!" 
00219                              << std::endl; exit(1);
00220                     }
00221     }
00222 
00223     return res;
00224 } // expx2m1
00225 
00226 
00227 
00228 //------------------------------------------------------------------
00229 
00230 real sqrt1px2(const real& x) throw()
00231 // sqrt(1+x^2); Blomquist 13.12.02;
00232 {
00233     if (expo(x) > 33) return abs(x);
00234     else return sqrt(1+x*x);
00235 }
00236 
00237 // Folgende Konstante sqrtp1m1_s wird gebraucht in 
00238 // real sqrtp1m1(const real& x) throw();  Blomquist, 05.08.03
00239 const real sqrtp1m1_s = 9007199254740984.0 / 1125899906842624.0;
00240 
00241 real sqrtp1m1(const real& x) throw()
00242 // sqrtp1m1(x) = sqrt(x+1)-1;
00243 // Blomquist, 05.08.03; 
00244 {
00245     real y = x;
00246     int ex = expo(x);
00247     if (ex<=-50) times2pown(y,-1);  // |x|<2^(-50); fast division by 2 
00248     else if (ex>=105) y = sqrt(x);  // x >= 2^(+104) = 2.02824...e+31
00249     else if (ex>=53) y = sqrt(x)-1; // x >= 2^(+52)  = 4.50359...e+15
00250     else if (x>-0.5234375 && x<=sqrtp1m1_s ) y = x / (sqrt(x+1) + 1);
00251     else y = sqrt(x+1)-1;
00252     return y;
00253 }
00254 
00255 
00256 //------------------------------------------------------------------
00257 
00258 real sqrtx2y2(const real& x, const real& y) throw()
00259 // calculating sqrt(x^2 + y^2) in high accuracy. Blomquist 01.12.02
00260 {
00261     real a,b,r;
00262     dotprecision dot;
00263     int exa,exb,ex;
00264     a = x;   b = y;
00265     exa = expo(a);  exb = expo(b);
00266     if (exb > exa)
00267     { // Permutation of a,b:
00268         r = a;  a = b;  b = r;
00269         ex = exa;  exa = exb;  exb = ex;
00270     }  // |a| >= |b|
00271     // a = |a| >= |b| > 0:
00272     ex = 511 - exa;   // scaling a with 2^ex --> expo(a) == 511, --> a*a and
00273     times2pown(a,ex); // b*b without overflow. An underflow of b*b will not
00274     times2pown(b,ex); // affect the accuracy of the result!
00275     dot = 0;
00276     accumulate(dot,a,a);
00277     accumulate(dot,b,b);
00278     r = rnd(dot);
00279     r = sqrt(r); // sqrt(...) declared in rmath.hpp
00280     times2pown(r,-ex); // back-scaling
00281     return r;
00282 } // sqrtx2y2
00283 
00284 //------------------------------------------------------------------
00285 
00286 real sqrtx2m1(const real& x) 
00287 // sqrt(x^2-1);  rel. Fehlerschranke:  eps = 2.221305E-16 = e(f)
00288 // Blomquist, 13.04.04;
00289 {    const real c1 = 1.000732421875, 
00290                c2 = 44000.0,
00291                c3 = 1024.0;  // c1,c2,c3 werden exakt gespeichert!
00292     real res,ep,ep2,s1,s2,x1,x2,arg=x;
00293     if (sign(arg)<0) arg = -arg; // arg = |x| >= 0
00294     if (arg <= c1) { // x = 1+ep; x^2-1 = 2*ep + ep^2;
00295         ep = x - 1; // Differenz rundungsfehlerfrei!
00296         ep2 = ep*ep;  // ep2 i.a. fehlerbehaftet!
00297         times2pown(ep,1);  // ep = 2*ep;
00298         res = sqrt(ep+ep2); // res=y0: Startwert;
00299         // x - y0^2 = (2*eps - s1^2) + [eps^2 - s2*(y0 + s1)]
00300         s1 = Cut26(res);  s2 = res - s1; // Startwert y0 = s1 + s2;
00301         arg = ep - s1*s1;  // arg = 2*eps - s1^2;
00302         arg += (ep2 - s2*(res+s1));  // arg = x - y0^2
00303         if (sign(arg)>0) {
00304         arg = arg / res;
00305         times2pown(arg,-1);
00306         res += arg;  // 1. Newton-Schritt beendet; eps = 2.221261E-16
00307         }
00308     } else if (arg<c2) { 
00309         // x-y0^2 = [(x1^2-1)-s1^2] + [x2*(x+x1)-s2*(y0+s1)]
00310         x1 = Cut26(arg);  x2 = arg - x1;  // arg = x = x1 + x2;
00311         ep2 = x2*(arg+x1);  // ep2 ist fehlerbehaftet
00312         x2 = x1*x1;  ep = x2-1;
00313         res = sqrt(ep+ep2); // res ist Startwert für Newton-Verfahren
00314         s1 = Cut26(res);  s2 = res - s1;  // Startwert = s1 + s2;
00315         ep2 = ep2 - s2 * (res+s1); // ep2 = [x2*(x+x1)-s2*(y0+s1)]
00316         if (arg<c3) ep -= s1*s1;   // ep = (x1^2-1) - s1^2;
00317         else {
00318             x2 -= s1*s1;  // x2 = x1^2-s1^2
00319             ep = x2 - 1; }         // ep = (x1^2-s1^2) - 1;
00320         ep += ep2;        // ep = x - y0^2;
00321         ep /= res;
00322         times2pown(ep,-1);
00323         res = res + ep;  // 1. Newton-Schritt in hoher Genauigkeit
00324                          // beendet;  eps = 2.221305E-16
00325     } else { // arg = |x| >= 44000;
00326         res = -1/arg;
00327         times2pown(res,-1); // Multiplikation mit 0.5;
00328         res += arg;  // res = x - 1/(2*x);  eps = 2.221114E-16
00329     }
00330     return res;
00331 } // sqrtx2m1 (Punktargumente)
00332 
00333 //------------------------------------------------------------------
00334 
00335 real sqrt1mx2(const real& x) throw(STD_FKT_OUT_OF_DEF)
00336 // sqrt(1-x2);  rel. Fehlerschranke:  eps = 3.700747E-16 = e(f)
00337 // Blomquist, 19.06.04;
00338 { 
00339     real t=x,res;
00340     int ex;
00341     if (sign(t)<0) t = -t; // Argument t >=0;
00342     if (t>1) cxscthrow(STD_FKT_OUT_OF_DEF("real sqrt1mx2(const real&)"));
00343     // For argument t now it holds: 0 <= t <=1;
00344     ex = expo(t);
00345     if (ex<=-26) res = 1; // t < 2^(-26) --> res = 1
00346     else if (ex<=-15) {   // t < 2^(-15) --> res = 1-x2/2
00347         res = t*t;  
00348         times2pown(res,-1);
00349         res = 1 - res;
00350     } else {
00351         if (ex>=0) 
00352         {  // ex>=0 --> t>=0.5; 
00353             res = 1-t;   // res: delta = 1-t;
00354             t = res * res;
00355             times2pown(res,1); // res: 2*delta
00356             res = res - t;     // res: 1-x2 = 2*delta - delta2
00357         } else res = 1-t*t;    // res: Maschinenwert von 1-x2
00358         res = sqrt(res);       // res: Nullte Naeherung 
00359     }
00360     return res;
00361 } // sqrt1mx2
00362 
00363 
00364 //------------------------------------------------------------------
00365 
00366 int Interval_Nr(int* v, const int& n, const int& ex)
00367 // n>0 subintervals:   |...\|...\|...\ ..... |...|
00368 // subinterval Nr.:      0    1    2   .....  n-1  
00369 {
00370     int i=0,j=n,k;  // n>0:  Number of subintervals
00371     do {
00372         k = (i+j)/2;
00373         if (ex < v[k]) j = k-1;
00374         else i = k+1;
00375     } while(i<=j);
00376     return j;  // x with ex=expo(x) lies in the subinterval number j  
00377 }
00378 
00379 //------------------------------------------------------------------
00380 
00381 real ln_sqrtx2y2(const real& x, const real& y) 
00382                                 throw(STD_FKT_OUT_OF_DEF)
00383 // ln( sqrt(x^2+y^2) ) == 0.5*ln(x^2+y^2); Blomquist, 21.11.03;
00384 // Relative error bound: 5.160563E-016;
00385 // Absolute error bound: 2.225075E-308; if x=1 and 0<=y<=b0;
00386 {
00387     int j,N;
00388     real a,b,r,r1;
00389     dotprecision dot;
00390    
00391     a = sign(x)<0 ? -x : x;  // a = |x| >= 0;
00392     b = sign(y)<0 ? -y : y;  // b = |y| >= 0;
00393     int exa=expo(a), exb=expo(b), ex;
00394     if (b > a)
00395     {
00396         r = a;   a = b;   b = r;
00397         ex = exa;   exa = exb;   exb = ex;
00398     }
00399     // It holds now:   0 <= b <= a 
00400     if (sign(a)==0)
00401         cxscthrow(STD_FKT_OUT_OF_DEF
00402                     ("real ln_sqrtx2y2(const real&, const real&)"));
00403     if (exa>20) // to avoid overflow by calculating a^2 + b^2
00404     {  // a>=2^(20):
00405         j = Interval_Nr(B_lnx2y2_1,21,exa); // j: No. of subinterval
00406         N = B_lnx2y2_N1[j];    // N: Optimal int value
00407         if (exb-exa > -25)
00408         {   // For (exb-exa>-25) we use the complete term:  
00409             // N*ln(2) + [ln(2^(-N)*a)+0.5*ln(1+(b/a)^2)]
00410             b = b/a;  // a > 0
00411             b = lnp1(b*b);
00412             times2pown(b,-1);  // exact division by 2
00413             times2pown(a,-N);
00414             r = b + ln(a); // [ ... ] calculated!
00415             r += B_lnx2y2_c1[j];
00416         }
00417         else { // For (exb-exa<=-25) only two summands!: 
00418             times2pown(a,-N);
00419             r = ln(a) + B_lnx2y2_c1[j];
00420         }
00421     }
00422     else  // exa<=20 or a<2^(20):
00423     {     // Now calculation of a^2+b^2 without overflow:
00424         if (exa<=-20) // to avoid underflow by calculating a^2+b^2
00425             if (exa<=-1022) // a in the denormalized range
00426             {
00427                 r = b/a;
00428                 r = lnp1(r*r);  times2pown(r,-1); // r: 0.5*ln(1+..)
00429                 times2pown(a,1067);
00430                 r += ln(a);     // [ .... ] ready
00431                 r -= ln2_1067;  // rel. error = 2.459639e-16;
00432             }
00433             else // MinReal=2^(-1022) <= a < 2^(-20)
00434             {    // Calculating the number j of the subinterval:
00435                 j = 20 - Interval_Nr(B_lnx2y2_2,21,exa);
00436                 r = a;  times2pown(r,B_lnx2y2_N1[j]);
00437                 r = ln(r);  // r: ln(2^N*a);
00438                 if (exb-exa > -25) { // calculating the complete term
00439                     b = b/a;
00440                     a = lnp1(b*b);
00441                     times2pown(a,-1);
00442                     r += a;  // [ ... ] ready now
00443                 }
00444                 // We now have: exb-exa<=-25,  ==>  b/a <= 2^(-24);
00445                 r -= B_lnx2y2_c1[j]; // 0.5*ln(1+(b/a)^2) neglected!
00446                 // relative error = 4.524090e-16 in both cases;
00447             }
00448         else // calculation of a^2+b^2 without overflow or underflow:
00449         {   // exa>-20  respective  a>=2^(-20):
00450             dot = 0;
00451             accumulate(dot,a,a);
00452             accumulate(dot,b,b);  // dot = a^2+b^2, exact!
00453             real s = rnd(dot);    // s = a^2 + b^2, rounded!
00454             if (s>=0.25 && s<=1.75)
00455                 if (s>=0.828125 && s<=1.171875)
00456                 { // Series:
00457                     if (a==1 && exb<=-28)
00458                     {
00459                         r = b;  times2pown(r,-1);
00460                         r *= b;
00461                     }
00462                     else {
00463                         dot -= 1;
00464                         r = rnd(dot); // r = a^2+b^2-1 rounded!
00465                         r = lnp1(r);
00466                         times2pown(r,-1);
00467                     }
00468                 }
00469                 else { // Reading dot = a^2+b^2 twice:
00470                     r = rnd(dot);
00471                     dot -= r;
00472                     r1 = rnd(dot); // a^2+b^2 = r+r1, rounded!
00473                     r1 = lnp1(r1/r);
00474                     r = ln(r) + r1;
00475                     times2pown(r,-1); // exact division by 2
00476                 }
00477             else { // calculating straight from: 0.5*ln(x^2+y^2)
00478                 r = ln(s);
00479                 times2pown(r,-1);
00480             } 
00481         }
00482     }
00483     return r;
00484 } // ln_sqrtx2y2
00485 
00486 typedef union { double f; char intern[8]; } help_real;
00487 
00488 real Cut24(const real& x){
00489     // y = Cut24(x) liefert ein y, das mit den ersten 24 Mantissenbits
00490     // von x übereinstimmt, das hidden bit ist dabei mitgezählt!
00491     // Die restlichen 53-24=29 Mantissenbits werden auf Null gesetzt.
00492   help_real y;
00493   y.f = _double(x);
00494 #if INTEL
00495   y.intern[3] &= 224;
00496   y.intern[0] = y.intern[1] = y.intern[2] = 0;
00497 #else
00498   y.intern[4] &= 224;
00499   y.intern[7] = y.intern[6] = y.intern[5] = 0;
00500 #endif
00501   return real(y.f);
00502 }
00503 
00504 real Cut25(const real& x){
00505     // y = Cut25(x) liefert ein y, das mit den ersten 25 Mantissenbits
00506     // von x übereinstimmt, das hidden bit ist dabei mitgezählt!
00507     // Die restlichen 53-25=28 Mantissenbits werden auf Null gesetzt.
00508   help_real y;
00509   y.f = _double(x);
00510 #if INTEL
00511   y.intern[3] &= 240;
00512   y.intern[0] = y.intern[1] = y.intern[2] = 0;
00513 #else
00514   y.intern[4] &= 240;
00515   y.intern[7] = y.intern[6] = y.intern[5] = 0;
00516 #endif
00517   return real(y.f);
00518 }
00519 
00520 real Cut26(const real& x){
00521     // y = Cut26(x) liefert ein y, das mit den ersten 26 Mantissenbits
00522     // von x übereinstimmt, das hidden bit ist dabei mitgezählt!
00523     // Die restlichen 53-26=27 Mantissenbits werden auf Null gesetzt.
00524   help_real y;
00525   y.f = _double(x);
00526 #if INTEL
00527   y.intern[3] &= 248;
00528   y.intern[0] = y.intern[1] = y.intern[2] = 0;
00529 #else
00530   y.intern[4] &= 248;
00531   y.intern[7] = y.intern[6] = y.intern[5] = 0;
00532 #endif
00533   return real(y.f);
00534 }
00535 
00536 int Round(const real& x) throw()
00537 // y = Round(x) delivers the rounded value y of x.
00538 // For |x| < 2147483647.5 the assignment  y = Round(x)  delivers
00539 // the next integer value y.
00540 // Examples:
00541 // y = Round(-0.1);            --->  y = 0;
00542 // y = Round(-123.5);          --->  y = -124.0;
00543 // y = Round(+123.5);          --->  y = +124.0;
00544 // y = Round(+123.499);        --->  y = +123.0;                
00545 // y = Round(+2147483647.499); --->  y = +2147483647.0;
00546 // y = Round(+2147483647.5);   --->  y = -2147483648.0;
00547 // y = Round(-2147483647.499); --->  y = -2147483647.0;
00548 // y = Round(-2147483647.5);   --->  y = -2147483648.0;
00549 // y = Round(-2147483648.501); --->  y = -2147483648.0; wrong rounding!
00550 // Blomquist, 29.10.2008;
00551 {
00552         double dbl;
00553         
00554         dbl = _double(x);
00555         return dbl<0 ? int(dbl-0.5) : int(dbl+0.5);
00556 }
00557 
00558 int ceil(const real& x) throw()
00559 {
00560         real y(x); 
00561         bool neg(y<0);
00562         if (neg) y = -y; // y >= 0;
00563         int res = int( _double(y) );
00564         y = y - res;
00565         if (!neg && y>0) res += 1;
00566         if (neg) res = -res;
00567         return res;
00568 }
00569 
00570 int ifloor(const real& x) throw()
00571 {
00572         real y(x); 
00573         bool neg(y<0);
00574         if (neg) y = -y; // y >= 0;
00575         int res = int( _double(y) );
00576         y = y - res;
00577         if (neg) 
00578         {
00579                 res = -res;
00580                 if (y>0) res -= 1;
00581         }
00582         return res;
00583 }
00584 
00585 static real q_acoshp1[5] = // Polynomial coefficients of Q_4(x) 
00586            // roundet to nearest. acosh(1+x) = sqrt(2*x)*Q_4(x)
00587 { 1.0 / 1.0,                                    // q_acoshp1[0]
00588   -6004799503160661.0 / 72057594037927936.0,
00589   +5404319552844595.0 / 288230376151711744.0,
00590   -6433713753386423.0 / 1152921504606846976.0,
00591   +8756999275442631.0 / 4611686018427387904.0   // q_acoshp1[4]
00592 };
00593 
00594 static const real c_ln2_B = 6243314768165359.0 / 9007199254740992.0;
00595 // c_ln2_B < ln(2) is the nearest machine number for ln(2) with an
00596 // absolute error < 2.3190469E-17;
00597 
00598 real acoshp1(const real& x) throw()
00599 // acoshp1(x) = acosh(1+x);  rel. error: eps = 7.792706E-16 = e(f)
00600 // Ausfuehrlich getestet;  Blomquist, 27.03.05;
00601 { 
00602     real res;
00603     int ex(expo(x));
00604     if (x<0) 
00605         cxscthrow(STD_FKT_OUT_OF_DEF("real acoshp1(const real&)"));
00606     // For argument x now it holds: 0 <= x <= MaxReal;
00607     if (ex<=-50) res = sqrt(2*x); // 0<=x<2^(-50): acoshp1(x)=sqrt(2x)
00608     else if (ex<=-9) // 2^(-50)<=x<2^{-9}: acoshp1(x)=sqrt(2x)*Q_4(x)
00609       res = sqrt(2*x)*((((q_acoshp1[4]*x+q_acoshp1[3])*x+q_acoshp1[2])
00610                             *x+q_acoshp1[1])*x + q_acoshp1[0]);
00611     else if (ex<=0) res = lnp1(x+sqrt(2*x+x*x));     // range A_3
00612     else if (ex<=50) res = lnp1(x*(1+sqrt(1+2/x)));  // range A_4
00613     else if (ex<=1022) res = ln(2*x);                // range A_5
00614     else res = ln(x) + c_ln2_B;                      // range A_6
00615 
00616     return res;
00617 } // acoshp1
00618 
00619 
00620 // *****************************************************************************
00621 //                               sin(Pi*x)/Pi                                  *
00622 // ********************************************************************************
00623 //   Sinpix_pi:           Teilpunkte des Intervalls [0,0.5]:
00624 // ********************************************************************************
00625 real a_sinpix_pi[8] = {0,
00626                        7.450580596923828125e-9, // 2^(-27)
00627                        0.05078125,
00628                        0.1015625,
00629                        0.15234375,
00630                        0.2578125,
00631                        0.3828125,
00632                        0.5 };
00633 // ********************************************************************************
00634 
00635 // ********************************************************************************
00636 //  Sinpix_pi:      a_k,b_k des Kettenbruchs K2 in A1 = [2^-27,0.05078125]:
00637 // ********************************************************************************
00638   const real q_sin1_a[3] = {0.0,
00639                             -7408124450506663.0 / 4503599627370496.0,
00640                             +4595816915236340.0 / 36028797018963968.0 };
00641 
00642   const real q_sin1_b[3] = {+1.0,
00643                             +8889749340277122.0 / 18014398509481984.0,
00644                             -6103596185201565.0 / 36028797018963968.0};
00645 // ********************************************************************************
00646 
00647 // ********************************************************************************
00648 //  Sinpix_pi:     a_k,b_k des Kettenbruchs K4 in A2 = [0.05078125,0.1015625]:
00649 // ********************************************************************************
00650    const real q_sin2_a[5] = {0.0,
00651                              -4829370543434630.0 / 18014398509481984.0,
00652                              +5229154385037560.0 / 140737488355328.0,
00653                              -6000119407941526.0 / 576460752303423488.0,
00654                              +4592842702840413.0 / 1125899906842624.0 };
00655 
00656   const real q_sin2_b[5] = {-6359670668682560.0 / 576460752303423488.0,
00657                             -6771299865719208.0 / 1125899906842624.0,
00658                             +6861350950372367.0 / 1125899906842624.0,
00659                             -4669728680615660.0 / 2251799813685248.0,
00660                             +4607746187351458.0 / 2251799813685248.0 };
00661 // ********************************************************************************
00662 
00663 // ********************************************************************************
00664 //   Sinpix_pi:    a_k,b_k des Kettenbruchs K4 in A3 = [0.1015625,0.15234375]:
00665 // ********************************************************************************
00666   const real q_sin3_a[5] = {0.0,
00667                                  -7954137925457403.0 / 18014398509481984.0,
00668                                                                     +7534721222201852.0 / 562949953421312.0,
00669                                                 -8481413201661564.0 / 288230376151711744.0,
00670                                                           +6484759126571114.0 / 4503599627370496.0};
00671   const real q_sin3_b[5] = {-8780869688499296.0 / 288230376151711744.0,
00672                                  -7929689932517363.0 / 2251799813685248.0,
00673                                                                     +8223190551048856.0 / 2251799813685248.0,
00674                                                 -5789057670884168.0 / 4503599627370496.0,
00675                                                           +5586166285064551.0 / 4503599627370496.0 };
00676 // ********************************************************************************
00677 
00678 // ********************************************************************************
00679 //   Sinpix_pi:    a_k,b_k des Kettenbruchs K5 in A4 = [0.15234375,0.2578125]:
00680 // ********************************************************************************
00681   const real q_sin4_a[6] = {+0.0,
00682                                  -5469389903424399.0 / 9007199254740992.0,
00683                                                                     +7704524171392577.0 / 1125899906842624.0,
00684                                                 -8513632445057187.0 / 144115188075855872.0,
00685                                                           +6447776473700967.0 / 9007199254740992.0,
00686                                                      -7362728796775514.0 / 288230376151711744.0 };
00687 
00688   const real q_sin4_b[6] = {-8529339040415872.0 / 144115188075855872.0,
00689                                  -5452397600582528.0 / 2251799813685248.0,
00690                                                                     +5848879160432416.0 / 2251799813685248.0,
00691                                                 -8609222354677976.0 / 9007199254740992.0,
00692                                                           +8031029418839465.0 / 9007199254740992.0,
00693                                                      -4904826237544497.0 / 9007199254740992.0 };
00694 // ********************************************************************************
00695 
00696 // ********************************************************************************
00697 //   Sinpix_pi:     a_k,b_k des Kettenbruchs K5 in A5 = [0.2578125,0.3828125]:
00698 // ********************************************************************************
00699   const real q_sin5_a[6] = {0.0,
00700                             4*(5276398025110506.0 / 9007199254740992.0),
00701                             +7169490078506431.0 / 1125899906842624.0,
00702                             -7877917296929161.0 / 36028797018963968.0,
00703                             +5470058796541645.0 / 9007199254740992.0,
00704                             -7233501679212276.0 / 144115188075855872.0 };
00705 
00706   const real q_sin5_b[6] = {+4598161403289824.0 / 144115188075855872.0,
00707                             +4893639363223553.0 / 2251799813685248.0,
00708                             -5525704702280927.0 / 2251799813685248.0,
00709                             +4608444394217766.0 / 4503599627370496.0,
00710                             -8018266172128683.0 / 9007199254740992.0,
00711                             +5822429161405618.0 / 9007199254740992.0 };
00712 // ********************************************************************************
00713 
00714 // ********************************************************************************
00715 //   Sinpix_pi:   a_k,b_k des Kettenbruchs K5 in A6 = [0.3828125,0.5]:
00716 // ********************************************************************************
00717   const real q_sin6_a[6] = {0.0,
00718                             +7461973733147728.0 / 36028797018963968.0,
00719                             +7979710108639590.0 / 140737488355328.0,
00720                             -6289496525317218.0 / 288230376151711744.0,
00721                             +6964694082726429.0 / 1125899906842624.0,
00722                             -5522978779702360.0 / 1152921504606846976.0 };
00723 
00724   const real q_sin6_b[6] = {+5609829449006976.0 / 18014398509481984.0,
00725                             +8354019229473476.0 / 1125899906842624.0,
00726                             -8475200820952730.0 / 1125899906842624.0,
00727                             +5829377160898302.0 / 2251799813685248.0,
00728                             -5714036688191727.0 / 2251799813685248.0,
00729                             +7900583259891052.0 / 4503599627370496.0 };
00730 // ********************************************************************************
00731 
00732 // interne Funktion, wird fuer Gammafunktion benoetigt
00733 int int_no(real *a,
00734            const int n,    // n: Anzahl der Elemente des Feldes a
00735            const real& x)  // x: Eine real-Zahl
00736 // Ein Intervall [A,B] wird durch ein Feld a mit n Elementen in
00737 // n-1 Teilintervalle unterteilt. Fuer ein x aus [A,B] wird eine
00738 // Intervall-Nummer zurueckgegeben, wobei das erste Teilintervall
00739 // die Nr. 0 und das letzte Teilintervall die Nr. n-2 erhaelt.
00740 // Fuer  x<A  ist Nr. = -1, und fuer  x>=B  gilt  Nr. = n-1; 
00741 {
00742     int i,j,k;
00743     i = 0;
00744     j = n-1;
00745     do 
00746     {
00747                 k = (i+j)/2;
00748                 if (x<a[k]) j = k-1;
00749                 else i = k+1;
00750     } 
00751     while (i <= j);
00752 
00753     return j;
00754 }
00755 // *****************************************************************************
00756 //    Sinpix_pi:         Approximation in A1 = [2^-27,0.05078125]:
00757 // *****************************************************************************
00758 
00759 real sinpi_A1(const real& x)
00760 {
00761         real y,v;
00762         
00763         v = 1/(x*x);
00764         // Approximation:  sin(pi*x)/pi = x*K2(v);
00765         y = q_sin1_a[2] / (  v + q_sin1_b[2]);
00766         y = q_sin1_a[1] / ( (v + q_sin1_b[1]) + y);
00767         y *= x;
00768         y += x;
00769         
00770         return y;
00771 }
00772 
00773 // ********************************************************************************
00774 //    Sinpix_pi:        Approximation in A2 = [0.05078125,0.1015625]:
00775 // ********************************************************************************
00776 real sinpi_A2(const real& x)
00777 {
00778         real y,v;
00779         
00780         if (x == 0.08203125) 
00781                 y = q_sin2_b[0];
00782         else
00783         {
00784                 v = 1/(x-0.08203125);
00785                 // Approximation:  sin(pi*x)/pi = x*K4(v);
00786                 y = q_sin2_a[4] / (  v + q_sin2_b[4]);
00787                 y = q_sin2_a[3] / ( (v + q_sin2_b[3]) + y);
00788                 y = q_sin2_a[2] / ( (v + q_sin2_b[2]) + y);
00789                 y = q_sin2_a[1] / ( (v + q_sin2_b[1]) + y) + q_sin2_b[0];
00790         }
00791         y *= x;
00792         y += x;
00793         
00794         return y;
00795 }
00796 
00797 // ********************************************************************************
00798 //    Sinpix_pi:        Approximation in A3 = [0.1015625,0.15234375]:
00799 // ********************************************************************************
00800 real sinpi_A3(const real& x)
00801 {
00802         real y,v;
00803         
00804         if (x == 0.13671875)
00805                 y = q_sin3_b[0];
00806         else
00807         {
00808                 v = 1/(x-0.13671875);
00809                 // Approximation:  sin(pi*x)/pi = x*K4(v);
00810                 y = q_sin3_a[4] / (  v + q_sin3_b[4]);
00811                 y = q_sin3_a[3] / ( (v + q_sin3_b[3]) + y);
00812                 y = q_sin3_a[2] / ( (v + q_sin3_b[2]) + y);
00813                 y = q_sin3_a[1] / ( (v + q_sin3_b[1]) + y) + q_sin3_b[0];
00814         }
00815         y *= x;
00816         y += x;
00817                 
00818         return y;
00819 }
00820 
00821 // ********************************************************************************
00822 //   Sinpix_pi:        Approximation in A4 = [0.15234375,0.2578125]:
00823 // ********************************************************************************
00824 real sinpi_A4(const real& x)
00825 {
00826         real y,v;
00827         
00828         if (x == 0.19140625)
00829                 y = q_sin4_b[0];
00830         else
00831         {
00832                 v = 1/(x-0.19140625);
00833                 // Approximation:  sin(pi*x)/pi = x*K5(v);
00834                 y = q_sin4_a[5] / (  v + q_sin4_b[5]);
00835                 y = q_sin4_a[4] / ( (v + q_sin4_b[4]) + y);
00836                 y = q_sin4_a[3] / ( (v + q_sin4_b[3]) + y);
00837                 y = q_sin4_a[2] / ( (v + q_sin4_b[2]) + y);
00838                 y = q_sin4_a[1] / ( (v + q_sin4_b[1]) + y) + q_sin4_b[0];
00839         }
00840         y *= x;
00841         y += x;
00842         
00843         return y;
00844 }
00845 
00846 // ********************************************************************************
00847 //   Sinpix_pi:        Approximation in A5 = [0.2578125,0.3828125]:
00848 // ********************************************************************************
00849 real sinpi_A5(const real& x)
00850 {
00851         real y,v;
00852         
00853         if (x == 0.30078125)
00854                 y = q_sin5_b[0];
00855         else
00856         {
00857                 v = 1/(x-0.30078125);
00858                 // Approximation:  sin(pi*x)/pi = K5(v);
00859                 y = q_sin5_a[5] / (  v + q_sin5_b[5]);
00860                 y = q_sin5_a[4] / ( (v + q_sin5_b[4]) + y);
00861                 y = q_sin5_a[3] / ( (v + q_sin5_b[3]) + y);
00862                 y = q_sin5_a[2] / ( (v + q_sin5_b[2]) + y);
00863                 y = q_sin5_a[1] / ( (v + q_sin5_b[1]) + y) + q_sin5_b[0];
00864         }
00865         y +=1.0;
00866         times2pown(y,-2); // Division durch 4
00867         
00868         return y;
00869 }
00870 
00871 // ********************************************************************************
00872 //    Sinpix_pi:         Approximation in A6 = [0.3828125,0.5]:
00873 // ********************************************************************************
00874 real sinpi_A6(const real& x)
00875 {
00876         real y,v;
00877         
00878         if (x == 0.43359375)
00879                 y = q_sin6_b[0];
00880         else
00881         {
00882                 v = 1/(x-0.43359375);
00883                 // Approximation:  sin(pi*x)/pi = K5(v);
00884                 y = q_sin6_a[5] / (  v + q_sin6_b[5]);
00885                 y = q_sin6_a[4] / ( (v + q_sin6_b[4]) + y);
00886                 y = q_sin6_a[3] / ( (v + q_sin6_b[3]) + y);
00887                 y = q_sin6_a[2] / ( (v + q_sin6_b[2]) + y);
00888                 y = q_sin6_a[1] / ( (v + q_sin6_b[1]) + y) + q_sin6_b[0];
00889         }
00890         
00891         return y;
00892 }
00893 
00894 real sinpix_pi(const real& x)
00895 {  // Sin(Pi*x)/Pi;  rel. Fehlerschranke: 4.870167e-16;
00896         real res(0),xr;
00897         int m,nr;
00898         bool neg;
00899 
00900         m = Round(x);
00901         if (m == -2147483648)
00902                 cxscthrow(STD_FKT_OUT_OF_DEF("real sinpix_pi(const real&)"));
00903         xr = x - m;  // xr: reduced argument, exactly calculated!
00904    // |xr| <= 0.5;
00905         neg = xr<0;
00906         if (neg) xr = -xr;
00907    // 0 <= xr <= 0.5;
00908         nr = int_no(a_sinpix_pi,8,xr);
00909         
00910         switch(nr)
00911         {
00912                 case 0: res = xr;           break;
00913                 case 1: res = sinpi_A1(xr); break;
00914                 case 2: res = sinpi_A2(xr); break;
00915                 case 3: res = sinpi_A3(xr); break;
00916                 case 4: res = sinpi_A4(xr); break;
00917                 case 5: res = sinpi_A5(xr); break;
00918                 case 6: res = sinpi_A6(xr); break;              
00919                 case 7: res = 5734161139222659.0 / 18014398509481984.0;
00920         }
00921         
00922         if (neg) res = -res;
00923         if (m%2 != 0) res = -res;
00924         
00925         return res;
00926 }
00927 
00928 // ********************************************************************************
00929 // ***************  Konstanten bez. Gamma(x) und 1/Gamma(x) ********************
00930 // ********************************************************************************
00931 
00932 // Interval bounds of the neighbour intervals with respect to the
00933 // intervals S_k for Gamma(x):
00934 real gam_f85[19] =
00935 {-0.5, 8.5, 16.5, 24.5, 35.5, 46.5, 57.5, 68.5, 79.5, 90.5,
00936  101.5, 112.5, 122.5, 132.5, 142.5, 150.5, 157.5, 164.5, 171.5 };
00937 
00938  // 1/Gamma(x): S0=[1.5,2.5], x0 = 2.0;
00939   const real q_gams0_a[8] = 
00940   {0.0, 
00941    -7616205496030159.0 / 18014398509481984.0,
00942    6808968414757234.0 / 9007199254740992.0,
00943   -7179656013820698.0 / 144115188075855872.0,
00944   7646522333236288.0 / 144115188075855872.0,
00945   -7301751514589296.0 / 144115188075855872.0,
00946   6062274112599236.0 / 288230376151711744.0,
00947   -7580146497393670.0 / 576460752303423488.0 };
00948          
00949   const real q_gams0_b[8] = 
00950   {1.0,
00951   -4965940208012024.0 / 9007199254740992.0,
00952   8627036189024519.0 / 9007199254740992.0,
00953   -7819065539096245.0 / 72057594037927936.0,
00954   7433219784613049.0 / 36028797018963968.0,
00955   7389378817674059.0 / 72057594037927936.0,
00956   -6060163955665826.0 / 144115188075855872.0,
00957   5769165265609039.0 / 18014398509481984.0 };
00958          
00959   // Gamma(x): S1=[10.5,11.5], x0 = 11.125;
00960   const real q_gams1_a[7] = 
00961   {0.0, 
00962   5262165366811426.0 / 72057594037927936.0,
00963   5574236285326272.0 / 9007199254740992.0,
00964   -8823729115230752.0 / 9223372036854775808.0,
00965   4639548867796070.0 / 72057594037927936.0,
00966   -4784760610310013.0 / 4611686018427387904.0,
00967   5900435828568799.0 / 288230376151711744.0 };
00968          
00969   const real q_gams1_b[7] = 
00970   {7108556377252296.0 / 36028797018963968.0,
00971   -7219014979973550.0 / 9007199254740992.0,
00972   7224885284258947.0 / 9007199254740992.0,
00973   -8569030245421939.0 / 36028797018963968.0,
00974   4838653997792226.0 / 18014398509481984.0,
00975   -4567015707194151.0 / 36028797018963968.0,
00976   5705105513288442.0 / 36028797018963968.0 };
00977          
00978   // Gamma(x): S2=[18.5,19.5], x0 = 18.96875;
00979   const real q_gams2_a[6] = 
00980   {0.0, 
00981   7322606177885925.0 / 72057594037927936.0,
00982   5856515731454725.0 / 144115188075855872.0,
00983   -5420981868254561.0 / 1152921504606846976.0,
00984   6209135328051357.0 / 1152921504606846976.0,
00985   -8261990134869242.0 / 2305843009213693952.0 };
00986          
00987   const real q_gams2_b[6] = 
00988   {-5267325780498998.0 / 18014398509481984.0,
00989   -4688643295042637.0 / 18014398509481984.0,
00990   6865435581764388.0 / 36028797018963968.0,
00991   -5147410677045000.0 / 144115188075855872.0,
00992   5878944809110092.0 / 144115188075855872.0,
00993   5315270486582075.0 / 576460752303423488.0 };
00994   
00995   // Gamma(x): S3=[29.5,30.5], x0 = 29.9375;
00996   const real q_gams3_a[6] = 
00997   {0.0, 
00998   4608894695736103.0 / 9007199254740992.0,
00999   4622056591672246.0 / 144115188075855872.0,
01000   7936321632956658.0 / 1152921504606846976.0,
01001   7533691379187725.0 / 2305843009213693952.0,
01002   8317076627520632.0 / 4611686018427387904.0 };
01003          
01004   const real q_gams3_b[6] = 
01005   {-5793090370845400.0 / 36028797018963968.0,
01006   -5993724737769479.0 / 18014398509481984.0,
01007   -7355055148590989.0 / 288230376151711744.0,
01008   -6811362569175008.0 / 288230376151711744.0,
01009   -5643865380181690.0 / 288230376151711744.0,
01010   -7257737463164532.0 / 288230376151711744.0 };
01011   
01012   // Gamma(x): S4=[40.5,41.5], x0 = 41.140625;
01013   const real q_gams4_a[6] = 
01014   {0.0, 
01015   -7504135817814391.0 / 18014398509481984.0,
01016   4990516833351222.0 / 576460752303423488.0,
01017   7409944094502278.0 / 4611686018427387904.0,
01018   8239437277348520.0 / 2305843009213693952.0,
01019   -5361398870955224.0 / 4611686018427387904.0 };
01020          
01021   const real q_gams4_b[6] = 
01022   {7405548010914168.0 / 18014398509481984.0,
01023   6819421126036941.0 / 36028797018963968.0,
01024   8474998846370093.0 / 288230376151711744.0,
01025   7733944242652532.0 / 144115188075855872.0,
01026   -6547033063052812.0 / 144115188075855872.0,
01027   -6636555235967309.0 / 576460752303423488.0 };
01028   
01029   // Gamma(x): S5=[51.5,52.5], x0 = 52.015625;
01030   const real q_gams5_a[6] = 
01031   {0.0, 
01032   -7282435522169731.0 / 144115188075855872.0,
01033   7812862759300002.0 / 288230376151711744.0,
01034   -7591789799457117.0 / 9223372036854775808.0,
01035   7290376444377792.0 / 2305843009213693952.0,
01036   -7051827298983325.0 / 9223372036854775808.0 };
01037          
01038   const real q_gams5_b[6] = 
01039   {-4692552597358020.0 / 36028797018963968.0,
01040   7065248316566075.0 / 36028797018963968.0,
01041   -5667667511562837.0 / 36028797018963968.0,
01042   7741409329322298.0 / 144115188075855872.0,
01043   -6371025742549657.0 / 144115188075855872.0,
01044   8045016330906311.0 / 288230376151711744.0 };
01045   
01046   // Gamma(x): S6=[62.5,63.5], x0 = 63.015625;
01047   const real q_gams6_a[6] = 
01048   {0.0, 
01049   6719861857699617.0 / 36028797018963968.0,
01050   6146108740657636.0 / 1152921504606846976.0,
01051   -8996491371873855.0 / 4611686018427387904.0,
01052   5082150623010284.0 / 2305843009213693952.0,
01053   -4650554434211955.0 / 18446744073709551616.0 };
01054          
01055   const real q_gams6_b[6] = 
01056   {6795471792542416.0 / 18014398509481984.0,
01057   -4567367130077106.0 / 36028797018963968.0,
01058   8403023343856889.0 / 288230376151711744.0,
01059   5083390457138636.0 / 144115188075855872.0,
01060   -4614092542204999.0 / 144115188075855872.0,
01061   7104477795873427.0 / 72057594037927936.0 };
01062   
01063   // Gamma(x): S7=[73.5,74.5], x0 = 74.16015625;
01064   const real q_gams7_a[6] = 
01065   {0.0 / 1.0,
01066   5372276754422009.0 / 18014398509481984.0,
01067   4663470139182266.0 / 576460752303423488.0,
01068   8173857739398457.0 / 4611686018427387904.0,
01069   4932343261664244.0 / 4611686018427387904.0,
01070   -5822870305854791.0 / 295147905179352825856.0 };
01071          
01072   const real q_gams7_b[6] = 
01073   {-4806355488125952.0 / 1152921504606846976.0,
01074   -6211400907258787.0 / 36028797018963968.0,
01075   -5429997192499743.0 / 288230376151711744.0,
01076   -5653057298064211.0 / 288230376151711744.0,
01077   -4746908442350821.0 / 1152921504606846976.0,
01078   8670733620601602.0 / 18014398509481984.0 };
01079   
01080   // Gamma(x): S8=[84.5,85.5], x0 = 85.1015625;
01081   const real q_gams8_a[5] = 
01082   { 0.0 / 1.0,
01083   -8754830070807807.0 / 72057594037927936.0,
01084   7931975738629914.0 / 2305843009213693952.0,
01085   6414800170661986.0 / 73786976294838206464.0,
01086   4558538137025832.0 / 18014398509481984.0};
01087          
01088   const real q_gams8_b[5] = 
01089   {-4924952929015874.0 / 18014398509481984.0,
01090   8571263173618757.0 / 72057594037927936.0,
01091   7912939311540494.0 / 576460752303423488.0,
01092   8975599623306617.0 / 18014398509481984.0,
01093   -4569152133381960.0 / 9007199254740992.0 };
01094   
01095   // Gamma(x): S9=[95.5,96.5], x0 = 95.984375;
01096   const real q_gams9_a[5] = 
01097   {0.0 / 1.0,
01098   -6140046099729716.0 / 144115188075855872.0,
01099   7278997066304035.0 / 576460752303423488.0,
01100   -4759432412103962.0 / 9223372036854775808.0,
01101   6912711986126027.0 / 4611686018427387904.0 };
01102          
01103   const real q_gams9_b[5] = 
01104   {-5611185402650416.0 / 72057594037927936.0,
01105   4915638659703209.0 / 36028797018963968.0,
01106   -7692469721238990.0 / 72057594037927936.0,
01107   5045754589592299.0 / 144115188075855872.0,
01108   -8290188163752846.0 / 288230376151711744.0 };
01109   
01110   // Gamma(x): S10=[106.5,107.5], x0 = 107.078125;
01111   const real q_gams10_a[5] = 
01112   {0.0 / 1.0,
01113   4717576521494070.0 / 72057594037927936.0,
01114   6906626643101502.0 / 1152921504606846976.0,
01115   -8147816895338065.0 / 9223372036854775808.0,
01116   7960901442281121.0 / 9223372036854775808.0 };
01117          
01118   const real q_gams10_b[5] = 
01119   {7952058244493952.0 / 288230376151711744.0,
01120   -7601355771801470.0 / 72057594037927936.0,
01121   4923655207606594.0 / 72057594037927936.0,
01122   -7238499638430702.0 / 576460752303423488.0,
01123   5852764550899526.0 / 576460752303423488.0 };
01124   
01125   // Gamma(x): S11=[117.5,118.5], x0 = 117.8671875;
01126   const real q_gams11_a[5] = 
01127   {0.0 / 1.0,
01128   5000157748740957.0 / 36028797018963968.0,
01129   6733781271808263.0 / 2305843009213693952.0,
01130   4912448110743899.0 / 18446744073709551616.0,
01131   6917678469221762.0 / 576460752303423488.0 };
01132          
01133   const real q_gams11_b[5] = 
01134   {-4805173503400964.0 / 36028797018963968.0,
01135   -7686561799456867.0 / 72057594037927936.0,
01136   -6649023763933472.0 / 576460752303423488.0,
01137   -7324214555543306.0 / 72057594037927936.0,
01138   8284786467509721.0 / 72057594037927936.0 };
01139   
01140   // Gamma(x): Stuetzinterval S12=[126.5,127.5], x0 = 126.7421875;
01141   const real q_gams12_a[5] = 
01142   {0.0 / 1.0,
01143   5226845396890227.0 / 36028797018963968.0,
01144   5602232597459763.0 / 1152921504606846976.0,
01145   4920783552645014.0 / 4611686018427387904.0,
01146   5570041179715558.0 / 9223372036854775808.0 };
01147          
01148   const real q_gams12_b[5] = 
01149   {-6799658092196962.0 / 18014398509481984.0,
01150   -4810318281953418.0 / 36028797018963968.0,
01151   -8340064868499919.0 / 576460752303423488.0,
01152   -8375346434191481.0 / 576460752303423488.0,
01153   -6231299816839781.0 / 1152921504606846976.0 };
01154   
01155   // Gamma(x): S13=[137.5,138.5], x0 = 138.0390625;
01156   const real q_gams13_a[5] = 
01157   {0.0 / 1.0,
01158   5077602289066213.0 / 18014398509481984.0,
01159   4971389438544026.0 / 576460752303423488.0,
01160   8326326060413143.0 / 4611686018427387904.0,
01161   7582214946963028.0 / 9223372036854775808.0 };
01162          
01163   const real q_gams13_b[5] = 
01164   {-8336661118182000.0 / 72057594037927936.0,
01165   -6152827173580884.0 / 36028797018963968.0,
01166   -6306415240260295.0 / 576460752303423488.0,
01167   -5949724688158565.0 / 576460752303423488.0,
01168   -5625482009313801.0 / 576460752303423488.0 };
01169   
01170   // Gamma(x): S14=[146.5,147.5], x0 = 146.94921875;
01171   const real q_gams14_a[6] = 
01172   {0.0 / 1.0,
01173   8624518348171320.0 / 18014398509481984.0,
01174   7049908374954842.0 / 576460752303423488.0,
01175   5769548496333642.0 / 2305843009213693952.0,
01176   5110826386599631.0 / 4611686018427387904.0,
01177   5902099580937297.0 / 9223372036854775808.0 };
01178          
01179   const real q_gams14_b[6] = 
01180   {4591842870321392.0 / 18014398509481984.0,
01181   -7195103997153274.0 / 36028797018963968.0,
01182   -5062133254875404.0 / 576460752303423488.0,
01183   -4899776496053167.0 / 576460752303423488.0,
01184   -4703735605739871.0 / 576460752303423488.0,
01185   -8717622510435264.0 / 1152921504606846976.0 };
01186   
01187   // Gamma(x): Stuetzinterval S15=[153.5,154.5], x0 = 153.90234375;
01188   const real q_gams15_a[6] = 
01189   {0.0 / 1.0,
01190   5047093325633351.0 / 9007199254740992.0,
01191   8838568428828895.0 / 576460752303423488.0,
01192   7169875770591066.0 / 2305843009213693952.0,
01193   6277116142443993.0 / 4611686018427387904.0,
01194   7160380168510556.0 / 9223372036854775808.0 };
01195          
01196   const real q_gams15_b[6] = 
01197   {5575895624898616.0 / 18014398509481984.0,
01198   -7982725134478240.0 / 36028797018963968.0,
01199   -8683131154525203.0 / 1152921504606846976.0,
01200   -8504376590080638.0 / 1152921504606846976.0,
01201   -8272396567408614.0 / 1152921504606846976.0,
01202   -7370684458128734.0 / 1152921504606846976.0 }; 
01203    
01204   // Gamma(x): S16=[160.5,161.5], x0 = 161.08984375;
01205   const real q_gams16_a[6] = 
01206   {0.0 / 1.0,
01207   8928052213955455.0 / 18014398509481984.0,
01208   5405751921217612.0 / 288230376151711744.0,
01209   8725861134087760.0 / 2305843009213693952.0,
01210   7583718370137025.0 / 4611686018427387904.0,
01211   8577302812898546.0 / 9223372036854775808.0 };
01212          
01213   const real q_gams16_b[6] = 
01214   {6669445708872192.0 / 144115188075855872.0,
01215   -8769965967955460.0 / 36028797018963968.0,
01216   -7524042392635917.0 / 1152921504606846976.0,
01217   -7422547685779587.0 / 1152921504606846976.0,
01218   -7284571995868819.0 / 1152921504606846976.0,
01219   -7797018373646746.0 / 1152921504606846976.0 }; 
01220   
01221   // Gamma(x): S17=[167.5,168.5], x0 = 168.0;
01222   const real q_gams17_a[6] = 
01223   {0.0 / 1.0,
01224   4642907317603488.0 / 9007199254740992.0,
01225   6403634366308572.0 / 288230376151711744.0,
01226   5153510478021874.0 / 1152921504606846976.0,
01227   8919203526515161.0 / 4611686018427387904.0,
01228   5015558226948384.0 / 4611686018427387904.0 };
01229          
01230   const real q_gams17_b[6] = 
01231   {-6229620043098112.0 / 9223372036854775808.0,
01232   -4750296278401558.0 / 18014398509481984.0,
01233   -6639785613322219.0 / 1152921504606846976.0,
01234   -6577910457226093.0 / 1152921504606846976.0,
01235   -6491061443859352.0 / 1152921504606846976.0,
01236   -6372399542597081.0 / 1152921504606846976.0 };                                                                         
01237 
01238 // ****************************************************************************
01239 // ******************** Gamma(x), 1/Gamma(x) **********************************
01240 // ****************************************************************************
01241 
01242  real gam_S0(const real& x)
01243 // Calculating approximations for 1/Gamma(x) in S0 = [1.5,2.5];
01244 // Rel. error bound:  3.930138e-16;
01245 // Blomquist, 19.05.2009;
01246 {
01247         real y(0),v;
01248 
01249          // Continued fraction:  K_7(v), v = 1/(x-x0), x0=2;
01250         if (x==2)
01251                 y = q_gams0_b[0];
01252         else
01253         {
01254                 v = 1/(x-2);
01255                 
01256                 y = q_gams0_a[7] / (  v + q_gams0_b[7]);
01257                 y = q_gams0_a[6] / ( (v + q_gams0_b[6]) + y );
01258                 y = q_gams0_a[5] / ( (v + q_gams0_b[5]) + y );  
01259                 y = q_gams0_a[4] / ( (v + q_gams0_b[4]) + y );
01260                 y = q_gams0_a[3] / ( (v + q_gams0_b[3]) + y );
01261                 y = q_gams0_a[2] / ( (v + q_gams0_b[2]) + y );
01262                 y = q_gams0_a[1] / ( (v + q_gams0_b[1]) + y ) + q_gams0_b[0];
01263         }       
01264         return y;
01265 }
01266  
01267  real gam_S0_n0(const real& x)
01268 // Calculating approximations for K_7(v),  v=1/x,  x in [-0.5,+0.5];
01269 // Rel. error bound:  3.930138e-16;
01270 // Blomquist, 19.05.2009;
01271 {
01272         real y(0),v;
01273 
01274          // Continued fraction:  K_7(v), v = 1/x, x0=0;
01275         if (x==0)
01276                 y = q_gams0_b[0];
01277         else
01278         {
01279                 v = 1/x;
01280                  
01281                 y = q_gams0_a[7] / (  v + q_gams0_b[7]);
01282                 y = q_gams0_a[6] / ( (v + q_gams0_b[6]) + y );
01283                 y = q_gams0_a[5] / ( (v + q_gams0_b[5]) + y );  
01284                 y = q_gams0_a[4] / ( (v + q_gams0_b[4]) + y );
01285                 y = q_gams0_a[3] / ( (v + q_gams0_b[3]) + y );
01286                 y = q_gams0_a[2] / ( (v + q_gams0_b[2]) + y );
01287                 y = q_gams0_a[1] / ( (v + q_gams0_b[1]) + y ) + q_gams0_b[0];
01288         }       
01289                 
01290         return y;
01291 }
01292 
01293 real gam_S0_n1(const real& x)
01294 // Calculating approximations for K_7(v),  v=1/(x-1),  x in [0.5,1.5];
01295 // Rel. error bound:  3.930138e-16;
01296 // Blomquist, 21.05.2009;
01297 {
01298         real y(0),v;
01299 
01300          // Continued fraction:  K_7(v), v = 1/(x-1), x0=1;
01301         if (x==1)
01302                 y = q_gams0_b[0];
01303         else
01304         {
01305                 v = 1/(x-1);
01306                  
01307                 y = q_gams0_a[7] / (  v + q_gams0_b[7]);
01308                 y = q_gams0_a[6] / ( (v + q_gams0_b[6]) + y );
01309                 y = q_gams0_a[5] / ( (v + q_gams0_b[5]) + y );  
01310                 y = q_gams0_a[4] / ( (v + q_gams0_b[4]) + y );
01311                 y = q_gams0_a[3] / ( (v + q_gams0_b[3]) + y );
01312                 y = q_gams0_a[2] / ( (v + q_gams0_b[2]) + y );
01313                 y = q_gams0_a[1] / ( (v + q_gams0_b[1]) + y ) + q_gams0_b[0];
01314         }       
01315                 
01316         return y;
01317 }
01318 
01319 real gammar_S0(const real& x)
01320 // Calculating approximations for 1/Gamma(x), x in (-0.5,+8.5];
01321 // eps = 1.725284e-15;
01322 // Blomquist, 21.05.2009;
01323 {
01324         real y,Ne;
01325         int n,p;
01326         
01327         n = Round(x);
01328         
01329         switch(n)
01330         {
01331                 case 0: if (expo(x)<=-52) y = x;
01332                 else y = x*(x+1)*gam_S0_n0(x); break; // x in (-0.5,+0.5);
01333                 case 1: y = x*gam_S0_n1(x);            break; // x in [0.5,1.5);
01334                 case 2: y = gam_S0(x);                 break; // x in [1.5,2.5);
01335 
01336                 default: // n=3,4,...,8; x in [2.5,8.5];
01337                         p = n-2; Ne = x-1;
01338                         for (int k=2; k<=p; k++) Ne *= (x-k);
01339                         y = gam_S0(x-p) / Ne;
01340         }
01341         
01342         return y;
01343 }
01344 
01345 real gamma_S0(const real& x)
01346 // Calculating approximations for Gamma(x) in S0 = (-0.5,+8.5];
01347 // eps = 1.947330E-15;
01348 // Blomquist, 21.05.2009;
01349 {
01350         real y = 1/gammar_S0(x);
01351         return y;
01352 }
01353 
01354 real gam_S1(const real& x)
01355 // Calculating approximations for Gamma(x) in  S1 = [10.5,11.5];
01356 // Rel. error bound:  7.696082e-16;
01357 // Blomquist, 22.05.2009;               
01358 {
01359         real y(0),v;
01360 
01361         // Continued fraction:  K_6(v), v = 1/(x-x0), x0=11.125;
01362         if (x==11.125)
01363                 y = q_gams1_b[0];
01364         else
01365         {
01366                 v = 1/(x-11.125);
01367                 
01368                 y = q_gams1_a[6] / (  v + q_gams1_b[6]);
01369                 y = q_gams1_a[5] / ( (v + q_gams1_b[5]) + y );  
01370                 y = q_gams1_a[4] / ( (v + q_gams1_b[4]) + y );
01371                 y = q_gams1_a[3] / ( (v + q_gams1_b[3]) + y );
01372                 y = q_gams1_a[2] / ( (v + q_gams1_b[2]) + y );
01373                 y = q_gams1_a[1] / ( (v + q_gams1_b[1]) + y ) + q_gams1_b[0];
01374         }       
01375         y += 1;
01376         y *= q_ex10(x);
01377         times2pown(y,-15);
01378                 
01379         return y;
01380 }       
01381         
01382 real gamma_S1(const real& x)
01383 // Calculating approximations for Gamma(x), x in [8.5,16.5];
01384 // eps = 1.879833e-15;
01385 // Blomquist, 24.05.2009;
01386 {
01387         real y,Pr;
01388         int n,p;
01389         
01390         n = Round(x);
01391         p = n-11;
01392 
01393         if (n>11) // neighbour intervals to the right
01394         {
01395                 Pr = x-1;
01396                 for (int k=2; k<=p; k++)
01397                         Pr *= x-k;
01398                 y = Pr*gam_S1(x-p);
01399         }
01400         else // neighbour intervals to the left and S1 itself
01401         {
01402                 p = -p;  Pr = x;
01403                 for (int k=1; k<=p-1; k++) 
01404                         Pr *= x+k;
01405                 y = (p==0)? gam_S1(x) : gam_S1(x+p)/Pr;         
01406         }
01407         
01408         return y;
01409 }       
01410 
01411 real gam_S2(const real& x)
01412 // Calculating approximations for Gamma(x) in S2 = [18.5,19.5];
01413 // Rel. error bound:  7.411454e-16;
01414 // Blomquist, 25.05.2009;               
01415 {
01416         real y(0),v;
01417 
01418         // Continued fraction:  K_5(v), v = 1/(x-x0), x0=18.96875;
01419         if (x==18.96875)
01420                 y = q_gams2_b[0];
01421         else
01422         {
01423                 v = 1/(x-18.96875);
01424 
01425                 y = q_gams2_a[5] / (  v + q_gams2_b[5]);
01426                 y = q_gams2_a[4] / ( (v + q_gams2_b[4]) + y );
01427                 y = q_gams2_a[3] / ( (v + q_gams2_b[3]) + y );
01428                 y = q_gams2_a[2] / ( (v + q_gams2_b[2]) + y );
01429                 y = q_gams2_a[1] / ( (v + q_gams2_b[1]) + y ) + q_gams2_b[0];
01430         }       
01431         y += 1;
01432         y *= q_exp2(4*x);
01433         times2pown(y,-23);
01434                 
01435         return y;
01436 }       
01437 
01438 real gamma_S2(const real& x)
01439 // Calculating approximations for Gamma(x), x in [16.5,24.5];
01440 // eps = 1.851370e-15;
01441 // Blomquist, 25.05.2009;
01442 {
01443         real y,Pr;
01444         int n,p;
01445         
01446         n = Round(x);
01447         p = n-19;
01448 
01449         if (n>19) // neighbour intervals to the right
01450         {
01451                 Pr = x-1;
01452                 for (int k=2; k<=p; k++)
01453                         Pr *= x-k;
01454                 y = Pr*gam_S2(x-p);
01455         }
01456         else // neighbour intervals to the left and S2 itself
01457         {
01458                 p = -p;  Pr = x;
01459                 for (int k=1; k<=p-1; k++) 
01460                         Pr *= x+k;
01461                 y = (p==0)? gam_S2(x) : gam_S2(x+p)/Pr;         
01462         }
01463         
01464         return y;
01465 }       
01466 
01467 real gam_S3(const real& x)
01468 // Calculating approximations for Gamma(x) in  S3 = [29.5,30.5];
01469 // Rel. error bound:  9.724029e-16;
01470 // Blomquist, 26.05.2009;               
01471 {
01472         real y(0),v;
01473 
01474         // Continued fraction:  K_5(v), v = 1/(x-x0), x0=29.9375;
01475         if (x==29.9375)
01476                 y = q_gams3_b[0];
01477         else
01478         {
01479                 v = 1/(x-29.9375);
01480 
01481                 y = q_gams3_a[5] / (  v + q_gams3_b[5]);
01482                 y = q_gams3_a[4] / ( (v + q_gams3_b[4]) + y );
01483                 y = q_gams3_a[3] / ( (v + q_gams3_b[3]) + y );
01484                 y = q_gams3_a[2] / ( (v + q_gams3_b[2]) + y );
01485                 y = q_gams3_a[1] / ( (v + q_gams3_b[1]) + y ) + q_gams3_b[0];
01486         }       
01487         y += 1;
01488         y *= q_exp2(4*x);
01489         times2pown(y,-17);
01490                 
01491         return y;
01492 }       
01493 
01494 real gamma_S3(const real& x)
01495 // Calculating approximations for Gamma(x), x in [24.5,35.5];
01496 // eps = 2.082627e-15;
01497 // Blomquist, 26.05.2009;
01498 {
01499         real y,Pr;
01500         int n,p;
01501         
01502         n = Round(x);
01503         p = n-30;
01504 
01505         if (n>30) // neighbour intervals to the right
01506         {
01507                 Pr = x-1;
01508                 for (int k=2; k<=p; k++)
01509                         Pr *= x-k;
01510                 y = Pr*gam_S3(x-p);
01511         }
01512         else // neighbour intervals to the left and S3 itself
01513         {
01514                 p = -p;  Pr = x;
01515                 for (int k=1; k<=p-1; k++) 
01516                         Pr *= x+k;
01517                 y = (p==0)? gam_S3(x) : gam_S3(x+p)/Pr;         
01518         }
01519         
01520         return y;
01521 }       
01522 
01523 real gam_S4(const real& x)
01524 // Calculating approximations for Gamma(x) in  S4 = [40.5,41.5];
01525 // Rel. error bound:  8.170628e-16;
01526 // Blomquist, 26.05.2009;               
01527 {
01528         real y(0),v;
01529 
01530         // Continued fraction:  K_5(v), v = 1/(x-x0), x0=41.140625;
01531         if (x==41.140625)
01532                 y = q_gams4_b[0];
01533         else
01534         {
01535                 v = 1/(x-41.140625);
01536 
01537                 y = q_gams4_a[5] / (  v + q_gams4_b[5]);
01538                 y = q_gams4_a[4] / ( (v + q_gams4_b[4]) + y );
01539                 y = q_gams4_a[3] / ( (v + q_gams4_b[3]) + y );
01540                 y = q_gams4_a[2] / ( (v + q_gams4_b[2]) + y );
01541                 y = q_gams4_a[1] / ( (v + q_gams4_b[1]) + y ) + q_gams4_b[0];
01542         }       
01543         y += 1;
01544         y *= exp(4*x);
01545         times2pown(y,-78);
01546                 
01547         return y;
01548 }       
01549 
01550 real gamma_S4(const real& x)
01551 // Calculating approximations for Gamma(x), x in [35.5,46.5];
01552 // eps = 1.927286e-15;
01553 // Blomquist, 28.05.2009;
01554 {
01555         real y,Pr;
01556         int n,p;
01557         
01558         n = Round(x);
01559         p = n-41;
01560 
01561         if (n>41) // neighbour intervals to the right
01562         {
01563                 Pr = x-1;
01564                 for (int k=2; k<=p; k++)
01565                         Pr *= x-k;
01566                 y = Pr*gam_S4(x-p);
01567         }
01568         else // neighbour intervals to the left and S4 itself
01569         {
01570                 p = -p;  Pr = x;
01571                 for (int k=1; k<=p-1; k++) 
01572                         Pr *= x+k;
01573                 y = (p==0)? gam_S4(x) : gam_S4(x+p)/Pr;         
01574         }
01575         
01576         return y;
01577 }       
01578 
01579 real gam_S5(const real& x)
01580 // Calculating approximations for Gamma(x) in S5 = [51.5,52.5];
01581 // Rel. error bound:  6.351369e-16;
01582 // Blomquist, 28.05.2009;               
01583 {
01584         real y(0),v;
01585 
01586         // Continued fraction:  K_5(v), v = 1/(x-x0), x0=52.015625;
01587         if (x==52.015625)
01588                 y = q_gams5_b[0];
01589         else
01590         {
01591                 v = 1/(x-52.015625);
01592 
01593                 y = q_gams5_a[5] / (  v + q_gams5_b[5]);
01594                 y = q_gams5_a[4] / ( (v + q_gams5_b[4]) + y );
01595                 y = q_gams5_a[3] / ( (v + q_gams5_b[3]) + y );
01596                 y = q_gams5_a[2] / ( (v + q_gams5_b[2]) + y );
01597                 y = q_gams5_a[1] / ( (v + q_gams5_b[1]) + y ) + q_gams5_b[0];
01598         }       
01599         y += 1;
01600         y *= exp(4*x);
01601         times2pown(y,-80);
01602                 
01603         return y;
01604 }       
01605 
01606 real gamma_S5(const real& x)
01607 // Calculating approximations for Gamma(x), x in [46.5,57.5];
01608 // eps = 1.745360e-15;
01609 // Blomquist, 29.05.2009;
01610 {
01611         real y,Pr;
01612         int n,p;
01613         
01614         n = Round(x);
01615         p = n-52;
01616 
01617         if (n>52) // neighbour intervals to the right
01618         {
01619                 Pr = x-1;
01620                 for (int k=2; k<=p; k++)
01621                         Pr *= x-k;
01622                 y = Pr*gam_S5(x-p);
01623         }
01624         else // neighbour intervals to the left and S5 itself
01625         {
01626                 p = -p;  Pr = x;
01627                 for (int k=1; k<=p-1; k++) 
01628                         Pr *= x+k;
01629                 y = (p==0)? gam_S5(x) : gam_S5(x+p)/Pr;         
01630         }
01631         
01632         return y;
01633 }       
01634 
01635 real gam_S6(const real& x)
01636 // Calculating approximations for Gamma(x) in S6 = [62.5,63.5];
01637 // Rel. error bound:  8.310104e-16;
01638 // Blomquist, 29.05.2009;               
01639 {
01640         real y(0),v;
01641 
01642         // Continued fraction:  K_5(v), v = 1/(x-x0), x0=63.015625;
01643         if (x==63.015625)
01644                 y = q_gams6_b[0];
01645         else
01646         {
01647                 v = 1/(x-63.015625);
01648 
01649                 y = q_gams6_a[5] / (  v + q_gams6_b[5]);
01650                 y = q_gams6_a[4] / ( (v + q_gams6_b[4]) + y );
01651                 y = q_gams6_a[3] / ( (v + q_gams6_b[3]) + y );
01652                 y = q_gams6_a[2] / ( (v + q_gams6_b[2]) + y );
01653                 y = q_gams6_a[1] / ( (v + q_gams6_b[1]) + y ) + q_gams6_b[0];
01654         }       
01655         y += 1;
01656         y *= exp(4*x);
01657         times2pown(y,-80);
01658                 
01659         return y;
01660 }       
01661 
01662 real gamma_S6(const real& x)
01663 // Calculating approximations for Gamma(x), x in [57.5,68.5];
01664 // eps = 1.941234e-15;
01665 // Blomquist, 06.06.2009;
01666 {
01667         real y,Pr;
01668         int n,p;
01669         
01670         n = Round(x);
01671         p = n-63;
01672 
01673         if (n>63) // neighbour intervals to the right
01674         {
01675                 Pr = x-1;
01676                 for (int k=2; k<=p; k++)
01677                         Pr *= x-k;
01678                 y = Pr*gam_S6(x-p);
01679         }
01680         else // neighbour intervals to the left and S6 itself
01681         {
01682                 p = -p;  Pr = x;
01683                 for (int k=1; k<=p-1; k++) 
01684                         Pr *= x+k;
01685                 y = (p==0)? gam_S6(x) : gam_S6(x+p)/Pr;         
01686         }
01687         
01688         return y;
01689 }       
01690 
01691 
01692 real gam_S7(const real& x)
01693 // Calculating approximations for Gamma(x) in S7 = [73.5,74.5];
01694 // Rel. error bound:  7.999142e-16;
01695 // Blomquist, 05.06.2009;               
01696 {
01697         real y(0),v;
01698 
01699         // Continued fraction:  K_5(v), v = 1/(x-x0), x0=74.16015625;
01700         if (x==74.16015625)
01701                 y = q_gams7_b[0];
01702         else
01703         {
01704                 v = 1/(x-74.16015625);
01705 
01706                 y = q_gams7_a[5] / (  v + q_gams7_b[5]);
01707                 y = q_gams7_a[4] / ( (v + q_gams7_b[4]) + y );
01708                 y = q_gams7_a[3] / ( (v + q_gams7_b[3]) + y );
01709                 y = q_gams7_a[2] / ( (v + q_gams7_b[2]) + y );
01710                 y = q_gams7_a[1] / ( (v + q_gams7_b[1]) + y ) + q_gams7_b[0];
01711         }       
01712         y += 1;
01713         y *= exp(4*x);
01714         times2pown(y,-76);
01715                 
01716         return y;
01717 }       
01718 
01719 real gamma_S7(const real& x)
01720 // Calculating approximations for Gamma(x), x in [68.5,79.5];
01721 // eps = 1.910138e-15;
01722 // Blomquist, 08.06.2009;
01723 {
01724         real y,Pr;
01725         int n,p;
01726         
01727         n = Round(x);
01728         p = n-74;
01729 
01730         if (n>74) // neighbour intervals to the right
01731         {
01732                 Pr = x-1;
01733                 for (int k=2; k<=p; k++)
01734                         Pr *= x-k;
01735                 y = Pr*gam_S7(x-p);
01736         }
01737         else // neighbour intervals to the left and S7 itself 
01738         {
01739                 p = -p;  Pr = x;
01740                 for (int k=1; k<=p-1; k++) 
01741                         Pr *= x+k;
01742                 y = (p==0)? gam_S7(x) : gam_S7(x+p)/Pr;         
01743         }
01744         
01745         return y;
01746 }       
01747 
01748 real gam_S8(const real& x)
01749 // Calculating approximations for Gamma(x) in  S8 = [84.5,85.5];
01750 // Rel. error bound:  7.192334e-16;
01751 // Blomquist, 08.06.2009;               
01752 {
01753         real y(0),v;
01754 
01755         // Continued fraction:  K_4(v), v = 1/(x-x0), x0=85.1015625;
01756         if (x==85.1015625)
01757                 y = q_gams8_b[0];
01758         else
01759         {
01760                 v = 1/(x-85.1015625);
01761 
01762                 y = q_gams8_a[4] / (  v + q_gams8_b[4]);
01763                 y = q_gams8_a[3] / ( (v + q_gams8_b[3]) + y );
01764                 y = q_gams8_a[2] / ( (v + q_gams8_b[2]) + y );
01765                 y = q_gams8_a[1] / ( (v + q_gams8_b[1]) + y ) + q_gams8_b[0];
01766         }       
01767         y += 1;
01768         y *= q_ex10(2*x);
01769         times2pown(y,-144);
01770                 
01771         return y;
01772 }       
01773 
01774 real gamma_S8(const real& x)
01775 // Calculating approximations for Gamma(x), x in [79.5,90.5];
01776 // eps = 1.829457e-15;
01777 // Blomquist, 09.06.2009;
01778 {
01779         real y,Pr;
01780         int n,p;
01781         
01782         n = Round(x);
01783         p = n-85;
01784 
01785         if (n>85) // neighbour intervals to the right
01786         {
01787                 Pr = x-1;
01788                 for (int k=2; k<=p; k++)
01789                         Pr *= x-k;
01790                 y = Pr*gam_S8(x-p);
01791         }
01792         else // neighbour intervals to the left and S8 itself
01793         {
01794                 p = -p;  Pr = x;
01795                 for (int k=1; k<=p-1; k++) 
01796                         Pr *= x+k;
01797                 y = (p==0)? gam_S8(x) : gam_S8(x+p)/Pr;         
01798         }
01799         
01800         return y;
01801 }       
01802 
01803 
01804 real gam_S9(const real& x)
01805 // Calculating approximations for Gamma(x) in S9 = [95.5,96.5];
01806 // Rel. error bound:  6.622145e-16;
01807 // Blomquist, 09.06.2009;               
01808 {
01809         real y(0),v;
01810 
01811         // Continued fraction:  K_4(v), v = 1/(x-x0), x0=95.984375;
01812         if (x==95.984375)
01813                 y = q_gams9_b[0];
01814         else
01815         {
01816                 v = 1/(x-95.984375);
01817 
01818                 y = q_gams9_a[4] / (  v + q_gams9_b[4]);
01819                 y = q_gams9_a[3] / ( (v + q_gams9_b[3]) + y );
01820                 y = q_gams9_a[2] / ( (v + q_gams9_b[2]) + y );
01821                 y = q_gams9_a[1] / ( (v + q_gams9_b[1]) + y ) + q_gams9_b[0];
01822         }       
01823         y += 1;
01824         y *= q_ex10(2*x);
01825         times2pown(y,-146);
01826                 
01827         return y;
01828 }       
01829 
01830 real gamma_S9(const real& x)
01831 // Calculating approximations for Gamma(x), x in [90.5,101.5];
01832 // eps = 1.772438e-15;
01833 // Blomquist, 10.06.2009;
01834 {
01835         real y,Pr;
01836         int n,p;
01837         
01838         n = Round(x);
01839         p = n-96;
01840 
01841         if (n>96) // neighbour intervals to the right
01842         {
01843                 Pr = x-1;
01844                 for (int k=2; k<=p; k++)
01845                         Pr *= x-k;
01846                 y = Pr*gam_S9(x-p);
01847         }
01848         else // neighbour intervals to the left and S9 itself
01849         {
01850                 p = -p;  Pr = x;
01851                 for (int k=1; k<=p-1; k++) 
01852                         Pr *= x+k;
01853                 y = (p==0)? gam_S9(x) : gam_S9(x+p)/Pr;         
01854         }
01855         
01856         return y;
01857 }       
01858 
01859 
01860 real gam_S10(const real& x)
01861 // Calculating approximations for Gamma(x) in S10 = [106.5,107.5];
01862 // Rel. error bound:  8.254965e-16;
01863 // Blomquist, 10.06.2009;               
01864 {
01865         real y(0),v;
01866 
01867         // Continued fraction:  K_4(v), v = 1/(x-x0), x0=107.078125;
01868         if (x==107.078125)
01869                 y = q_gams10_b[0];
01870         else
01871         {
01872                 v = 1/(x-107.078125);
01873 
01874                 y = q_gams10_a[4] / (  v + q_gams10_b[4]);
01875                 y = q_gams10_a[3] / ( (v + q_gams10_b[3]) + y );
01876                 y = q_gams10_a[2] / ( (v + q_gams10_b[2]) + y );
01877                 y = q_gams10_a[1] / ( (v + q_gams10_b[1]) + y ) + q_gams10_b[0];
01878         }       
01879         y += 1;
01880         y *= q_ex10(2*x);
01881         times2pown(y,-146);
01882                 
01883         return y;
01884 }       
01885 
01886 real gamma_S10(const real& x)
01887 // Calculating approximations for Gamma(x), x in [101.5,112.5];
01888 // eps = 1.935720e-15;
01889 // Blomquist, 12.06.2009;
01890 {
01891         real y,Pr;
01892         int n,p;
01893         
01894         n = Round(x);
01895         p = n-107;
01896 
01897         if (n>107) // neighbour intervals to the right
01898         {
01899                 Pr = x-1;
01900                 for (int k=2; k<=p; k++)
01901                         Pr *= x-k;
01902                 y = Pr*gam_S10(x-p);
01903         }
01904         else // neighbour intervals to the left and S10 itself
01905         {
01906                 p = -p;  Pr = x;
01907                 for (int k=1; k<=p-1; k++) 
01908                         Pr *= x+k;
01909                 y = (p==0)? gam_S10(x) : gam_S10(x+p)/Pr;               
01910         }
01911         
01912         return y;
01913 }       
01914 
01915 real gam_S11(const real& x)
01916 // Calculating approximations for Gamma(x) in S11 = [117.5,118.5];
01917 // Rel. error bound:  6.766734e-16;
01918 // Blomquist, 11.06.2009;               
01919 {
01920         real y(0),v;
01921 
01922         // Continued fraction:  K_4(v), v = 1/(x-x0), x0=117.8671875;
01923         if (x==117.8671875)
01924                 y = q_gams11_b[0];
01925         else
01926         {
01927                 v = 1/(x-117.8671875);
01928 
01929                 y = q_gams11_a[4] / (  v + q_gams11_b[4]);
01930                 y = q_gams11_a[3] / ( (v + q_gams11_b[3]) + y );
01931                 y = q_gams11_a[2] / ( (v + q_gams11_b[2]) + y );
01932                 y = q_gams11_a[1] / ( (v + q_gams11_b[1]) + y ) + q_gams11_b[0];
01933         }       
01934         y += 1;
01935         y *= q_ex10(2*x);
01936         times2pown(y,-144);
01937                 
01938         return y;
01939 }       
01940 
01941 real gamma_S11(const real& x)
01942 // Calculating approximations for Gamma(x) in [112.5,122.5];
01943 // eps = 1.786897e-15;
01944 // Blomquist, 14.06.2009;
01945 {
01946         real y,Pr;
01947         int n,p;
01948         
01949         n = Round(x);
01950         p = n-118;
01951         
01952         if (n>118) // neighbour intervals to the right
01953         {
01954                 Pr = x-1;
01955                 for (int k=2; k<=p; k++)
01956                         Pr *= x-k;
01957                 y = Pr*gam_S11(x-p);
01958         }
01959         else // left neighbour intervals and S11 itself
01960         {
01961                 p = -p;  Pr = x;
01962                 for (int k=1; k<=p-1; k++) 
01963                         Pr *= x+k;
01964                 y = (p==0)? gam_S11(x) : gam_S11(x+p)/Pr;               
01965         }
01966         
01967         return y;
01968 }       
01969 
01970 real gam_S12(const real& x)
01971 // Calculating approximations for Gamma(x) in S12 = [126.5,127.5];
01972 // Rel. error bound:  8.175777e-16;
01973 // Blomquist, 13.06.2009;               
01974 {
01975         real y(0),v;
01976 
01977         // Continued fraction:  K_4(v), v = 1/(x-x0), x0=126.7421875;
01978         if (x==126.7421875)
01979                 y = q_gams12_b[0];
01980         else
01981         {
01982                 v = 1/(x-126.7421875);
01983 
01984                 y = q_gams12_a[4] / (  v + q_gams12_b[4]);
01985                 y = q_gams12_a[3] / ( (v + q_gams12_b[3]) + y );
01986                 y = q_gams12_a[2] / ( (v + q_gams12_b[2]) + y );
01987                 y = q_gams12_a[1] / ( (v + q_gams12_b[1]) + y ) + q_gams12_b[0];
01988         }       
01989         y += 1;
01990         y *= q_ex10(2*x);
01991         times2pown(y,-141);
01992                 
01993         return y;
01994 }       
01995 
01996 real gamma_S12(const real& x)
01997 // Calculating approximations for Gamma(x) in [122.5,132.5];
01998 // eps = 1.927801e-15;
01999 // Blomquist, 15.06.2009;
02000 {
02001         real y,Pr;
02002         int n,p;
02003         
02004         n = Round(x);
02005         p = n-127;
02006         if (n>127) // neighbour intervals to the right
02007         {
02008                 Pr = x-1;
02009                 for (int k=2; k<=p; k++)
02010                         Pr *= x-k;
02011                 y = Pr*gam_S12(x-p);
02012         }
02013         else // left neighbour intervals and S12 itself
02014         {
02015                 p = -p;  Pr = x;
02016                 for (int k=1; k<=p-1; k++) 
02017                         Pr *= x+k;
02018                 y = (p==0)? gam_S12(x) : gam_S12(x+p)/Pr;               
02019         }
02020         
02021         return y;
02022 }       
02023 
02024 real gam_S13(const real& x)
02025 // Calculating approximations for Gamma(x) in S13 = [137.5,138.5];
02026 // Rel. error bound:  8.563188e-16;
02027 // Blomquist, 14.06.2009;               
02028 {
02029         real y(0),v;
02030 
02031         // Continued fraction:  K_4(v), v = 1/(x-x0), x0=138.0390625;
02032         if (x==138.0390625)
02033                 y = q_gams13_b[0];
02034         else
02035         {
02036                 v = 1/(x-138.0390625);
02037 
02038                 y = q_gams13_a[4] / (  v + q_gams13_b[4]);
02039                 y = q_gams13_a[3] / ( (v + q_gams13_b[3]) + y );
02040                 y = q_gams13_a[2] / ( (v + q_gams13_b[2]) + y );
02041                 y = q_gams13_a[1] / ( (v + q_gams13_b[1]) + y ) + q_gams13_b[0];
02042         }       
02043         y += 1;
02044         y *= q_ex10(2*x);
02045         times2pown(y,-137);
02046                 
02047         return y;
02048 }       
02049 
02050 real gamma_S13(const real& x)
02051 // Calculating approximations for Gamma(x) in [132.5,142.5];
02052 // eps = 1.966542e-15;
02053 // Blomquist, 15.06.2009;
02054 {
02055         real y,Pr;
02056         int n,p;
02057         
02058         n = Round(x);
02059         p = n-138;
02060         if (n>138) // neighbour intervals to the right
02061         {
02062                 Pr = x-1;
02063                 for (int k=2; k<=p; k++)
02064                         Pr *= x-k;
02065                 y = Pr*gam_S13(x-p);
02066         }
02067         else // left neighbour intervals and S13 itself
02068         {
02069                 p = -p;  Pr = x;
02070                 for (int k=1; k<=p-1; k++) 
02071                         Pr *= x+k;
02072                 y = (p==0)? gam_S13(x) : gam_S13(x+p)/Pr;               
02073         }
02074         
02075         return y;
02076 }       
02077 
02078 real gam_S14(const real& x)
02079 // Calculating approximations for Gamma(x) in S14 = [146.5,147.5];
02080 // Rel. error bound:  8.570174e-16;
02081 // Blomquist, 15.06.2009;               
02082 {
02083         real y(0),v;
02084 
02085         // Continued fraction:  K_5(v), v = 1/(x-x0), x0=146.94921875;
02086         if (x==146.94921875)
02087                 y = q_gams14_b[0];
02088         else
02089         {
02090                 v = 1/(x-146.94921875);
02091                 
02092                 y = q_gams14_a[5] / (  v + q_gams14_b[5]);
02093                 y = q_gams14_a[4] / ( (v + q_gams14_b[4]) + y );
02094                 y = q_gams14_a[3] / ( (v + q_gams14_b[3]) + y );
02095                 y = q_gams14_a[2] / ( (v + q_gams14_b[2]) + y );
02096                 y = q_gams14_a[1] / ( (v + q_gams14_b[1]) + y ) + q_gams14_b[0];
02097         }       
02098         y += 1;
02099         y *= q_ex10(2*x);
02100         times2pown(y,-133);
02101                 
02102         return y;
02103 }       
02104 
02105 real gamma_S14(const real& x)
02106 // Calculating approximations for Gamma(x) in [142.5,150.5];
02107 // eps = 1.745196e-15;
02108 // Blomquist, 16.06.2009;
02109 {
02110         real y,Pr;
02111         int n,p;
02112         
02113         n = Round(x);
02114         p = n-147;
02115         if (n>147) // neighbour intervals to the right:
02116         {
02117                 Pr = x-1;
02118                 for (int k=2; k<=p; k++)
02119                         Pr *= x-k;
02120                 y = Pr*gam_S14(x-p);
02121         }
02122         else // left neighbour intervals and S14 itself:
02123         {
02124                 p = -p;  Pr = x;
02125                 for (int k=1; k<=p-1; k++) 
02126                         Pr *= x+k;
02127                 y = (p==0)? gam_S14(x) : gam_S14(x+p)/Pr;               
02128         }
02129         
02130         return y;
02131 }       
02132 
02133 real gam_S15(const real& x)
02134 // Calculating approximations for Gamma(x) in S15 = [153.5,154.5];
02135 // Rel. error bound:  1.279047e-15;
02136 // Blomquist, 16.06.2009;               
02137 {
02138         real y(0),v;
02139 
02140         // Continued fraction:  K_5(v), v = 1/(x-x0), x0=153.90234375;
02141         if (x==153.90234375)
02142                 y = q_gams15_b[0];
02143         else
02144         {
02145                 v = 1/(x-153.90234375);
02146                 
02147                 y = q_gams15_a[5] / (  v + q_gams15_b[5]);
02148                 y = q_gams15_a[4] / ( (v + q_gams15_b[4]) + y );
02149                 y = q_gams15_a[3] / ( (v + q_gams15_b[3]) + y );
02150                 y = q_gams15_a[2] / ( (v + q_gams15_b[2]) + y );
02151                 y = q_gams15_a[1] / ( (v + q_gams15_b[1]) + y ) + q_gams15_b[0];
02152         }       
02153         y += 1;
02154         v = q_ex10(x);
02155         times2pown(y,-129);
02156         y *= v;
02157         y *= v;
02158 
02159         return y;
02160 }
02161 
02162 real gamma_S15(const real& x)
02163 // Calculating approximations for Gamma(x) in [150.5,157.5];
02164 // eps = 1.945181e-15;
02165 // Blomquist, 19.06.2009;
02166 {
02167         real y,Pr;
02168         int n,p;
02169         
02170         n = Round(x);
02171         p = n-154;
02172         if (n>154) // neighbour intervals to the right
02173         {
02174                 Pr = x-1;
02175                 for (int k=2; k<=p; k++)
02176                         Pr *= x-k;
02177                 y = Pr*gam_S15(x-p);
02178         }
02179         else // left neighbour intervals and S15 itself
02180         {
02181                 p = -p;  Pr = x;
02182                 for (int k=1; k<=p-1; k++) 
02183                         Pr *= x+k;
02184                 y = (p==0)? gam_S15(x) : gam_S15(x+p)/Pr;               
02185         }
02186         
02187         return y;
02188 }       
02189 
02190 real gam_S16(const real& x)
02191 // Calculating approximations for Gamma(x) in S16 = [160.5,161.5];
02192 // Rel. error bound:  1.381209e-15;
02193 // Blomquist, 19.06.2009;               
02194 {
02195         real y(0),v;
02196 
02197         // Continued fraction:  K_5(v), v = 1/(x-x0), x0=161.08984375;
02198         if (x==161.08984375)
02199                 y = q_gams16_b[0];
02200         else
02201         {
02202                 v = 1/(x-161.08984375);
02203                 
02204                 y = q_gams16_a[5] / (  v + q_gams16_b[5]);
02205                 y = q_gams16_a[4] / ( (v + q_gams16_b[4]) + y );
02206                 y = q_gams16_a[3] / ( (v + q_gams16_b[3]) + y );
02207                 y = q_gams16_a[2] / ( (v + q_gams16_b[2]) + y );
02208                 y = q_gams16_a[1] / ( (v + q_gams16_b[1]) + y ) + q_gams16_b[0];
02209         }       
02210         y += 1;
02211         v = q_ex10(x);
02212         times2pown(v,-62);
02213         y *= sqr(v);
02214         return y;
02215 }
02216 
02217 real gamma_S16(const real& x)
02218 // Calculating approximations for Gamma(x) in [157.5,164.5];
02219 // eps = 2.047343e-15;
02220 // Blomquist, 20.06.2009;
02221 {
02222         real y,Pr;
02223         int n,p;
02224         
02225         n = Round(x);
02226         p = n-161;
02227         if (n>161) // neighbour intervals to the right
02228         {
02229                 Pr = x-1;
02230                 for (int k=2; k<=p; k++)
02231                         Pr *= x-k;
02232                 y = Pr*gam_S16(x-p);
02233         }
02234         else  // left neighbour intervals and S16 itself
02235         {
02236                 p = -p;  Pr = x;
02237                 for (int k=1; k<=p-1; k++) 
02238                         Pr *= x+k;
02239                 y = (p==0)? gam_S16(x) : gam_S16(x+p)/Pr;               
02240         }
02241         
02242         return y;
02243 }       
02244 
02245 real gam_S17(const real& x)
02246 // Calculating approximations for Gamma(x) in:  S17 = [167.5,168.5];
02247 // Rel. error bound:  1.390239e-15;
02248 // Blomquist, 20.06.2009;               
02249 {
02250         real y(0),v;
02251 
02252         // Continued fraction:  K_5(v), v = 1/(x-x0), x0=168.0;
02253         if (x==168.0)
02254                 y = q_gams17_b[0];
02255         else
02256         {
02257                 v = 1/(x-168.0);
02258                 
02259                 y = q_gams17_a[5] / (  v + q_gams17_b[5]);
02260                 y = q_gams17_a[4] / ( (v + q_gams17_b[4]) + y );
02261                 y = q_gams17_a[3] / ( (v + q_gams17_b[3]) + y );
02262                 y = q_gams17_a[2] / ( (v + q_gams17_b[2]) + y );
02263                 y = q_gams17_a[1] / ( (v + q_gams17_b[1]) + y ) + q_gams17_b[0];
02264         }       
02265         y += 1;
02266         v = q_ex10(x);
02267         times2pown(y,-119);
02268         y *= v;
02269         y *= v; 
02270         return y;
02271 }
02272 
02273 real gamma_S17(const real& x)
02274 // Calculating approximations for Gamma(x), x in [164.5,171.5];
02275 // eps = 2.056373e-15;
02276 // Blomquist, 20.06.2009;
02277 {
02278         real y,Pr;
02279         int n,p;
02280         
02281         n = Round(x);
02282         p = n-168;
02283         if (n>168) // neighbour intervals to the right
02284         {
02285                 Pr = x-1;
02286                 for (int k=2; k<=p; k++)
02287                         Pr *= x-k;
02288                 y = Pr*gam_S17(x-p);
02289         }
02290         else // left neighbour intervals and S17 itself
02291         {
02292                 p = -p;  Pr = x;
02293                 for (int k=1; k<=p-1; k++) 
02294                         Pr *= x+k;
02295                 y = (p==0)? gam_S17(x) : gam_S17(x+p)/Pr;               
02296         }
02297         
02298         return y;
02299 }       
02300 
02301 real gamma_05(const real& x)
02302 // Calculating Gamma(x), x in (-0.5,+171.5]
02303 // eps = 2.082627e-15;
02304 // Blomquist, 21.06.2009;
02305 {
02306         real y(0);
02307         int Nr;
02308         
02309         Nr = int_no(gam_f85,19,x);
02310         
02311         switch(Nr)
02312         {
02313                 case 0:  y = gamma_S0(x);   break; // x in (-0.5,8.5);
02314                 case 1:  y = gamma_S1(x);   break; // x in [8.5,16.5);
02315                 case 2:  y = gamma_S2(x);   break; // x in [16.5,24.5);
02316                 case 3:  y = gamma_S3(x);   break; // x in [24.5,35.5);
02317                 case 4:  y = gamma_S4(x);   break; // x in [35.5,46.5);
02318                 case 5:  y = gamma_S5(x);   break; // x in [46.5,57.5);
02319                 case 6:  y = gamma_S6(x);   break; // x in [57.5,68.5);
02320                 case 7:  y = gamma_S7(x);   break; // x in [68.5,79.5);
02321                 case 8:  y = gamma_S8(x);   break; // x in [79.5,90.5);
02322                 case 9:  y = gamma_S9(x);   break; // x in [90.5,101.5);
02323                 case 10: y = gamma_S10(x);  break; // x in [101.5,112.5);
02324                 case 11: y = gamma_S11(x);  break; // x in [112.5,122.5);
02325                 case 12: y = gamma_S12(x);  break; // x in [122.5,132.5);
02326                 case 13: y = gamma_S13(x);  break; // x in [132.5,142.5);
02327                 case 14: y = gamma_S14(x);  break; // x in [142.5,150.5);
02328                 case 15: y = gamma_S15(x);  break; // x in [150.5,157.5);
02329                 case 16: y = gamma_S16(x);  break; // x in [157.5,164.5);
02330                 case 17: y = gamma_S17(x);  break; // x in [164.5,171.5);
02331 
02332                 default: // n=18; x in [171.5, ...];
02333                         y = gamma_S17(x);
02334         }
02335         
02336         return y;
02337 }
02338 
02339 real gammar(const real& x)
02340 // Calculating 1/Gamma(x) for x in [-170.0,+171.0];
02341 // eps = 2.866906e-15;
02342 // Blomquist, 21.06.2009;
02343 {
02344         real y;
02345         
02346         if (x<-170.0 || x>171.0) 
02347                 cxscthrow(STD_FKT_OUT_OF_DEF("real gammar(const real& x)"));
02348         
02349         if (x<=-0.5)
02350                 y = -sinpix_pi(x)*x*gamma_05(-x);
02351         else
02352                 if (x<=8.5)
02353                         y = gammar_S0(x);
02354         else y = 1/gamma_05(x);
02355         
02356         return y;
02357 }
02358 
02359 real gamma(const real& x)
02360 // Calculating Gamma(x) for x in [-170.0,+171.5]
02361 // eps = 3.088951e-15;
02362 // Blomquist, 21.06.2009;
02363 {
02364         real y;
02365         
02366         if (x>171.5 || x<-170.0) 
02367                 cxscthrow(STD_FKT_OUT_OF_DEF("real gamma(const real& x)"));
02368         
02369         if (x<=-0.5)
02370                 y = -1/( sinpix_pi(x)*x*gamma_05(-x) );
02371         else
02372                 y = gamma_05(x);
02373         
02374         return y;
02375 }
02376 
02377 
02378 extern "C" {
02379    void r_lfsr(void) {;} // Siehe real.hpp in real_ari...?!?!
02380 }
02381 
02382 } // namespace cxsc
02383