C-XSC - A C++ Class Library for Extended Scientific Computing
2.5.4
|
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