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: l_imath.cpp,v 1.38 2014/01/30 17:23:46 cxsc Exp $ */ 00025 00026 #include <math.h> 00027 #include "l_imath.hpp" 00028 #include "imath.hpp" 00029 #include "rmath.hpp" 00030 00031 namespace cxsc { 00032 00033 #define CXSC_One 1. 00034 #define CXSC_Zero 0. 00035 #define CXSC_MinusOne -1. 00036 00037 l_interval pow(const l_interval & x, const l_interval & e) throw(ERROR_LINTERVAL_STD_FKT_OUT_OF_DEF) 00038 { 00039 int stagsave = stagprec, 00040 stagmax = 19, 00041 intexp; 00042 00043 bool fertig; 00044 00045 l_interval y; 00046 interval dx = interval(x), 00047 de = interval(e), 00048 einfachgenau; 00049 real supabsde = Sup(abs(de)); 00050 00051 einfachgenau = pow(dx,de); 00052 00053 fertig = false; 00054 if (Inf(de) == Sup(de)) // && 00055 if (supabsde < 32768.0) 00056 { 00057 intexp = int(_double(real(Sup(e)))); 00058 if (real(intexp) == Sup(e)) 00059 { 00060 y = power(x,intexp); // Integer-Potenz wesentlich schneller 00061 fertig = true; 00062 } 00063 } 00064 00065 if (!fertig) 00066 { 00067 if (Inf(x) < 0.0) 00068 cxscthrow(ERROR_LINTERVAL_STD_FKT_OUT_OF_DEF("l_interval pow(const l_interval & x, const l_interval & e)")); 00069 else if (Inf(dx) == Sup(dx) && Sup(dx) == CXSC_One) 00070 y = x; 00071 else if (Inf(de) == Sup(de) && Sup(de) == CXSC_One) 00072 y = x; 00073 else if (Inf(de) == Sup(de) && Sup(de) == CXSC_Zero) 00074 y = 1.0; 00075 else 00076 { 00077 if (stagprec < stagmax) 00078 stagprec++; 00079 else 00080 stagprec = stagmax; 00081 y = exp(e*ln(x)); 00082 stagprec = stagsave; 00083 y = adjust(y); 00084 y = y & einfachgenau; 00085 } 00086 } 00087 00088 return y; 00089 } 00090 00091 l_interval power(const l_interval &x, int n) // Power(x,n) 00092 { 00093 int stagsave = stagprec, 00094 stagmax = 19; 00095 00096 bool neg = false; 00097 00098 long int zhi = 2; 00099 interval dx = interval(x), 00100 einfachgenau; 00101 l_interval y, neu; 00102 00103 einfachgenau = Power(dx,n); 00104 00105 if (Inf(dx) == Sup(dx) && Sup(dx) == CXSC_One) 00106 y = x; 00107 else if (n == 0) 00108 y = adjust(l_interval(1.0)); 00109 else 00110 { 00111 if (stagprec < stagmax) 00112 stagprec++; 00113 else 00114 stagprec = stagmax; 00115 00116 if (n == 1) 00117 y = x; 00118 else if (n == 2) 00119 y = sqr(x); 00120 else 00121 { 00122 if (n < 0) 00123 { 00124 neg = true; 00125 n = -n; 00126 } 00127 // Initialisierung 00128 if (n%2) 00129 y = x; 00130 else 00131 y = l_interval(1.0); // Praezision wird bei 1 Mult. auf 00132 // aktuellen Wert gesetzt; 00133 // Berechnugn durch binaere Darstellung der n 00134 neu = sqr(x); // neu = x*x; 00135 do { 00136 if ((n/zhi)%2) y *= neu; 00137 zhi += zhi; 00138 if (zhi <= n) // letzte Mult. entfaellt --> schneller 00139 neu *= neu; 00140 } while (zhi <= n); 00141 00142 if (neg) 00143 y = 1.0/(y); 00144 } 00145 stagprec = stagsave; 00146 y = adjust(y); 00147 y = y & einfachgenau; 00148 } 00149 00150 return y; 00151 } 00152 00153 l_interval sqr(const l_interval & x) // Sqr(x) 00154 { 00155 l_interval y; 00156 00157 if (Inf(x) >= 0.0) /* result = [x.inf*x.inf, x.sup*x.sup] */ 00158 y = x * x; 00159 else if (Sup(x)<= 0.0) 00160 { 00161 /* result = [x.sup*x.sup, x.inf*x.inf] */ 00162 y = l_interval (-Sup(x) , -Inf(x)); 00163 y = (y) * (y); 00164 } else 00165 { 00166 /* result = [0.0, max(x.sup*x.sup,x.inf*x.inf)] */ 00167 if (abs(Inf(x)) >= abs(Sup(x))) 00168 y = l_interval(0.0, abs(Inf(x))); 00169 else 00170 y = l_interval(0.0, abs(Sup(x))); 00171 y = (y) * (y); 00172 } 00173 return y; 00174 } 00175 00176 l_interval sqrt(const l_interval & x) 00177 throw(ERROR_LINTERVAL_STD_FKT_OUT_OF_DEF) // Sqrt(x) 00178 { // Blomquist: scaling with 2^ex is necessary if expo(Sup(dx)) is too small! 00179 int stagsave = stagprec, 00180 stagmax = 30, 00181 stagsave2; 00182 interval dx = interval(x), 00183 einfachgenau; 00184 l_interval a1,y,t,mt; 00185 bool Inf_Zero; 00186 00187 einfachgenau = sqrt(dx); 00188 00189 if (Inf(x) < 0.0) 00190 cxscthrow(ERROR_LINTERVAL_STD_FKT_OUT_OF_DEF("l_interval sqrt(const l_interval & x)")); 00191 else if (Inf(dx) == Sup(dx) && 00192 (Sup(dx) == CXSC_Zero || Sup(dx) == CXSC_One)) 00193 y = x; 00194 else 00195 { // scaling necessary if exponent ex < 0 00196 l_interval x1 = x; 00197 Inf_Zero = (Inf(dx)==0); 00198 if (Inf_Zero) x1 = Sup(x1); 00199 int ex = expo(Sup(dx)); 00200 if (ex>0) ex = 0; // scaling unnecessary if ex>0 00201 else 00202 { 00203 ex = -ex; // ex >= 0 00204 if (ex > 1023) ex = 1023; // ex==1023 is sufficient 00205 if (ex%2) ex--; // ex>=0 is even 00206 } 00207 if (ex) times2pown(x1,ex); // ex >= 0 ---> exact scaling! 00208 // Interval-Newton-methode: y = m(y)-f(m(y))/f'(y) 00209 t = sqrt(interval(x1)); 00210 if (stagprec < stagmax) 00211 stagsave2 = stagprec+1; 00212 else 00213 stagsave2 = stagmax; 00214 stagprec = 1; 00215 while (stagprec < stagsave2) 00216 { 00217 stagprec += stagprec; 00218 if (stagprec > stagmax) 00219 stagprec = stagmax; 00220 mt = mid(t); times2pown(t,1); 00221 t = mt-((mt*mt-x1)/t); 00222 } 00223 if (ex) times2pown(t,-ex/2); // ex!=0 --> backscaling with 2^(-ex/2) 00224 stagprec = stagsave; // restore the previous precision 00225 y = adjust(t); // matching to previous stagprec. 00226 if (Inf_Zero) SetInf(y,0.0); 00227 y = y & einfachgenau; // seeking optimal inclusion with intersection 00228 } 00229 return y; 00230 } 00231 00232 l_interval sqrt(const l_interval &x, int n) 00233 throw(ERROR_LINTERVAL_STD_FKT_OUT_OF_DEF) // Sqrt(x,n) 00234 // -1073741822 <= n <= +1073741823, sonst autom. Programm-Abbruch 00235 // sqrt(x,n) jetzt mit Skalierung --> Hohe Genauigkeit in allen Bereichen! 00236 // Blomquist, 28.12.03; 00237 { 00238 int stagsave = stagprec, 00239 stagshort, staglong, i, ex, N, 00240 stagmax = 19, 00241 max = 2; 00242 l_interval my, fy, corr, y, xx; 00243 interval dx = interval(x), einfachgenau; 00244 00245 einfachgenau = sqrt(dx,n); 00246 00247 if (Inf(x) < 0.0) 00248 cxscthrow(ERROR_LINTERVAL_STD_FKT_OUT_OF_DEF 00249 ("l_interval sqrt(const l_interval &x, int n)")); 00250 else if (stagprec == 1) 00251 y = pow(dx, interval(1.0)/n); 00252 else if (Inf(dx) == Sup(dx) && (Sup(dx) == CXSC_Zero || 00253 Sup(dx) == CXSC_One)) 00254 y = x; 00255 else 00256 { 00257 if (stagprec < stagmax) 00258 stagprec++; 00259 else stagprec = stagmax; 00260 00261 while (max < stagprec) 00262 max += max; // quadratische Konvergenz kann 00263 // nicht eingehalten werden, 00264 max += max; // deshalb eine Schleife mehr 00265 00266 xx = x; 00267 ex = expo(Sup(dx)); 00268 N = -ex; 00269 // Skalierung mit 2^N, so dass dx,xx etwa = 1; ---------------- 00270 times2pown(dx,N); 00271 if (N>1023) { 00272 times2pown(xx,1023); 00273 times2pown(xx,N-1023); 00274 } 00275 else times2pown(xx,N); 00276 // Skalierung beendet, Blomquist 28.12.03 --------------------- 00277 00278 y = pow(dx, interval(1.0/n)); 00279 stagprec = 1; 00280 // Intervall-Newton-Verfahren 00281 for (i = 2; i <= max; i += i) 00282 { 00283 // Verdoppelung der Genauigkeit: 00284 stagshort = stagprec; 00285 stagprec += stagprec; // Verdoppelung hier 00286 if (stagprec > stagmax) 00287 stagprec = stagmax; 00288 staglong = stagprec; 00289 my = l_interval(mid(y)); 00290 fy = power(my, n)-xx; 00291 // Fehlerauswertung nur in halber Genauigkeit notwendig! 00292 stagprec = stagshort; 00293 corr = fy/(real(n)*power(y, n-1)); 00294 stagprec = staglong; 00295 // Fehlerkorrektur in normaler Genauigkeit: 00296 y = my-corr; 00297 } 00298 // Rueckskalierung mit dem Faktor 2^(-N/n): 00299 fy = l_interval(-N)/n; 00300 y *= exp(fy*li_ln2()); // li_ln2() = ln(2) 00301 00302 stagprec = stagsave; 00303 y = adjust(y); // Anpassung an Ausgangs-stagprec = stagsave. 00304 y = y & einfachgenau; // Falls y breiter ist als einfachgenau 00305 } 00306 00307 return y; 00308 } // sqrt(x,n) 00309 00310 l_interval sqrt1px2(const l_interval& x) throw() 00311 // Calculation of an optimal inclusion of sqrt(1+x^2); Blomquist, 13.12.02; 00312 // With stagmax=19 we get about 16*19=304 exact decimal digits. 00313 { // First step: Calculation of sqrt(1+x*x) in simple type interval: 00314 interval einfach = sqrt1px2( interval(x) ); // Only interval types 00315 int stagsave, stagmax=19; 00316 stagsave = stagprec; 00317 if (stagprec > stagmax) stagprec = stagmax; 00318 // Second step: Inclusion of sqrt(1+x^2) with stagprec <= 19 00319 const int exmax=512; 00320 l_interval y = abs(x); 00321 l_real lr = Sup( 1 + l_interval(Sup(y)) ); 00322 interval z = interval(Sup(x)); 00323 int ex = expo(Sup(z)); 00324 if (ex > exmax) 00325 { // scaling to avoid overflow by sqr(y): 00326 ex = exmax - ex; // ex = 512 - ex; 00327 times2pown(y,ex); // scaling to avoid overflow by sqr(y) 00328 y = sqrt( comp(0.5,2*ex+1) + sqr(y) ); // sqr(y) without overflow! 00329 times2pown(y,-ex); // backscaling of the result with 2^(-ex) 00330 } else 00331 { // no scaling necessary: 00332 y = sqrt(1+sqr(x)); 00333 } 00334 if (Inf(y)<1.0) SetInf(y,1.0); // improvement, Blomquist, 25.02.07; 00335 if (Sup(y)>lr) SetSup(y,lr); // improvement, Blomquist, 26.02.07; 00336 stagprec = stagsave; // restore the old stagprec value 00337 y = adjust(y); // y gets the previous precision 00338 y = einfach & y; // This intersection delivers for intervals x with large 00339 // diameters the optimal result inclusion 00340 return y; 00341 } // sqrt1px2(...) 00342 00343 l_interval sqrtx2y2(const l_interval& x, const l_interval& y) throw() 00344 // Inclusion of sqrt(x^2+y^2); Blomquist, 14.12.02; 00345 { 00346 interval ia = abs(interval(x)), ib = abs(interval(y)), einfach; 00347 einfach = sqrtx2y2(ia,ib); // Inclusion only with type interval. 00348 if (!einfach) return l_interval(0.0); 00349 00350 int stagsave=stagprec, stagmax=19; 00351 if (stagprec>stagmax) stagprec = stagmax; 00352 00353 l_interval a=abs(x), b=abs(y), r; 00354 int exa=expo(Sup(ia)), exb=expo(Sup(ib)), ex; 00355 if (exb > exa) 00356 { // Permutation of a,b: 00357 r = a; a = b; b = r; 00358 ex = exa; exa = exb; exb = ex; 00359 } 00360 ex = 511 - exa; exa = 0; 00361 if (ex>1022) 00362 { exa = ex - 1022; ex = 1022; } // ex > 1022 --> scaling in two steps 00363 times2pown(a,ex); // is necessary! 00364 if (exa) times2pown(a,exa); 00365 times2pown(b,ex); // First step: scaling b with 2^ex; 00366 if (exa) times2pown(b,exa); // Second step: scaling b with 2^exa; 00367 r = sqrt(a*a + b*b); 00368 times2pown(r,-ex); // Backscaling, first step 00369 if (exa) times2pown(r,-exa); // Backscaling, second step 00370 stagprec = stagsave; 00371 r = adjust(r); 00372 r = einfach & r; 00373 return r; 00374 } // sqrtx2y2 00375 00376 l_interval sqrtp1m1(const l_interval& x) throw(STD_FKT_OUT_OF_DEF) 00377 // sqrtp1m1(x) calculates an inclusion of sqrt(x+1)-1; 00378 // Blomquist, 05.08.03; 00379 { 00380 int stagsave=stagprec, stagmax=19; 00381 stagprec++; 00382 if (stagprec>stagmax) stagprec = stagmax; 00383 l_interval y,tmp; 00384 interval z = interval(x); // z is an inclusion of x 00385 real r = Inf(z); 00386 if (r < -1) 00387 cxscthrow(STD_FKT_OUT_OF_DEF("l_interval sqrtp1m1(const l_interval&)")); 00388 const real c = 1e-10; 00389 tmp = x+1; 00390 y = x<=interval(-c,c) ? x / (sqrt(tmp)+1) : sqrt(tmp)-1; 00391 stagprec = stagsave; 00392 y = adjust(y); 00393 return y; 00394 } // sqrtp1m1 00395 00396 00397 l_interval sin(const l_interval & x) throw(ERROR_LINTERVAL_FAK_OVERFLOW) // Sin(x) 00398 { 00399 int stagsave = stagprec, 00400 stagmax = 19; 00401 l_interval pihalbe, 00402 y; 00403 interval dx = interval(x), 00404 einfachgenau; 00405 00406 einfachgenau = sin(dx); 00407 00408 if (stagprec == 1) 00409 y = sin(dx); 00410 else if (Inf(dx) == Sup(dx) && Sup(dx) == CXSC_Zero) 00411 y = adjust(l_interval(0.0)); 00412 else 00413 { 00414 if (stagprec < stagmax) 00415 stagprec++; 00416 else 00417 stagprec = stagmax; 00418 00419 // pihalbe = 2.0*atan(l_interval(1.0)); 00420 pihalbe = li_pi4(); 00421 times2pown(pihalbe,1); // Blomquist, 05.12.03; 00422 y = x-pihalbe; 00423 00424 try { 00425 y = cos(y); 00426 } 00427 catch(const ERROR_LINTERVAL_FAK_OVERFLOW &) // Damit Funktionsname stimmt 00428 { 00429 cxscthrow(ERROR_LINTERVAL_FAK_OVERFLOW("l_interval sin(const l_interval & x)")); 00430 } 00431 00432 00433 stagprec = stagsave; 00434 y = adjust(y); 00435 y = y & einfachgenau; 00436 } 00437 00438 return y; 00439 } 00440 00441 l_interval cos(const l_interval & x) throw(ERROR_LINTERVAL_FAK_OVERFLOW) // Cos(x) 00442 { 00443 long int mm = 6; 00444 int stagsave = stagprec, 00445 stagmax = 19, 00446 n = 0, 00447 digits = 53, 00448 sign = 0, 00449 degree, k, 00450 m2; 00451 00452 bool xinf = false, 00453 xsup = false, 00454 fertig = false; 00455 00456 real abst, m, eps, 00457 // bas, tn, t4, 00458 lneps, lnt, 00459 zk, 00460 zhn = 1.0, 00461 lnb = 0.69314718, 00462 fak = 720.0; // 6! 00463 interval dx = interval(x), 00464 einfachgenau, 00465 dt, extr, error; 00466 l_interval zwopi, t, t2, p, ph, test, 00467 y; 00468 00469 einfachgenau = cos(dx); 00470 if (stagprec == 1) 00471 y = cos(dx); 00472 else if (Inf(dx) == Sup(dx) && Sup(dx) == CXSC_Zero) 00473 y = adjust(l_interval(1.0)); 00474 else 00475 { 00476 if (stagprec < stagmax) 00477 stagprec++; 00478 else 00479 stagprec = stagmax; 00480 00481 // stagprec++; 00482 // zwopi = 8.0*atan(l_interval(1.0)); 00483 zwopi = li_pi4(); 00484 times2pown(zwopi,3); // Blomquist, 05.12.03; 00485 // stagprec--; 00486 00487 // erste Argumentreduktion ( cos(x) = cos(x+2k*pi) ) 00488 if (Sup(zwopi) < Sup(abs(x))) 00489 { 00490 m = floor(_double(real(Sup((x/zwopi))))); // floor rundet zur naechsten 00491 t = x-m*zwopi; // ganzen Zahl kleiner x 00492 if (Sup(zwopi) < Sup(abs(t))) 00493 { // das ganze nochmal, falls m uebergelaufen ist! 00494 // ohne Sup wird Inclusion geprueft! 00495 m = floor(_double(real(Sup((t/zwopi)))));// rundet zur naechsten Zahl < x 00496 t = t-m*zwopi; 00497 } 00498 } else 00499 t = x; 00500 00501 // ueberpruefen, ob Maximum oder Minimum im Inneren des Intervalls 00502 extr = interval(2.0/zwopi*t); 00503 m2 = int(double(floor(_double(Sup(extr))))); 00504 if (interval(real(m2)) <= extr) 00505 { 00506 if (interval(real(m2-1),real(m2)) <= extr) 00507 { 00508 y = l_interval(-1.0,1.0); 00509 fertig = true; 00510 } else 00511 if (m2%2) 00512 xinf = TRUE; 00513 else 00514 xsup = TRUE; 00515 } 00516 00517 if (!(fertig)) 00518 { 00519 // zweite Argumentreduktion 00520 dt = interval(t); 00521 // eps = 0.01/(stagprec*stagprec*stagprec); 00522 eps = 0.01; 00523 while (Sup(abs(dt))/zhn >= eps) 00524 { 00525 n++; 00526 zhn += zhn; 00527 } 00528 t /= zhn; 00529 00530 // Abschaetzung der Genauigkeit 00531 00532 t2 = t*t; 00533 abst = real(Sup(abs(t))); 00534 if (abst < MinReal) 00535 abst = MinReal; // nur zur Sicherheit 00536 lnt = ln(abst); 00537 lneps = (1.0-digits*stagprec)*lnb; 00538 while (lneps-(real(mm)*lnt+ln(2.0/fak)) <= 0.0) 00539 { 00540 mm += 4; 00541 if (mm > 170) 00542 { 00543 // 170! = 7.26*10^306 00544 cxscthrow(ERROR_LINTERVAL_FAK_OVERFLOW("l_interval cos(const l_interval & x)")); 00545 mm = 170; 00546 break; 00547 } 00548 fak *= mm*(mm-1.0)*(mm-2.0)*(mm-3.0); 00549 } 00550 /* 00551 bas = real(Inf(power(l_interval(2.0),(1-digits*stagprec)))); 00552 tn = abs(real(Sup(power(t,6)))); 00553 t4 = real(mid(power(t,4))); 00554 while (bas-2.0*tn/fak <= 0.0) 00555 { 00556 mm += 4; 00557 if (mm > 170) 00558 { // 170! = 7.26*10^306 00559 errmon (ERROR_LINTERVAL(FAKOVERFLOW)); 00560 mm = 170; 00561 break; 00562 } 00563 tn *= t4; 00564 fak *= mm*(mm-1.0)*(mm-2.0)*(mm-3.0); 00565 } 00566 */ 00567 00568 degree = mm-2; // Achtung mm := 2n+2 ! 00569 00570 // Polynomauswertung 00571 00572 sign = (degree/2)%2; 00573 zk = real(degree)*real(degree-1); 00574 if (sign) 00575 p = -t2/zk; 00576 else 00577 p = t2/zk; 00578 for (k = degree-2; k >= 2; k -= 2) 00579 { 00580 sign = 1-sign; 00581 if (sign) 00582 p -= 1.0; 00583 else 00584 p += 1.0; 00585 zk = real(k)*real(k-1); 00586 p *= t2/zk; 00587 } 00588 00589 error = pow(interval(2.0), interval(1.0-digits*stagprec)) 00590 * interval(-0.5,0.5); 00591 p = l_interval(1.0)+error+p; 00592 00593 // Rueckgaengigmachung der zweiten Argumentreduktion 00594 for (int i = 0; i < n; i++) 00595 p = 2.0*p*p-1.0; 00596 00597 stagprec = stagsave; 00598 y = adjust(p); 00599 if (Inf(y) < -1.0) 00600 SetInf(y,-1.0); 00601 if (Sup(y) > 1.0) 00602 SetSup(y,1.0); 00603 if (xinf) 00604 SetInf(y,-1.0); 00605 if (xsup) 00606 SetSup(y,1.0); 00607 } 00608 y = y & einfachgenau; 00609 } 00610 return y; 00611 } 00612 00613 l_interval tan(const l_interval & x) throw(ERROR_LINTERVAL_FAK_OVERFLOW,ERROR_LINTERVAL_STD_FKT_OUT_OF_DEF) // Tan(x) 00614 { 00615 interval dx = interval(x), 00616 einfachgenau; 00617 l_interval s, c, y; 00618 00619 einfachgenau = tan(dx); 00620 00621 if (stagprec == 1) 00622 y = tan(dx); 00623 else if (Inf(dx) == Sup(dx) && Sup(dx) == CXSC_Zero) 00624 y = adjust(l_interval(0.0)); 00625 else 00626 { 00627 try 00628 { 00629 c = cos(x); 00630 } 00631 catch(const ERROR_LINTERVAL_FAK_OVERFLOW &) // Damit Funktionsname stimmt 00632 { 00633 cxscthrow(ERROR_LINTERVAL_FAK_OVERFLOW("l_interval tan(const l_interval &x)")); 00634 } 00635 00636 if (interval(0.0) <= c) 00637 { 00638 cxscthrow(ERROR_LINTERVAL_STD_FKT_OUT_OF_DEF("l_interval tan(const l_interval &x)")); 00639 } 00640 try 00641 { 00642 s = sin(x); 00643 } 00644 catch(const ERROR_LINTERVAL_FAK_OVERFLOW &) 00645 { 00646 cxscthrow(ERROR_LINTERVAL_FAK_OVERFLOW("l_interval tan(const l_interval &x)")); 00647 } 00648 stagprec++; 00649 y = s/c; 00650 stagprec--; 00651 y = adjust(y); 00652 y = y & einfachgenau; 00653 } 00654 return y; 00655 } 00656 00657 l_interval cot(const l_interval & x) throw(ERROR_LINTERVAL_FAK_OVERFLOW,ERROR_LINTERVAL_STD_FKT_OUT_OF_DEF) // Cot(x) 00658 { 00659 interval dx = interval(x), 00660 einfachgenau; 00661 l_interval s, c, y; 00662 00663 einfachgenau = cot(dx); 00664 00665 if (stagprec == 1) 00666 y = tan(dx); 00667 else if (Inf(dx) == Sup(dx) && Sup(dx) == CXSC_Zero) 00668 y = adjust(l_interval(0.0)); 00669 else 00670 { 00671 try 00672 { 00673 s = sin(x); 00674 } 00675 catch(const ERROR_LINTERVAL_FAK_OVERFLOW &) // Damit Funktionsname stimmt 00676 { 00677 cxscthrow(ERROR_LINTERVAL_FAK_OVERFLOW("l_interval cot(const l_interval &x)")); 00678 } 00679 00680 if (interval(0.0) <= s) 00681 { 00682 cxscthrow(ERROR_LINTERVAL_STD_FKT_OUT_OF_DEF("l_interval cot(const l_interval &x)")); 00683 } 00684 00685 try 00686 { 00687 c = cos(x); 00688 } 00689 catch(const ERROR_LINTERVAL_FAK_OVERFLOW &) // Damit Funktionsname stimmt 00690 { 00691 cxscthrow(ERROR_LINTERVAL_FAK_OVERFLOW("l_interval cot(const l_interval &x)")); 00692 } 00693 00694 c = cos(x); 00695 stagprec++; 00696 y = c/s; 00697 stagprec--; 00698 y = adjust(y); 00699 y = y & einfachgenau; 00700 } 00701 00702 return y; 00703 } 00704 00705 l_interval asin(const l_interval & x) throw(ERROR_LINTERVAL_STD_FKT_OUT_OF_DEF) // ASin(x) 00706 { 00707 l_interval t, ta, u, pihalbe, 00708 y; 00709 interval dx = interval(x), 00710 einfachgenau; 00711 real supabsdx = Sup(abs(dx)), 00712 infdx = Inf(dx), 00713 supdx = Sup(dx); 00714 00715 einfachgenau = asin(dx); 00716 00717 stagprec++; 00718 // pihalbe = 2.0*atan(l_interval(1.0)); 00719 pihalbe = li_pi4(); 00720 times2pown(pihalbe,1); // Blomquist, 05.12.03; 00721 stagprec--; 00722 00723 if (Inf(x) < CXSC_MinusOne || Sup(x) > CXSC_One) 00724 cxscthrow(ERROR_LINTERVAL_STD_FKT_OUT_OF_DEF("l_interval asin(const l_interval & x)")); 00725 else if (stagprec == 1) 00726 y = asin(dx); 00727 else if (infdx == supdx && supdx == CXSC_Zero) 00728 y = adjust(l_interval(0.0)); 00729 else if (infdx == supdx && supabsdx == CXSC_One) 00730 { 00731 if (supdx == 1.0) 00732 y = pihalbe; 00733 else 00734 y = -pihalbe; 00735 } else 00736 { 00737 stagprec++; 00738 try 00739 { 00740 if (supabsdx <= 0.75) 00741 u = x; 00742 else 00743 u = 2.0*x*sqrt((1.0-x)*(1.0+x)); 00744 t = u/sqrt((1.0-u)*(1.0+u)); 00745 stagprec--; 00746 ta = atan(t); 00747 } 00748 catch(const ERROR_LINTERVAL_STD_FKT_OUT_OF_DEF &) 00749 { 00750 cxscthrow(ERROR_LINTERVAL_STD_FKT_OUT_OF_DEF("l_interval asin(const l_interval & x)")); 00751 } 00752 stagprec++; 00753 if (supabsdx <= 0.75) 00754 y = ta; 00755 else if (Sup(t) >= 0.0) 00756 y = pihalbe-0.5*ta; 00757 else 00758 y = -pihalbe-0.5*ta; 00759 stagprec--; 00760 y = adjust(y); 00761 y = y & einfachgenau; 00762 } 00763 00764 return y; 00765 } 00766 00767 l_interval acos(const l_interval & x) throw(ERROR_LINTERVAL_STD_FKT_OUT_OF_DEF) // ACos(x) 00768 { 00769 bool neg=false; 00770 l_interval pi, s, y; 00771 interval dx = interval(x), 00772 einfachgenau; 00773 real supabsdx = Sup(abs(dx)), 00774 infdx = Inf(dx), 00775 supdx = Sup(dx); 00776 00777 try { 00778 00779 einfachgenau = acos(dx); 00780 00781 stagprec++; 00782 // pi = 4.0*atan(l_interval(1.0)); 00783 pi = li_pi4(); 00784 times2pown(pi,2); // Blomquist, 05.12.03; 00785 stagprec--; 00786 if (Inf(x) < CXSC_MinusOne || Sup(x) > CXSC_One) 00787 cxscthrow(ERROR_LINTERVAL_STD_FKT_OUT_OF_DEF()); 00788 else if (stagprec == 1) 00789 y = acos(dx); 00790 else if (infdx == supdx && supabsdx == CXSC_One) 00791 { 00792 if (supdx == 1.0) 00793 y = adjust(l_interval(0.0)); 00794 else 00795 y = adjust(pi); 00796 } else 00797 { 00798 if (supdx < 0.0) 00799 { 00800 y = -x; 00801 neg = true; 00802 } else 00803 { 00804 y = x; 00805 neg = false; 00806 } 00807 stagprec++; 00808 if (supabsdx > 0.75) 00809 { 00810 y = (1.0-(y))*(1.0+(y)); 00811 // stagprec--; 00812 y = asin(sqrt(y)); 00813 // stagprec++; 00814 if (neg) 00815 y = pi-(y); 00816 } else 00817 { 00818 // stagprec--; 00819 s = asin(y); 00820 // stagprec++; 00821 if (neg) 00822 y = 0.5*pi+s; 00823 else 00824 y = 0.5*pi-s; 00825 } 00826 00827 // Fehler in der Systemumgebung, deshalb die Aufblaehung des Intervalls 00828 real err1, err2; 00829 interval error; 00830 err1 = 5.0*real(diam(y)); 00831 err2 = 5.0*power(10.0,-16*(stagprec-1)-1); 00832 if (err1 > err2) 00833 error = interval(-err1, err1); 00834 else 00835 error = interval(-err2, err2); 00836 y += error; 00837 00838 stagprec--; 00839 y = adjust(y); 00840 y = y & einfachgenau; 00841 } 00842 } 00843 catch(const ERROR_LINTERVAL_STD_FKT_OUT_OF_DEF &) 00844 { 00845 cxscthrow(ERROR_LINTERVAL_STD_FKT_OUT_OF_DEF("l_interval acos(const l_interval & x)")); 00846 } 00847 return y; 00848 } 00849 00850 static real CXSC_ln2[21]; // CXSC_ln2[0], ... CXSC_ln2[20] 00851 static bool CXSC_ln2_initialized = false; 00852 00853 // l_interval li_ln2() throw() 00854 l_interval Ln2_l_interval() throw() 00855 // Inclusion of ln(2), Blomquist, 04.12.03; 00856 { 00857 l_interval y; 00858 int stagsave = stagprec, 00859 stagmax = 20; 00860 00861 if (!CXSC_ln2_initialized) 00862 { 00863 std::string str; 00864 std::cout << SaveOpt; 00865 std::cout << Hex; 00866 str = "+162E42FEFA39EFe3FE"; 00867 str >> CXSC_ln2[0]; 00868 str = "+1ABC9E3B39803Fe3C7"; 00869 str >> CXSC_ln2[1]; 00870 str = "+17B57A079A1934e390"; 00871 str >> CXSC_ln2[2]; 00872 str = "-1ACE93A4EBE5D1e35A"; 00873 str >> CXSC_ln2[3]; 00874 str = "-123A2A82EA0C24e324"; 00875 str >> CXSC_ln2[4]; 00876 str = "+1D881B7AEB2615e2ED"; 00877 str >> CXSC_ln2[5]; 00878 str = "+19552FB4AFA1B1e2B7"; 00879 str >> CXSC_ln2[6]; 00880 str = "+1DA5D5C6B82704e27E"; 00881 str >> CXSC_ln2[7]; 00882 str = "+14427573B29117e247"; 00883 str >> CXSC_ln2[8]; 00884 str = "-191F6B05A4D7A7e211"; 00885 str >> CXSC_ln2[9]; 00886 str = "-1DB5173AE53426e1DB"; 00887 str >> CXSC_ln2[10]; 00888 str = "+11317C387EB9EBe1A3"; 00889 str >> CXSC_ln2[11]; 00890 str = "-190F13B267F137e16D"; 00891 str >> CXSC_ln2[12]; 00892 str = "+16FA0EC7657F75e137"; 00893 str >> CXSC_ln2[13]; 00894 str = "-1234C5E1398A6Be101"; 00895 str >> CXSC_ln2[14]; 00896 str = "+1195EBBF4D7A70e0CA"; 00897 str >> CXSC_ln2[15]; 00898 str = "+18192432AFD0C4e094"; 00899 str >> CXSC_ln2[16]; 00900 str = "-1A1BE38BA4BA4De05E"; 00901 str >> CXSC_ln2[17]; 00902 str = "-1D7860151CFC06e024"; 00903 str >> CXSC_ln2[18]; 00904 str = "+1000032847ED6Fe000"; 00905 str >> CXSC_ln2[19]; 00906 str = "+1000032847ED70e000"; 00907 str >> CXSC_ln2[20]; 00908 00909 CXSC_ln2_initialized = true; 00910 std::cout << RestoreOpt; 00911 } 00912 stagprec = stagmax; 00913 y = adjust(l_interval(0.0)); 00914 for (int i=0; i <= stagmax; i++) 00915 y.data[i] = CXSC_ln2[i]; 00916 stagprec = stagsave; 00917 y = adjust(y); 00918 return y; 00919 } 00920 00921 static real CXSC_ln10[21]; // CXSC_ln10[0], ... CXSC_ln10[20] 00922 static bool CXSC_ln10_initialized = false; 00923 00924 // l_interval li_ln10() throw() 00925 l_interval Ln10_l_interval() throw() 00926 { 00927 l_interval y; 00928 int stagsave = stagprec, 00929 stagmax = 20; 00930 00931 if (!CXSC_ln10_initialized) 00932 { 00933 std::string str; 00934 std::cout << SaveOpt; 00935 std::cout << Hex; 00936 str = "+126BB1BBB55516e400"; 00937 str >> CXSC_ln10[0]; 00938 str = "-1F48AD494EA3E9e3CA"; 00939 str >> CXSC_ln10[1]; 00940 str = "-19EBAE3AE0260Ce394"; 00941 str >> CXSC_ln10[2]; 00942 str = "-12D10378BE1CF1e35E"; 00943 str >> CXSC_ln10[3]; 00944 str = "+10403E05AE52C6e328"; 00945 str >> CXSC_ln10[4]; 00946 str = "-1FA509CAFDF466e2F0"; 00947 str >> CXSC_ln10[5]; 00948 str = "-1C79A1FE9D0795e2BA"; 00949 str >> CXSC_ln10[6]; 00950 str = "+1058C448308218e284"; 00951 str >> CXSC_ln10[7]; 00952 str = "-1D250470877BFDe24D"; 00953 str >> CXSC_ln10[8]; 00954 str = "-1AE92987D3075De215"; 00955 str >> CXSC_ln10[9]; 00956 str = "-1D5CDBB8626956e1DF"; 00957 str >> CXSC_ln10[10]; 00958 str = "-13C4F27CE0410Ae1A9"; 00959 str >> CXSC_ln10[11]; 00960 str = "+1B3AC12ACF1BE9e173"; 00961 str >> CXSC_ln10[12]; 00962 str = "+1161BB49D219C8e13D"; 00963 str >> CXSC_ln10[13]; 00964 str = "-110D6613293728e107"; 00965 str >> CXSC_ln10[14]; 00966 str = "+142163A4CDA351e0CF"; 00967 str >> CXSC_ln10[15]; 00968 str = "+1E2713D6C22C16e097"; 00969 str >> CXSC_ln10[16]; 00970 str = "-15090EF85CB0ADe05E"; 00971 str >> CXSC_ln10[17]; 00972 str = "-1C5B3E859F876Ee027"; 00973 str >> CXSC_ln10[18]; 00974 str = "-10000703552C52e000"; 00975 str >> CXSC_ln10[19]; 00976 str = "-10000703552C51e000"; 00977 str >> CXSC_ln10[20]; 00978 00979 CXSC_ln10_initialized = true; 00980 std::cout << RestoreOpt; 00981 } 00982 stagprec = stagmax; 00983 y = adjust(l_interval(0.0)); 00984 for (int i=0; i <= stagmax; i++) 00985 y.data[i] = CXSC_ln10[i]; 00986 stagprec = stagsave; 00987 y = adjust(y); 00988 return y; 00989 } 00990 00991 static real CXSC_Rln10[21]; // CXSC_Rln10[0], ... CXSC_Rln10[20] 00992 static bool CXSC_Rln10_initialized = false; 00993 00994 // l_interval li_Rln10() throw() 00995 l_interval Ln10r_l_interval() throw() 00996 { 00997 l_interval y; 00998 int stagsave = stagprec, 00999 stagmax = 20; 01000 01001 if (!CXSC_Rln10_initialized) 01002 { 01003 std::string str; 01004 std::cout << SaveOpt; 01005 std::cout << Hex; 01006 str = "+1BCB7B1526E50Ee3FD"; 01007 str >> CXSC_Rln10[0]; 01008 str = "+195355BAAAFAD3e3C6"; 01009 str >> CXSC_Rln10[1]; 01010 str = "+1EE191F71A3012e38F"; 01011 str >> CXSC_Rln10[2]; 01012 str = "+17268808E8FCB5e358"; 01013 str >> CXSC_Rln10[3]; 01014 str = "+13DE3A94F1D509e320"; 01015 str >> CXSC_Rln10[4]; 01016 str = "+1DF42805E7E524e2E9"; 01017 str >> CXSC_Rln10[5]; 01018 str = "+11AAC96323250Be2B3"; 01019 str >> CXSC_Rln10[6]; 01020 str = "-1CE63884C058E4e27D"; 01021 str >> CXSC_Rln10[7]; 01022 str = "-1A1C82EA3969BAe247"; 01023 str >> CXSC_Rln10[8]; 01024 str = "+1B4F6686AD7A33e211"; 01025 str >> CXSC_Rln10[9]; 01026 str = "-1B97C8035FFC70e1DB"; 01027 str >> CXSC_Rln10[10]; 01028 str = "+1630771369962Ee1A0"; 01029 str >> CXSC_Rln10[11]; 01030 str = "-1E15BD37B295AFe16A"; 01031 str >> CXSC_Rln10[12]; 01032 str = "-132484B432318Be134"; 01033 str >> CXSC_Rln10[13]; 01034 str = "+15430212AE68C0e0FE"; 01035 str >> CXSC_Rln10[14]; 01036 str = "+1351923B322731e0C8"; 01037 str >> CXSC_Rln10[15]; 01038 str = "+11F934D794D64Fe092"; 01039 str >> CXSC_Rln10[16]; 01040 str = "+13E4B475D9FF20e05B"; 01041 str >> CXSC_Rln10[17]; 01042 str = "+185D9B63ED9A24e025"; 01043 str >> CXSC_Rln10[18]; 01044 str = "+1000035B8CA18Ce000"; 01045 str >> CXSC_Rln10[19]; 01046 str = "+1000035B8CA18De000"; 01047 str >> CXSC_Rln10[20]; 01048 01049 CXSC_Rln10_initialized = true; 01050 std::cout << RestoreOpt; 01051 } 01052 stagprec = stagmax; 01053 y = adjust(l_interval(0.0)); 01054 for (int i=0; i <= stagmax; i++) 01055 y.data[i] = CXSC_Rln10[i]; 01056 stagprec = stagsave; 01057 y = adjust(y); 01058 return y; 01059 } 01060 01061 static real CXSC_pi4[21]; 01062 static bool CXSC_pi4_initialized = false; 01063 01064 // l_interval li_pi4() throw() 01065 l_interval Pid4_l_interval() throw() 01066 { 01067 l_interval y; 01068 int stagsave = stagprec, 01069 stagmax = 20; 01070 01071 if (!CXSC_pi4_initialized) 01072 { 01073 std::string str; 01074 std::cout << SaveOpt; 01075 std::cout << Hex; 01076 str = "+1921FB54442D18e3FE"; 01077 str >> CXSC_pi4[0]; 01078 str = "+11A62633145C06e3C8"; 01079 str >> CXSC_pi4[1]; 01080 str = "+1C1CD129024E08e393"; 01081 str >> CXSC_pi4[2]; 01082 str = "+114CF98E804178e35E"; 01083 str >> CXSC_pi4[3]; 01084 str = "-159C4EC64DDAECe327"; 01085 str >> CXSC_pi4[4]; 01086 str = "+1410F31C6809BCe2F2"; 01087 str >> CXSC_pi4[5]; 01088 str = "-106AE64C32C5BCe2BB"; 01089 str >> CXSC_pi4[6]; 01090 str = "-1C99FA9EB241B4e286"; 01091 str >> CXSC_pi4[7]; 01092 str = "-1D791603D95252e24E"; 01093 str >> CXSC_pi4[8]; 01094 str = "-1571EDD0DBD254e218"; 01095 str >> CXSC_pi4[9]; 01096 str = "-133B4302721768e1E2"; 01097 str >> CXSC_pi4[10]; 01098 str = "+10BA698DFB5AC2e1AD"; 01099 str >> CXSC_pi4[11]; 01100 str = "+1FFAE5B7A035C0e178"; 01101 str >> CXSC_pi4[12]; 01102 str = "-1211C79404A576e143"; 01103 str >> CXSC_pi4[13]; 01104 str = "-1816945836FBA0e10D"; 01105 str >> CXSC_pi4[14]; 01106 str = "-1DA700CDB6BCCEe0D8"; 01107 str >> CXSC_pi4[15]; 01108 str = "+11ECE45B3DC200e0A3"; 01109 str >> CXSC_pi4[16]; 01110 str = "+1F2E2858EFC166e06D"; 01111 str >> CXSC_pi4[17]; 01112 str = "+1B4906C38ABA73e036"; 01113 str >> CXSC_pi4[18]; 01114 str = "+19A458FEA3F493e000"; 01115 str >> CXSC_pi4[19]; 01116 str = "+19A458FEA3F494e000"; 01117 str >> CXSC_pi4[20]; 01118 01119 CXSC_pi4_initialized = true; 01120 std::cout << RestoreOpt; 01121 } 01122 stagprec = stagmax; 01123 y = adjust(l_interval(0.0)); 01124 for (int i=0; i <= stagmax; i++) 01125 y.data[i] = CXSC_pi4[i]; 01126 stagprec = stagsave; 01127 y = adjust(y); 01128 return y; 01129 } 01130 01131 static real CXSC_sqrt2[21]; 01132 static bool CXSC_sqrt2_initialized = false; 01133 01134 // l_interval li_sqrt2() throw() 01135 l_interval Sqrt2_l_interval() throw() 01136 { 01137 l_interval y; 01138 int stagsave = stagprec, 01139 stagmax = 20; 01140 01141 if (!CXSC_sqrt2_initialized) 01142 { 01143 std::string str; 01144 std::cout << SaveOpt; 01145 std::cout << Hex; 01146 str = "+16A09E667F3BCDe3FF"; 01147 str >> CXSC_sqrt2[0]; 01148 str = "-1BDD3413B26456e3C9"; 01149 str >> CXSC_sqrt2[1]; 01150 str = "+157D3E3ADEC175e393"; 01151 str >> CXSC_sqrt2[2]; 01152 str = "+12775099DA2F59e35B"; 01153 str >> CXSC_sqrt2[3]; 01154 str = "+160CCE64552BF2e322"; 01155 str >> CXSC_sqrt2[4]; 01156 str = "+1821D5C5161D46e2E9"; 01157 str >> CXSC_sqrt2[5]; 01158 str = "-1C032046F8498Ee2B3"; 01159 str >> CXSC_sqrt2[6]; 01160 str = "+1EE950BC8738F7e27B"; 01161 str >> CXSC_sqrt2[7]; 01162 str = "-1AC3FDBC64E103e245"; 01163 str >> CXSC_sqrt2[8]; 01164 str = "+13B469101743A1e20D"; 01165 str >> CXSC_sqrt2[9]; 01166 str = "+15E3E9CA60B38Ce1D7"; 01167 str >> CXSC_sqrt2[10]; 01168 str = "+11BC337BCAB1BDe19C"; 01169 str >> CXSC_sqrt2[11]; 01170 str = "-1BBA5DEE9D6E7De166"; 01171 str >> CXSC_sqrt2[12]; 01172 str = "-1438DD083B1CC4e130"; 01173 str >> CXSC_sqrt2[13]; 01174 str = "+1B56A28E2EDFA7e0FA"; 01175 str >> CXSC_sqrt2[14]; 01176 str = "+1CCB2A634331F4e0C4"; 01177 str >> CXSC_sqrt2[15]; 01178 str = "-1BD9056876F83Ee08D"; 01179 str >> CXSC_sqrt2[16]; 01180 str = "-1234FA22AB6BEFe057"; 01181 str >> CXSC_sqrt2[17]; 01182 str = "+19040CA4A81395e020"; 01183 str >> CXSC_sqrt2[18]; 01184 str = "-1000002A493818e000"; 01185 str >> CXSC_sqrt2[19]; 01186 str = "-1000002A493817e000"; 01187 str >> CXSC_sqrt2[20]; 01188 01189 CXSC_sqrt2_initialized = true; 01190 std::cout << RestoreOpt; 01191 } 01192 stagprec = stagmax; 01193 y = adjust(l_interval(0.0)); 01194 for (int i=0; i <= stagmax; i++) 01195 y.data[i] = CXSC_sqrt2[i]; 01196 stagprec = stagsave; 01197 y = adjust(y); 01198 return y; 01199 } 01200 01201 static real CXSC_sqrt5[21]; // CXSC_sqrt5[0], ... CXSC_sqrt5[20] 01202 static bool CXSC_sqrt5_initialized = false; 01203 01204 l_interval Sqrt5_l_interval() throw() 01205 // Inclusion of sqrt(5), Blomquist, 30.11.2008; 01206 { 01207 l_interval y; 01208 int stagsave = stagprec, 01209 stagmax = 20; 01210 01211 if (!CXSC_sqrt5_initialized) 01212 { 01213 std::string str; 01214 std::cout << SaveOpt; 01215 std::cout << Hex; 01216 str = "+11E3779B97F4A8e400"; 01217 str >> CXSC_sqrt5[0]; 01218 str = "-1F506319FCFD19e3C9"; 01219 str >> CXSC_sqrt5[1]; 01220 str = "+1B906821044ED8e393"; 01221 str >> CXSC_sqrt5[2]; 01222 str = "-18BB1B5C0F272Ce35B"; 01223 str >> CXSC_sqrt5[3]; 01224 str = "+11D0C18E952768e324"; 01225 str >> CXSC_sqrt5[4]; 01226 str = "-1E9D585B0901F9e2EB"; 01227 str >> CXSC_sqrt5[5]; 01228 str = "-1C7DD252073EC0e2B5"; 01229 str >> CXSC_sqrt5[6]; 01230 str = "-1FCEF21EDAF7FAe27F"; 01231 str >> CXSC_sqrt5[7]; 01232 str = "+160EB25D20799Be241"; 01233 str >> CXSC_sqrt5[8]; 01234 str = "-1C90F95285168Fe208"; 01235 str >> CXSC_sqrt5[9]; 01236 str = "+1E1DFA160E75BCe1D2"; 01237 str >> CXSC_sqrt5[10]; 01238 str = "-10A08E66CB368Ce196"; 01239 str >> CXSC_sqrt5[11]; 01240 str = "+1C5371682CADD1e160"; 01241 str >> CXSC_sqrt5[12]; 01242 str = "-1998100220F4EDe129"; 01243 str >> CXSC_sqrt5[13]; 01244 str = "+1C6771A0968663e0F3"; 01245 str >> CXSC_sqrt5[14]; 01246 str = "+1DFB9E3C86CA7Ce0BD"; 01247 str >> CXSC_sqrt5[15]; 01248 str = "-18AE38ED5304B1e086"; 01249 str >> CXSC_sqrt5[16]; 01250 str = "+182A5FEC507706e050"; 01251 str >> CXSC_sqrt5[17]; 01252 str = "-1B5191A18C5647e018"; 01253 str >> CXSC_sqrt5[18]; 01254 str = "+100000000F9D52e000"; 01255 str >> CXSC_sqrt5[19]; 01256 str = "+100000000F9D53e000"; 01257 str >> CXSC_sqrt5[20]; 01258 01259 CXSC_sqrt5_initialized = true; 01260 std::cout << RestoreOpt; 01261 } 01262 stagprec = stagmax; 01263 y = adjust(l_interval(0.0)); 01264 for (int i=0; i <= stagmax; i++) 01265 y.data[i] = CXSC_sqrt5[i]; 01266 stagprec = stagsave; 01267 y = adjust(y); 01268 return y; 01269 } 01270 01271 // ************************************************************************ 01272 01273 static real CXSC_sqrt7[21]; // CXSC_sqrt7[0], ... CXSC_sqrt7[20] 01274 static bool CXSC_sqrt7_initialized = false; 01275 01276 l_interval Sqrt7_l_interval() throw() 01277 // Inclusion of sqrt(7), Blomquist, 30.11.2008; 01278 { 01279 l_interval y; 01280 int stagsave = stagprec, 01281 stagmax = 20; 01282 01283 if (!CXSC_sqrt7_initialized) 01284 { 01285 std::string str; 01286 std::cout << SaveOpt; 01287 std::cout << Hex; 01288 str = "+152A7FA9D2F8EAe400"; 01289 str >> CXSC_sqrt7[0]; 01290 str = "-121C62B033C079e3CA"; 01291 str >> CXSC_sqrt7[1]; 01292 str = "-177CAAD6200612e391"; 01293 str >> CXSC_sqrt7[2]; 01294 str = "-1EFA880DC72D64e359"; 01295 str >> CXSC_sqrt7[3]; 01296 str = "-171D206D5B1A4Ce31F"; 01297 str >> CXSC_sqrt7[4]; 01298 str = "+119392FA9B0494e2E6"; 01299 str >> CXSC_sqrt7[5]; 01300 str = "+17BB8A64890057e2AD"; 01301 str >> CXSC_sqrt7[6]; 01302 str = "-17E89300383DDEe277"; 01303 str >> CXSC_sqrt7[7]; 01304 str = "+130FB7AF68A6FBe241"; 01305 str >> CXSC_sqrt7[8]; 01306 str = "+1322281D303D36e209"; 01307 str >> CXSC_sqrt7[9]; 01308 str = "+1996109A16D3B1e1D3"; 01309 str >> CXSC_sqrt7[10]; 01310 str = "+1F239C301DFBB4e19C"; 01311 str >> CXSC_sqrt7[11]; 01312 str = "-1B5CA40AB771A2e163"; 01313 str >> CXSC_sqrt7[12]; 01314 str = "-1675711487FEAAe12A"; 01315 str >> CXSC_sqrt7[13]; 01316 str = "+122CB7FA26ABA5e0F4"; 01317 str >> CXSC_sqrt7[14]; 01318 str = "+1059211B7D5398e0BD"; 01319 str >> CXSC_sqrt7[15]; 01320 str = "-10F15BFA46EB7Fe087"; 01321 str >> CXSC_sqrt7[16]; 01322 str = "+15AB71566CE72Be051"; 01323 str >> CXSC_sqrt7[17]; 01324 str = "-1386BDCA3845C7e01A"; 01325 str >> CXSC_sqrt7[18]; 01326 str = "+10000000AC4BC7e000"; 01327 str >> CXSC_sqrt7[19]; 01328 str = "+10000000AC4BC8e000"; 01329 str >> CXSC_sqrt7[20]; 01330 01331 CXSC_sqrt7_initialized = true; 01332 std::cout << RestoreOpt; 01333 } 01334 stagprec = stagmax; 01335 y = adjust(l_interval(0.0)); 01336 for (int i=0; i <= stagmax; i++) 01337 y.data[i] = CXSC_sqrt7[i]; 01338 stagprec = stagsave; 01339 y = adjust(y); 01340 return y; 01341 } 01342 01343 01344 // ************************************************************************ 01345 01346 static real CXSC_ln2r[21]; // CXSC_ln2r[0], ... CXSC_ln2r[20] 01347 static bool CXSC_ln2r_initialized = false; 01348 01349 l_interval Ln2r_l_interval() throw() 01350 // Inclusion of 1/ln(2), Blomquist, 12.09.06; 01351 { 01352 l_interval y; 01353 int stagsave = stagprec, 01354 stagmax = 20; 01355 01356 if (!CXSC_ln2r_initialized) 01357 { 01358 std::string str; 01359 std::cout << SaveOpt; 01360 std::cout << Hex; 01361 str = "+171547652B82FEe3FF"; 01362 str >> CXSC_ln2r[0]; 01363 str = "+1777D0FFDA0D24e3C7"; 01364 str >> CXSC_ln2r[1]; 01365 str = "-160BB8A5442AB9e391"; 01366 str >> CXSC_ln2r[2]; 01367 str = "-14B52D3BA6D74De359"; 01368 str >> CXSC_ln2r[3]; 01369 str = "+19A342648FBC39e323"; 01370 str >> CXSC_ln2r[4]; 01371 str = "-1E0455744994EEe2ED"; 01372 str >> CXSC_ln2r[5]; 01373 str = "+1B25EEB82D7C16e2B7"; 01374 str >> CXSC_ln2r[6]; 01375 str = "+1F5485CF306255e281"; 01376 str >> CXSC_ln2r[7]; 01377 str = "-1EC07680A1F958e24B"; 01378 str >> CXSC_ln2r[8]; 01379 str = "-106326680EB5B6e215"; 01380 str >> CXSC_ln2r[9]; 01381 str = "-1B3D04C549BC98e1DF"; 01382 str >> CXSC_ln2r[10]; 01383 str = "+1EABCEAD10305Be1A9"; 01384 str >> CXSC_ln2r[11]; 01385 str = "-14440C57D7AB97e170"; 01386 str >> CXSC_ln2r[12]; 01387 str = "-17185D42A4E6D6e139"; 01388 str >> CXSC_ln2r[13]; 01389 str = "-1F332B5BE48526e101"; 01390 str >> CXSC_ln2r[14]; 01391 str = "+12CE4F199E108De0CB"; 01392 str >> CXSC_ln2r[15]; 01393 str = "-18DAFCC6077F2Ae092"; 01394 str >> CXSC_ln2r[16]; 01395 str = "+19ABB71EC25E12e05B"; 01396 str >> CXSC_ln2r[17]; 01397 str = "-11473D7A3366BDe022"; 01398 str >> CXSC_ln2r[18]; 01399 str = "-1000004977D38Be000"; 01400 str >> CXSC_ln2r[19]; 01401 str = "-1000004977D38Ae000"; 01402 str >> CXSC_ln2r[20]; 01403 01404 CXSC_ln2r_initialized = true; 01405 std::cout << RestoreOpt; 01406 } 01407 stagprec = stagmax; 01408 y = adjust(l_interval(0.0)); 01409 for (int i=0; i <= stagmax; i++) 01410 y.data[i] = CXSC_ln2r[i]; 01411 // y[i+1] = CXSC_ln2r[i]; 01412 stagprec = stagsave; 01413 y = adjust(y); 01414 return y; 01415 } 01416 01417 01418 // ********************************************************************** 01419 01420 static real CXSC_Pi[21]; // CXSC_Pi[0], ... CXSC_Pi[20] 01421 static bool CXSC_Pi_initialized = false; 01422 01423 l_interval Pi_l_interval() throw() 01424 // Inclusion of Pi, Blomquist, 12.09.06; 01425 { 01426 l_interval y; 01427 int stagsave = stagprec, 01428 stagmax = 20; 01429 01430 if (!CXSC_Pi_initialized) 01431 { 01432 std::string str; 01433 std::cout << SaveOpt; 01434 std::cout << Hex; 01435 str = "+1921FB54442D18e400"; 01436 str >> CXSC_Pi[0]; 01437 str = "+11A62633145C07e3CA"; 01438 str >> CXSC_Pi[1]; 01439 str = "-1F1976B7ED8FBCe392"; 01440 str >> CXSC_Pi[2]; 01441 str = "+14CF98E804177De35C"; 01442 str >> CXSC_Pi[3]; 01443 str = "+131D89CD9128A5e326"; 01444 str >> CXSC_Pi[4]; 01445 str = "+10F31C6809BBDFe2EC"; 01446 str >> CXSC_Pi[5]; 01447 str = "+1519B3CD3A431Be2B5"; 01448 str >> CXSC_Pi[6]; 01449 str = "+18158536F92F8Ae27E"; 01450 str >> CXSC_Pi[7]; 01451 str = "+1BA7F09AB6B6A9e246"; 01452 str >> CXSC_Pi[8]; 01453 str = "-1EDD0DBD2544CFe20E"; 01454 str >> CXSC_Pi[9]; 01455 str = "+179FB1BD1310BAe1D7"; 01456 str >> CXSC_Pi[10]; 01457 str = "+1A637ED6B0BFF6e1A1"; 01458 str >> CXSC_Pi[11]; 01459 str = "-1A485FCA40908Ee16A"; 01460 str >> CXSC_Pi[12]; 01461 str = "-1E501295D98169e133"; 01462 str >> CXSC_Pi[13]; 01463 str = "-1160DBEE83B4E0e0FD"; 01464 str >> CXSC_Pi[14]; 01465 str = "-19B6D799AE131Ce0C5"; 01466 str >> CXSC_Pi[15]; 01467 str = "+16CF70801F2E28e08F"; 01468 str >> CXSC_Pi[16]; 01469 str = "+163BF0598DA483e059"; 01470 str >> CXSC_Pi[17]; 01471 str = "+1871574E69A459e023"; 01472 str >> CXSC_Pi[18]; 01473 str = "-10000005702DB4e000"; 01474 str >> CXSC_Pi[19]; 01475 str = "-10000005702DB3e000"; 01476 str >> CXSC_Pi[20]; 01477 01478 CXSC_Pi_initialized = true; 01479 std::cout << RestoreOpt; 01480 } 01481 stagprec = stagmax; 01482 y = adjust(l_interval(0.0)); 01483 for (int i=0; i <= stagmax; i++) 01484 y.data[i] = CXSC_Pi[i]; 01485 // y[i+1] = CXSC_Pi[i]; 01486 stagprec = stagsave; 01487 y = adjust(y); 01488 return y; 01489 } 01490 01491 // ********************************************************************** 01492 01493 static real CXSC_Pid2[21]; // CXSC_Pid2[0], ... CXSC_Pid2[20] 01494 static bool CXSC_Pid2_initialized = false; 01495 01496 l_interval Pid2_l_interval() throw() 01497 // Inclusion of Pi/2, Blomquist, 12.09.06; 01498 { 01499 l_interval y; 01500 int stagsave = stagprec, 01501 stagmax = 20; 01502 01503 if (!CXSC_Pid2_initialized) 01504 { 01505 std::string str; 01506 std::cout << SaveOpt; 01507 std::cout << Hex; 01508 str = "+1921FB54442D18e3FF"; 01509 str >> CXSC_Pid2[0]; 01510 str = "+11A62633145C07e3C9"; 01511 str >> CXSC_Pid2[1]; 01512 str = "-1F1976B7ED8FBCe391"; 01513 str >> CXSC_Pid2[2]; 01514 str = "+14CF98E804177De35B"; 01515 str >> CXSC_Pid2[3]; 01516 str = "+131D89CD9128A5e325"; 01517 str >> CXSC_Pid2[4]; 01518 str = "+10F31C6809BBDFe2EB"; 01519 str >> CXSC_Pid2[5]; 01520 str = "+1519B3CD3A431Be2B4"; 01521 str >> CXSC_Pid2[6]; 01522 str = "+18158536F92F8Ae27D"; 01523 str >> CXSC_Pid2[7]; 01524 str = "+1BA7F09AB6B6A9e245"; 01525 str >> CXSC_Pid2[8]; 01526 str = "-1EDD0DBD2544CFe20D"; 01527 str >> CXSC_Pid2[9]; 01528 str = "+179FB1BD1310BAe1D6"; 01529 str >> CXSC_Pid2[10]; 01530 str = "+1A637ED6B0BFF6e1A0"; 01531 str >> CXSC_Pid2[11]; 01532 str = "-1A485FCA40908Ee169"; 01533 str >> CXSC_Pid2[12]; 01534 str = "-1E501295D98169e132"; 01535 str >> CXSC_Pid2[13]; 01536 str = "-1160DBEE83B4E0e0FC"; 01537 str >> CXSC_Pid2[14]; 01538 str = "-19B6D799AE131Ce0C4"; 01539 str >> CXSC_Pid2[15]; 01540 str = "+16CF70801F2E28e08E"; 01541 str >> CXSC_Pid2[16]; 01542 str = "+163BF0598DA483e058"; 01543 str >> CXSC_Pid2[17]; 01544 str = "+1871574E69A459e022"; 01545 str >> CXSC_Pid2[18]; 01546 str = "-10000002B816DAe000"; 01547 str >> CXSC_Pid2[19]; 01548 str = "-10000002B816D0e000"; 01549 str >> CXSC_Pid2[20]; 01550 01551 CXSC_Pid2_initialized = true; 01552 std::cout << RestoreOpt; 01553 } 01554 stagprec = stagmax; 01555 y = adjust(l_interval(0.0)); 01556 for (int i=0; i <= stagmax; i++) 01557 y.data[i] = CXSC_Pid2[i]; 01558 // y[i+1] = CXSC_Pid2[i]; 01559 stagprec = stagsave; 01560 y = adjust(y); 01561 return y; 01562 } 01563 01564 // ********************************************************************** 01565 01566 static real CXSC_Pi2[21]; // CXSC_Pi2[0], ... CXSC_Pi2[20] 01567 static bool CXSC_Pi2_initialized = false; 01568 01569 l_interval Pi2_l_interval() throw() 01570 // Inclusion of 2*Pi, Blomquist, 12.09.06; 01571 { 01572 l_interval y; 01573 int stagsave = stagprec, 01574 stagmax = 20; 01575 01576 if (!CXSC_Pi2_initialized) 01577 { 01578 std::string str; 01579 std::cout << SaveOpt; 01580 std::cout << Hex; 01581 str = "+1921FB54442D18e401"; 01582 str >> CXSC_Pi2[0]; 01583 str = "+11A62633145C07e3CB"; 01584 str >> CXSC_Pi2[1]; 01585 str = "-1F1976B7ED8FBCe393"; 01586 str >> CXSC_Pi2[2]; 01587 str = "+14CF98E804177De35D"; 01588 str >> CXSC_Pi2[3]; 01589 str = "+131D89CD9128A5e327"; 01590 str >> CXSC_Pi2[4]; 01591 str = "+10F31C6809BBDFe2ED"; 01592 str >> CXSC_Pi2[5]; 01593 str = "+1519B3CD3A431Be2B6"; 01594 str >> CXSC_Pi2[6]; 01595 str = "+18158536F92F8Ae27F"; 01596 str >> CXSC_Pi2[7]; 01597 str = "+1BA7F09AB6B6A9e247"; 01598 str >> CXSC_Pi2[8]; 01599 str = "-1EDD0DBD2544CFe20F"; 01600 str >> CXSC_Pi2[9]; 01601 str = "+179FB1BD1310BAe1D8"; 01602 str >> CXSC_Pi2[10]; 01603 str = "+1A637ED6B0BFF6e1A2"; 01604 str >> CXSC_Pi2[11]; 01605 str = "-1A485FCA40908Ee16B"; 01606 str >> CXSC_Pi2[12]; 01607 str = "-1E501295D98169e134"; 01608 str >> CXSC_Pi2[13]; 01609 str = "-1160DBEE83B4E0e0FE"; 01610 str >> CXSC_Pi2[14]; 01611 str = "-19B6D799AE131Ce0C6"; 01612 str >> CXSC_Pi2[15]; 01613 str = "+16CF70801F2E28e090"; 01614 str >> CXSC_Pi2[16]; 01615 str = "+163BF0598DA483e05A"; 01616 str >> CXSC_Pi2[17]; 01617 str = "+1871574E69A459e024"; 01618 str >> CXSC_Pi2[18]; 01619 str = "-1000000AE05B67e000"; 01620 str >> CXSC_Pi2[19]; 01621 str = "-1000000AE05B66e000"; 01622 str >> CXSC_Pi2[20]; 01623 01624 CXSC_Pi2_initialized = true; 01625 std::cout << RestoreOpt; 01626 } 01627 stagprec = stagmax; 01628 y = adjust(l_interval(0.0)); 01629 for (int i=0; i <= stagmax; i++) 01630 y.data[i] = CXSC_Pi2[i]; 01631 // y[i+1] = CXSC_Pi2[i]; 01632 stagprec = stagsave; 01633 y = adjust(y); 01634 return y; 01635 } 01636 01637 // ********************************************************************** 01638 01639 static real CXSC_Pid3[21]; // CXSC_Pid3[0], ... CXSC_Pid3[20] 01640 static bool CXSC_Pid3_initialized = false; 01641 01642 l_interval Pid3_l_interval() throw() 01643 // Inclusion of Pi/3, Blomquist, 12.09.06; 01644 { 01645 l_interval y; 01646 int stagsave = stagprec, 01647 stagmax = 20; 01648 01649 if (!CXSC_Pid3_initialized) 01650 { 01651 std::string str; 01652 std::cout << SaveOpt; 01653 std::cout << Hex; 01654 str = "+10C152382D7366e3FF"; 01655 str >> CXSC_Pid3[0]; 01656 str = "-1EE6913347C2A6e3C9"; 01657 str >> CXSC_Pid3[1]; 01658 str = "-14BBA47A9E5FD2e391"; 01659 str >> CXSC_Pid3[2]; 01660 str = "-1CCAEF65529B02e35B"; 01661 str >> CXSC_Pid3[3]; 01662 str = "+197CB7BCC18B87e324"; 01663 str >> CXSC_Pid3[4]; 01664 str = "-13EBBDA1FF3058e2EE"; 01665 str >> CXSC_Pid3[5]; 01666 str = "-11D10CB320F4D1e2B6"; 01667 str >> CXSC_Pid3[6]; 01668 str = "+1958EB892987ECe27F"; 01669 str >> CXSC_Pid3[7]; 01670 str = "+167C54B11CF247e249"; 01671 str >> CXSC_Pid3[8]; 01672 str = "+12C2E985923A44e210"; 01673 str >> CXSC_Pid3[9]; 01674 str = "+1945484A2DD81Fe1D8"; 01675 str >> CXSC_Pid3[10]; 01676 str = "+1197A9E475D54Fe1A0"; 01677 str >> CXSC_Pid3[11]; 01678 str = "-1E181FEE158585e16A"; 01679 str >> CXSC_Pid3[12]; 01680 str = "+1047FCE7066A6Ee134"; 01681 str >> CXSC_Pid3[13]; 01682 str = "+1D1A8602EA0C85e0FE"; 01683 str >> CXSC_Pid3[14]; 01684 str = "+14430C5998BF34e0C8"; 01685 str >> CXSC_Pid3[15]; 01686 str = "+173BF40AAD43D9e091"; 01687 str >> CXSC_Pid3[16]; 01688 str = "-137B014DDEDCF5e05B"; 01689 str >> CXSC_Pid3[17]; 01690 str = "-1A5F1B210EE7C5e022"; 01691 str >> CXSC_Pid3[18]; 01692 str = "+100000A8DA9B6Ee000"; 01693 str >> CXSC_Pid3[19]; 01694 str = "+100000A8DA9B6Fe000"; 01695 str >> CXSC_Pid3[20]; 01696 01697 CXSC_Pid3_initialized = true; 01698 std::cout << RestoreOpt; 01699 } 01700 stagprec = stagmax; 01701 y = adjust(l_interval(0.0)); 01702 for (int i=0; i <= stagmax; i++) 01703 y.data[i] = CXSC_Pid3[i]; 01704 // y[i+1] = CXSC_Pid3[i]; 01705 stagprec = stagsave; 01706 y = adjust(y); 01707 return y; 01708 } 01709 01710 01711 // ********************************************************************** 01712 01713 static real CXSC_Pir[21]; // CXSC_Pir[0], ... CXSC_Pir[20] 01714 static bool CXSC_Pir_initialized = false; 01715 01716 l_interval Pir_l_interval() throw() 01717 // Inclusion of 1/Pi, Blomquist, 12.09.06; 01718 { 01719 l_interval y; 01720 int stagsave = stagprec, 01721 stagmax = 20; 01722 01723 if (!CXSC_Pir_initialized) 01724 { 01725 std::string str; 01726 std::cout << SaveOpt; 01727 std::cout << Hex; 01728 str = "+145F306DC9C883e3FD"; 01729 str >> CXSC_Pir[0]; 01730 str = "-16B01EC5417056e3C7"; 01731 str >> CXSC_Pir[1]; 01732 str = "-16447E493AD4CEe391"; 01733 str >> CXSC_Pir[2]; 01734 str = "+1E21C820FF28B2e35B"; 01735 str >> CXSC_Pir[3]; 01736 str = "-1508510EA79237e324"; 01737 str >> CXSC_Pir[4]; 01738 str = "+1B8E909374B802e2EC"; 01739 str >> CXSC_Pir[5]; 01740 str = "-1B6D115F62E6DEe2B6"; 01741 str >> CXSC_Pir[6]; 01742 str = "-180F10A71A76B3e27F"; 01743 str >> CXSC_Pir[7]; 01744 str = "+1CFBA208D7D4BBe248"; 01745 str >> CXSC_Pir[8]; 01746 str = "-12EDEC598E3F65e210"; 01747 str >> CXSC_Pir[9]; 01748 str = "-1741037D8CDC54e1D9"; 01749 str >> CXSC_Pir[10]; 01750 str = "+1CC1A99CFA4E42e1A3"; 01751 str >> CXSC_Pir[11]; 01752 str = "+17E2EF7E4A0EC8e16C"; 01753 str >> CXSC_Pir[12]; 01754 str = "-1DA00087E99FC0e130"; 01755 str >> CXSC_Pir[13]; 01756 str = "-10D0EE74A5F593e0FA"; 01757 str >> CXSC_Pir[14]; 01758 str = "+1F6D367ECF27CBe0C2"; 01759 str >> CXSC_Pir[15]; 01760 str = "+136E9E8C7ECD3De089"; 01761 str >> CXSC_Pir[16]; 01762 str = "-100AE9456C229Ce053"; 01763 str >> CXSC_Pir[17]; 01764 str = "-141A0E84C2F8C6e01A"; 01765 str >> CXSC_Pir[18]; 01766 str = "-1000000010EB5Be000"; 01767 str >> CXSC_Pir[19]; 01768 str = "-1000000010EB5Ae000"; 01769 str >> CXSC_Pir[20]; 01770 01771 CXSC_Pir_initialized = true; 01772 std::cout << RestoreOpt; 01773 } 01774 stagprec = stagmax; 01775 y = adjust(l_interval(0.0)); 01776 for (int i=0; i <= stagmax; i++) 01777 y.data[i] = CXSC_Pir[i]; 01778 // y[i+1] = CXSC_Pir[i]; 01779 stagprec = stagsave; 01780 y = adjust(y); 01781 return y; 01782 } 01783 01784 // ********************************************************************** 01785 01786 static real CXSC_Pi2r[21]; // CXSC_Pi2r[0], ... CXSC_Pi2r[20] 01787 static bool CXSC_Pi2r_initialized = false; 01788 01789 l_interval Pi2r_l_interval() throw() 01790 // Inclusion of 1/(2*Pi), Blomquist, 12.09.06; 01791 { 01792 l_interval y; 01793 int stagsave = stagprec, 01794 stagmax = 20; 01795 01796 if (!CXSC_Pi2r_initialized) 01797 { 01798 std::string str; 01799 std::cout << SaveOpt; 01800 std::cout << Hex; 01801 str = "+145F306DC9C883e3FC"; 01802 str >> CXSC_Pi2r[0]; 01803 str = "-16B01EC5417056e3C6"; 01804 str >> CXSC_Pi2r[1]; 01805 str = "-16447E493AD4CEe390"; 01806 str >> CXSC_Pi2r[2]; 01807 str = "+1E21C820FF28B2e35A"; 01808 str >> CXSC_Pi2r[3]; 01809 str = "-1508510EA79237e323"; 01810 str >> CXSC_Pi2r[4]; 01811 str = "+1B8E909374B802e2EB"; 01812 str >> CXSC_Pi2r[5]; 01813 str = "-1B6D115F62E6DEe2B5"; 01814 str >> CXSC_Pi2r[6]; 01815 str = "-180F10A71A76B3e27E"; 01816 str >> CXSC_Pi2r[7]; 01817 str = "+1CFBA208D7D4BBe247"; 01818 str >> CXSC_Pi2r[8]; 01819 str = "-12EDEC598E3F65e20F"; 01820 str >> CXSC_Pi2r[9]; 01821 str = "-1741037D8CDC54e1D8"; 01822 str >> CXSC_Pi2r[10]; 01823 str = "+1CC1A99CFA4E42e1A2"; 01824 str >> CXSC_Pi2r[11]; 01825 str = "+17E2EF7E4A0EC8e16B"; 01826 str >> CXSC_Pi2r[12]; 01827 str = "-1DA00087E99FC0e12F"; 01828 str >> CXSC_Pi2r[13]; 01829 str = "-10D0EE74A5F593e0F9"; 01830 str >> CXSC_Pi2r[14]; 01831 str = "+1F6D367ECF27CBe0C1"; 01832 str >> CXSC_Pi2r[15]; 01833 str = "+136E9E8C7ECD3De088"; 01834 str >> CXSC_Pi2r[16]; 01835 str = "-100AE9456C229Ce052"; 01836 str >> CXSC_Pi2r[17]; 01837 str = "-141A0E84C2F8C6e019"; 01838 str >> CXSC_Pi2r[18]; 01839 str = "-100000000875AEe000"; 01840 str >> CXSC_Pi2r[19]; 01841 str = "-100000000875ADe000"; 01842 str >> CXSC_Pi2r[20]; 01843 01844 CXSC_Pi2r_initialized = true; 01845 std::cout << RestoreOpt; 01846 } 01847 stagprec = stagmax; 01848 y = adjust(l_interval(0.0)); 01849 for (int i=0; i <= stagmax; i++) 01850 y.data[i] = CXSC_Pi2r[i]; 01851 // y[i+1] = CXSC_Pi2r[i]; 01852 stagprec = stagsave; 01853 y = adjust(y); 01854 return y; 01855 } 01856 01857 01858 // ********************************************************************** 01859 01860 static real CXSC_SqrtPi[21]; // CXSC_SqrtPi[0], ... CXSC_SqrtPi[20] 01861 static bool CXSC_SqrtPi_initialized = false; 01862 01863 l_interval SqrtPi_l_interval() throw() 01864 // Inclusion of Sqrt(Pi), Blomquist, 12.09.06; 01865 { 01866 l_interval y; 01867 int stagsave = stagprec, 01868 stagmax = 20; 01869 01870 if (!CXSC_SqrtPi_initialized) 01871 { 01872 std::string str; 01873 std::cout << SaveOpt; 01874 std::cout << Hex; 01875 str = "+1C5BF891B4EF6Be3FF"; 01876 str >> CXSC_SqrtPi[0]; 01877 str = "-1618F13EB7CA89e3C9"; 01878 str >> CXSC_SqrtPi[1]; 01879 str = "-1B1F0071B7AAE4e391"; 01880 str >> CXSC_SqrtPi[2]; 01881 str = "-1389B5A46BDFE8e35A"; 01882 str >> CXSC_SqrtPi[3]; 01883 str = "-160AF5C5C89448e324"; 01884 str >> CXSC_SqrtPi[4]; 01885 str = "-14835F07122994e2E8"; 01886 str >> CXSC_SqrtPi[5]; 01887 str = "+1CEC283C18EE8Fe2B2"; 01888 str >> CXSC_SqrtPi[6]; 01889 str = "-13ADEBB9223CA8e27B"; 01890 str >> CXSC_SqrtPi[7]; 01891 str = "+1454912430D291e245"; 01892 str >> CXSC_SqrtPi[8]; 01893 str = "-1E8B2345020EF6e20F"; 01894 str >> CXSC_SqrtPi[9]; 01895 str = "-17262982556291e1D8"; 01896 str >> CXSC_SqrtPi[10]; 01897 str = "+1196FA9B140CABe1A1"; 01898 str >> CXSC_SqrtPi[11]; 01899 str = "-175EEE59D91D39e16B"; 01900 str >> CXSC_SqrtPi[12]; 01901 str = "+1789268B7D9D48e130"; 01902 str >> CXSC_SqrtPi[13]; 01903 str = "+17162E2F06B89Ce0FA"; 01904 str >> CXSC_SqrtPi[14]; 01905 str = "+1EC9C08F40A3DBe0C3"; 01906 str >> CXSC_SqrtPi[15]; 01907 str = "+1B6048DD0729E2e08D"; 01908 str >> CXSC_SqrtPi[16]; 01909 str = "+1471CF4C33FF6Be056"; 01910 str >> CXSC_SqrtPi[17]; 01911 str = "+1D75FBD8B36F94e020"; 01912 str >> CXSC_SqrtPi[18]; 01913 str = "+1000002D74B3A2e000"; 01914 str >> CXSC_SqrtPi[19]; 01915 str = "+1000002D74B3A3e000"; 01916 str >> CXSC_SqrtPi[20]; 01917 01918 CXSC_SqrtPi_initialized = true; 01919 std::cout << RestoreOpt; 01920 } 01921 stagprec = stagmax; 01922 y = adjust(l_interval(0.0)); 01923 for (int i=0; i <= stagmax; i++) 01924 y.data[i] = CXSC_SqrtPi[i]; 01925 // y[i+1] = CXSC_SqrtPi[i]; 01926 stagprec = stagsave; 01927 y = adjust(y); 01928 return y; 01929 } 01930 01931 01932 // ********************************************************************** 01933 01934 static real CXSC_Sqrt2Pi[21]; // CXSC_Sqrt2Pi[0], ... CXSC_Sqrt2Pi[20] 01935 static bool CXSC_Sqrt2Pi_initialized = false; 01936 01937 l_interval Sqrt2Pi_l_interval() throw() 01938 // Inclusion of sqrt(2*Pi), Blomquist, 12.09.06; 01939 { 01940 l_interval y; 01941 int stagsave = stagprec, 01942 stagmax = 20; 01943 01944 if (!CXSC_Sqrt2Pi_initialized) 01945 { 01946 std::string str; 01947 std::cout << SaveOpt; 01948 std::cout << Hex; 01949 str = "+140D931FF62706e400"; 01950 str >> CXSC_Sqrt2Pi[0]; 01951 str = "-1A6A0D6F814637e3CA"; 01952 str >> CXSC_Sqrt2Pi[1]; 01953 str = "-1311D073060ACEe394"; 01954 str >> CXSC_Sqrt2Pi[2]; 01955 str = "+16000B50DC2F41e35B"; 01956 str >> CXSC_Sqrt2Pi[3]; 01957 str = "+16EF75CA45A834e324"; 01958 str >> CXSC_Sqrt2Pi[4]; 01959 str = "+19BDB2B4C39342e2EC"; 01960 str >> CXSC_Sqrt2Pi[5]; 01961 str = "+1F5582E2063EE6e2B5"; 01962 str >> CXSC_Sqrt2Pi[6]; 01963 str = "+183F879BEA150Ce27C"; 01964 str >> CXSC_Sqrt2Pi[7]; 01965 str = "-1F1EA3CA289B00e244"; 01966 str >> CXSC_Sqrt2Pi[8]; 01967 str = "-1699CDA77736F9e20D"; 01968 str >> CXSC_Sqrt2Pi[9]; 01969 str = "-11A379D298B55Ee1D4"; 01970 str >> CXSC_Sqrt2Pi[10]; 01971 str = "-1A6DDB0152BA94e19E"; 01972 str >> CXSC_Sqrt2Pi[11]; 01973 str = "-1957E2E58A02FEe167"; 01974 str >> CXSC_Sqrt2Pi[12]; 01975 str = "-1D6160F18E604De131"; 01976 str >> CXSC_Sqrt2Pi[13]; 01977 str = "+1311860CDF7215e0F8"; 01978 str >> CXSC_Sqrt2Pi[14]; 01979 str = "+12271F44C50274e0C1"; 01980 str >> CXSC_Sqrt2Pi[15]; 01981 str = "-100BF5C5497A21e08A"; 01982 str >> CXSC_Sqrt2Pi[16]; 01983 str = "+1E94B6E6AD51E2e052"; 01984 str >> CXSC_Sqrt2Pi[17]; 01985 str = "-1C910B5F3D27CEe019"; 01986 str >> CXSC_Sqrt2Pi[18]; 01987 str = "+100000007C99B0e000"; 01988 str >> CXSC_Sqrt2Pi[19]; 01989 str = "+100000007C99B1e000"; 01990 str >> CXSC_Sqrt2Pi[20]; 01991 01992 CXSC_Sqrt2Pi_initialized = true; 01993 std::cout << RestoreOpt; 01994 } 01995 stagprec = stagmax; 01996 y = adjust(l_interval(0.0)); 01997 for (int i=0; i <= stagmax; i++) 01998 y.data[i] = CXSC_Sqrt2Pi[i]; 01999 // y[i+1] = CXSC_Sqrt2Pi[i]; 02000 stagprec = stagsave; 02001 y = adjust(y); 02002 return y; 02003 } 02004 02005 02006 // ********************************************************************** 02007 02008 static real CXSC_SqrtPir[21]; // CXSC_SqrtPir[0], ... CXSC_SqrtPir[20] 02009 static bool CXSC_SqrtPir_initialized = false; 02010 02011 l_interval SqrtPir_l_interval() throw() 02012 // Inclusion of 1/sqrt(Pi), Blomquist, 12.09.06; 02013 { 02014 l_interval y; 02015 int stagsave = stagprec, 02016 stagmax = 20; 02017 02018 if (!CXSC_SqrtPir_initialized) 02019 { 02020 std::string str; 02021 std::cout << SaveOpt; 02022 std::cout << Hex; 02023 str = "+120DD750429B6De3FE"; 02024 str >> CXSC_SqrtPir[0]; 02025 str = "+11AE3A914FED80e3C6"; 02026 str >> CXSC_SqrtPir[1]; 02027 str = "-13CBBEBF65F145e38F"; 02028 str >> CXSC_SqrtPir[2]; 02029 str = "-1E0C574632F53Ee358"; 02030 str >> CXSC_SqrtPir[3]; 02031 str = "-1E6633BE9E7F15e322"; 02032 str >> CXSC_SqrtPir[4]; 02033 str = "+1CF859270F1141e2EB"; 02034 str >> CXSC_SqrtPir[5]; 02035 str = "-1FE4FB499C328Ae2B4"; 02036 str >> CXSC_SqrtPir[6]; 02037 str = "-10B82C446DC78De27D"; 02038 str >> CXSC_SqrtPir[7]; 02039 str = "-1878B089078800e247"; 02040 str >> CXSC_SqrtPir[8]; 02041 str = "-13DAEADA9E233Ee20F"; 02042 str >> CXSC_SqrtPir[9]; 02043 str = "+1137197A708BD2e1D9"; 02044 str >> CXSC_SqrtPir[10]; 02045 str = "-109009506D5BA2e19E"; 02046 str >> CXSC_SqrtPir[11]; 02047 str = "+17C9F0B5951E94e168"; 02048 str >> CXSC_SqrtPir[12]; 02049 str = "-1735F4949633A4e131"; 02050 str >> CXSC_SqrtPir[13]; 02051 str = "-146014DBC90D0Ee0FB"; 02052 str >> CXSC_SqrtPir[14]; 02053 str = "+1CAB0B222EEEA0e0C5"; 02054 str >> CXSC_SqrtPir[15]; 02055 str = "+1B1C750754B40Ae08F"; 02056 str >> CXSC_SqrtPir[16]; 02057 str = "-16B2CD2E72C16Ee057"; 02058 str >> CXSC_SqrtPir[17]; 02059 str = "-148C024FF194B2e021"; 02060 str >> CXSC_SqrtPir[18]; 02061 str = "+10000073E19B74e000"; 02062 str >> CXSC_SqrtPir[19]; 02063 str = "+10000073E19B75e000"; 02064 str >> CXSC_SqrtPir[20]; 02065 02066 CXSC_SqrtPir_initialized = true; 02067 std::cout << RestoreOpt; 02068 } 02069 stagprec = stagmax; 02070 y = adjust(l_interval(0.0)); 02071 for (int i=0; i <= stagmax; i++) 02072 y.data[i] = CXSC_SqrtPir[i]; 02073 // y[i+1] = CXSC_SqrtPir[i]; 02074 stagprec = stagsave; 02075 y = adjust(y); 02076 return y; 02077 } 02078 02079 02080 // ********************************************************************** 02081 02082 static real CXSC_Sqrt2Pir[21]; // CXSC_Sqrt2Pir[0], ... CXSC_Sqrt2Pir[20] 02083 static bool CXSC_Sqrt2Pir_initialized = false; 02084 02085 l_interval Sqrt2Pir_l_interval() throw() 02086 // Inclusion of 1/sqrt(2*Pi), Blomquist, 12.09.06; 02087 { 02088 l_interval y; 02089 int stagsave = stagprec, 02090 stagmax = 20; 02091 02092 if (!CXSC_Sqrt2Pir_initialized) 02093 { 02094 std::string str; 02095 std::cout << SaveOpt; 02096 std::cout << Hex; 02097 str = "+19884533D43651e3FD"; 02098 str >> CXSC_Sqrt2Pir[0]; 02099 str = "-1CBC0D30EBFD15e3C7"; 02100 str >> CXSC_Sqrt2Pir[1]; 02101 str = "-1C7402C7D60CFBe38F"; 02102 str >> CXSC_Sqrt2Pir[2]; 02103 str = "+12706D8C0471B5e357"; 02104 str >> CXSC_Sqrt2Pir[3]; 02105 str = "-1FF6718B45881De321"; 02106 str >> CXSC_Sqrt2Pir[4]; 02107 str = "-13AABB82C248DCe2EB"; 02108 str >> CXSC_Sqrt2Pir[5]; 02109 str = "-1458A899162EE4e2B2"; 02110 str >> CXSC_Sqrt2Pir[6]; 02111 str = "-14EBD8868F41EBe27B"; 02112 str >> CXSC_Sqrt2Pir[7]; 02113 str = "+13278E993445F1e243"; 02114 str >> CXSC_Sqrt2Pir[8]; 02115 str = "-1CC019F5F4780Ae20D"; 02116 str >> CXSC_Sqrt2Pir[9]; 02117 str = "+147CE4B4ECDBD7e1D7"; 02118 str >> CXSC_Sqrt2Pir[10]; 02119 str = "-19A3DCC6A3534Be19F"; 02120 str >> CXSC_Sqrt2Pir[11]; 02121 str = "+11379A7BA8CB0Ae169"; 02122 str >> CXSC_Sqrt2Pir[12]; 02123 str = "-12D909C875312Ee132"; 02124 str >> CXSC_Sqrt2Pir[13]; 02125 str = "+1C1CEC4882C77Be0FB"; 02126 str >> CXSC_Sqrt2Pir[14]; 02127 str = "-14C4078263DF36e0C5"; 02128 str >> CXSC_Sqrt2Pir[15]; 02129 str = "+1AB3FC8D2AB243e08F"; 02130 str >> CXSC_Sqrt2Pir[16]; 02131 str = "+17B9172454310Ae059"; 02132 str >> CXSC_Sqrt2Pir[17]; 02133 str = "-1444B6B781B7F2e023"; 02134 str >> CXSC_Sqrt2Pir[18]; 02135 str = "-100001DB5C6774e000"; 02136 str >> CXSC_Sqrt2Pir[19]; 02137 str = "-100001DB5C6773e000"; 02138 str >> CXSC_Sqrt2Pir[20]; 02139 02140 CXSC_Sqrt2Pir_initialized = true; 02141 std::cout << RestoreOpt; 02142 } 02143 stagprec = stagmax; 02144 y = adjust(l_interval(0.0)); 02145 for (int i=0; i <= stagmax; i++) 02146 y.data[i] = CXSC_Sqrt2Pir[i]; 02147 // y[i+1] = CXSC_Sqrt2Pir[i]; 02148 stagprec = stagsave; 02149 y = adjust(y); 02150 return y; 02151 } 02152 02153 02154 // ********************************************************************** 02155 02156 static real CXSC_Pip2[21]; // CXSC_Pip2[0], ... CXSC_Pip2[20] 02157 static bool CXSC_Pip2_initialized = false; 02158 02159 l_interval Pip2_l_interval() throw() 02160 // Inclusion of Pi^2, Blomquist, 12.09.06; 02161 { 02162 l_interval y; 02163 int stagsave = stagprec, 02164 stagmax = 20; 02165 02166 if (!CXSC_Pip2_initialized) 02167 { 02168 std::string str; 02169 std::cout << SaveOpt; 02170 std::cout << Hex; 02171 str = "+13BD3CC9BE45DEe402"; 02172 str >> CXSC_Pip2[0]; 02173 str = "+1692B71366CC04e3CC"; 02174 str >> CXSC_Pip2[1]; 02175 str = "+18358E10ACD480e396"; 02176 str >> CXSC_Pip2[2]; 02177 str = "-1F2F5DD7997DDFe35F"; 02178 str >> CXSC_Pip2[3]; 02179 str = "+129E39B47B884Ee324"; 02180 str >> CXSC_Pip2[4]; 02181 str = "-12CF7459DD5DAFe2EE"; 02182 str >> CXSC_Pip2[5]; 02183 str = "-11842F87B5FE0Fe2B8"; 02184 str >> CXSC_Pip2[6]; 02185 str = "+1FFD8A79616A21e282"; 02186 str >> CXSC_Pip2[7]; 02187 str = "+12492A6663E899e24C"; 02188 str >> CXSC_Pip2[8]; 02189 str = "-1A15F4352CC511e215"; 02190 str >> CXSC_Pip2[9]; 02191 str = "-1301AA1792FF3Ce1DE"; 02192 str >> CXSC_Pip2[10]; 02193 str = "+122B6F31626EFEe1A8"; 02194 str >> CXSC_Pip2[11]; 02195 str = "+1B317FA13BDD8Fe172"; 02196 str >> CXSC_Pip2[12]; 02197 str = "+16F83B49040075e13C"; 02198 str >> CXSC_Pip2[13]; 02199 str = "-1B1890A945FE17e106"; 02200 str >> CXSC_Pip2[14]; 02201 str = "+12DCD389B96CDBe0D0"; 02202 str >> CXSC_Pip2[15]; 02203 str = "-1743F5DDE2F157e097"; 02204 str >> CXSC_Pip2[16]; 02205 str = "-153F96FFD4AEB5e060"; 02206 str >> CXSC_Pip2[17]; 02207 str = "+13CD6F5847D569e028"; 02208 str >> CXSC_Pip2[18]; 02209 str = "+10001471E79A7Be000"; 02210 str >> CXSC_Pip2[19]; 02211 str = "+10001471E79A8Be000"; 02212 str >> CXSC_Pip2[20]; 02213 02214 CXSC_Pip2_initialized = true; 02215 std::cout << RestoreOpt; 02216 } 02217 stagprec = stagmax; 02218 y = adjust(l_interval(0.0)); 02219 for (int i=0; i <= stagmax; i++) 02220 y.data[i] = CXSC_Pip2[i]; 02221 // y[i+1] = CXSC_Pip2[i]; 02222 stagprec = stagsave; 02223 y = adjust(y); 02224 return y; 02225 } 02226 02227 02228 // ********************************************************************** 02229 02230 static real CXSC_Sqrt2r[21]; // CXSC_Sqrt2r[0], ... CXSC_Sqrt2r[20] 02231 static bool CXSC_Sqrt2r_initialized = false; 02232 02233 l_interval Sqrt2r_l_interval() throw() 02234 // Inclusion of 1/sqrt(2), Blomquist, 12.09.06; 02235 { 02236 l_interval y; 02237 int stagsave = stagprec, 02238 stagmax = 20; 02239 02240 if (!CXSC_Sqrt2r_initialized) 02241 { 02242 std::string str; 02243 std::cout << SaveOpt; 02244 std::cout << Hex; 02245 str = "+16A09E667F3BCDe3FE"; 02246 str >> CXSC_Sqrt2r[0]; 02247 str = "-1BDD3413B26456e3C8"; 02248 str >> CXSC_Sqrt2r[1]; 02249 str = "+157D3E3ADEC175e392"; 02250 str >> CXSC_Sqrt2r[2]; 02251 str = "+12775099DA2F59e35A"; 02252 str >> CXSC_Sqrt2r[3]; 02253 str = "+160CCE64552BF2e321"; 02254 str >> CXSC_Sqrt2r[4]; 02255 str = "+1821D5C5161D46e2E8"; 02256 str >> CXSC_Sqrt2r[5]; 02257 str = "-1C032046F8498Ee2B2"; 02258 str >> CXSC_Sqrt2r[6]; 02259 str = "+1EE950BC8738F7e27A"; 02260 str >> CXSC_Sqrt2r[7]; 02261 str = "-1AC3FDBC64E103e244"; 02262 str >> CXSC_Sqrt2r[8]; 02263 str = "+13B469101743A1e20C"; 02264 str >> CXSC_Sqrt2r[9]; 02265 str = "+15E3E9CA60B38Ce1D6"; 02266 str >> CXSC_Sqrt2r[10]; 02267 str = "+11BC337BCAB1BDe19B"; 02268 str >> CXSC_Sqrt2r[11]; 02269 str = "-1BBA5DEE9D6E7De165"; 02270 str >> CXSC_Sqrt2r[12]; 02271 str = "-1438DD083B1CC4e12F"; 02272 str >> CXSC_Sqrt2r[13]; 02273 str = "+1B56A28E2EDFA7e0F9"; 02274 str >> CXSC_Sqrt2r[14]; 02275 str = "+1CCB2A634331F4e0C3"; 02276 str >> CXSC_Sqrt2r[15]; 02277 str = "-1BD9056876F83Ee08C"; 02278 str >> CXSC_Sqrt2r[16]; 02279 str = "-1234FA22AB6BEFe056"; 02280 str >> CXSC_Sqrt2r[17]; 02281 str = "+19040CA4A81395e01F"; 02282 str >> CXSC_Sqrt2r[18]; 02283 str = "-10000015249C0Ce000"; 02284 str >> CXSC_Sqrt2r[19]; 02285 str = "-10000015249C0Be000"; 02286 str >> CXSC_Sqrt2r[20]; 02287 02288 CXSC_Sqrt2r_initialized = true; 02289 std::cout << RestoreOpt; 02290 } 02291 stagprec = stagmax; 02292 y = adjust(l_interval(0.0)); 02293 for (int i=0; i <= stagmax; i++) 02294 y.data[i] = CXSC_Sqrt2r[i]; 02295 // y[i+1] = CXSC_Sqrt2r[i]; 02296 stagprec = stagsave; 02297 y = adjust(y); 02298 return y; 02299 } 02300 02301 02302 // ********************************************************************** 02303 02304 static real CXSC_Sqrt3[21]; // CXSC_Sqrt3[0], ... CXSC_Sqrt3[20] 02305 static bool CXSC_Sqrt3_initialized = false; 02306 02307 l_interval Sqrt3_l_interval() throw() 02308 // Inclusion of sqrt(3), Blomquist, 12.09.06; 02309 { 02310 l_interval y; 02311 int stagsave = stagprec, 02312 stagmax = 20; 02313 02314 if (!CXSC_Sqrt3_initialized) 02315 { 02316 std::string str; 02317 std::cout << SaveOpt; 02318 std::cout << Hex; 02319 str = "+1BB67AE8584CAAe3FF"; 02320 str >> CXSC_Sqrt3[0]; 02321 str = "+1CEC95D0B5C1E3e3C9"; 02322 str >> CXSC_Sqrt3[1]; 02323 str = "-1F11DB689F2CCFe391"; 02324 str >> CXSC_Sqrt3[2]; 02325 str = "+13DA4798C720A6e35B"; 02326 str >> CXSC_Sqrt3[3]; 02327 str = "+121B9169B89243e325"; 02328 str >> CXSC_Sqrt3[4]; 02329 str = "-1813508751212Be2EC"; 02330 str >> CXSC_Sqrt3[5]; 02331 str = "-1B3D547B775C1Ee2B5"; 02332 str >> CXSC_Sqrt3[6]; 02333 str = "-19D986D92E2F0Ae27C"; 02334 str >> CXSC_Sqrt3[7]; 02335 str = "+1A34334CE806B6e245"; 02336 str >> CXSC_Sqrt3[8]; 02337 str = "+1A383B9E122E61e20F"; 02338 str >> CXSC_Sqrt3[9]; 02339 str = "+1C61D736F2F6F2e1D8"; 02340 str >> CXSC_Sqrt3[10]; 02341 str = "-10AF49233F9250e1A1"; 02342 str >> CXSC_Sqrt3[11]; 02343 str = "-1558A109EC0523e16A"; 02344 str >> CXSC_Sqrt3[12]; 02345 str = "+1F799D4D4FF2BCe134"; 02346 str >> CXSC_Sqrt3[13]; 02347 str = "-1AD7B219E34EDBe0FE"; 02348 str >> CXSC_Sqrt3[14]; 02349 str = "+15AB940B6677E3e0C8"; 02350 str >> CXSC_Sqrt3[15]; 02351 str = "-1D9B2A8203B8F0e091"; 02352 str >> CXSC_Sqrt3[16]; 02353 str = "-1DB0C8975A3834e05B"; 02354 str >> CXSC_Sqrt3[17]; 02355 str = "-1BCAAB3F6BE884e025"; 02356 str >> CXSC_Sqrt3[18]; 02357 str = "+100000531C2B6Ce000"; 02358 str >> CXSC_Sqrt3[19]; 02359 str = "+100000531C2B6De000"; 02360 str >> CXSC_Sqrt3[20]; 02361 02362 CXSC_Sqrt3_initialized = true; 02363 std::cout << RestoreOpt; 02364 } 02365 stagprec = stagmax; 02366 y = adjust(l_interval(0.0)); 02367 for (int i=0; i <= stagmax; i++) 02368 y.data[i] = CXSC_Sqrt3[i]; 02369 // y[i+1] = CXSC_Sqrt3[i]; 02370 stagprec = stagsave; 02371 y = adjust(y); 02372 return y; 02373 } 02374 02375 02376 // ********************************************************************** 02377 02378 static real CXSC_Sqrt3d2[21]; // CXSC_Sqrt3d2[0], ... CXSC_Sqrt3d2[20] 02379 static bool CXSC_Sqrt3d2_initialized = false; 02380 02381 l_interval Sqrt3d2_l_interval() throw() 02382 // Inclusion of Sqrt(3)/2, Blomquist, 12.09.06; 02383 { 02384 l_interval y; 02385 int stagsave = stagprec, 02386 stagmax = 20; 02387 02388 if (!CXSC_Sqrt3d2_initialized) 02389 { 02390 std::string str; 02391 std::cout << SaveOpt; 02392 std::cout << Hex; 02393 str = "+1BB67AE8584CAAe3FE"; 02394 str >> CXSC_Sqrt3d2[0]; 02395 str = "+1CEC95D0B5C1E3e3C8"; 02396 str >> CXSC_Sqrt3d2[1]; 02397 str = "-1F11DB689F2CCFe390"; 02398 str >> CXSC_Sqrt3d2[2]; 02399 str = " +13DA4798C720A6e35A"; 02400 str >> CXSC_Sqrt3d2[3]; 02401 str = "+121B9169B89243e324"; 02402 str >> CXSC_Sqrt3d2[4]; 02403 str = " -1813508751212Be2EB"; 02404 str >> CXSC_Sqrt3d2[5]; 02405 str = "-1B3D547B775C1Ee2B4"; 02406 str >> CXSC_Sqrt3d2[6]; 02407 str = "-19D986D92E2F0Ae27B"; 02408 str >> CXSC_Sqrt3d2[7]; 02409 str = "+1A34334CE806B6e244"; 02410 str >> CXSC_Sqrt3d2[8]; 02411 str = "+1A383B9E122E61e20E"; 02412 str >> CXSC_Sqrt3d2[9]; 02413 str = "+1C61D736F2F6F2e1D7"; 02414 str >> CXSC_Sqrt3d2[10]; 02415 str = "-10AF49233F9250e1A0"; 02416 str >> CXSC_Sqrt3d2[11]; 02417 str = " -1558A109EC0523e169"; 02418 str >> CXSC_Sqrt3d2[12]; 02419 str = "+1F799D4D4FF2BCe133"; 02420 str >> CXSC_Sqrt3d2[13]; 02421 str = "-1AD7B219E34EDBe0FD"; 02422 str >> CXSC_Sqrt3d2[14]; 02423 str = "+15AB940B6677E3e0C7"; 02424 str >> CXSC_Sqrt3d2[15]; 02425 str = "-1D9B2A8203B8F0e090"; 02426 str >> CXSC_Sqrt3d2[16]; 02427 str = "-1DB0C8975A3834e05A"; 02428 str >> CXSC_Sqrt3d2[17]; 02429 str = "-1BCAAB3F6BE884e024"; 02430 str >> CXSC_Sqrt3d2[18]; 02431 str = "+100000298E15B6e000"; 02432 str >> CXSC_Sqrt3d2[19]; 02433 str = "+100000298E15B7e000"; 02434 str >> CXSC_Sqrt3d2[20]; 02435 02436 CXSC_Sqrt3d2_initialized = true; 02437 std::cout << RestoreOpt; 02438 } 02439 stagprec = stagmax; 02440 y = adjust(l_interval(0.0)); 02441 for (int i=0; i <= stagmax; i++) 02442 y.data[i] = CXSC_Sqrt3d2[i]; 02443 // y[i+1] = CXSC_Sqrt3d2[i]; 02444 stagprec = stagsave; 02445 y = adjust(y); 02446 return y; 02447 } 02448 02449 02450 // ********************************************************************** 02451 02452 static real CXSC_Sqrt3r[21]; // CXSC_Sqrt3r[0], ... CXSC_Sqrt3r[20] 02453 static bool CXSC_Sqrt3r_initialized = false; 02454 02455 l_interval Sqrt3r_l_interval() throw() 02456 // Inclusion of 1/Sqrt(3), Blomquist, 12.09.06; 02457 { 02458 l_interval y; 02459 int stagsave = stagprec, 02460 stagmax = 20; 02461 02462 if (!CXSC_Sqrt3r_initialized) 02463 { 02464 std::string str; 02465 std::cout << SaveOpt; 02466 std::cout << Hex; 02467 str = "+1279A74590331Ce3FE"; 02468 str >> CXSC_Sqrt3r[0]; 02469 str = "+134863E0792BEDe3C8"; 02470 str >> CXSC_Sqrt3r[1]; 02471 str = "-1A82F9E6C53222e392"; 02472 str >> CXSC_Sqrt3r[2]; 02473 str = "-1CB0F41134253Ae35C"; 02474 str >> CXSC_Sqrt3r[3]; 02475 str = "+1859ED919EC30Be326"; 02476 str >> CXSC_Sqrt3r[4]; 02477 str = "+1454874FB1F3F4e2EF"; 02478 str >> CXSC_Sqrt3r[5]; 02479 str = "-1DE69C6D3D2741e2B9"; 02480 str >> CXSC_Sqrt3r[6]; 02481 str = "+17EEC450C48BE1e283"; 02482 str >> CXSC_Sqrt3r[7]; 02483 str = "-16F743EEE65D53e24D"; 02484 str >> CXSC_Sqrt3r[8]; 02485 str = "-1887B505D7E7C2e215"; 02486 str >> CXSC_Sqrt3r[9]; 02487 str = "-1484D2E10C1161e1DE"; 02488 str >> CXSC_Sqrt3r[10]; 02489 str = "-1A0B1F86177FB7e1A8"; 02490 str >> CXSC_Sqrt3r[11]; 02491 str = "+1FE389D3F2C54Ee170"; 02492 str >> CXSC_Sqrt3r[12]; 02493 str = "+1F29F77C671544e13A"; 02494 str >> CXSC_Sqrt3r[13]; 02495 str = "-16CE74ED77D9BEe104"; 02496 str >> CXSC_Sqrt3r[14]; 02497 str = "-1E38708FF0CCB5e0CE"; 02498 str >> CXSC_Sqrt3r[15]; 02499 str = "-1F13BCC70157D1e098"; 02500 str >> CXSC_Sqrt3r[16]; 02501 str = "+17EC34CF9B1930e062"; 02502 str >> CXSC_Sqrt3r[17]; 02503 str = "-117A638EFF3A8Be02B"; 02504 str >> CXSC_Sqrt3r[18]; 02505 str = "-10016A8EF69C32e000"; 02506 str >> CXSC_Sqrt3r[19]; 02507 str = "-10016A8EF69C31e000"; 02508 str >> CXSC_Sqrt3r[20]; 02509 02510 CXSC_Sqrt3r_initialized = true; 02511 std::cout << RestoreOpt; 02512 } 02513 stagprec = stagmax; 02514 y = adjust(l_interval(0.0)); 02515 for (int i=0; i <= stagmax; i++) 02516 y.data[i] = CXSC_Sqrt3r[i]; 02517 // y[i+1] = CXSC_Sqrt3r[i]; 02518 stagprec = stagsave; 02519 y = adjust(y); 02520 return y; 02521 } 02522 02523 02524 // ********************************************************************** 02525 02526 static real CXSC_LnPi[21]; // CXSC_LnPi[0], ... CXSC_LnPi[20] 02527 static bool CXSC_LnPi_initialized = false; 02528 02529 l_interval LnPi_l_interval() throw() 02530 // Inclusion of ln(Pi), Blomquist, 12.09.06; 02531 { 02532 l_interval y; 02533 int stagsave = stagprec, 02534 stagmax = 20; 02535 02536 if (!CXSC_LnPi_initialized) 02537 { 02538 std::string str; 02539 std::cout << SaveOpt; 02540 std::cout << Hex; 02541 str = "+1250D048E7A1BDe3FF"; 02542 str >> CXSC_LnPi[0]; 02543 str = "+17ABF2AD8D5088e3C6"; 02544 str >> CXSC_LnPi[1]; 02545 str = "-16CCF43244818Ae38E"; 02546 str >> CXSC_LnPi[2]; 02547 str = "+1F9303719C0176e358"; 02548 str >> CXSC_LnPi[3]; 02549 str = "+15DF52611CB54Ee322"; 02550 str >> CXSC_LnPi[4]; 02551 str = "-1D9056E74F8C97e2EC"; 02552 str >> CXSC_LnPi[5]; 02553 str = "+100B095B6C2E1Ae2B5"; 02554 str >> CXSC_LnPi[6]; 02555 str = "-18C7557878A9E7e27F"; 02556 str >> CXSC_LnPi[7]; 02557 str = "+1B9BBBB4F4CEE7e248"; 02558 str >> CXSC_LnPi[8]; 02559 str = "+1B477FCC702F86e212"; 02560 str >> CXSC_LnPi[9]; 02561 str = "+141F1344A31799e1DC"; 02562 str >> CXSC_LnPi[10]; 02563 str = "+1B6740BE95CD58e1A6"; 02564 str >> CXSC_LnPi[11]; 02565 str = "-1F2C63904D27DBe16E"; 02566 str >> CXSC_LnPi[12]; 02567 str = "+1426F00B933976e136"; 02568 str >> CXSC_LnPi[13]; 02569 str = "+125703BE5FAA20e100"; 02570 str >> CXSC_LnPi[14]; 02571 str = "-1DADAE5397F95Be0C9"; 02572 str >> CXSC_LnPi[15]; 02573 str = "+17C9D110381543e091"; 02574 str >> CXSC_LnPi[16]; 02575 str = "-1259230E627FCAe05B"; 02576 str >> CXSC_LnPi[17]; 02577 str = "+191CEAB6B13A33e024"; 02578 str >> CXSC_LnPi[18]; 02579 str = "+10000109D49A14e000"; 02580 str >> CXSC_LnPi[19]; 02581 str = "+10000109D49A15e000"; 02582 str >> CXSC_LnPi[20]; 02583 02584 CXSC_LnPi_initialized = true; 02585 std::cout << RestoreOpt; 02586 } 02587 stagprec = stagmax; 02588 y = adjust(l_interval(0.0)); 02589 for (int i=0; i <= stagmax; i++) 02590 y.data[i] = CXSC_LnPi[i]; 02591 // y[i+1] = CXSC_LnPi[i]; 02592 stagprec = stagsave; 02593 y = adjust(y); 02594 return y; 02595 } 02596 02597 02598 // ********************************************************************** 02599 02600 static real CXSC_Ln2Pi[21]; // CXSC_Ln2Pi[0], ... CXSC_Ln2Pi[20] 02601 static bool CXSC_Ln2Pi_initialized = false; 02602 02603 l_interval Ln2Pi_l_interval() throw() 02604 // Inclusion of ln(2*Pi), Blomquist, 12.09.06; 02605 { 02606 l_interval y; 02607 int stagsave = stagprec, 02608 stagmax = 20; 02609 02610 if (!CXSC_Ln2Pi_initialized) 02611 { 02612 std::string str; 02613 std::cout << SaveOpt; 02614 std::cout << Hex; 02615 str = "+1D67F1C864BEB5e3FF"; 02616 str >> CXSC_Ln2Pi[0]; 02617 str = "-165B5A1B7FF5DFe3C9"; 02618 str >> CXSC_Ln2Pi[1]; 02619 str = "-1B7F70C13DC1CCe392"; 02620 str >> CXSC_Ln2Pi[2]; 02621 str = "+13458B4DDEC6A3e35C"; 02622 str >> CXSC_Ln2Pi[3]; 02623 str = "+133DAA155D2130e324"; 02624 str >> CXSC_Ln2Pi[4]; 02625 str = "-18A007FC5E501Be2EE"; 02626 str >> CXSC_Ln2Pi[5]; 02627 str = "-15406FA3AA9644e2B4"; 02628 str >> CXSC_Ln2Pi[6]; 02629 str = "-13E8D52A392CC9e27E"; 02630 str >> CXSC_Ln2Pi[7]; 02631 str = "-1A43099131E88De248"; 02632 str >> CXSC_Ln2Pi[8]; 02633 str = "-114835B6623C4De212"; 02634 str >> CXSC_Ln2Pi[9]; 02635 str = "-1ABB7858CF827Ae1DC"; 02636 str >> CXSC_Ln2Pi[10]; 02637 str = "+1D8D7045A5A495e1A6"; 02638 str >> CXSC_Ln2Pi[11]; 02639 str = "+1A26094B3F6FC5e16F"; 02640 str >> CXSC_Ln2Pi[12]; 02641 str = "-1EF27932D0E3D0e137"; 02642 str >> CXSC_Ln2Pi[13]; 02643 str = "-12128804136AB6e100"; 02644 str >> CXSC_Ln2Pi[14]; 02645 str = "+15F8A4AC0BEE17e0C7"; 02646 str >> CXSC_Ln2Pi[15]; 02647 str = "+1892F2A5B69B5Fe091"; 02648 str >> CXSC_Ln2Pi[16]; 02649 str = "+1CC7C09477ADCEe05B"; 02650 str >> CXSC_Ln2Pi[17]; 02651 str = "-116DD579AF074Ae022"; 02652 str >> CXSC_Ln2Pi[18]; 02653 str = "+100000321C8783e000"; 02654 str >> CXSC_Ln2Pi[19]; 02655 str = "+100000321C8784e000"; 02656 str >> CXSC_Ln2Pi[20]; 02657 02658 CXSC_Ln2Pi_initialized = true; 02659 std::cout << RestoreOpt; 02660 } 02661 stagprec = stagmax; 02662 y = adjust(l_interval(0.0)); 02663 for (int i=0; i <= stagmax; i++) 02664 y.data[i] = CXSC_Ln2Pi[i]; 02665 // y[i+1] = CXSC_Ln2Pi[i]; 02666 stagprec = stagsave; 02667 y = adjust(y); 02668 return y; 02669 } 02670 02671 02672 // ********************************************************************** 02673 02674 static real CXSC_E[21]; // CXSC_E[0], ... CXSC_E[20] 02675 static bool CXSC_E_initialized = false; 02676 02677 l_interval E_l_interval() throw() 02678 // Inclusion of e=exp(1), Blomquist, 12.09.06; 02679 { 02680 l_interval y; 02681 int stagsave = stagprec, 02682 stagmax = 20; 02683 02684 if (!CXSC_E_initialized) 02685 { 02686 std::string str; 02687 std::cout << SaveOpt; 02688 std::cout << Hex; 02689 str = "+15BF0A8B145769e400"; 02690 str >> CXSC_E[0]; 02691 str = "+14D57EE2B1013Ae3CA"; 02692 str >> CXSC_E[1]; 02693 str = "-1618713A31D3E2e392"; 02694 str >> CXSC_E[2]; 02695 str = "+1C5A6D2B53C26De35C"; 02696 str >> CXSC_E[3]; 02697 str = "-1F75CDE60219B6e326"; 02698 str >> CXSC_E[4]; 02699 str = "-188C76D93041A1e2EF"; 02700 str >> CXSC_E[5]; 02701 str = "+12FE363630C75Ee2B9"; 02702 str >> CXSC_E[6]; 02703 str = "-1C25F937F544EEe283"; 02704 str >> CXSC_E[7]; 02705 str = "-1E852C20E12A2Ae24D"; 02706 str >> CXSC_E[8]; 02707 str = "-14D4F6DE605705e212"; 02708 str >> CXSC_E[9]; 02709 str = "-1F3225EF539355e1D8"; 02710 str >> CXSC_E[10]; 02711 str = "-16109728625547e1A2"; 02712 str >> CXSC_E[11]; 02713 str = "-194301506D94CFe16C"; 02714 str >> CXSC_E[12]; 02715 str = "-1879C78F8CBA44e136"; 02716 str >> CXSC_E[13]; 02717 str = "-1D5976250C1018e0FD"; 02718 str >> CXSC_E[14]; 02719 str = "+1C877C56284DABe0C7"; 02720 str >> CXSC_E[15]; 02721 str = "+1E73530ACCA4F5e091"; 02722 str >> CXSC_E[16]; 02723 str = "-1F161A150FD53Ae05B"; 02724 str >> CXSC_E[17]; 02725 str = "+159927DB0E8845e022"; 02726 str >> CXSC_E[18]; 02727 str = "+10000094BB2C8Ee000"; 02728 str >> CXSC_E[19]; 02729 str = "+10000094BB2C8Fe000"; 02730 str >> CXSC_E[20]; 02731 02732 CXSC_E_initialized = true; 02733 std::cout << RestoreOpt; 02734 } 02735 stagprec = stagmax; 02736 y = adjust(l_interval(0.0)); 02737 for (int i=0; i <= stagmax; i++) 02738 y.data[i] = CXSC_E[i]; 02739 // y[i+1] = CXSC_E[i]; 02740 stagprec = stagsave; 02741 y = adjust(y); 02742 return y; 02743 } 02744 02745 02746 // ********************************************************************** 02747 02748 static real CXSC_Er[21]; // CXSC_Er[0], ... CXSC_Er[20] 02749 static bool CXSC_Er_initialized = false; 02750 02751 l_interval Er_l_interval() throw() 02752 // Inclusion of 1/e, Blomquist, 12.09.06; 02753 { 02754 l_interval y; 02755 int stagsave = stagprec, 02756 stagmax = 20; 02757 02758 if (!CXSC_Er_initialized) 02759 { 02760 std::string str; 02761 std::cout << SaveOpt; 02762 std::cout << Hex; 02763 str = "+178B56362CEF38e3FD"; 02764 str >> CXSC_Er[0]; 02765 str = "-1CA8A4270FADF5e3C6"; 02766 str >> CXSC_Er[1]; 02767 str = "-1837912B3FD2AAe390"; 02768 str >> CXSC_Er[2]; 02769 str = "-152711999FB68Ce35A"; 02770 str >> CXSC_Er[3]; 02771 str = "-17AD7C1289274Ee324"; 02772 str >> CXSC_Er[4]; 02773 str = "+17E8E56842B705e2E6"; 02774 str >> CXSC_Er[5]; 02775 str = "-1D24CB13796C2De2B0"; 02776 str >> CXSC_Er[6]; 02777 str = "-1456AABDA5C8F2e279"; 02778 str >> CXSC_Er[7]; 02779 str = "+1229F03C6276DDe243"; 02780 str >> CXSC_Er[8]; 02781 str = "-1569CFC4F53109e20D"; 02782 str >> CXSC_Er[9]; 02783 str = "-155B63C9B68091e1D5"; 02784 str >> CXSC_Er[10]; 02785 str = "+1580CF14DC087Ce19F"; 02786 str >> CXSC_Er[11]; 02787 str = "+1F9FF222313669e168"; 02788 str >> CXSC_Er[12]; 02789 str = "+15BC9CB1A22487e132"; 02790 str >> CXSC_Er[13]; 02791 str = "-1857E415C89B13e0FB"; 02792 str >> CXSC_Er[14]; 02793 str = "+13DF75706E3643e0C5"; 02794 str >> CXSC_Er[15]; 02795 str = "+13BDF5B7646234e08D"; 02796 str >> CXSC_Er[16]; 02797 str = "+1C956A5A3BE55De057"; 02798 str >> CXSC_Er[17]; 02799 str = "-167243FE9CD95Ee020"; 02800 str >> CXSC_Er[18]; 02801 str = "+1000002F30CCDBe000"; 02802 str >> CXSC_Er[19]; 02803 str = "+1000002F30CCDCe000"; 02804 str >> CXSC_Er[20]; 02805 02806 CXSC_Er_initialized = true; 02807 std::cout << RestoreOpt; 02808 } 02809 stagprec = stagmax; 02810 y = adjust(l_interval(0.0)); 02811 for (int i=0; i <= stagmax; i++) 02812 y.data[i] = CXSC_Er[i]; 02813 // y[i+1] = CXSC_Er[i]; 02814 stagprec = stagsave; 02815 y = adjust(y); 02816 return y; 02817 } 02818 02819 02820 // ********************************************************************** 02821 02822 static real CXSC_Ep2[21]; // CXSC_Ep2[0], ... CXSC_Ep2[20] 02823 static bool CXSC_Ep2_initialized = false; 02824 02825 l_interval Ep2_l_interval() throw() 02826 // Inclusion of e^2, Blomquist, 12.09.06; 02827 { 02828 l_interval y; 02829 int stagsave = stagprec, 02830 stagmax = 20; 02831 02832 if (!CXSC_Ep2_initialized) 02833 { 02834 std::string str; 02835 std::cout << SaveOpt; 02836 std::cout << Hex; 02837 str = "+1D8E64B8D4DDAEe401"; 02838 str >> CXSC_Ep2[0]; 02839 str = "-19E62E22EFCA4Ce3CA"; 02840 str >> CXSC_Ep2[1]; 02841 str = "+1577508F5CF5EDe394"; 02842 str >> CXSC_Ep2[2]; 02843 str = "-186EF0294C2511e35E"; 02844 str >> CXSC_Ep2[3]; 02845 str = "+177D109F148782e327"; 02846 str >> CXSC_Ep2[4]; 02847 str = "+166BBC354AB700e2F0"; 02848 str >> CXSC_Ep2[5]; 02849 str = "-1273AEC0115969e2BA"; 02850 str >> CXSC_Ep2[6]; 02851 str = "-1C5AE00D3BEEF1e284"; 02852 str >> CXSC_Ep2[7]; 02853 str = "+15ACA3FDC9595Fe24C"; 02854 str >> CXSC_Ep2[8]; 02855 str = "-113FCDFE2B1F0Ce215"; 02856 str >> CXSC_Ep2[9]; 02857 str = "+10EEDFD1AE90C9e1DF"; 02858 str >> CXSC_Ep2[10]; 02859 str = "+1D2CB8EDC7078Be1A9"; 02860 str >> CXSC_Ep2[11]; 02861 str = "+11827A19F175F8e173"; 02862 str >> CXSC_Ep2[12]; 02863 str = "-10267512A9BFB2e13C"; 02864 str >> CXSC_Ep2[13]; 02865 str = "-19A1E2FC413AE3e105"; 02866 str >> CXSC_Ep2[14]; 02867 str = "+1170C7A5981ADBe0CF"; 02868 str >> CXSC_Ep2[15]; 02869 str = "-1FC991480067CFe099"; 02870 str >> CXSC_Ep2[16]; 02871 str = "-12E9A54CF5CFB5e062"; 02872 str >> CXSC_Ep2[17]; 02873 str = "-166FA6C468910Ae02A"; 02874 str >> CXSC_Ep2[18]; 02875 str = "+100043EA6DC142e000"; 02876 str >> CXSC_Ep2[19]; 02877 str = "+100043EA6DC143e000"; 02878 str >> CXSC_Ep2[20]; 02879 02880 CXSC_Ep2_initialized = true; 02881 std::cout << RestoreOpt; 02882 } 02883 stagprec = stagmax; 02884 y = adjust(l_interval(0.0)); 02885 for (int i=0; i <= stagmax; i++) 02886 y.data[i] = CXSC_Ep2[i]; 02887 // y[i+1] = CXSC_Ep2[i]; 02888 stagprec = stagsave; 02889 y = adjust(y); 02890 return y; 02891 } 02892 02893 02894 // ********************************************************************** 02895 02896 static real CXSC_Ep2r[21]; // CXSC_Ep2r[0], ... CXSC_Ep2r[20] 02897 static bool CXSC_Ep2r_initialized = false; 02898 02899 l_interval Ep2r_l_interval() throw() 02900 // Inclusion of 1/e^2, Blomquist, 12.09.06; 02901 { 02902 l_interval y; 02903 int stagsave = stagprec, 02904 stagmax = 20; 02905 02906 if (!CXSC_Ep2r_initialized) 02907 { 02908 std::string str; 02909 std::cout << SaveOpt; 02910 std::cout << Hex; 02911 str = "+1152AAA3BF81CCe3FC"; 02912 str >> CXSC_Ep2r[0]; 02913 str = "-1809224547B4BFe3C6"; 02914 str >> CXSC_Ep2r[1]; 02915 str = "-16A8E079134F13e390"; 02916 str >> CXSC_Ep2r[2]; 02917 str = "+14564CACF0994Ee358"; 02918 str >> CXSC_Ep2r[3]; 02919 str = "+1B796438129AF8e322"; 02920 str >> CXSC_Ep2r[4]; 02921 str = "-1ACFED57EF2AE5e2EC"; 02922 str >> CXSC_Ep2r[5]; 02923 str = "-1A968CBDBB5D9De2B5"; 02924 str >> CXSC_Ep2r[6]; 02925 str = "+1A7238CBD97B71e27C"; 02926 str >> CXSC_Ep2r[7]; 02927 str = "-146C53DB77BB01e245"; 02928 str >> CXSC_Ep2r[8]; 02929 str = "-1EEC161C3EBBD7e20C"; 02930 str >> CXSC_Ep2r[9]; 02931 str = "-12D084DC157ACEe1D5"; 02932 str >> CXSC_Ep2r[10]; 02933 str = "+12A61F46883347e19F"; 02934 str >> CXSC_Ep2r[11]; 02935 str = "+1993BAF10CAE0Be164"; 02936 str >> CXSC_Ep2r[12]; 02937 str = "+1F9224351178FFe12E"; 02938 str >> CXSC_Ep2r[13]; 02939 str = "-1C366D1C7BA64Ae0F7"; 02940 str >> CXSC_Ep2r[14]; 02941 str = "-17D9938EFA4657e0C0"; 02942 str >> CXSC_Ep2r[15]; 02943 str = "+1B6668DF0C1286e08A"; 02944 str >> CXSC_Ep2r[16]; 02945 str = "+1F7A4FFC9B48C6e050"; 02946 str >> CXSC_Ep2r[17]; 02947 str = "+1F3E3AF6F17591e01A"; 02948 str >> CXSC_Ep2r[18]; 02949 str = "+100000006C7831e000"; 02950 str >> CXSC_Ep2r[19]; 02951 str = "+100000006C7832e000"; 02952 str >> CXSC_Ep2r[20]; 02953 02954 CXSC_Ep2r_initialized = true; 02955 std::cout << RestoreOpt; 02956 } 02957 stagprec = stagmax; 02958 y = adjust(l_interval(0.0)); 02959 for (int i=0; i <= stagmax; i++) 02960 y.data[i] = CXSC_Ep2r[i]; 02961 // y[i+1] = CXSC_Ep2r[i]; 02962 stagprec = stagsave; 02963 y = adjust(y); 02964 return y; 02965 } 02966 02967 02968 // ********************************************************************** 02969 02970 static real CXSC_EpPi[21]; // CXSC_EpPi[0], ... CXSC_EpPi[20] 02971 static bool CXSC_EpPi_initialized = false; 02972 02973 l_interval EpPi_l_interval() throw() 02974 // Inclusion of e^(Pi), Blomquist, 12.09.06; 02975 { 02976 l_interval y; 02977 int stagsave = stagprec, 02978 stagmax = 20; 02979 02980 if (!CXSC_EpPi_initialized) 02981 { 02982 std::string str; 02983 std::cout << SaveOpt; 02984 std::cout << Hex; 02985 str = "+1724046EB0933Ae403"; 02986 str >> CXSC_EpPi[0]; 02987 str = "-184C962DD81952e3CD"; 02988 str >> CXSC_EpPi[1]; 02989 str = "-12D659C0BCD22Ee396"; 02990 str >> CXSC_EpPi[2]; 02991 str = "+117496B8A92F91e360"; 02992 str >> CXSC_EpPi[3]; 02993 str = "+16A8C4203E5FCDe32A"; 02994 str >> CXSC_EpPi[4]; 02995 str = "-166B11F99A663Be2F4"; 02996 str >> CXSC_EpPi[5]; 02997 str = "-118EC2076DABB1e2BE"; 02998 str >> CXSC_EpPi[6]; 02999 str = "+19776E5BEB18A5e288"; 03000 str >> CXSC_EpPi[7]; 03001 str = "+1AD4091E84B051e252"; 03002 str >> CXSC_EpPi[8]; 03003 str = "+1E89AA12909B40e21C"; 03004 str >> CXSC_EpPi[9]; 03005 str = "+1ACE3C0DDBB994e1E3"; 03006 str >> CXSC_EpPi[10]; 03007 str = "+141EC9379CBBFEe1AC"; 03008 str >> CXSC_EpPi[11]; 03009 str = "+1FC4E78D00A016e173"; 03010 str >> CXSC_EpPi[12]; 03011 str = "+1608BE35B9A409e13D"; 03012 str >> CXSC_EpPi[13]; 03013 str = "-1A0D8AA90EB6B9e103"; 03014 str >> CXSC_EpPi[14]; 03015 str = "+106FE8AFD21ACFe0CD"; 03016 str >> CXSC_EpPi[15]; 03017 str = "+1C072FEA1BFCAFe095"; 03018 str >> CXSC_EpPi[16]; 03019 str = "+1915B9F352EC68e05B"; 03020 str >> CXSC_EpPi[17]; 03021 str = "-13FA07C37897E9e024"; 03022 str >> CXSC_EpPi[18]; 03023 str = "-100003D8039138e000"; 03024 str >> CXSC_EpPi[19]; 03025 str = "-100003D8039137e000"; 03026 str >> CXSC_EpPi[20]; 03027 03028 CXSC_EpPi_initialized = true; 03029 std::cout << RestoreOpt; 03030 } 03031 stagprec = stagmax; 03032 y = adjust(l_interval(0.0)); 03033 for (int i=0; i <= stagmax; i++) 03034 y.data[i] = CXSC_EpPi[i]; 03035 // y[i+1] = CXSC_EpPi[i]; 03036 stagprec = stagsave; 03037 y = adjust(y); 03038 return y; 03039 } 03040 03041 03042 // ********************************************************************** 03043 03044 static real CXSC_Ep2Pi[21]; // CXSC_Ep2Pi[0], ... CXSC_Ep2Pi[20] 03045 static bool CXSC_Ep2Pi_initialized = false; 03046 03047 l_interval Ep2Pi_l_interval() throw() 03048 // Inclusion of e^(2*Pi), Blomquist, 12.09.06; 03049 { 03050 l_interval y; 03051 int stagsave = stagprec, 03052 stagmax = 20; 03053 03054 if (!CXSC_Ep2Pi_initialized) 03055 { 03056 std::string str; 03057 std::cout << SaveOpt; 03058 std::cout << Hex; 03059 str = "+10BBEEE9177E19e408"; 03060 str >> CXSC_Ep2Pi[0]; 03061 str = "+1C7DD9272526B1e3D0"; 03062 str >> CXSC_Ep2Pi[1]; 03063 str = "+15200F57AB89EDe39A"; 03064 str >> CXSC_Ep2Pi[2]; 03065 str = "-1FCCB6EDBE9C36e363"; 03066 str >> CXSC_Ep2Pi[3]; 03067 str = "+1BEA0BF179A589e32D"; 03068 str >> CXSC_Ep2Pi[4]; 03069 str = "-1F3AD5A6B77F9Ee2F7"; 03070 str >> CXSC_Ep2Pi[5]; 03071 str = "-1622F702B57637e2C0"; 03072 str >> CXSC_Ep2Pi[6]; 03073 str = "-100C09AE818734e287"; 03074 str >> CXSC_Ep2Pi[7]; 03075 str = "+10DA7ADA79EFE6e24D"; 03076 str >> CXSC_Ep2Pi[8]; 03077 str = "+1FF9BF48B72959e216"; 03078 str >> CXSC_Ep2Pi[9]; 03079 str = "-17AD7A3F6D2A14e1E0"; 03080 str >> CXSC_Ep2Pi[10]; 03081 str = "+1FCD4B0FA971E4e1A9"; 03082 str >> CXSC_Ep2Pi[11]; 03083 str = "+193A2CDC04526Be172"; 03084 str >> CXSC_Ep2Pi[12]; 03085 str = "-18CBE5FDFAF25Fe13C"; 03086 str >> CXSC_Ep2Pi[13]; 03087 str = "+1D47EEE171DA93e105"; 03088 str >> CXSC_Ep2Pi[14]; 03089 str = "-15B0F8DA29DB32e0CF"; 03090 str >> CXSC_Ep2Pi[15]; 03091 str = "-19207AD7E637D8e097"; 03092 str >> CXSC_Ep2Pi[16]; 03093 str = "+191CA743F265A6e061"; 03094 str >> CXSC_Ep2Pi[17]; 03095 str = "+1A15069182EF28e02A"; 03096 str >> CXSC_Ep2Pi[18]; 03097 str = "-1000EAC5FC05A9e000"; 03098 str >> CXSC_Ep2Pi[19]; 03099 str = "-1000EAC5FC05A8e000"; 03100 str >> CXSC_Ep2Pi[20]; 03101 03102 CXSC_Ep2Pi_initialized = true; 03103 std::cout << RestoreOpt; 03104 } 03105 stagprec = stagmax; 03106 y = adjust(l_interval(0.0)); 03107 for (int i=0; i <= stagmax; i++) 03108 y.data[i] = CXSC_Ep2Pi[i]; 03109 // y[i+1] = CXSC_Ep2Pi[i]; 03110 stagprec = stagsave; 03111 y = adjust(y); 03112 return y; 03113 } 03114 03115 03116 // ********************************************************************** 03117 03118 static real CXSC_EpPid2[21]; // CXSC_EpPid2[0], ... CXSC_EpPid2[20] 03119 static bool CXSC_EpPid2_initialized = false; 03120 03121 l_interval EpPid2_l_interval() throw() 03122 // Inclusion of e^(Pi/2), Blomquist, 12.09.06; 03123 { 03124 l_interval y; 03125 int stagsave = stagprec, 03126 stagmax = 20; 03127 03128 if (!CXSC_EpPid2_initialized) 03129 { 03130 std::string str; 03131 std::cout << SaveOpt; 03132 std::cout << Hex; 03133 str = "+133DEDC855935Fe401"; 03134 str >> CXSC_EpPid2[0]; 03135 str = "+13E45A768FB73Ce3CB"; 03136 str >> CXSC_EpPid2[1]; 03137 str = "-1FB31CF300FF3Ce395"; 03138 str >> CXSC_EpPid2[2]; 03139 str = "-1E80D8BEB83F79e35F"; 03140 str >> CXSC_EpPid2[3]; 03141 str = "-14A3DE039142DDe326"; 03142 str >> CXSC_EpPid2[4]; 03143 str = "-18792D7A37282Be2EB"; 03144 str >> CXSC_EpPid2[5]; 03145 str = "-19DF43A5980C28e2B5"; 03146 str >> CXSC_EpPid2[6]; 03147 str = "-1C6F0F641C0D67e27F"; 03148 str >> CXSC_EpPid2[7]; 03149 str = "-1779C86C2DB5ACe249"; 03150 str >> CXSC_EpPid2[8]; 03151 str = "+168521EE91B16Fe211"; 03152 str >> CXSC_EpPid2[9]; 03153 str = "+12530F905D97BDe1DB"; 03154 str >> CXSC_EpPid2[10]; 03155 str = "+13498112CB7585e1A5"; 03156 str >> CXSC_EpPid2[11]; 03157 str = "+1BA4546B13A434e16F"; 03158 str >> CXSC_EpPid2[12]; 03159 str = "+14FF791C56421Ce138"; 03160 str >> CXSC_EpPid2[13]; 03161 str = "-1F375C223A2152e102"; 03162 str >> CXSC_EpPid2[14]; 03163 str = "-126AB0C8C77412e0CC"; 03164 str >> CXSC_EpPid2[15]; 03165 str = "-1B39C9C0B8C54Ae094"; 03166 str >> CXSC_EpPid2[16]; 03167 str = "-167741414E31E3e05D"; 03168 str >> CXSC_EpPid2[17]; 03169 str = "+1DEFB4462546C1e025"; 03170 str >> CXSC_EpPid2[18]; 03171 str = "-1000010F7B89CDe000"; 03172 str >> CXSC_EpPid2[19]; 03173 str = "-1000010F7B89CCe000"; 03174 str >> CXSC_EpPid2[20]; 03175 03176 CXSC_EpPid2_initialized = true; 03177 std::cout << RestoreOpt; 03178 } 03179 stagprec = stagmax; 03180 y = adjust(l_interval(0.0)); 03181 for (int i=0; i <= stagmax; i++) 03182 y.data[i] = CXSC_EpPid2[i]; 03183 // y[i+1] = CXSC_EpPid2[i]; 03184 stagprec = stagsave; 03185 y = adjust(y); 03186 return y; 03187 } 03188 03189 03190 // ********************************************************************** 03191 03192 static real CXSC_EpPid4[21]; // CXSC_EpPid4[0], ... CXSC_EpPid4[20] 03193 static bool CXSC_EpPid4_initialized = false; 03194 03195 l_interval EpPid4_l_interval() throw() 03196 // Inclusion of e^(Pi/4), Blomquist, 12.09.06; 03197 { 03198 l_interval y; 03199 int stagsave = stagprec, 03200 stagmax = 20; 03201 03202 if (!CXSC_EpPid4_initialized) 03203 { 03204 std::string str; 03205 std::cout << SaveOpt; 03206 std::cout << Hex; 03207 str = "+118BD669471CAAe400"; 03208 str >> CXSC_EpPid4[0]; 03209 str = "+1F0ED609715756e3CA"; 03210 str >> CXSC_EpPid4[1]; 03211 str = "-1B9C7B871FE1DBe394"; 03212 str >> CXSC_EpPid4[2]; 03213 str = "+15C0FECE98F209e35D"; 03214 str >> CXSC_EpPid4[3]; 03215 str = "+18C9FACC5DF3CEe327"; 03216 str >> CXSC_EpPid4[4]; 03217 str = "+15EDE838B4A399e2EF"; 03218 str >> CXSC_EpPid4[5]; 03219 str = "-1C7EFACA363051e2B9"; 03220 str >> CXSC_EpPid4[6]; 03221 str = "-1A1EBEA1646411e283"; 03222 str >> CXSC_EpPid4[7]; 03223 str = "+1AEF54E68CE03Be24C"; 03224 str >> CXSC_EpPid4[8]; 03225 str = "-11250CB97FDDBFe212"; 03226 str >> CXSC_EpPid4[9]; 03227 str = "-169ADC0E65B8A7e1DB"; 03228 str >> CXSC_EpPid4[10]; 03229 str = "+198A501DB90EDDe1A5"; 03230 str >> CXSC_EpPid4[11]; 03231 str = "-1586909A3F6365e16E"; 03232 str >> CXSC_EpPid4[12]; 03233 str = "+1BE542410F8CE7e138"; 03234 str >> CXSC_EpPid4[13]; 03235 str = "+1E7EEC51889EECe102"; 03236 str >> CXSC_EpPid4[14]; 03237 str = "-1913C9FC19333Ce0CC"; 03238 str >> CXSC_EpPid4[15]; 03239 str = "+1112C71EA1E6F0e095"; 03240 str >> CXSC_EpPid4[16]; 03241 str = "-1C4CCF0F5D1E14e05E"; 03242 str >> CXSC_EpPid4[17]; 03243 str = "+1AC4A72310FA27e028"; 03244 str >> CXSC_EpPid4[18]; 03245 str = "-100013EC6A07AEe000"; 03246 str >> CXSC_EpPid4[19]; 03247 str = "-100013EC6A07ADe000"; 03248 str >> CXSC_EpPid4[20]; 03249 03250 CXSC_EpPid4_initialized = true; 03251 std::cout << RestoreOpt; 03252 } 03253 stagprec = stagmax; 03254 y = adjust(l_interval(0.0)); 03255 for (int i=0; i <= stagmax; i++) 03256 y.data[i] = CXSC_EpPid4[i]; 03257 // y[i+1] = CXSC_EpPid4[i]; 03258 stagprec = stagsave; 03259 y = adjust(y); 03260 return y; 03261 } 03262 03263 03264 // ********************************************************************** 03265 03266 static real CXSC_EulerGa[21]; // CXSC_EulerGa[0], ... CXSC_EulerGa[20] 03267 static bool CXSC_EulerGa_initialized = false; 03268 03269 l_interval EulerGa_l_interval() throw() 03270 // Inclusion of EulerGamma, Blomquist, 12.09.06; 03271 { 03272 l_interval y; 03273 int stagsave = stagprec, 03274 stagmax = 20; 03275 03276 if (!CXSC_EulerGa_initialized) 03277 { 03278 std::string str; 03279 std::cout << SaveOpt; 03280 std::cout << Hex; 03281 str = "+12788CFC6FB619e3FE"; 03282 str >> CXSC_EulerGa[0]; 03283 str = "-16CB90701FBFABe3C5"; 03284 str >> CXSC_EulerGa[1]; 03285 str = "-134A95E3133C51e38F"; 03286 str >> CXSC_EulerGa[2]; 03287 str = "+19730064300F7De359"; 03288 str >> CXSC_EulerGa[3]; 03289 str = "-171ECA0084E369e322"; 03290 str >> CXSC_EulerGa[4]; 03291 str = "-1302FE2B078898e2EC"; 03292 str >> CXSC_EulerGa[5]; 03293 str = "+192732D88415F4e2B5"; 03294 str >> CXSC_EulerGa[6]; 03295 str = "+11056AE9132136e27F"; 03296 str >> CXSC_EulerGa[7]; 03297 str = "-17DC6F12E630A3e249"; 03298 str >> CXSC_EulerGa[8]; 03299 str = "+175FD4B1BD70F2e212"; 03300 str >> CXSC_EulerGa[9]; 03301 str = "-19BC9466120C20e1DC"; 03302 str >> CXSC_EulerGa[10]; 03303 str = "-18FD5699260EADe1A6"; 03304 str >> CXSC_EulerGa[11]; 03305 str = "-12EA987665551Fe16F"; 03306 str >> CXSC_EulerGa[12]; 03307 str = "-1FB159BA4A423De138"; 03308 str >> CXSC_EulerGa[13]; 03309 str = "+1FA543D43BCC60e102"; 03310 str >> CXSC_EulerGa[14]; 03311 str = "-1E6F04E0F639F6e0C9"; 03312 str >> CXSC_EulerGa[15]; 03313 str = "-1A23768654F43De091"; 03314 str >> CXSC_EulerGa[16]; 03315 str = "-14F1C5CB4F55EBe058"; 03316 str >> CXSC_EulerGa[17]; 03317 str = "+1E71DF52EDAA7Fe020"; 03318 str >> CXSC_EulerGa[18]; 03319 str = "+1000001C398F9Be000"; 03320 str >> CXSC_EulerGa[19]; 03321 str = "+1000001C398F9Ce000"; 03322 str >> CXSC_EulerGa[20]; 03323 03324 CXSC_EulerGa_initialized = true; 03325 std::cout << RestoreOpt; 03326 } 03327 stagprec = stagmax; 03328 y = adjust(l_interval(0.0)); 03329 for (int i=0; i <= stagmax; i++) 03330 y.data[i] = CXSC_EulerGa[i]; 03331 // y[i+1] = CXSC_EulerGa[i]; 03332 stagprec = stagsave; 03333 y = adjust(y); 03334 return y; 03335 } 03336 03337 03338 // ********************************************************************** 03339 03340 static real CXSC_Catalan[21]; // CXSC_Catalan[0], ... CXSC_Catalan[20] 03341 static bool CXSC_Catalan_initialized = false; 03342 03343 l_interval Catalan_l_interval() throw() 03344 // Inclusion of Catalan-Constant, Blomquist, 12.09.06; 03345 { 03346 l_interval y; 03347 int stagsave = stagprec, 03348 stagmax = 20; 03349 03350 if (!CXSC_Catalan_initialized) 03351 { 03352 std::string str; 03353 std::cout << SaveOpt; 03354 std::cout << Hex; 03355 str = "+1D4F9713E8135De3FE"; 03356 str >> CXSC_Catalan[0]; 03357 str = "+11485608B8DF4De3C5"; 03358 str >> CXSC_Catalan[1]; 03359 str = "-12F39C13BC1EC8e38F"; 03360 str >> CXSC_Catalan[2]; 03361 str = "+1C2FF8094A263Ee357"; 03362 str >> CXSC_Catalan[3]; 03363 str = "+168F335DBE5370e321"; 03364 str >> CXSC_Catalan[4]; 03365 str = "+16291BBB16163Ee2E9"; 03366 str >> CXSC_Catalan[5]; 03367 str = "+124D663F739C43e2B3"; 03368 str >> CXSC_Catalan[6]; 03369 str = "-136A0725ED0E94e27B"; 03370 str >> CXSC_Catalan[7]; 03371 str = "-1D3A26F9C06FCEe240"; 03372 str >> CXSC_Catalan[8]; 03373 str = "-164E42486BFCD2e209"; 03374 str >> CXSC_Catalan[9]; 03375 str = "14F358CFDEC843e1D3"; 03376 str >> CXSC_Catalan[10]; 03377 str = "-11EB82210976ABe19D"; 03378 str >> CXSC_Catalan[11]; 03379 str = "-17D31F6DF5E801e167"; 03380 str >> CXSC_Catalan[12]; 03381 str = "+13FD19CE3E396Ae131"; 03382 str >> CXSC_Catalan[13]; 03383 str = "-1C8CBB3852FF3Fe0F8"; 03384 str >> CXSC_Catalan[14]; 03385 str = "+1A86EB34EAD01Ae0C2"; 03386 str >> CXSC_Catalan[15]; 03387 str = "+1C68C37800513Be087"; 03388 str >> CXSC_Catalan[16]; 03389 str = "+1D46EBB334D7C9e050"; 03390 str >> CXSC_Catalan[17]; 03391 str = "-1944C5E2711625e019"; 03392 str >> CXSC_Catalan[18]; 03393 str = "-100000005E2172e000"; 03394 str >> CXSC_Catalan[19]; 03395 str = "-100000005E2171e000"; 03396 str >> CXSC_Catalan[20]; 03397 03398 CXSC_Catalan_initialized = true; 03399 std::cout << RestoreOpt; 03400 } 03401 stagprec = stagmax; 03402 y = adjust(l_interval(0.0)); 03403 for (int i=0; i <= stagmax; i++) 03404 y.data[i] = CXSC_Catalan[i]; 03405 // y[i+1] = CXSC_Catalan[i]; 03406 stagprec = stagsave; 03407 y = adjust(y); 03408 return y; 03409 } 03410 03411 03412 l_interval atan(const l_interval & x) throw() // aTan(x) 03413 { 03414 int stagsave = stagprec, 03415 stagmax = 19, 03416 digits = 53, 03417 m = 1, 03418 stagsave2, 03419 degree, sign, zk; 03420 real lnt2, lneps, ln3m, eps, 03421 // bas, tn, t3, 03422 zhn = 1.0, 03423 lnb = 0.69314718, // ln(2.0) 03424 ln3 = 1.098612289; // ln(3.0), Blomquist 27.12.03 03425 l_interval t, t2, p, 03426 y; 03427 interval dx = interval(x), 03428 einfachgenau, 03429 err; 03430 03431 einfachgenau = atan(dx); 03432 03433 if (stagprec == 1) 03434 y = atan(dx); 03435 else if (Inf(dx) == Sup(dx) && Sup(dx) == CXSC_Zero) 03436 y = adjust(l_interval(0.0)); 03437 else 03438 { 03439 if (dx == interval(1.0)) y = li_pi4(); // y = Pi/4 03440 else 03441 { 03442 if (stagprec < stagmax) 03443 stagprec++; 03444 else 03445 stagprec = stagmax; 03446 03447 // Argumentreduktion 03448 t = x; 03449 eps = 0.01/stagprec; // vorher 0.1, aber zu gross 03450 while (Sup(abs(interval(t))) > eps) 03451 { // t = t/(1.0+sqrt(1.0+t*t)); 03452 t = t/(1.0+sqrt1px2(t)); // Blomquist, 13.12.02; 03453 zhn += zhn; 03454 } 03455 t2 = t*t; 03456 03457 // Bestimmung des Approximationsgrades 03458 // Beginn Blomquist, 27.12.03 ------------------------------ 03459 err = interval(abs(t2)); 03460 if (expo(Sup(err)) < -300) m = 4; // Vermeidet alte Fehlermeldg. 03461 else 03462 { 03463 // lnt2 = Sup(ln(err)); // Erzeugt alten Fehler nicht mehr!! ALT!! 03464 lnt2 = ln(Sup(err)); // Neu, 24.08.12; 03465 ln3m = ln(3.0-real(Sup(t2))); 03466 lneps = (1.0-digits*stagprec)*lnb; 03467 do { 03468 m += 3; 03469 } while (lneps-ln3-real(m+1)*lnt2+ln(2.0*m+3.0)+ln3m <= 0.0); 03470 } 03471 // Ende Blomquist, 27.12.03 --------------------------------- 03472 03473 degree = m; 03474 03475 // Approximation 03476 sign = (degree%2) ? -1 : 1; 03477 zk = sign*(2*degree+1); 03478 p = l_interval(1.0)/real(zk); 03479 for (int k = degree-1; k >= 0; k--) 03480 { 03481 sign = -sign; 03482 zk = sign*(2*k+1); 03483 p = p*t2+l_interval(1.0)/real(zk); 03484 } // letztes t wird beim Fehler mit einmultipliziert 03485 03486 // Fehler bestimmen 03487 stagsave2 = stagprec; 03488 stagprec = 2; 03489 l_interval relerr; 03490 stagprec = stagsave2; 03491 err = pow(interval(2.0), interval(1.0-digits*stagprec)) 03492 * interval(-1.0, 1.0); // von 2^(2-d*s) ge\x{201E}ndert ! 03493 relerr[1] = 1.0; 03494 relerr[2] = Inf(err); 03495 relerr[3] = Sup(err); 03496 03497 // Argumentreduktion rueckgaengig machen, mit t mult. 03498 // und Fehler einbauen 03499 y = zhn*t*p*relerr; 03500 03501 stagprec = stagsave; 03502 y = adjust(y); 03503 y = y & einfachgenau; 03504 } 03505 } 03506 03507 return y; 03508 } 03509 03510 l_interval acot(const l_interval &x) throw() // ACot(x) 03511 { 03512 l_interval pihalbe, y; 03513 interval dx = interval(x), 03514 einfachgenau; 03515 03516 einfachgenau = acot(dx); 03517 stagprec++; 03518 // pihalbe = 2.0*atan(l_interval(1.0)); 03519 pihalbe = li_pi4(); 03520 times2pown(pihalbe,1); // Blomquist, 05.12.03; 03521 stagprec--; 03522 03523 if (stagprec == 1) 03524 y = acot(dx); 03525 else if (Inf(dx) == Sup(dx) && Sup(dx) == CXSC_Zero) 03526 y = adjust(pihalbe); 03527 else 03528 { 03529 y = pihalbe-atan(x); 03530 // y = atan(1.0/x); andere M�glichkeit der Berechnung 03531 y = adjust(y); 03532 y = y & einfachgenau; 03533 } 03534 03535 return y; 03536 } 03537 03538 l_interval exp(const l_interval & x) throw(ERROR_LINTERVAL_FAK_OVERFLOW) // exp(x) 03539 { 03540 long int n = 2; 03541 int stagsave = stagprec, 03542 stagmax = 19, 03543 degree, 03544 rednum = 0, 03545 digits = 53; 03546 real fak = 2.0, 03547 zhn = 1.0, // kann groesser 32768 werden 03548 tmp, lny, lneps, eps, 03549 lnb = 0.69314718; 03550 l_interval p, t, 03551 y; 03552 interval dx = interval(x), 03553 einfachgenau, 03554 dt, error; 03555 // ilnb =interval(0.6931471, 0.6931572), 03556 // ilny, ilneps; 03557 03558 einfachgenau = exp(dx); 03559 03560 // gueltigen Bereich pruefen 03561 if (stagprec == 1) 03562 y = exp(dx); 03563 else if (Inf(dx) == Sup(dx) && Inf(dx) == CXSC_Zero) 03564 y = adjust(l_interval(1.0)); 03565 else if (Inf(dx)<-708.396418532264) y = einfachgenau; // Blomquist,12.07.04 03566 else 03567 { 03568 if (stagprec < stagmax) 03569 stagprec++; 03570 else 03571 stagprec = stagmax; 03572 03573 if (Sup(dx) <= 0.0) 03574 t = -x; // Werte unter -1e308 nicht auswertbar, 03575 else 03576 t = x; // dafuer andere Ergebnisse exakter 03577 03578 dt = interval(t); 03579 03580 // Argumentreduktion 03581 // eps = 0.000001*real(Sup(t))/stagprec; 03582 // eps = 1.0/pow10(stagprec/4); 03583 eps = 0.1; 03584 while(Sup(dt)/zhn > eps) 03585 { 03586 rednum++; 03587 zhn += zhn; 03588 } 03589 t /= l_interval(zhn); 03590 03591 // Taylor Approximation 03592 dt = interval(t); 03593 03594 // Anzahl der Durchlaeufe bei der Approximation 03595 tmp = Sup(abs(dt)); 03596 if (MinReal > tmp) 03597 tmp = MinReal; 03598 lny = ln(tmp); 03599 lneps = (1.0-digits*stagprec)*lnb; 03600 03601 while(lneps-lnb+ln(fak)-real(n)*lny <= 0.0) 03602 { 03603 n += 3; 03604 if (n > 170) 03605 { // 170! = 7.26*10306 03606 cxscthrow(ERROR_LINTERVAL_FAK_OVERFLOW("l_interval exp(const l_interval & x)")); 03607 n = 170; 03608 break; 03609 } 03610 fak = fak*real(n)*(n-1.0)*(n-2.0); 03611 } 03612 degree = n; 03613 03614 p = t/real(degree); 03615 int i; 03616 for (i=degree-1; i >= 1; i--) 03617 p = (p+1.0)*t/real(i); 03618 03619 // Fehlerabschaetzung 03620 error = interval(-2.0, 2.0)*pow(dt, interval(real(n)))/fak; 03621 03622 // Fehler einbauen 03623 p += 1.0+l_interval(error); 03624 03625 // Argumentreduktion wird rueckgaengig gemacht 03626 for (i = 1; i <= rednum; i++) 03627 p *= p; 03628 03629 if (Sup(dx) <= 0.0) 03630 p = 1.0/p; 03631 03632 stagprec = stagsave; 03633 y = adjust(p); 03634 y = y & einfachgenau; 03635 } 03636 03637 return y; 03638 } 03639 03640 l_interval exp2(const l_interval & x) // 2^x 03641 { 03642 int stagsave = stagprec, 03643 stagmax = 19; 03644 if (stagprec>stagmax) 03645 stagprec = stagmax; 03646 l_interval y; 03647 03648 y = exp(x*Ln2_l_interval()); 03649 03650 stagprec = stagsave; 03651 y = adjust(y); 03652 03653 return y; 03654 } 03655 03656 l_interval exp10(const l_interval & x) // 10^x 03657 { 03658 int stagsave = stagprec, 03659 stagmax = 19; 03660 if (stagprec>stagmax) 03661 stagprec = stagmax; 03662 l_interval y; 03663 03664 y = exp(x*Ln10_l_interval()); 03665 03666 stagprec = stagsave; 03667 y = adjust(y); 03668 03669 return y; 03670 } 03671 03672 l_interval expm1(const l_interval & x) throw() 03673 // exp(x)-1; 03674 { 03675 l_interval y(0.0); 03676 real B,S_absdx,Sbru,abserr; 03677 interval dx = interval(x), 03678 einfachgenau,sx,bru; 03679 int stagsave = stagprec,ex,N, 03680 stagmax = 19; // from exponential function 03681 einfachgenau = expm1(dx); // produces all necessary error messages! 03682 03683 if (stagprec == 1) y = einfachgenau; 03684 else 03685 { 03686 if (stagprec>stagmax) stagprec = stagmax; 03687 S_absdx = Sup( abs(dx) ); 03688 ex = expo(S_absdx); 03689 if (ex<-49) { // Taylor series 03690 // P_N(x) := x + x^2/2! + ... + x^N/N! 03691 sx = interval(S_absdx); // sx: point interval 03692 N = ex-53*stagprec; // This N is here not the polynomial degree! 03693 if (N<-1022) N = -1022; 03694 if (N > 0) goto Fertig; // because S_absdx = 0; 03695 B = comp(0.5,N); // B: rough bound of the absolute error 03696 // The absolute error should be smaller than B. 03697 N=0; bru = sx; // N is now the polynomial degree! 03698 // Calculation of the polynomia degree N: 03699 // Calculation of the absolute error: 03700 // |absolute error| := |(exp(x)-1) - P_N(x)| <= AE; 03701 // AE = S_absdx^(N+1)/[(N+1)!*(1-S_absdx)] 03702 do 03703 { 03704 N++; 03705 bru = (bru * S_absdx)/(N+1); 03706 Sbru = Sup(bru); 03707 } 03708 while(Sbru > B); 03709 bru = bru * 1.00000001; // for 1/(1-S_absdx) 03710 abserr = Sup(bru); 03711 // |absolute error| <= abserr; 03712 // Caculating an inclusion y of P_N(x): 03713 y = x/N; 03714 for (int i=N-1; i>=1; i--) 03715 y = (y+1)*x/i; 03716 // x + x^2/2! + ... + x^N/N! = P_N(x) included by y 03717 // Conserning the absolute error: 03718 y = y + interval(-abserr,abserr); 03719 03720 } else { 03721 stagprec++; 03722 if (stagprec>stagmax) stagprec = stagmax; 03723 y = exp(x) - 1; 03724 } 03725 Fertig: 03726 stagprec = stagsave; // restore old stagprec value 03727 y = adjust(y); // adjust precision of y to old stagprec value 03728 y = y & einfachgenau; 03729 } 03730 return y; 03731 } // expm1 03732 03733 l_interval expmx2(const l_interval& x) 03734 // e^(-x^2); Blomquist, 13.04.04; 03735 { 03736 int stagsave = stagprec, 03737 stagmax = 19; 03738 l_interval y,z = abs(x); 03739 l_real r1,r2; 03740 interval dx = interval(z), 03741 einfachgenau; 03742 einfachgenau = expmx2(dx); 03743 if (stagprec>stagmax) stagprec = stagmax; 03744 03745 if (stagprec == 1) y = expmx2(dx); // simple precision 03746 else if (Inf(z)>30) y = einfachgenau; // verhindert Overflow bei z*z! 03747 else y = exp(-z*z); 03748 03749 stagprec = stagsave; 03750 y = adjust(y); 03751 y = y & einfachgenau; 03752 return y; 03753 } // expmx2 03754 03755 03756 03757 03758 int cxsc_zweihoch(int) throw(); 03759 03760 l_interval ln(const l_interval & x) throw(ERROR_LINTERVAL_STD_FKT_OUT_OF_DEF) // Ln(x) 03761 { 03762 int stagsave = stagprec, 03763 stagmax = 19, 03764 k, 03765 m = 0, 03766 n = 0, 03767 digits = 53, // Anzahl der Mantissenstellen 03768 cmp, 03769 zhn1 = 1, 03770 zhn2 = 1; 03771 bool zwei=false; 03772 real mx, tmp, 03773 // zkp1, 03774 ln2 = 0.69314718, 03775 lnb = 0.69314718, 03776 lny, lneps; 03777 l_interval y, t, t2, p; 03778 interval dx = interval(x), 03779 einfachgenau, 03780 dy, error; 03781 03782 einfachgenau = ln(dx); 03783 // --------------------------------------------------------------------- 03784 if (Sup(dx) > succ(succ(Inf(dx)))) { // Dieser Teil vermeidet eine 03785 y = einfachgenau; y = adjust(y); // Fehlermeldung der pow(..)-Fkt, 03786 goto fertig; // wenn x zu breit ist und die 1 03787 } // enthaelt. Bei zu breiten Inter- 03788 // vallen wird y daher mit einfachgenau berechnet. Blomquist,05.12.03; 03789 // --------------------------------------------------------------------- 03790 if (Inf(x) <= 0.0) 03791 cxscthrow(ERROR_LINTERVAL_STD_FKT_OUT_OF_DEF("l_interval ln(const l_interval & x)")); 03792 else if (stagprec == 1) 03793 y = ln(dx); 03794 else if (Inf(dx) == Sup(dx) && Sup(dx) == CXSC_One) 03795 y = adjust(l_interval(0.0)); 03796 else 03797 { 03798 if (stagprec < stagmax) 03799 stagprec++; 03800 else 03801 stagprec = stagmax; 03802 03803 // erste Argumentreduktion 03804 /* 03805 if (Inf(dx) > 1.0) 03806 { 03807 y = 1.0/x; 03808 dx = interval(x); 03809 neg = true; 03810 } else 03811 y = x; 03812 */ 03813 y = x; 03814 03815 // zweite Argumentreduktion 03816 mx = Sup(dx); 03817 tmp = 1.0/Sup(dx); 03818 if (tmp > mx) 03819 mx=tmp; 03820 tmp = 1.0/Inf(dx); 03821 if (tmp > mx) 03822 mx = tmp; 03823 if (mx > 1.1) 03824 { 03825 cmp = int(_double(ln(ln(mx)/ln(1.1))/ln2)); 03826 if (cmp > 0) 03827 { 03828 if (cmp > 14) 03829 n = 14; 03830 else 03831 n = cmp; 03832 zhn1 = cxsc_zweihoch(n); 03833 } 03834 } 03835 // y = sqrt(y, zhn1); // Old slow version 03836 if (zhn1 != 1) 03837 for (int k=1; k<=n; k++) 03838 y = sqrt(y); // Blomquist, 26.06.2007; 03839 mx = abs(real(Sup(y))); 03840 if (mx > 1.1) 03841 { 03842 cmp = int(_double(ln(ln(mx)/ln(1.1))/ln2)); 03843 if (cmp > 0) 03844 { 03845 if (cmp > 14) 03846 n = 14; 03847 else 03848 n = cmp; 03849 zhn2 = cxsc_zweihoch(n); 03850 } 03851 // y = sqrt(y, zhn2); // Old slow version 03852 if (zhn2 != 1) 03853 for (int k=1; k<=n; k++) 03854 y = sqrt(y); // Blomquist, 26.06.2007; 03855 zwei = TRUE; 03856 } 03857 03858 // dritte Argumentreduktion 03859 t = (y-1.0)/(y+1.0); 03860 03861 // Abschaetzung der Genauigkeit 03862 t2 = t*t; 03863 tmp = Sup(interval(t2)); 03864 if (tmp == 0.0) 03865 tmp = MinReal; 03866 dy = tmp; 03867 lny = Sup(ln(dy)); 03868 lneps = (1.0-digits*stagprec)*lnb; 03869 03870 do { 03871 m += 2; 03872 } while (lneps-lnb+ln(2.0*m+3.0)-real(m+1)*lny <= 0.0); 03873 03874 // Polynomauswertung 03875 p = 0.0; 03876 for (k = 2*m+1; k >= 1; k -= 2) 03877 p = p*t2+l_interval(2.0)/real(k); 03878 p *= t; 03879 03880 // Fehlereinschliessung 03881 error = interval(-4.0, 4.0)*pow(interval(t2), interval(real(m+1)))/(2.0*m+3.0); 03882 03883 // Fehler einbauen 03884 p = p+error; 03885 03886 // Argumentreduktion 2 wird rueckgaengig gemacht 03887 y = real(zhn1)*p; 03888 if (zwei) 03889 y *= real(zhn2); 03890 03891 // Argumentreduktion 1 wird rueckgaengig gemacht 03892 /* if (neg) y = -(y);*/ 03893 03894 stagprec = stagsave; 03895 y = adjust(y); 03896 y = y & einfachgenau; 03897 03898 // Zusaetzliche Fehlerbehandlung noetig, da Underflow-Fehler 03899 // nicht abgefangen werden kann 03900 /* 03901 error = interval(-4.0*pow10(-stagprec*16-1), 4.0*pow10(-stagprec*16-1)); 03902 y += error; 03903 */ 03904 03905 } 03906 fertig: // Blomquist,05.12.03; 03907 return y; 03908 } 03909 03910 int cxsc_zweihoch(int n) throw() 03911 { 03912 // liefert 2^n 03913 03914 int res = 1; 03915 switch (n) 03916 { 03917 case 0: 03918 return 1; 03919 case 1: 03920 return 2; 03921 case 2: 03922 return 4; 03923 default: 03924 { 03925 int y = 1, 03926 x = 2, 03927 zhi = 2; 03928 03929 if (n%2) 03930 res = 2; 03931 y = x*x; 03932 do { 03933 if ((n/zhi)%2) 03934 res *= y; 03935 zhi += zhi; 03936 if (zhi <= n) 03937 y *= y; 03938 } while (zhi <= n); 03939 } 03940 } 03941 return res; 03942 } 03943 03944 l_interval log2(const l_interval & x) // log2(x) 03945 { 03946 int stagsave = stagprec, 03947 stagmax = 19; 03948 if (stagprec>stagmax) 03949 stagprec = stagmax; 03950 l_interval y; 03951 03952 y = ln(x)/Ln2_l_interval(); 03953 03954 stagprec = stagsave; 03955 y = adjust(y); 03956 03957 return y; 03958 } 03959 03960 l_interval log10(const l_interval & x) // log10(x) 03961 { 03962 int stagsave = stagprec, 03963 stagmax = 19; 03964 if (stagprec>stagmax) 03965 stagprec = stagmax; 03966 l_interval y; 03967 03968 y = ln(x)/Ln10_l_interval(); 03969 03970 stagprec = stagsave; 03971 y = adjust(y); 03972 03973 return y; 03974 } 03975 03976 l_interval sinh(const l_interval & x) throw(ERROR_LINTERVAL_FAK_OVERFLOW) // Sinh(x) 03977 { 03978 long int n = 1; 03979 int stagsave = stagprec, 03980 stagmax = 19, 03981 sign = 1, 03982 digits = 53, 03983 degree, stagsave2; 03984 real tmp, 03985 lnt, lneps, 03986 fak = 1.0, 03987 lnb = 0.69314718; 03988 l_interval y, 03989 t, t2, p, pot; 03990 interval dy, 03991 dx = interval(x), 03992 einfachgenau, 03993 ibase = interval(2.0), 03994 // ilnb = interval(0.6931471,0.6931472), 03995 // ilnx, ilneps, in, 03996 err; 03997 // ifak = interval(fak); 03998 03999 einfachgenau = sinh(dx); 04000 04001 if (stagprec == 1) 04002 y = sinh(dx); 04003 else if (Inf(dx) == Sup(dx) && Sup(dx) == CXSC_Zero) 04004 y = x; 04005 else 04006 { 04007 if (stagprec < stagmax) 04008 stagprec++; 04009 else 04010 stagprec = stagmax; 04011 04012 // Ausnuztung der Punktsymmetrie 04013 if (Sup(dx) < 0.0) 04014 { 04015 y = -x; 04016 sign = -1; 04017 } else 04018 y = x; 04019 04020 dy = interval(y); 04021 04022 // Bei Werten �ber 0.5 -- Auswertung �ber e-Funktion 04023 if (Sup(dy) > 0.5) 04024 { 04025 try { 04026 t = exp(y); 04027 } 04028 catch(const ERROR_LINTERVAL_FAK_OVERFLOW &) 04029 { 04030 cxscthrow( ERROR_LINTERVAL_FAK_OVERFLOW("l_interval sinh(const l_interval & x)")); 04031 } 04032 y = sign*0.5*(t-1.0/t); 04033 } 04034 // Auswertung �ber Potenzreihenentwicklung 04035 else 04036 { 04037 t = y; 04038 t2 = t*t; 04039 04040 // Absch�tzung der Genauigkeit 04041 tmp = real(Sup(abs(t))); // bei Andreas und Christian "t2" 04042 if (tmp < MinReal) 04043 tmp = MinReal; 04044 lnt = ln(tmp); 04045 lneps = (1.0-digits*stagprec)*lnb; 04046 do { 04047 n += 3; 04048 if (n > 170) 04049 { // 170! = 7.26*10^306 04050 cxscthrow(ERROR_LINTERVAL_FAK_OVERFLOW("l_interval sinh(const l_interval & x)")); 04051 n = 170; 04052 fak *= n*(n-1.0); 04053 break; 04054 } 04055 fak *= n*(n-1.0)*(n-2.0); 04056 } while(lneps+ln(fak)-real(n)*lnt-lnb <= 0.0); 04057 /* 04058 real bas, tn, t2d; 04059 bas = double(real(power(l_real(2.0),1-digits*stagprec))); 04060 tn = abs(double(real(mid(t)))); 04061 t2d = double(real(mid(t2))); 04062 do { 04063 n += 2.0; 04064 if (n > 170) 04065 { // 170! = 7.26*10^306 04066 errmon(ERROR_LINTERVAL(FAKOVERFLOW)); 04067 n = 170; 04068 fak *= n; 04069 break; 04070 } 04071 tn *= t2d; 04072 fak *= n*(n-1.0); 04073 } while(bas-tn/fak <= 0.0); 04074 */ 04075 degree = 2*(n/2+1); // degree ist gerade 04076 // Achtung degree = 2n+2! 04077 04078 // Auswertung mit Horner-Schema 04079 p = 1.0; 04080 for (int i = degree; i >= 2; i -= 2) 04081 { 04082 pot = interval((i+1.0)*i); 04083 p = (p*(t2))/pot+1.0; 04084 } 04085 p *= t; 04086 04087 // Negative Werte zur�ckwandeln 04088 if (sign == -1) 04089 p = -p; 04090 04091 // Fehlerauswertung 04092 stagsave2 = stagprec; 04093 stagprec = 2; 04094 l_interval relerr; 04095 stagprec = stagsave2; 04096 err = pow(ibase , interval(1.0-digits*stagprec)) * interval(-1.0, 1.0); 04097 relerr[1] = 1.0; 04098 relerr[2] = Inf(err); 04099 relerr[3] = Sup(err); 04100 04101 // Fehler einbauen 04102 y = p*relerr; 04103 } 04104 stagprec = stagsave; 04105 y = adjust(y); 04106 y = y & einfachgenau; 04107 } 04108 return y; 04109 } 04110 04111 l_interval cosh(const l_interval & x) throw(ERROR_LINTERVAL_FAK_OVERFLOW) // Cosh(x) 04112 { 04113 int stagsave = stagprec, 04114 stagmax = 19; 04115 l_interval y, t; 04116 interval dx = interval(x), 04117 einfachgenau; 04118 04119 einfachgenau = cosh(dx); 04120 04121 if (stagprec == 1) 04122 y = cosh(dx); 04123 else if (Inf(dx) == Sup(dx) && Sup(dx) == CXSC_Zero) 04124 y = adjust(l_interval(1.0)); 04125 else 04126 { 04127 if (stagprec < stagmax) 04128 stagprec++; 04129 else 04130 stagprec = stagmax; 04131 if (Sup(dx) <= 0.0) 04132 y = -x; 04133 else 04134 y = x; 04135 04136 try 04137 { 04138 t = exp(y); 04139 } 04140 catch(const ERROR_LINTERVAL_FAK_OVERFLOW &) 04141 { 04142 cxscthrow(ERROR_LINTERVAL_FAK_OVERFLOW("l_interval cosh(const l_interval & x)")); 04143 } 04144 04145 y = 0.5*(t+1.0/t); 04146 if (interval(0.0) <= dx) 04147 SetInf(y,1.0); 04148 stagprec = stagsave; 04149 y = adjust(y); 04150 y = y & einfachgenau; 04151 } 04152 04153 return y; 04154 } 04155 04156 // Fields for tanh-function 04157 int ex_tanh[8] = { 1,-1,-2,-4,-5,-6,-8,-9 }; 04158 int c_tanhN[8] = { 1,-1,2,-17,62,-1382,21844,-929569 }; 04159 int c_tanhD[8] = { 1,3,15,315,2835,155925,6081075,638512875 }; 04160 04161 l_interval tanh(const l_interval & x) throw() 04162 // tanh(x) with Taylor-series and x --> MaxReal; Blomquist 24.05.04; 04163 { 04164 l_interval s,t,y; 04165 interval dx = interval(x), 04166 einfachgenau,r,r2; 04167 int stagsave = stagprec,ex,m,N,k, 04168 stagmax = 19; 04169 bool neg; 04170 real Sdx; 04171 einfachgenau = tanh(dx); 04172 if (stagprec>stagmax) stagprec = stagmax; 04173 if (stagprec == 1) y = tanh(dx); 04174 else if (dx==0) y = 0; 04175 else if (0<=dx) y = einfachgenau; 04176 else // 0 not in x and not in dx: 04177 { 04178 t = x; 04179 neg = Sup(dx)<0; Sdx = Sup(dx); 04180 if (neg) {t = -t; dx = -dx; Sdx = Sup(dx);} // Inf(dx) > 0: 04181 ex = expo(Sdx); 04182 if (ex < -70) { 04183 m = ex - 53*stagprec; // m >= -1074 necessary !! 04184 if (m < -1074) m = -1074; 04185 r = interval(Sdx); 04186 r2 = r*r; 04187 N = 0; 04188 do { 04189 N++; 04190 r = r*r2; 04191 k = expo(Sup(r)) + ex_tanh[N]; 04192 } while (k>m); 04193 r2 = interval(c_tanhN[N])/c_tanhD[N]; 04194 r2 = r * abs(r2); // r2: inclusion of the absolute error 04195 N--; // N: Polynomial degree 04196 y = l_interval(c_tanhN[N])/c_tanhD[N]; 04197 s = t*t; 04198 for (k=N-1; k>=0; k--) // Horner-scheme 04199 y = y*s + l_interval(c_tanhN[k])/c_tanhD[k]; 04200 y = y * t; 04201 y = y + interval(-Sup(r2),Sup(r2)); // Approxim. error 04202 } else if (ex<-4) { 04203 s = sinh(t); y = cosh(t); 04204 if (stagprec<stagmax) stagprec++; 04205 y = s/y; 04206 } else if (Sdx<352) { 04207 times2pown(t,1); 04208 y = 1 - 2/(exp(t)+1); // --> Sup(y) <= 1 !! 04209 } else { 04210 if (Inf(dx)<352) { 04211 t = Inf(t); times2pown(t,1); 04212 y = 1 - 2/(exp(t)+1); 04213 } else y = l_interval(1) - comp(0.5,-1013); 04214 SetSup(y,1); 04215 } 04216 if (neg) y = -y; 04217 } // else 04218 04219 stagprec = stagsave; // restore old stagprec value 04220 y = adjust(y); // adjust precision of y to old stagprec value 04221 y = y & einfachgenau; // optimal inclusion in case of too 04222 // large relative diameters of y; 04223 return y; 04224 } // tanh 04225 04226 // ************************************************************************ 04227 // Fields for coth-function 04228 int ex_coth[8] = { -1,-5,-8,-12,-15,-18,-22,-25 }; 04229 int c_cothN[8] = { 1,-1,2,-1,2,-1382,4,-3617 }; 04230 real c_cothD[8] = { 1,15,315,1575,31185,212837625, 04231 6081075,54273594375.0 }; // components exactly stored! 04232 // ************************************************************************ 04233 04234 l_interval coth(const l_interval & x) throw() 04235 // coth(x); Blomquist 17.04.04; 04236 { 04237 l_interval t,s,c,y; 04238 interval dx = interval(x), 04239 einfachgenau,r,r2; 04240 int stagsave = stagprec,ex,m,N,k, 04241 stagmax = 19; 04242 bool neg; 04243 einfachgenau = coth(dx); // produces all necessary error messages! 04244 04245 if (stagprec == 1) y = coth(dx); 04246 else 04247 { 04248 if (stagprec>stagmax) stagprec = stagmax; 04249 neg = Sup(dx)<0.0; 04250 t = x; 04251 if (neg) { t = -t; dx = -dx; } // Inf(t),Inf(dx) > 0; 04252 ex = expo(Inf(dx)); 04253 if (ex<-66) { // Laurent series 04254 m = -ex - 53*stagprec; 04255 r = interval(Sup(dx)); r2 = r*r; 04256 N = 0; 04257 do { 04258 N++; 04259 if (N>7) { y = einfachgenau; goto Fertig; } 04260 r = r*r2; 04261 k = expo(Sup(r)) + ex_coth[N]; 04262 } while (k>m); 04263 r2 = interval(c_cothN[N])/c_cothD[N]; 04264 r2 = r*abs(r2)/3; // r2: inclusion of the absolute error 04265 N--; // Polynomial degree 04266 y = l_interval(c_cothN[N])/c_cothD[N]; 04267 s = t*t; 04268 for (k=N-1; k>=0; k--) 04269 y = y*s + l_interval(c_cothN[k])/c_cothD[k]; 04270 y = (y*t)/3 + 1/t; 04271 y = y + interval(-Sup(r2),Sup(r2)); // with approx. error 04272 } else if (ex < 2) { 04273 if (stagprec<stagmax) stagprec++; // stagprec <= 19 04274 c = cosh(t); 04275 s = sinh(t); 04276 y = c/s; 04277 } else if (Inf(dx)<353.0) { 04278 if (stagprec<stagmax) stagprec++; // stagprec <= 19 04279 times2pown(t,1); 04280 y = 1.0 + 2.0/(exp(t)-1.0); 04281 } else { // Inf(dx) >= 353.0 04282 y = l_interval(1.0) + comp(0.5,-1016); 04283 SetInf(y,1.0); 04284 } 04285 if (interval(1.0) <= y) SetInf(y,1.0); 04286 if (neg) y = -y; 04287 Fertig: stagprec = stagsave; // restore old stagprec value 04288 y = adjust(y); // adjust precision of y to old stagprec value 04289 y = y & einfachgenau; 04290 } 04291 return y; 04292 } // coth 04293 04294 04295 l_interval acosh(const l_interval & x) throw() 04296 // acosh(x); Blomquist 14.04.04; 04297 { 04298 int stagsave = stagprec,ex1,ex2, 04299 stagmax = 19; 04300 interval dx = interval(x), 04301 einfachgenau; 04302 l_interval y,t; 04303 04304 einfachgenau = acosh(dx); // false definition range stops program 04305 04306 if (stagprec == 1) y = acosh(dx); 04307 else if (Inf(dx) == Sup(dx) && Sup(dx) == 1.0) 04308 y = adjust(l_interval(0.0)); 04309 else 04310 { 04311 if (stagprec < stagmax) 04312 stagprec++; 04313 else stagprec = stagmax; 04314 ex1 = expo(Inf(dx)); ex2 = expo(Sup(dx)); 04315 if (ex1>500) { // acosh(x) approx ln(2) + ln(x): 04316 y = li_ln2() + ln(x); 04317 // Absolute approximation error: 04318 t = 1.0/Inf(x); 04319 y += l_interval(-Sup(t*t),0); // approxim. error realized here 04320 } else if (ex2<2) { 04321 t = x-1; 04322 y = lnp1( t+sqrt(t*(2+t)) ); 04323 } else y = ln(x+sqrtx2m1(x)); 04324 04325 stagprec = stagsave; 04326 y = adjust(y); 04327 y = y & einfachgenau; 04328 } 04329 return y; 04330 } 04331 04332 l_interval asinh(const l_interval & x) 04333 throw(ERROR_LINTERVAL_STD_FKT_OUT_OF_DEF,ERROR_LINTERVAL_FAK_OVERFLOW) 04334 // ASinh(x) hier ohne zeitaufwendiges Newton-Verfahren, da die Funktion 04335 // sqrt1px2(x) = sqrt(1+x*x) zur Verfuegung steht, die Overflow vermeidet! 04336 // Blomquist, 28.12.03; 04337 { 04338 int stagsave = stagprec, 04339 stagmax = 19; 04340 l_interval y, my; 04341 interval dx = interval(x), einfachgenau, 04342 error = abs(dx); 04343 real absSupdx = Sup(error), r; 04344 04345 try 04346 { 04347 einfachgenau = asinh(dx); 04348 04349 if (stagprec == 1) 04350 y = asinh(dx); 04351 else if (Inf(dx) == Sup(dx) && Sup(dx) == CXSC_Zero) 04352 y = x; 04353 else // if (absSupdx < MaxReal) 04354 { 04355 if (stagprec < stagmax) stagprec++; 04356 else stagprec = stagmax; // Erhoete Praezision: 04357 04358 if (absSupdx < 2E-108) { // Taylorreihe von ln(x+sqrt(1+x2)) 04359 // ln(x+sqrt(1+x2)) = x - x3/6 +- ..... 04360 y = x; 04361 dx = interval(absSupdx); 04362 error = dx*dx*dx/6; 04363 r = Sup(error); // r: Oberschranke des absoluten Fehlers 04364 // Jetzt folgt die Addition einer garantierten Einschliessung 04365 // des abs. Fehlers: 04366 y += l_interval(-r,r); 04367 } else 04368 if (Sup(x) < 0.0) 04369 if (absSupdx < 1E+10) y = lnp1( -x+sqrtp1m1(x*x) ); 04370 else y = -ln(-x+sqrt1px2(x)); 04371 else 04372 if (absSupdx < 1E+10) y = lnp1( x+sqrtp1m1(x*x) ); 04373 else y = ln(x+sqrt1px2(x)); 04374 04375 stagprec = stagsave; 04376 y = adjust(y); 04377 y = y & einfachgenau; 04378 } 04379 } // try 04380 catch(const ERROR_LINTERVAL_STD_FKT_OUT_OF_DEF &) 04381 { 04382 cxscthrow(ERROR_LINTERVAL_STD_FKT_OUT_OF_DEF("l_interval asinh(const l_interval & x)")); 04383 } 04384 catch(const ERROR_LINTERVAL_FAK_OVERFLOW &) 04385 { 04386 cxscthrow(ERROR_LINTERVAL_FAK_OVERFLOW("l_interval asinh(const l_interval & x)")); 04387 } 04388 return y; 04389 } // asinh(x) 04390 04391 l_interval atanh(const l_interval & x) 04392 throw(ERROR_LINTERVAL_STD_FKT_OUT_OF_DEF,ERROR_LINTERVAL_FAK_OVERFLOW) 04393 // ATanh(x), Blomquist 29.12.03; 04394 { 04395 int stagsave = stagprec, 04396 stagmax = 19; 04397 l_interval y, my; 04398 interval dx = interval(x), 04399 error = abs(dx), 04400 einfachgenau; 04401 real absSupdx = Sup(error); 04402 04403 einfachgenau = atanh(dx); 04404 04405 if (Inf(x) <= CXSC_MinusOne || Sup(x) >= CXSC_One) 04406 cxscthrow(ERROR_LINTERVAL_STD_FKT_OUT_OF_DEF("l_interval atanh(const l_interval & x)")); 04407 else if (stagprec == 1) 04408 y = atanh(dx); 04409 else if (Inf(dx) == Sup(dx) && Sup(dx) == CXSC_Zero) 04410 y = x; 04411 else 04412 { 04413 if (stagprec < stagmax) stagprec++; 04414 else stagprec = stagmax; // Ab jetzt hoehere Praezision 04415 if (absSupdx < 0.125) { 04416 y = x / (1-x); 04417 times2pown(y,1); // Schnelle Multiplikation mit 2 04418 y = lnp1(y); 04419 times2pown(y,-1); // Schnelle Division durch 2 04420 } else { 04421 y = ln((1.0+x)/(1.0-x)); 04422 times2pown(y,-1); 04423 } 04424 stagprec = stagsave; 04425 y = adjust(y); 04426 y = y & einfachgenau; 04427 } 04428 04429 return y; 04430 } // atanh 04431 04432 l_interval acoth(const l_interval & x) 04433 throw(ERROR_LINTERVAL_STD_FKT_OUT_OF_DEF,ERROR_LINTERVAL_FAK_OVERFLOW) 04434 // Acoth(x) hier ohne das langsame Newtonverfahren; Fue grosse x gilt: 04435 // Acoth = 0.5*ln(1+2/(x-1)); Blomquist, 28.12.03; 04436 { 04437 int stagsave = stagprec, 04438 stagmax = 19; 04439 l_interval y, my; 04440 interval dx = interval(x), 04441 einfachgenau, 04442 error = abs(dx); 04443 real absSupdx = Sup(error); 04444 04445 einfachgenau = acoth(dx); 04446 04447 if ((l_interval(Inf(x))) < l_interval(CXSC_MinusOne,CXSC_One) 04448 || (l_interval(Sup(x))) < l_interval(CXSC_MinusOne,CXSC_One)) 04449 cxscthrow(ERROR_LINTERVAL_STD_FKT_OUT_OF_DEF 04450 ("l_interval acoth(const l_interval & x)")); 04451 else if (stagprec == 1) 04452 y = acoth(dx); 04453 else 04454 { 04455 if (stagprec < stagmax) stagprec++; 04456 else stagprec = stagmax; // Rechg. in erhoeter Praezision: 04457 04458 if (absSupdx > 1E10) { 04459 y = l_interval(2)/(x-1); 04460 y = lnp1(y); 04461 times2pown(y,-1); 04462 } else { 04463 y = ln((x+1.0)/(x-1.0)); 04464 times2pown(y,-1); 04465 } 04466 04467 stagprec = stagsave; 04468 y = adjust(y); 04469 y = y & einfachgenau; 04470 } 04471 04472 return y; 04473 } // acoth 04474 04475 l_interval lnp1(const l_interval& x) throw() 04476 // ln(1+x) = zeta * P(zeta); zeta := x/(2+x); 04477 // P(zeta) := 2 + (2/3)zeta^2 + (2/5)zeta^4 + ... 04478 // Approximation of P(zeta) with P_N(zeta), defined by 04479 // P_N(zeta) := 2 + (2/3)zeta^(2*1) + 04480 // + (2/5)zeta^(2*2) + ... + (2/2N+1)zeta^(2*N); 04481 // Absolute approximation error Delta := P(zeta) - P_N(zeta); 04482 // With z:=zeta^2 it holds: Delta <= (2*z^(N+1))/{(2N+3)*(1-z)}; 04483 { 04484 const real c1(0.693147181); 04485 int stagsave(stagprec), 04486 stagmax(19); // For ln-function 04487 if (stagprec>stagmax) stagprec = stagmax; // stagprec now <= 19; 04488 interval dx(x), // dx of type interval inclusing x. 04489 einfachgenau(lnp1(dx)), // einfachgenau inclusing lnp1(x). 04490 ax( abs(dx) ); // ax inclusing the absolute values of dx. 04491 real gr(Sup(ax)), // gr: maximum of the absolute values from x. 04492 lngr,u; 04493 l_interval t(x); 04494 if (gr < 1E-8) // Series expansion, gr = 1E-8 is Ok! 04495 if (sign(gr)==0) t=0; 04496 else { // Calculation of N for the approximation of P(zeta) 04497 // using the polynomial P_N(zeta) of degree N: 04498 int N(0),TwoNp3(3),k; 04499 k = expo(gr); 04500 if (k>-1019) { 04501 k = k - 1 - 53*stagprec; 04502 lngr = ln(gr/2); 04503 if (k < -1074) k = -1074; 04504 k--; 04505 u = k*c1; // u = (k-1)*ln(2); 04506 do { 04507 N++; 04508 TwoNp3 += 2; 04509 } while(TwoNp3*lngr - ln(TwoNp3) > u); 04510 } 04511 // Polynomial degree 0<= N <=17 is now calculated! 04512 // Calculation of the N+1 polynomial coefficients: 04513 l_interval a_lnp1[18], // Declaration of the coeffits. 04514 zeta,z2; 04515 a_lnp1[0] = 2; 04516 for (int i=1; i<=N; ++i) { 04517 a_lnp1[i] = a_lnp1[0]/(2*i + 1); 04518 } 04519 // Polyn.-coeff. a_lnp1[i] are now calculated; 0<=i<=N; 04520 // The calculation of P_N(zeta) follows now: 04521 zeta = t/(2+t); 04522 z2 = sqr(zeta); 04523 dx = z2; // is used for the approximation error! 04524 t = a_lnp1[N]; // Horner-scheme: 04525 for (int i=N-1; i>=0; --i) { 04526 t = t*z2 + a_lnp1[i]; 04527 } // t = P_N(zeta) is now calculated 04528 // Calculating now the approximation error 04529 // with dx of type interval: 04530 dx = interval(Sup(dx)); 04531 ax = Power(dx,N+1); 04532 times2pown(ax,1); // Multiplication with 2; 04533 dx = ax / ((2*N+3)*(1-dx)); 04534 t = t + l_interval(0.0,Sup(dx)); 04535 // Approximation error implemented now 04536 t = zeta * t; 04537 } 04538 else 04539 if (gr<1) { // Calculation in higher accuracy: 04540 stagprec++; 04541 if (stagprec>stagmax) stagprec = stagmax; 04542 t = ln(1+t); 04543 } else t = ln(1+t); // Normal accuracy: 04544 stagprec = stagsave; 04545 t = adjust(t); 04546 t = t & einfachgenau; // intersection of t and einfachgenau; 04547 // If x is too wide and therefore t because of the inevitable 04548 // interval overestimations, t is the best possible inclusion. 04549 return t; 04550 } 04551 04552 l_interval sqrtx2m1(const l_interval& x) 04553 // sqrt(x^2-1); Blomquist, 13.04.04; 04554 { 04555 int stagsave = stagprec, 04556 stagmax = 19; 04557 l_interval y,z = abs(x); 04558 l_real r1,r2; 04559 interval dx = interval(z), 04560 einfachgenau; 04561 einfachgenau = sqrtx2m1(dx); 04562 if (stagprec>stagmax) stagprec = stagmax; 04563 04564 if (stagprec == 1) y = sqrtx2m1(dx); // simple interval 04565 else if (Inf(z)==1) { 04566 r1=0; z = Sup(z); z = sqrt(z*z-1); 04567 r2 = Sup(z); 04568 y = l_interval(r1,r2); 04569 } else if (expo(Sup(dx))<500) y = sqrt(z*z-1.0); // no overflow! 04570 else { // avoiding overflow using: x-1/(2x) < sqrt(x^2-1) < x 04571 r2 = Sup(z); 04572 z = Inf(z); // z: Point interval 04573 y = 1.0/z; 04574 times2pown(y,-1); 04575 y = z - y; 04576 r1 = Inf(y); 04577 y = l_interval(r1,r2); 04578 } 04579 stagprec = stagsave; 04580 y = adjust(y); 04581 y = y & einfachgenau; 04582 return y; 04583 } // sqrtx2m1 04584 04585 l_interval sqrt1mx2(const l_interval& x) 04586 // sqrt(1-x^2); Blomquist, 13.04.04; 04587 { 04588 int stagsave = stagprec, 04589 stagmax = 19; 04590 l_interval y,z = abs(x); 04591 l_real r1,r2; 04592 interval dx = interval(z), 04593 einfachgenau; 04594 einfachgenau = sqrt1mx2(dx); 04595 if (stagprec>stagmax) stagprec = stagmax; 04596 04597 if (stagprec == 1) y = sqrt1mx2(dx); // simple interval 04598 else { 04599 y = comp(0.5,1023); // y = 2^(+1022) 04600 times2pown(z,511); 04601 y = sqrt(y-z*z); 04602 times2pown(y,-511); 04603 } 04604 stagprec = stagsave; 04605 y = adjust(y); 04606 y = y & einfachgenau; 04607 return y; 04608 } // sqrt1mx2 04609 04610 04611 l_interval ln_sqrtx2y2(const l_interval& x, const l_interval& y) throw() 04612 { // Inclusion of ln(sqrt{x^2+y^2}) in staggered arithmetic 04613 int stagsave = stagprec; 04614 // stagmax = 19; 04615 interval dx = interval(x), 04616 dy = interval(y), 04617 einfachgenau = ln_sqrtx2y2(dx,dy); 04618 dx = abs(dx); dy = abs(dy); 04619 l_interval ax(abs(x)),ay(abs(y)),ar; 04620 int ex_x(expo(Sup(dx))), ex_y(expo(Sup(dy))), 04621 N,N1; 04622 if (ex_y>ex_x) ex_x = ex_y; 04623 if (ex_x>508) { // To avoid overflow: 04624 N = ex_x-500; 04625 times2pown(ax,-N); times2pown(ay,-N); 04626 ar = ax*ax + ay*ay; 04627 ar = ln(ar); 04628 times2pown(ar,-1); 04629 ar += N*li_ln2(); 04630 } else 04631 if (ex_x<-20) { // To avoid underflow: 04632 N = 500 - ex_x; 04633 if (N>1023) { // Avoiding an range error with times2pown(...) 04634 N1 = N-1023; 04635 times2pown(ax,1023); times2pown(ax,N1); 04636 times2pown(ay,1023); times2pown(ay,N1); 04637 } else { 04638 times2pown(ax,N); 04639 times2pown(ay,N); 04640 } 04641 ar = ax*ax + ay*ay; 04642 ar = ln(ar); 04643 times2pown(ar,-1); 04644 ar -= N*li_ln2(); // ar = ar - N*ln(2) 04645 } else { // Calculation with function lnp1: 04646 ar = sqr(ax) + sqr(ay) - 1; 04647 ar = lnp1(ar); 04648 times2pown(ar,-1); 04649 } 04650 stagprec = stagsave; 04651 ar = adjust(ar); 04652 ar = ar & einfachgenau; 04653 04654 return ar; 04655 } 04656 04657 l_interval acoshp1(const l_interval& x) 04658 // acoshp1(x) = acosh(1+x); Blomquist, 20.04.05; 04659 { 04660 int stagsave = stagprec,ex, 04661 stagmax = 19; 04662 l_interval y,t; 04663 l_real lr; 04664 interval dx = interval(x), 04665 einfachgenau; 04666 einfachgenau = acoshp1(dx); 04667 if (stagprec>stagmax) stagprec = stagmax; 04668 04669 if (stagprec == 1) y = acoshp1(dx); // simple interval 04670 else { 04671 ex = expo(Sup(dx)); 04672 if (ex<=-1016) { 04673 t = x; 04674 times2pown(t,1); // fast multiplication with 2 = 2^1; 04675 t = sqrt(t); // t = sqrt(2x); 04676 lr = Inf(t); // Calculating the infimum (Leibniz-series) 04677 y = l_interval(lr)*(1.0 - l_interval(lr)/12.0); 04678 y = l_interval(Inf(y),Sup(t)); // using alternating Leibniz-series 04679 } 04680 else 04681 if (ex<-400) y = lnp1(x*(1.0+sqrt(1+2.0/x))); 04682 else y = acosh(1+x); // x = MaxReal not allowed 04683 } 04684 stagprec = stagsave; 04685 y = adjust(y); 04686 y = y & einfachgenau; 04687 return y; 04688 } // acoshp1 04689 04690 04691 } // namespace cxsc 04692