C-XSC - A C++ Class Library for Extended Scientific Computing  2.5.4
dot.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: dot.cpp,v 1.28 2014/01/30 17:23:44 cxsc Exp $ */
00025 
00026 #include "dot.hpp"
00027 #include "real.hpp"
00028 #include "interval.hpp"
00029 
00030 #include "RtsTyp.h" 
00031 #include "RtsFunc.h"
00032 #include "ioflags.hpp"
00033 
00034 // Fuer memcpy... Kann auch <string.h> sein.
00035 #include <memory.h> 
00036 // #include <sstream>
00037 
00038 namespace cxsc {
00039 
00040 #ifdef CXSC_USE_TLS_PREC
00041 
00042 #ifdef _WIN32
00043 __declspec(thread) unsigned int opdotprec = 0;
00044 #else
00045 __thread unsigned int opdotprec = 0;
00046 #endif
00047 
00048 #else
00049 
00050 unsigned int opdotprec = 0;
00051 
00052 #endif
00053 
00054 //----------------------------------------------------------------------
00055 // global verf�gbare Dotprecision Variablen
00056 //
00057 //  dotakku[0..3] stehen fuer Matrix, Langzahl u.a. Pakete zur
00058 //                Verfuegung
00059 //  dotakku[4]    wird in den skalaren Paketen intern verwendet
00060 //dotprecision dotakku[MAXDOTAKKU];
00061 
00062 dotprecision::dotprecision(void) throw() : akku(new a_btyp[A_LENGTH]), err(0.0), k(0)
00063 {
00064    memset(akku,0,BUFFERSIZE); // d_clr(&akku);
00065    // d_clr(&akku); // Da d_init calloc verwendet, nicht noetig!?!?
00066    // d_clr  ist definiert in p88rts.h, also RtsTyp.h
00067 }
00068 
00069 dotprecision::dotprecision(const dotprecision &from) throw()
00070                                  : akku(new a_btyp[A_LENGTH]), err(from.err), k(from.k)
00071 {
00072    memcpy((void *)akku,(void *)from.akku,BUFFERSIZE);
00073 }
00074 
00075 dotprecision::dotprecision(const real &r) throw()
00076                                           : akku(new a_btyp[A_LENGTH]), err(0.0), k(0)
00077 {
00078    real dummy(r);
00079    memset(akku,0,BUFFERSIZE); // d_clr(&akku);
00080 
00081    d_radd(&akku,*((a_real *)&dummy));   
00082 }
00083 
00084 dotprecision & dotprecision::operator =(const dotprecision &from) throw()
00085 {
00086    if(&from==this) // Handle self-assignment
00087       return *this; // Hier kein new, bzw delete, deswegen diese Methode
00088 
00089    memcpy((void *)akku,(void *)from.akku,BUFFERSIZE);
00090    err = from.err;
00091    //k = from.k;
00092 
00093    return *this;
00094 }
00095 
00096 dotprecision & dotprecision::operator =(const real &r) throw()
00097 {
00098    real dummy(r);
00099    memset(akku,0,BUFFERSIZE); 
00100    d_radd(&akku,*((a_real *)&dummy));   // KRANK!
00101    err = 0.0;
00102    return *this;
00103 }
00104 
00105 dotprecision::~dotprecision()
00106 {
00107    delete [] akku;
00108 }
00109 
00110 dotprecision operator +(const dotprecision &a,const dotprecision &b) throw()
00111 {
00112    dotprecision tmp(a); 
00113    d_dadd(&(tmp.akku),(Dotprecision)b.akku);
00114    tmp.err = addu(a.err,b.err);
00115    //tmp.k = 0;
00116    return tmp;
00117 }
00118 
00119 dotprecision operator -(const dotprecision &a,const dotprecision &b) throw()
00120 {
00121    dotprecision tmp(b);
00122    tmp.negdot();
00123    d_dadd(&(tmp.akku),(Dotprecision)a.akku);
00124    tmp.err = addu(a.err,b.err);
00125    //tmp.k = 0;
00126    return tmp;
00127 }
00128 
00129 dotprecision operator +(const dotprecision &d,const real &r) throw()
00130 {
00131    dotprecision erg(d);
00132    d_radd(&(erg.akku),*((a_real *)&r));
00133    return erg;
00134 }
00135 
00136 dotprecision operator +(const real &r,const dotprecision &d) throw()
00137 {
00138    dotprecision erg(d);
00139    d_radd(&(erg.akku),*((a_real *)&r));
00140    return erg;
00141 }
00142 
00143 dotprecision operator -(const dotprecision &d,const real &r) throw()
00144 {
00145    dotprecision erg(d);
00146    d_radd(&(erg.akku),-*((a_real *)&r));
00147    return erg;
00148 }
00149 
00150 dotprecision operator -(const real &r,const dotprecision &d) throw()
00151 {
00152    dotprecision erg(d);
00153    erg.negdot();
00154    d_radd(&(erg.akku),*((a_real *)&r));
00155    return erg;
00156 }
00157 
00158 dotprecision & operator +=(dotprecision &d,const real &r) throw()
00159 {
00160    d_radd(&d.akku,*((a_real *)&r));
00161    return d;
00162 }
00163 
00164 dotprecision & operator +=(dotprecision &d,const dotprecision &d2) throw()
00165 {
00166    d_dadd(&(d.akku),(Dotprecision)(d2.akku));
00167    d.err = addu(d.err,d2.err);
00168    return d;
00169 }
00170 
00171 dotprecision & operator -=(dotprecision &d,const real &r) throw()
00172 {
00173    d_radd(&d.akku,-*((a_real *)&r));
00174    return d;
00175 }
00176 
00177 dotprecision & operator -=(dotprecision &d,const dotprecision &d2) throw()
00178 {
00179    dotprecision tmp(d2);
00180    tmp.negdot();
00181    d_dadd(&(d.akku),(Dotprecision)(tmp.akku));
00182    d.err = addu(d.err,d2.err);
00183    return d;
00184 }
00185 
00186 dotprecision operator -(const dotprecision &a) throw()
00187 {
00188    dotprecision tmp(a);
00189    tmp.negdot();
00190    return tmp;
00191 }
00192 
00193 dotprecision operator +(const dotprecision &a) throw()
00194 {
00195    return a;
00196 }
00197 
00198 // ---------------------------------------------------------------------------
00199 // Alle Vergleichsoperatoren werden auf die folgenden zwei (==, <=) zurueckgefuehrt!
00200 
00201 bool operator ==(const dotprecision &a,const dotprecision &b)  throw()  // Uebernommen aus dot.cpp "testequal"
00202 {
00203    int res = true;
00204 
00205    int leftsign = sign(a);
00206    if (leftsign != sign(b) )
00207       res = false;
00208    else 
00209    {
00210       int leftstart  = (a_intg)a.akku[A_BEGIN];
00211       int leftend    = (a_intg)a.akku[A_END];
00212       int rightstart = (a_intg)b.akku[A_BEGIN];
00213       int rightend   = (a_intg)b.akku[A_END];
00214 
00215       if (leftend < rightstart || rightend < leftstart) 
00216          res = false;
00217       else 
00218       {
00219          res = true;
00220          while (res && leftstart < rightstart)
00221             res = (a.akku[leftstart++] == ZERO);
00222 
00223          while (res && rightstart < leftstart)
00224             res = (b.akku[rightstart++] == ZERO);
00225 
00226          while (res && leftstart <= leftend && leftstart <= rightend) 
00227          {
00228             res = (a.akku[leftstart++] == b.akku[rightstart++]);
00229          }
00230 
00231          while (res && leftstart <= leftend)
00232             res = (a.akku[leftstart++] == ZERO);
00233          while (res && rightstart <= rightend)
00234             res = (b.akku[rightstart++] == ZERO);
00235       }
00236    } 
00237 
00238    return res && (a.err == b.err);
00239 }
00240 
00241 bool operator <=(const dotprecision &x,const dotprecision &y) throw()   // Uebernommen aus dot.cpp "testlessequal"
00242 {
00243    int res = true, cont;
00244    // Dotprecision dotakku = ((dotprecision*) &dot)->akku; b.akku
00245 
00246    dotprecision a = x + x.err;
00247    dotprecision b = y - y.err;
00248 
00249    int leftsign = sign(a);
00250    int rightsign = sign(b);
00251    if (leftsign != rightsign) 
00252       res = (leftsign < rightsign);
00253    else if (leftsign == 0) 
00254       res = true;
00255    else 
00256    {
00257       int leftstart  = (a_intg)a.akku[A_BEGIN];
00258       int leftend    = (a_intg)a.akku[A_END];
00259       int rightstart = (a_intg)b.akku[A_BEGIN];
00260       int rightend   = (a_intg)b.akku[A_END];
00261 
00262       if (leftend < rightstart) 
00263          res = (leftsign == -1);
00264       else if (rightend < leftstart) 
00265          res = (leftsign != -1);
00266       else 
00267       {
00268          cont = true;
00269 
00270          // ---------------------------------------------------------
00271          //  hat a Mantissenelemente vor b ==> false
00272 
00273          while (cont && leftstart < rightstart) 
00274          {
00275             cont = (a.akku[leftstart++] == ZERO);
00276             if (!cont) 
00277                res = false;
00278          }
00279          // ---------------------------------------------------------
00280          //  hat b Mantissenelemente vor a ==> true
00281 
00282          while (cont && rightstart < leftstart) 
00283          {
00284             cont = (b.akku[rightstart++] == ZERO);
00285             if (!cont) 
00286                res = true;
00287          }
00288 
00289          // ---------------------------------------------------------
00290          //  ein gemeinsames Mantissenelement a <= b  ?
00291 
00292          while (cont && leftstart <= leftend && leftstart <= rightend) 
00293          {
00294             cont = (a.akku[leftstart] == b.akku[rightstart]);
00295             if (!cont) 
00296             {
00297                res = (a.akku[leftstart] <= b.akku[rightstart]);
00298             }
00299             leftstart++,rightstart++;
00300          }
00301 
00302          // ---------------------------------------------------------
00303          //  hat a weitere Mantissenelemente ==> false
00304 
00305          while (cont && leftstart <= leftend) 
00306          {
00307             cont = (a.akku[leftstart++] == ZERO);
00308             if (!cont) 
00309                res = false;
00310          }
00311          // ---------------------------------------------------------
00312          //  hat b weitere Mantissenelemente ==> true
00313 
00314          while (cont && rightstart <= rightend) 
00315          {
00316             cont = (b.akku[rightstart++] == ZERO);
00317             if (!cont) 
00318                res = true;
00319          }
00320 
00321          if (cont)            // --------------------------------
00322             res = true;       // Mantissen waren gleich
00323          else                 // --------------------------------
00324             if (leftsign == -1)// Mantissen waren unterschiedlich
00325                res = !res;    // negiere Vergleichsresultat falls
00326                               // Akku's negativ waren
00327       }
00328    }
00329 
00330    return res;
00331 }
00332 
00333 bool operator !=(const dotprecision &a,const dotprecision &b) throw() { return !(a==b); }
00334 bool operator  <(const dotprecision &a,const dotprecision &b) throw() { return !(b<=a); }
00335 bool operator  >(const dotprecision &a,const dotprecision &b) throw() { return !(a<=b); }
00336 bool operator >=(const dotprecision &a,const dotprecision &b) throw() { return (b<=a);  }
00337 
00338 bool operator ==(const real &r,const dotprecision &d) throw() { return dotprecision(r)==d; }
00339 bool operator !=(const real &r,const dotprecision &d) throw() { return dotprecision(r)!=d; }
00340 bool operator  <(const real &r,const dotprecision &d) throw() { return dotprecision(r)<d;  }  
00341 bool operator  >(const real &r,const dotprecision &d) throw() { return dotprecision(r)>d;  }
00342 bool operator <=(const real &r,const dotprecision &d) throw() { return dotprecision(r)<=d; }
00343 bool operator >=(const real &r,const dotprecision &d) throw() { return dotprecision(r)>=d; }
00344 
00345 bool operator ==(const dotprecision &d,const real &r) throw() { return d==dotprecision(r); }
00346 bool operator !=(const dotprecision &d,const real &r) throw() { return d!=dotprecision(r); }
00347 bool operator  <(const dotprecision &d,const real &r) throw() { return d<dotprecision(r);  }
00348 bool operator  >(const dotprecision &d,const real &r) throw() { return d>dotprecision(r);  }
00349 bool operator <=(const dotprecision &d,const real &r) throw() { return d<=dotprecision(r); }
00350 bool operator >=(const dotprecision &d,const real &r) throw() { return d>=dotprecision(r); }
00351 
00352 bool operator  !(const dotprecision &a) throw()               { return sign(a)==0; }
00353 
00354 
00355 void rnd (const dotprecision& d, real& r, rndtype type) throw() 
00356 {
00357    if(type==RND_NEXT) {
00358       *((a_real *)(&r))=d_stan((Dotprecision)d.akku);
00359    } else if(type==RND_UP) {
00360       *((a_real *)(&r))=d_stau((Dotprecision)d.akku);
00361       r = addu(r,d.err);
00362    } else {
00363       *((a_real *)(&r))=d_stad((Dotprecision)d.akku); 
00364       r = subd(r,d.err);
00365    }
00366 }
00367 
00368 void rnd (const dotprecision& d, real& l, real& u) throw() 
00369 {
00370    *((a_real *)(&l))=d_stad(d.akku);
00371    *((a_real *)(&u))=d_stau(d.akku);
00372 
00373    l = subd(l,d.err);
00374    u = addu(u,d.err);
00375 }
00376 
00377 void rnd (const dotprecision& d, interval& x) throw() 
00378 {
00379     real a,b;
00380     rnd(d,a,b);
00381     x = interval(a,b);
00382 }
00383 
00384 real rnd (const dotprecision& d, rndtype type) throw() 
00385 {
00386     real r;
00387    if(type==RND_NEXT) {
00388       *((a_real *)(&r))=d_stan((Dotprecision)d.akku);
00389    } else if(type==RND_UP) {
00390       *((a_real *)(&r))=d_stau((Dotprecision)d.akku);
00391       r = addu(r,d.err);
00392    } else {
00393       *((a_real *)(&r))=d_stad((Dotprecision)d.akku); 
00394       r = subd(r,d.err);
00395    }
00396    return r;
00397 }
00398 
00399 /*std::string  & operator << (std::string  & s, const dotprecision & d) throw() 
00400 {
00401    //std::ostringstream o(s);
00402    
00403    if (ioflags.isset(IOFlags::realformat))
00404    {   
00405       
00406       real rl,ru; 
00407       rnd (d, rl,ru);    
00408       s += "dot(";
00409       //s << SaveOpt << RndDown;
00410       // o << rl;
00411       //sprintf (&sh[strlen(sh)] << rl, ", "); sh << RndUp;
00412       //
00413       s+=",";
00414       //      sprintf (&sh[strlen(sh)] << ru, ")");  sh << RestoreOpt;
00415       // o << ru;
00416       
00417       s+=")";
00418       //          strcpy (s, sh);
00419    } else
00420    {
00421       //real r=rnd(d);
00422 
00423       //o << r;
00424       s+="<zahl>";
00425              
00426       / *
00427       rndtype rnd;
00428       
00429       unsigned long length, formatflag, addblanks;
00430       unsigned long FracDigits = dotdigits;
00431       
00432       if (ioflags.isset(ioflags::rndup)) rnd = RND_UP;
00433       else if (ioflags.isset(ioflags::rnddown)) rnd = RND_DOWN;
00434       else rnd = RND_NEXT;
00435       
00436       if (ioflags.isset(IOFlags::variable))
00437          formatflag = dotwidth; 
00438       else if (ioflags.isset(IOFlags::varfixwidth))
00439          formatflag = dotwidth, FracDigits = -FracDigits;
00440       else
00441          formatflag = (ioflags.isset(IOFlags::fixed)) ? 0 : -1;
00442          
00443       // d_outp (str = dm, this->akku, formatflag, digits, rnd, &length);
00444       {
00445          long   dexpo,bdp,len;
00446          long   expo,i,digits,IntDigits,DecPlaces,vz;
00447          long   HoldWidth = (FracDigits < 0);
00448          
00449          if (HoldWidth) FracDigits = -FracDigits;
00450              
00451          // Kehre noetigenfalls Rundungsrichtung um
00452          if (d.sign() < 0) 
00453          {
00454               rnd = -rnd;
00455          }
00456          
00457          bdp = 2+A_I_DIGITS;
00458          len = bdp+1;
00459          if (c[A_END] > A_D_P) len += B_LENGTH * ((a_intg)c[A_END] - A_D_P);
00460 
00461       
00462       }
00463       * /
00464           
00465    }
00466    //s=o.str();                             
00467    return s;
00468 }
00469 
00470 
00471 std::ostream & operator << (std::ostream & o, const dotprecision & d) throw() 
00472 {
00473    std::string buffer;
00474    buffer << d;
00475    o << d;
00476    return o;
00477    
00478 / *
00479    // nur real-Ausgabe   
00480    real rl,ru;
00481    *((a_real *)(&rl))=d_stad(d.akku);
00482    *((a_real *)(&ru))=d_stau(d.akku);
00483 
00484    o << "dot(" << rl << "," << ru << ") " << ru-rl;
00485 * /
00486 }
00487 
00488 std::istream & operator >> (std::istream & i,dotprecision & d) throw() 
00489 {
00490    // Dies ist noch _schlecht_!
00491    double b;
00492    i >> b;
00493    d=b;
00494    return i;
00495 }
00496 */
00497 dotprecision & dotprecision::negdot(void) throw() 
00498 {
00499    this->akku[A_SIGN] = 1 - this->akku[A_SIGN];
00500    return *this;
00501 }
00502 
00503 int sign(const dotprecision &a) throw() 
00504 {
00505    if(a.akku[A_BEGIN]==ZERO)
00506       return 0;
00507    return (a.akku[A_SIGN] ? -1 : 1);
00508 }
00509 
00510 dotprecision abs(const dotprecision &a) throw() 
00511 {
00512    if(sign(a)>=0)
00513       return a;
00514    dotprecision tmp(a);
00515    tmp.negdot();
00516    return tmp;
00517 }
00518 
00519 dotprecision & accumulate (dotprecision& dot, const real& a, const real& b) throw() 
00520 {
00521    // nur dann accumulieren, wenn noetig
00522    if (!!a && !!b)
00523       d_padd (&dot.akku,*(a_real*)&a,*(a_real*)&b);
00524    return dot;  
00525 }
00526 
00527 } // namespace cxsc
00528