C-XSC - A C++ Class Library for Extended Scientific Computing  2.5.4
l_rmath.cpp
00001 /*
00002 **  CXSC is a C++ library for eXtended Scientific Computing (V 2.5.4)
00003 **
00004 **  Copyright (C) 1990-2000 Institut fuer Angewandte Mathematik,
00005 **                          Universitaet Karlsruhe, Germany
00006 **            (C) 2000-2014 Wiss. Rechnen/Softwaretechnologie
00007 **                          Universitaet Wuppertal, Germany   
00008 **
00009 **  This library is free software; you can redistribute it and/or
00010 **  modify it under the terms of the GNU Library General Public
00011 **  License as published by the Free Software Foundation; either
00012 **  version 2 of the License, or (at your option) any later version.
00013 **
00014 **  This library is distributed in the hope that it will be useful,
00015 **  but WITHOUT ANY WARRANTY; without even the implied warranty of
00016 **  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00017 **  Library General Public License for more details.
00018 **
00019 **  You should have received a copy of the GNU Library General Public
00020 **  License along with this library; if not, write to the Free
00021 **  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00022 */
00023 
00024 /* CVS $Id: l_rmath.cpp,v 1.28 2014/01/30 17:23:46 cxsc Exp $ */
00025 
00026 #include <l_rmath.hpp>
00027 #include <rmath.hpp>
00028 #include <l_interval.hpp>
00029 
00030 namespace cxsc {
00031 
00032 l_real sqrt(const l_real &x) throw(ERROR_LREAL_STD_FKT_OUT_OF_DEF)
00033 // Blomquist, additional scaling, 10.12.02
00034 {
00035    int stagsave = stagprec, stagmax = 19, stagcalc;
00036    l_real y;
00037    if (sign(x[1])<0) // Blomquist: is faster than  x < 0.0;
00038    cxscthrow(ERROR_LREAL_STD_FKT_OUT_OF_DEF("l_real sqrt(const l_real &x)"));
00039    else if (zero_(x) || x == 1.0 ) // Blomquist, faster than x == 0;
00040       y = x;
00041    else 
00042    {
00043        l_real x1 = x;
00044        int ex = expo(x1[1]);
00045        ex = 1021-ex;        // ex: optimal scaling factor
00046        if (ex%2) ex--;      // ex is now even;
00047        times2pown(x1,ex);   // scaling x1 with 2^ex
00048 
00049        // Einsatz des Newton-Verfahrens y = y-f(y)/f'(y)
00050        if (stagprec < stagmax) stagcalc = stagprec+1;  
00051        else stagcalc = stagmax+1;  // Blomquist with +1; 30.11.02;
00052        y = sqrt(_real(x1));
00053        stagprec = 1; 
00054        while (stagprec < stagcalc) 
00055        {
00056            stagprec += stagprec; // now calculation in double precision:
00057            if (stagprec > stagcalc) stagprec = stagcalc;
00058            // y = 0.5*((x1/y)+y);
00059            y += x1/y;  // Blomquist, 30.11.02;
00060            times2pown(y,-1); // Blomquist, this multiplication is faster!!
00061        }
00062        times2pown(y,-ex/2); // backscaling
00063        stagprec = stagsave; // restore old stagprec
00064        adjust(y);
00065    }
00066    return y;
00067 }
00068 
00069 l_real sqrtx2y2(const l_real& x, 
00070                        const l_real& y) throw() // Blomquist 6.12.02
00071 // Calculation of an approximation of sqrt(x^2+y^2).
00072 // In general the maximum precision is stagprec=19, predifined by the used
00073 // sqrt-function declared in l_rmath.hpp.
00074 // If the difference |exa-exb| of the exponents (base 2) is sufficient high,
00075 // precision and accuracy can be choosen greater 19.
00076 { 
00077     int stagsave = stagprec, stagcalc, stagmax = 19;
00078     l_real a,b,r,r1;
00079     int exa,exb,ex;
00080     a = x;  b = y;
00081     exa = expo(a[1]);
00082     exb = expo(b[1]);
00083     if (exb > exa)
00084     {  // Permutation of a,b:
00085         r = a;  a = b;  b = r;
00086         ex = exa;  exa = exb;  exb = ex;
00087     } // |a| >= |b|
00088     if (sign(a[1]) < 0) a = -a;  // a == abs(a);
00089     if (sign(b[1]==0)) return a;
00090 
00091     // |a| = a >= |b| > 0:
00092     ex = 0;  // initialization
00093     if (6*exb < 5*exa-1071)
00094     {   // approximation with a + 0.5*b^2/a - b^4/(8*a^3), if the next
00095         // Taylor term b^6/(16*a^5) is less than 2^-1074==minreal;
00096         // minreal is the smallest positive number in IEEE-format.
00097         r1 = (b/a);
00098         r = r1*r1;  
00099         times2pown(r,-2);   // r = r/4;
00100         r = 1 - r;
00101         r1 *= b;
00102         times2pown(r1,-1);  // r = r/2;
00103         r *= r1;
00104         r += a; // r = a + 0.5*b^2/a - b^4/(8*a^3);
00105     } else
00106     {   // Scaling ubwards or downwards to achiev optimal accuracy and 
00107         // to avoid overflow by accu-reading:
00108         stagcalc = stagprec;  
00109         if (stagcalc > stagmax) stagcalc = stagmax; // using sqrt(...) a 
00110         stagprec = stagcalc; // higher precision than stagmax makes no sense!
00111 
00112         ex = 511 - exa; // scaling factor to avoide overflow by accu-reading
00113         if (ex < 0)
00114         {   // sqrt(a^2+b^2) = a*sqrt( 1+(b/a)^2 )
00115             r = b/a;
00116             r = a*sqrt(1+r*r);  // sqrt(...) declared in l_rmath.hpp
00117         } else
00118         {   
00119             times2pown(a,ex);   // exact scaling with ex >= 0
00120             times2pown(b,ex);   // exact scaling eith ex >= 0
00121             dotprecision dot(0.0);
00122             accumulate(dot,a,a);
00123             accumulate(dot,b,b);
00124             r = dot; // r with no overflow!
00125             r = sqrt(r);  // sqrt(...) declared in l_rmath.hpp
00126             times2pown(r,-ex);  // back-scaling 
00127         }
00128         stagprec = stagsave;  // restore old stagprec
00129     }
00130     return r;
00131 } // sqrtx2y2(...)
00132 
00133 l_real sqrt1px2(const l_real& x) throw()
00134 // Inclusion of sqrt(1+x^2); Blomquist, 13.12.02;
00135 // With stagmax=19 we get about 16*19=304 exact decimal digits.
00136 {
00137     l_real y,t;
00138     int stagsave, stagmax=19;
00139     stagsave = stagprec;
00140     if (stagprec > stagmax) stagprec = stagmax;
00141     if (expo(x[1]) > 260)
00142     {  // sqrt(1+x^2) = |x| + 1/(2*|x|)
00143         y = abs(x);
00144         t = 1/y;
00145         times2pown(t,-1);
00146         y += t;
00147     } else y = sqrt(1+x*x);
00148     stagprec = stagsave;
00149     y = adjust(y);
00150     return y;
00151 }
00152 
00153 l_real power(const l_real& x, int n)
00154 {
00155    int       stagsave = stagprec, stagmax = 19;
00156    long int  zhi = 2;
00157    l_real    y, neu;
00158    bool      neg=false;
00159 
00160    if (x == 1.0)
00161       y = x;
00162    else if (n == 0)
00163       y = adjust(_l_real(1.0));
00164    else 
00165    {
00166       if (stagprec < stagmax) 
00167          stagprec++;
00168       else
00169          stagprec = stagmax;
00170 
00171       if (n == 1)
00172          y = x;
00173       else if (n == 2)
00174          y = x*x;
00175       else 
00176       {
00177          if (n < 0) 
00178          {
00179             neg = true;
00180             n = -n;
00181          }
00182 
00183          // Initialisierung
00184          if (n%2)   
00185             y = x;
00186          else      
00187             y = 1.0;  // Praezision wird bei 1 Mult. auf
00188                       // aktuellen Wert gesetzt;
00189                       // Berechnung durch binaere Darstellung der n
00190          neu = x*x;
00191          do {
00192             if ((n/zhi)%2)
00193                y *= neu;
00194             zhi += zhi;
00195             if (zhi <= n)  // letzte Mult. entfaellt --> schneller
00196                neu *= neu;
00197          } while (zhi <= n);
00198          if (neg) 
00199             y = 1/(y);
00200       }
00201       stagprec = stagsave;
00202       y = adjust(y);
00203    }
00204    return y;
00205 }
00206 
00207 // real staggered constants (the same as in l_interval.hpp):
00208 l_real Ln2_l_real()   throw()   // ln(2) 
00209 { return mid( Ln2_l_interval() ); }
00210 l_real Ln10_l_real()  throw()   // ln(10)
00211 { return mid( Ln10_l_interval() ); }
00212 l_real Ln10r_l_real() throw()   // 1/ln(10)
00213 { return mid( Ln10r_l_interval()); }
00214 l_real Pid4_l_real()  throw()   // Pi/4
00215 { return mid( Pid4_l_interval() ); }
00216 l_real Sqrt2_l_real() throw()   // sqrt(2)
00217 { return mid( Sqrt2_l_interval() ); }
00218 l_real Sqrt5_l_real() throw()   // sqrt(5)
00219 { return mid( Sqrt5_l_interval() ); }
00220 l_real Sqrt7_l_real() throw()   // sqrt(7)
00221 { return mid( Sqrt7_l_interval() ); }
00222 l_real Ln2r_l_real() throw()     // 1/ln(2)
00223 { return mid( Ln2r_l_interval() ); }
00224 l_real Pi_l_real() throw()       // Pi
00225 { return mid( Pi_l_interval() ); }
00226 l_real Pid2_l_real() throw()     // Pi/2
00227 { return mid( Pid2_l_interval() ); }
00228 l_real Pi2_l_real() throw()      // 2*Pi
00229 { return mid( Pi2_l_interval() ); }
00230 l_real Pid3_l_real() throw()     // Pi/3
00231 { return mid( Pid3_l_interval() ); }
00232 l_real Pir_l_real() throw()      // 1/Pi
00233 { return mid( Pir_l_interval() ); }
00234 l_real Pi2r_l_real() throw()     // 1/(2*Pi)
00235 { return mid( Pi2r_l_interval() ); }
00236 l_real SqrtPi_l_real() throw()   // sqrt(Pi)
00237 { return mid( SqrtPi_l_interval() ); }
00238 l_real Sqrt2Pi_l_real() throw()  // sqrt(2*Pi)
00239 { return mid( Sqrt2Pi_l_interval() ); }
00240 l_real SqrtPir_l_real() throw()  // 1/sqrt(Pi)
00241 { return mid( SqrtPir_l_interval() ); }
00242 l_real Sqrt2Pir_l_real() throw() // 1/sqrt(2*Pi)
00243 { return mid( Sqrt2Pir_l_interval() ); }
00244 l_real Pip2_l_real() throw()     // Pi^2
00245 { return mid( Pip2_l_interval() ); }
00246 l_real Sqrt2r_l_real() throw()   // 1/sqrt(2)
00247 { return mid( Sqrt2r_l_interval() ); }
00248 l_real Sqrt3_l_real() throw()    // sqrt(3)
00249 { return mid( Sqrt3_l_interval() ); }
00250 l_real Sqrt3d2_l_real() throw()  // sqrt(3)/2
00251 { return mid( Sqrt3d2_l_interval() ); }
00252 l_real Sqrt3r_l_real() throw()   // 1/sqrt(3)
00253 { return mid( Sqrt3r_l_interval() ); }
00254 l_real LnPi_l_real() throw()     // ln(Pi)
00255 { return mid( LnPi_l_interval() ); }
00256 l_real Ln2Pi_l_real() throw()    // ln(2*Pi)
00257 { return mid( Ln2Pi_l_interval() ); }
00258 l_real E_l_real() throw()        // e = exp(1)
00259 { return mid( E_l_interval() ); }
00260 l_real Er_l_real() throw()       // 1/e
00261 { return mid( Er_l_interval() ); }
00262 l_real Ep2_l_real() throw()      // e^2
00263 { return mid( Ep2_l_interval() ); }
00264 l_real Ep2r_l_real() throw()     // 1/e^2
00265 { return mid( Ep2r_l_interval() ); }
00266 l_real EpPi_l_real() throw()     // e^Pi
00267 { return mid( EpPi_l_interval() ); }
00268 l_real Ep2Pi_l_real() throw()    // e^(2*Pi)
00269 { return mid( Ep2Pi_l_interval() ); }
00270 l_real EpPid2_l_real() throw()   // e^(Pi/2)
00271 { return mid( EpPid2_l_interval() ); }
00272 l_real EpPid4_l_real() throw()   // e^(Pi/4)
00273 { return mid( EpPid4_l_interval() ); }
00274 l_real EulerGa_l_real() throw()  // EulerGamma
00275 { return mid(EulerGa_l_interval()  ); }
00276 l_real Catalan_l_real() throw()  // Catalan
00277 { return mid( Catalan_l_interval() ); }
00278 
00279 } // namespace cxsc
00280