C-XSC - A C++ Class Library for Extended Scientific Computing  2.5.4
l_imath.cpp
00001 /*
00002 **  CXSC is a C++ library for eXtended Scientific Computing (V 2.5.4)
00003 **
00004 **  Copyright (C) 1990-2000 Institut fuer Angewandte Mathematik,
00005 **                          Universitaet Karlsruhe, Germany
00006 **            (C) 2000-2014 Wiss. Rechnen/Softwaretechnologie
00007 **                          Universitaet Wuppertal, Germany   
00008 **
00009 **  This library is free software; you can redistribute it and/or
00010 **  modify it under the terms of the GNU Library General Public
00011 **  License as published by the Free Software Foundation; either
00012 **  version 2 of the License, or (at your option) any later version.
00013 **
00014 **  This library is distributed in the hope that it will be useful,
00015 **  but WITHOUT ANY WARRANTY; without even the implied warranty of
00016 **  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00017 **  Library General Public License for more details.
00018 **
00019 **  You should have received a copy of the GNU Library General Public
00020 **  License along with this library; if not, write to the Free
00021 **  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00022 */
00023 
00024 /* CVS $Id: 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