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