C-XSC - A C++ Class Library for Extended Scientific Computing: l

00001 /*
00002 **  CXSC is a C++ library for eXtended Scientific Computing (V 2.5.1)
00003 **
00004 **  Copyright (C) 1990-2000 Institut fuer Angewandte Mathematik,
00005 **                          Universitaet Karlsruhe, Germany
00006 **            (C) 2000-2011 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_real.hpp,v 1.39 2011/06/07 15:17:40 cxsc Exp $ */
00025 
00026 #ifndef _CXSC_L_REAL_HPP_INCLUDED
00027 #define _CXSC_L_REAL_HPP_INCLUDED
00028 
00029 #include <iostream>
00030 #include <string>
00031 #include <cstdlib>
00032 
00033 #include "dot.hpp"
00034 #include "real.hpp"
00035 #include "interval.hpp" // Blomquist 03.10.02; for _interval(t) in times2pown
00036 #include "except.hpp"
00037 
00038 namespace cxsc {
00039 
00040 class l_interval;
00041 class interval;
00042 
00043 #ifdef CXSC_USE_TLS_PREC
00044 
00045 #ifdef _WIN32
00046 extern __declspec(thread) int stagprec;
00047 #else
00048 extern __thread int stagprec;
00049 #endif
00050 
00051 #else
00052 
00053 extern int stagprec;
00054 
00055 #endif
00056 
00057 
00058 
00059 
00061 
00077 class l_real
00078 {
00079    private:
00080       // ---- Data elements ----
00081       int   prec;
00082       real  *data;
00083 
00084    public:
00085       // ---- Constructors ----
00087       l_real(void) throw();
00089       l_real(const l_real &) throw();
00090       
00091 
00093       l_real & operator =(const l_real &) throw();
00095       l_real & operator =(const real &) throw();
00097       l_real & operator =(const dotprecision &) throw();
00099       l_real & operator =(const lx_real &) throw(); // Blomquist, 12.11.2008;
00100 
00101       // ---- Destructors    ----
00102       ~l_real(void) throw();
00103 
00104       // ---- Type casting ----
00106       explicit l_real(int)          throw();
00108       explicit l_real(long)          throw();
00110       explicit l_real(const real &) throw();
00112       explicit l_real(const dotprecision &) throw();
00114       explicit l_real(const double &) throw(); // Blomquist 10.09.02. {l_real.cpp}
00115       
00116       friend real::real(const l_real &) throw();
00117       friend dotprecision::dotprecision(const l_real &) throw();
00118       friend dotprecision & dotprecision::operator =(const l_real &) throw();
00119 
00120 //      friend real _real(const l_real & a) throw() { return real(a); }
00121 //      friend l_real _l_real(const real & a) throw()  { return l_real(a); }
00122 //      friend dotprecision _dotprecision(const l_real & a) throw()
00123 //                                               { return dotprecision(a); }
00124 //      friend l_real _l_real(const dotprecision & a) throw() 
00125 //                                               { return l_real(a); }
00126       friend l_interval _unchecked_l_interval(const l_real &, 
00127                                               const l_real &) throw();
00128 
00129       // The following are defined in the specific vector, matrix-files
00130 #if(CXSC_INDEX_CHECK) 
00131 
00132       explicit INLINE l_real(const l_rvector &)       throw (ERROR_LRVECTOR_TYPE_CAST_OF_THICK_OBJ,ERROR_LRVECTOR_USE_OF_UNINITIALIZED_OBJ);
00134       explicit INLINE l_real(const l_rvector_slice &) throw (ERROR_LRVECTOR_TYPE_CAST_OF_THICK_OBJ,ERROR_LRVECTOR_USE_OF_UNINITIALIZED_OBJ);
00136       explicit INLINE l_real(const l_rmatrix &)       throw (ERROR_LRMATRIX_TYPE_CAST_OF_THICK_OBJ,ERROR_LRMATRIX_USE_OF_UNINITIALIZED_OBJ);
00138       explicit INLINE l_real(const l_rmatrix_slice &) throw (ERROR_LRMATRIX_TYPE_CAST_OF_THICK_OBJ,ERROR_LRMATRIX_USE_OF_UNINITIALIZED_OBJ);
00140 
00145       friend INLINE real _l_real(const l_rvector &)       throw (ERROR_LRVECTOR_TYPE_CAST_OF_THICK_OBJ,ERROR_LRVECTOR_USE_OF_UNINITIALIZED_OBJ);
00147 
00152       friend INLINE real _l_real(const l_rvector_slice &) throw (ERROR_LRVECTOR_TYPE_CAST_OF_THICK_OBJ,ERROR_LRVECTOR_USE_OF_UNINITIALIZED_OBJ);
00154 
00159       friend INLINE real _l_real(const l_rmatrix &)       throw (ERROR_LRMATRIX_TYPE_CAST_OF_THICK_OBJ,ERROR_LRMATRIX_USE_OF_UNINITIALIZED_OBJ);
00161 
00166       friend INLINE real _l_real(const l_rmatrix_slice &) throw (ERROR_LRMATRIX_TYPE_CAST_OF_THICK_OBJ,ERROR_LRMATRIX_USE_OF_UNINITIALIZED_OBJ);
00167 #else
00168 
00169       explicit INLINE l_real(const l_rvector &)       throw ();
00171       explicit INLINE l_real(const l_rvector_slice &) throw ();
00173       explicit INLINE l_real(const l_rmatrix &)       throw ();
00175       explicit INLINE l_real(const l_rmatrix_slice &) throw ();
00177 
00182       friend INLINE real _l_real(const l_rvector &)       throw ();
00184 
00189       friend INLINE real _l_real(const l_rvector_slice &) throw ();
00191 
00196       friend INLINE real _l_real(const l_rmatrix &)       throw ();
00198 
00203       friend INLINE real _l_real(const l_rmatrix_slice &) throw ();
00204 #endif
00205 
00206 
00207       // ---- Output functions ----
00209       friend std::ostream & operator <<(std::ostream &,const l_real &) throw();
00211       friend std::istream & operator >>(std::istream &,l_real &)       throw();
00213       friend std::string & operator <<(std::string &,const l_real &)   throw();
00215       friend std::string & operator >>(std::string &,l_real &)         throw();
00217       friend void          operator >>(const std::string &,l_real &)   throw();
00219       friend void          operator >>(const char *,l_real &)          throw();
00220 
00221       // ---- Standard functions ---- (arithmetic operators)
00223       real&             operator[](int) const throw();
00224 
00226       friend     l_real operator -(const l_real& lr1) throw();
00228       friend     l_real operator +(const l_real& lr1) throw();
00229 
00231       friend     l_real operator +(const l_real &,const l_real &) throw();
00233       friend     l_real operator -(const l_real &,const l_real &) throw();
00235       friend     l_real operator *(const l_real &,const l_real &) throw();
00237       friend     l_real operator /(const l_real &,const l_real &) throw(DIV_BY_ZERO);
00239       friend inline l_interval operator |(const l_real &,const l_real &) throw();
00240 
00242       friend     l_real operator +(const l_real &,const real &) throw();
00244       friend     l_real operator +(const real &,const l_real &) throw();
00246       friend     l_real operator -(const l_real &,const real &) throw();
00248       friend     l_real operator -(const real &,const l_real &) throw();
00250       friend     l_real operator *(const l_real &,const real &) throw();
00252       friend     l_real operator *(const real &,const l_real &) throw();
00254       friend     l_real operator /(const l_real &,const real &) throw();
00256       friend     l_real operator /(const real &,const l_real &) throw();
00258       friend inline l_interval operator |(const real &,const l_real &) throw();
00260       friend inline l_interval operator |(const l_real &,const real &) throw();
00261 
00263       friend  dotprecision operator +(const l_real &,const dotprecision &) throw();
00265       friend  dotprecision operator +(const dotprecision &,const l_real &) throw();
00267       friend  dotprecision operator -(const l_real &,const dotprecision &) throw();
00269       friend  dotprecision operator -(const dotprecision &,const l_real &) throw();
00271       friend inline idotprecision operator |(const dotprecision &,const l_real &) throw();
00273       friend inline idotprecision operator |(const l_real &,const dotprecision &) throw();
00274 
00276       friend     l_real & operator +=(l_real &,const l_real &) throw();
00278       friend     l_real & operator -=(l_real &,const l_real &) throw();
00280       friend     l_real & operator *=(l_real &,const l_real &) throw();
00282       friend     l_real & operator /=(l_real &,const l_real &) throw();
00283 
00285       friend     l_real & operator +=(l_real &,const real &) throw();      
00287       friend     l_real & operator -=(l_real &,const real &) throw();
00289       friend     l_real & operator *=(l_real &,const real &) throw();      
00291       friend     l_real & operator /=(l_real &,const real &) throw();
00292 
00294       friend     real   & operator +=(real &,const l_real &) throw();
00296       friend     real   & operator -=(real &,const l_real &) throw();
00298       friend     real   & operator *=(real &,const l_real &) throw();
00300       friend     real   & operator /=(real &,const l_real &) throw();
00301       
00303       friend     inline dotprecision & operator +=(dotprecision &d,const l_real &lr) throw() { lr._akku_add(d); return d; }
00305       friend     inline dotprecision & operator -=(dotprecision &d,const l_real &lr) throw() { lr._akku_sub(d); return d; }
00306 
00307       // ---- Compare operators ----
00309       friend bool operator ==(const l_real &,const l_real &) throw();
00311       friend bool operator !=(const l_real &,const l_real &) throw();
00313       friend bool operator  <(const l_real &,const l_real &) throw();
00315       friend bool operator  >(const l_real &,const l_real &) throw();
00317       friend bool operator <=(const l_real &,const l_real &) throw();
00319       friend bool operator >=(const l_real &,const l_real &) throw();
00320 
00322       friend bool operator ==(const real &,const l_real &) throw();
00324       friend bool operator !=(const real &,const l_real &) throw();
00326       friend bool operator  <(const real &,const l_real &) throw();
00328       friend bool operator  >(const real &,const l_real &) throw();
00330       friend bool operator <=(const real &,const l_real &) throw();
00332       friend bool operator >=(const real &,const l_real &) throw();
00333 
00335       friend bool operator ==(const l_real &,const real &) throw();
00337       friend bool operator !=(const l_real &,const real &) throw();
00339       friend bool operator  <(const l_real &,const real &) throw();
00341       friend bool operator  >(const l_real &,const real &) throw();
00343       friend bool operator <=(const l_real &,const real &) throw();
00345       friend bool operator >=(const l_real &,const real &) throw();
00346 
00348       friend bool operator ==(const dotprecision &,const l_real &) throw();
00350       friend bool operator !=(const dotprecision &,const l_real &) throw();
00352       friend bool operator  <(const dotprecision &,const l_real &) throw();
00354       friend bool operator  >(const dotprecision &,const l_real &) throw();
00356       friend bool operator <=(const dotprecision &,const l_real &) throw();
00358       friend bool operator >=(const dotprecision &,const l_real &) throw();
00359 
00361       friend bool operator ==(const l_real &,const dotprecision &) throw();
00363       friend bool operator !=(const l_real &,const dotprecision &) throw();
00365       friend bool operator  <(const l_real &,const dotprecision &) throw();
00367       friend bool operator  >(const l_real &,const dotprecision &) throw();
00369       friend bool operator <=(const l_real &,const dotprecision &) throw();
00371       friend bool operator >=(const l_real &,const dotprecision &) throw();
00372 
00374       friend bool operator ==(const interval &,const l_real &) throw();
00376       friend bool operator !=(const interval &,const l_real &) throw();
00378       friend bool operator  <(const interval &,const l_real &) throw();
00380       friend bool operator  >(const interval &,const l_real &) throw();
00382       friend bool operator <=(const interval &,const l_real &) throw();
00384       friend bool operator >=(const interval &,const l_real &) throw();
00385 
00387       friend bool operator ==(const l_real &,const interval &) throw();
00389       friend bool operator !=(const l_real &,const interval &) throw();
00391       friend bool operator  <(const l_real &,const interval &) throw();
00393       friend bool operator  >(const l_real &,const interval &) throw();
00395       friend bool operator <=(const l_real &,const interval &) throw();
00397       friend bool operator >=(const l_real &,const interval &) throw();
00398 
00400       friend bool operator ==(const idotprecision &,const l_real &) throw();
00402       friend bool operator !=(const idotprecision &,const l_real &) throw();
00404       friend bool operator  <(const idotprecision &,const l_real &) throw();
00406       friend bool operator  >(const idotprecision &,const l_real &) throw();
00408       friend bool operator <=(const idotprecision &,const l_real &) throw();
00410       friend bool operator >=(const idotprecision &,const l_real &) throw();
00411 
00413       friend bool operator ==(const l_real &,const idotprecision &) throw();
00415       friend bool operator !=(const l_real &,const idotprecision &) throw();
00417       friend bool operator  <(const l_real &,const idotprecision &) throw();
00419       friend bool operator  >(const l_real &,const idotprecision &) throw();
00421       friend bool operator <=(const l_real &,const idotprecision &) throw();
00423       friend bool operator >=(const l_real &,const idotprecision &) throw();
00424 
00425 
00427       friend bool operator!(const l_real& lr) throw();
00428       
00429       // ---- functions ----
00431       friend void accumulate(dotprecision&, const real&, const l_real&) throw();
00433       friend void accumulate(dotprecision&, const l_real&, const real&) throw();
00435       friend void accumulate(dotprecision&, const l_real&, const l_real&) throw();
00436       
00438       friend void accumulate(idotprecision&,const real&, const l_real&) throw();
00440       friend void accumulate(idotprecision&,const l_real&,const real&) throw();
00442       friend void accumulate(idotprecision&,const l_real&, const l_real&) throw();
00443 
00445       friend l_real   abs  (const l_real&) throw();
00447       friend int      sign (const l_real&) throw();
00449       friend int      StagPrec(const l_real&) throw();
00451       friend l_real   adjust(const l_real&) throw();
00452 
00454       friend l_real   rnd_up(const dotprecision&);    // Blomquist, 20.11.2006;
00456       friend l_real   rnd_down(const dotprecision&);  // Blomquist, 20.11.2006;
00457 
00459     friend int expo_sm(const l_real&); // Blomquist, 25.03.2007;
00460     // Calculating expo(x[k]) of the smallest |x[k]| <> 0.
00461 
00463     friend int expo_gr(const l_real&); // Blomquist, 25.03.2007;
00464     // Calculating expo(x[k]) of the greatest |x[k]|.
00465 
00466       // ---- Friends      -----
00467       
00469       friend inline l_real Inf(const l_interval &) throw();
00471       friend inline l_real Sup(const l_interval &) throw();
00473       friend        l_real mid(const l_interval &) throw();
00475       friend inline bool zero_(const l_real &) throw(); // Blomquist,27.11.02
00476    private:
00477       void _clear(int) throw(); // filling a l_real number from element int p
00478                                 //  up to the end with zero.
00479       void _akku_out(const dotprecision&) throw(); // The dotprecision value is rounded to the
00480                         // activated l_real number in its own precision.
00481       void _akku_out_up(const dotprecision&) throw(); // The dotprecision value is rounded up to 
00482                       // the activated l_real number in its own precision.
00483       void _akku_out_down(const dotprecision&) throw(); // The dotprecision value is rounded down 
00484                       // to the activated l_real number in its own precision.
00485       void _akku_add(dotprecision&) const throw(); // adding the activated
00486                                      // l_real number to the accumulator d.
00487       void _akku_sub(dotprecision&) const throw(); // subtracting the 
00488       // activated l_real number to the accumulator d of type dotprecision.
00489       inline real& elem(int i) const throw() {  return data[i-1];  }
00490 };
00491 
00492 
00493 
00494       // ---- Compare operators ----
00495       bool operator ==(const l_real &,const l_real &) throw();
00496       bool operator !=(const l_real &,const l_real &) throw();
00497       bool operator  <(const l_real &,const l_real &) throw();
00498       bool operator  >(const l_real &,const l_real &) throw();
00499       bool operator <=(const l_real &,const l_real &) throw();
00500       bool operator >=(const l_real &,const l_real &) throw();
00501 
00502       bool operator ==(const real &,const l_real &) throw();
00503       bool operator !=(const real &,const l_real &) throw();
00504       bool operator  <(const real &,const l_real &) throw();
00505       bool operator  >(const real &,const l_real &) throw();
00506       bool operator <=(const real &,const l_real &) throw();
00507       bool operator >=(const real &,const l_real &) throw();
00508 
00509       bool operator ==(const l_real &,const real &) throw();
00510       bool operator !=(const l_real &,const real &) throw();
00511       bool operator  <(const l_real &,const real &) throw();
00512       bool operator  >(const l_real &,const real &) throw();
00513       bool operator <=(const l_real &,const real &) throw();
00514       bool operator >=(const l_real &,const real &) throw();
00515 
00516       bool operator ==(const dotprecision &,const l_real &) throw();
00517       bool operator !=(const dotprecision &,const l_real &) throw();
00518       bool operator  <(const dotprecision &,const l_real &) throw();
00519       bool operator  >(const dotprecision &,const l_real &) throw();
00520       bool operator <=(const dotprecision &,const l_real &) throw();
00521       bool operator >=(const dotprecision &,const l_real &) throw();
00522 
00523       bool operator ==(const l_real &,const dotprecision &) throw();
00524       bool operator !=(const l_real &,const dotprecision &) throw();
00525       bool operator  <(const l_real &,const dotprecision &) throw();
00526       bool operator  >(const l_real &,const dotprecision &) throw();
00527       bool operator <=(const l_real &,const dotprecision &) throw();
00528       bool operator >=(const l_real &,const dotprecision &) throw();
00529 
00530       bool operator ==(const interval &,const l_real &) throw();
00531       bool operator !=(const interval &,const l_real &) throw();
00532       bool operator  <(const interval &,const l_real &) throw();
00533       bool operator  >(const interval &,const l_real &) throw();
00534       bool operator <=(const interval &,const l_real &) throw();
00535       bool operator >=(const interval &,const l_real &) throw();
00536 
00537       bool operator ==(const l_real &,const interval &) throw();
00538       bool operator !=(const l_real &,const interval &) throw();
00539       bool operator  <(const l_real &,const interval &) throw();
00540       bool operator  >(const l_real &,const interval &) throw();
00541       bool operator <=(const l_real &,const interval &) throw();
00542       bool operator >=(const l_real &,const interval &) throw();
00543 
00544       bool operator ==(const idotprecision &,const l_real &) throw();
00545       bool operator !=(const idotprecision &,const l_real &) throw();
00546       bool operator  <(const idotprecision &,const l_real &) throw();
00547       bool operator  >(const idotprecision &,const l_real &) throw();
00548       bool operator <=(const idotprecision &,const l_real &) throw();
00549       bool operator >=(const idotprecision &,const l_real &) throw();
00550 
00551       bool operator ==(const l_real &,const idotprecision &) throw();
00552       bool operator !=(const l_real &,const idotprecision &) throw();
00553       bool operator  <(const l_real &,const idotprecision &) throw();
00554       bool operator  >(const l_real &,const idotprecision &) throw();
00555       bool operator <=(const l_real &,const idotprecision &) throw();
00556       bool operator >=(const l_real &,const idotprecision &) throw();
00557 
00558 
00559       bool operator!(const l_real& lr) throw();
00560 
00561 
00562 inline real _real(const l_real & a) throw() { return real(a); }
00563 inline l_real _l_real(const real & a) throw()  { return l_real(a); }
00564 inline dotprecision _dotprecision(const l_real & a) throw()
00565                                                { return dotprecision(a); }
00566 inline l_real _l_real(const dotprecision & a) throw() 
00567                                                { return l_real(a); }
00568 l_interval _unchecked_l_interval(const l_real &, const l_real &) throw();
00569 
00571 l_real   rnd_up(const dotprecision&);    // Blomquist, 20.11.2006;
00573 l_real   rnd_down(const dotprecision&);  // Blomquist, 20.11.2006;
00574 
00576     inline l_real l_pow2n(const int n) throw()
00577 {   // Fast and exact calculation of 2^n; -1074 <= n <= 1023;
00578     // Blomquist 01.10.02.
00579     return l_real( comp(0.5,n+1) );
00580 }
00581 
00583 inline void times2pown(l_real& lr, const int n) throw() // Blomquist 03.10.02
00584 { // lr is multiplied with 2^n; if lr[i]*2^n are all normalized, the result
00585   // is exact. if one of the lr[i]*2^n are denormalized, the result is not
00586   // exact in general.
00587     int k = StagPrec(lr);
00588     for (int i=1; i<=k; i++)
00589     {
00590         times2pown(lr[i],n);
00591     }
00592 }
00593 
00595 
00607 inline void times2pown(l_real& lr, interval& z, const int n) throw() 
00608 {  // Blomquist 03.10.02;
00609 
00610     if ( n<-1074 || n>1023 ) 
00611     { std::cerr << "Error in:  " 
00612            << "times2pown(l_real& lr, interval& z, const int n): " << std::endl
00613            << " -1074 <= n <= +1023 not fulfilled" << std::endl; exit(0); } 
00614     int k = StagPrec(lr);
00615     z = 0;
00616     real mt,t;
00617     real F = comp(0.5,n+1);
00618     for (int i=1; i<=k; i++)
00619     {
00620         mt = mant(lr[i]);
00621         t = lr[i];
00622         times2pown(lr[i],n);
00623         if ( mt != mant(lr[i]) ) 
00624         {
00625             lr[i] = 0;
00626             z += _interval(t) * F;
00627         }
00628     }
00629 
00630 }
00631 
00633 
00644 inline void Times2pown(l_real& a, interval& z, int n) throw()
00645 // If we denote the old value of a with y then with the
00646 // new calculated values of a and z it holds:
00647 //  -----  a+z is an inclusion of y*2^n;  ------
00648 // If z==0 then it holds a = y*2^n; (exact multiplication!).
00649 // Especially, if n>=0, the multiplication with 2^n is exact
00650 // if no oveflow occurs.
00651 {
00652     int fac,rest;
00653     interval z1;
00654     z=0;
00655     if (n>=0)
00656     {   // if (n>=0) a*2^n is exactly calculated
00657         // if no overflow occurs:
00658         fac = n/1023;  rest = n%1023;
00659         for (int k=1; k<=fac; k++)
00660             times2pown(a,1023);
00661         times2pown(a,rest);
00662     }
00663     else // n < 0:
00664         if (n<-2100)
00665         {
00666             if(a>0) z = interval(0,minreal);
00667             else 
00668                 if (a<0) z = interval(-minreal,0);
00669                 else z=0;
00670             a = 0;
00671         } else // -2100<=n<0
00672         {
00673             fac = n/-1074;  rest = n%-1074;
00674             for (int k=1; k<=fac; k++)
00675             {
00676                 times2pown(a,z1,-1074);
00677                 times2pown(z,-1074);
00678                 z += z1;
00679             }
00680             times2pown(a,z1,rest);
00681             times2pown(z,rest);
00682             z += z1;
00683         }
00684 } // void Times2pown(...)
00685 
00687 
00695 inline void Times2pown(l_real& a, const real& p) throw()
00696 // The first parameter delivers an approximation of a * 2^p;
00697 // For p in [-2100,+2100] p must be an integer value.
00698 // This condition is NOT tested in this function!
00699 // For p outside [-2100,+2100] an approximation of a * 2^p is
00700 // calculated for any p of type real, unless an overflow occurs.
00701 // If the function is activated with the second parameter of type int,
00702 // then the first parameter delivers approximations of a * 2^p,
00703 // unless an overflow occurs.
00704 // Blomquist, 04.11.2008;
00705 {
00706  const int c1 = -1000000, 
00707            c2 = 2100,
00708                           c3 = 1023,
00709                 c4 = -1074;
00710  int ex(expo_gr(a)),fac,rest,n;
00711  double dbl;
00712 
00713  if (ex > c1)
00714  {
00715          if (p>=0)
00716                  if (p>c2)
00717                          times2pown(a,c2); // Produces an error
00718          else // 0 <= p <= 2100
00719          {  // By too great p-values overflow is possible!
00720                  dbl = _double(p);
00721                  n = (int) dbl;
00722                  fac = n/c3;
00723                  rest = n%c3;
00724                  for (int k=1; k<=fac; k++)
00725                          times2pown(a,c3);
00726                  times2pown(a,rest);
00727          }
00728          else // p<0; No overflow or underflow!
00729                  if (p<-c2) a = 0.0;
00730          else // -2100 <= p < 0
00731          {
00732                  dbl = _double(p);
00733                  n = (int) dbl;
00734                  fac = n/c4;
00735                  rest = n%c4;
00736                  for (int k=1; k<=fac; k++)
00737                          times2pown(a,c4);
00738                  times2pown(a,rest);
00739          }
00740  }
00741 } // Times2pown(...)
00742 
00743 
00744 inline bool zero_(const l_real& lr) throw()
00745 {  // returns only true if all lr.elem(i) == 0; Blomquist, 27.11.02; 
00746     int i=1, p=StagPrec(lr);
00747     bool tmp = true;
00748     do
00749     {
00750         if (sign(lr.elem(i))!=0) tmp = false;
00751         i++;
00752     }  while(tmp && i <= p );
00753     return tmp;
00754 }
00755 
00756 // real staggered constants:
00758 l_real Ln2_l_real()   throw();   // ln(2) 
00760 l_real Ln10_l_real()  throw();   // ln(10)
00762 l_real Ln10r_l_real() throw();   // 1/ln(10)
00764 l_real Pid4_l_real()  throw();   // Pi/4
00766 l_real Sqrt2_l_real() throw();   // sqrt(2)
00768 l_real Sqrt5_l_real() throw();   // sqrt(5)
00770 l_real Sqrt7_l_real() throw();   // sqrt(7)
00772 l_real Ln2r_l_real() throw();     // 1/ln(2)
00774 l_real Pi_l_real() throw();       // Pi
00776 l_real Pid2_l_real() throw();     // Pi/2
00778 l_real Pi2_l_real() throw();      // 2*Pi
00780 l_real Pid3_l_real() throw();     // Pi/3
00782 l_real Pir_l_real() throw();      // 1/Pi
00784 l_real Pi2r_l_real() throw();     // 1/(2*Pi)
00786 l_real SqrtPi_l_real() throw();   // sqrt(Pi)
00788 l_real Sqrt2Pi_l_real() throw();  // sqrt(2*Pi)
00790 l_real SqrtPir_l_real() throw();  // 1/sqrt(Pi)
00792 l_real Sqrt2Pir_l_real() throw(); // 1/sqrt(2*Pi)
00794 l_real Pip2_l_real() throw();     // Pi^2
00796 l_real Sqrt2r_l_real() throw();   // 1/sqrt(2)
00798 l_real Sqrt3_l_real() throw();    // sqrt(3)
00800 l_real Sqrt3d2_l_real() throw();  // sqrt(3)/2
00802 l_real Sqrt3r_l_real() throw();   // 1/sqrt(3)
00804 l_real LnPi_l_real() throw();     // ln(Pi)
00806 l_real Ln2Pi_l_real() throw();    // ln(2*Pi)
00808 l_real E_l_real() throw();        // e = exp(1)
00810 l_real Er_l_real() throw();       // 1/e
00812 l_real Ep2_l_real() throw();      // e^2
00814 l_real Ep2r_l_real() throw();     // 1/e^2
00816 l_real EpPi_l_real() throw();     // e^Pi
00818 l_real Ep2Pi_l_real() throw();    // e^(2*Pi)
00820 l_real EpPid2_l_real() throw();   // e^(Pi/2)
00822 l_real EpPid4_l_real() throw();   // e^(Pi/4)
00824 l_real EulerGa_l_real() throw();  // EulerGamma
00826 l_real Catalan_l_real() throw();  // Catalan
00827 
00828 } // namespace cxsc 
00829 
00830 #endif // _CXSC_L_REAL_HPP_INCLUDED
00831