C-XSC - A C++ Class Library for Extended Scientific Computing  2.5.4
dot.cpp
1 /*
2 ** CXSC is a C++ library for eXtended Scientific Computing (V 2.5.4)
3 **
4 ** Copyright (C) 1990-2000 Institut fuer Angewandte Mathematik,
5 ** Universitaet Karlsruhe, Germany
6 ** (C) 2000-2014 Wiss. Rechnen/Softwaretechnologie
7 ** Universitaet Wuppertal, Germany
8 **
9 ** This library is free software; you can redistribute it and/or
10 ** modify it under the terms of the GNU Library General Public
11 ** License as published by the Free Software Foundation; either
12 ** version 2 of the License, or (at your option) any later version.
13 **
14 ** This library is distributed in the hope that it will be useful,
15 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
16 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 ** Library General Public License for more details.
18 **
19 ** You should have received a copy of the GNU Library General Public
20 ** License along with this library; if not, write to the Free
21 ** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
22 */
23 
24 /* CVS $Id: dot.cpp,v 1.28 2014/01/30 17:23:44 cxsc Exp $ */
25 
26 #include "dot.hpp"
27 #include "real.hpp"
28 #include "interval.hpp"
29 
30 #include "RtsTyp.h"
31 #include "RtsFunc.h"
32 #include "ioflags.hpp"
33 
34 // Fuer memcpy... Kann auch <string.h> sein.
35 #include <memory.h>
36 // #include <sstream>
37 
38 namespace cxsc {
39 
40 #ifdef CXSC_USE_TLS_PREC
41 
42 #ifdef _WIN32
43 __declspec(thread) unsigned int opdotprec = 0;
44 #else
45 __thread unsigned int opdotprec = 0;
46 #endif
47 
48 #else
49 
50 unsigned int opdotprec = 0;
51 
52 #endif
53 
54 //----------------------------------------------------------------------
55 // global verf�gbare Dotprecision Variablen
56 //
57 // dotakku[0..3] stehen fuer Matrix, Langzahl u.a. Pakete zur
58 // Verfuegung
59 // dotakku[4] wird in den skalaren Paketen intern verwendet
60 //dotprecision dotakku[MAXDOTAKKU];
61 
62 dotprecision::dotprecision(void) throw() : akku(new a_btyp[A_LENGTH]), err(0.0), k(0)
63 {
64  memset(akku,0,BUFFERSIZE); // d_clr(&akku);
65  // d_clr(&akku); // Da d_init calloc verwendet, nicht noetig!?!?
66  // d_clr ist definiert in p88rts.h, also RtsTyp.h
67 }
68 
70  : akku(new a_btyp[A_LENGTH]), err(from.err), k(from.k)
71 {
72  memcpy((void *)akku,(void *)from.akku,BUFFERSIZE);
73 }
74 
76  : akku(new a_btyp[A_LENGTH]), err(0.0), k(0)
77 {
78  real dummy(r);
79  memset(akku,0,BUFFERSIZE); // d_clr(&akku);
80 
81  d_radd(&akku,*((a_real *)&dummy));
82 }
83 
85 {
86  if(&from==this) // Handle self-assignment
87  return *this; // Hier kein new, bzw delete, deswegen diese Methode
88 
89  memcpy((void *)akku,(void *)from.akku,BUFFERSIZE);
90  err = from.err;
91  //k = from.k;
92 
93  return *this;
94 }
95 
97 {
98  real dummy(r);
99  memset(akku,0,BUFFERSIZE);
100  d_radd(&akku,*((a_real *)&dummy)); // KRANK!
101  err = 0.0;
102  return *this;
103 }
104 
105 dotprecision::~dotprecision()
106 {
107  delete [] akku;
108 }
109 
110 dotprecision operator +(const dotprecision &a,const dotprecision &b) throw()
111 {
112  dotprecision tmp(a);
113  d_dadd(&(tmp.akku),(Dotprecision)b.akku);
114  tmp.err = addu(a.err,b.err);
115  //tmp.k = 0;
116  return tmp;
117 }
118 
119 dotprecision operator -(const dotprecision &a,const dotprecision &b) throw()
120 {
121  dotprecision tmp(b);
122  tmp.negdot();
123  d_dadd(&(tmp.akku),(Dotprecision)a.akku);
124  tmp.err = addu(a.err,b.err);
125  //tmp.k = 0;
126  return tmp;
127 }
128 
129 dotprecision operator +(const dotprecision &d,const real &r) throw()
130 {
131  dotprecision erg(d);
132  d_radd(&(erg.akku),*((a_real *)&r));
133  return erg;
134 }
135 
136 dotprecision operator +(const real &r,const dotprecision &d) throw()
137 {
138  dotprecision erg(d);
139  d_radd(&(erg.akku),*((a_real *)&r));
140  return erg;
141 }
142 
143 dotprecision operator -(const dotprecision &d,const real &r) throw()
144 {
145  dotprecision erg(d);
146  d_radd(&(erg.akku),-*((a_real *)&r));
147  return erg;
148 }
149 
150 dotprecision operator -(const real &r,const dotprecision &d) throw()
151 {
152  dotprecision erg(d);
153  erg.negdot();
154  d_radd(&(erg.akku),*((a_real *)&r));
155  return erg;
156 }
157 
159 {
160  d_radd(&d.akku,*((a_real *)&r));
161  return d;
162 }
163 
165 {
166  d_dadd(&(d.akku),(Dotprecision)(d2.akku));
167  d.err = addu(d.err,d2.err);
168  return d;
169 }
170 
171 dotprecision & operator -=(dotprecision &d,const real &r) throw()
172 {
173  d_radd(&d.akku,-*((a_real *)&r));
174  return d;
175 }
176 
177 dotprecision & operator -=(dotprecision &d,const dotprecision &d2) throw()
178 {
179  dotprecision tmp(d2);
180  tmp.negdot();
181  d_dadd(&(d.akku),(Dotprecision)(tmp.akku));
182  d.err = addu(d.err,d2.err);
183  return d;
184 }
185 
186 dotprecision operator -(const dotprecision &a) throw()
187 {
188  dotprecision tmp(a);
189  tmp.negdot();
190  return tmp;
191 }
192 
193 dotprecision operator +(const dotprecision &a) throw()
194 {
195  return a;
196 }
197 
198 // ---------------------------------------------------------------------------
199 // Alle Vergleichsoperatoren werden auf die folgenden zwei (==, <=) zurueckgefuehrt!
200 
201 bool operator ==(const dotprecision &a,const dotprecision &b) throw() // Uebernommen aus dot.cpp "testequal"
202 {
203  int res = true;
204 
205  int leftsign = sign(a);
206  if (leftsign != sign(b) )
207  res = false;
208  else
209  {
210  int leftstart = (a_intg)a.akku[A_BEGIN];
211  int leftend = (a_intg)a.akku[A_END];
212  int rightstart = (a_intg)b.akku[A_BEGIN];
213  int rightend = (a_intg)b.akku[A_END];
214 
215  if (leftend < rightstart || rightend < leftstart)
216  res = false;
217  else
218  {
219  res = true;
220  while (res && leftstart < rightstart)
221  res = (a.akku[leftstart++] == ZERO);
222 
223  while (res && rightstart < leftstart)
224  res = (b.akku[rightstart++] == ZERO);
225 
226  while (res && leftstart <= leftend && leftstart <= rightend)
227  {
228  res = (a.akku[leftstart++] == b.akku[rightstart++]);
229  }
230 
231  while (res && leftstart <= leftend)
232  res = (a.akku[leftstart++] == ZERO);
233  while (res && rightstart <= rightend)
234  res = (b.akku[rightstart++] == ZERO);
235  }
236  }
237 
238  return res && (a.err == b.err);
239 }
240 
241 bool operator <=(const dotprecision &x,const dotprecision &y) throw() // Uebernommen aus dot.cpp "testlessequal"
242 {
243  int res = true, cont;
244  // Dotprecision dotakku = ((dotprecision*) &dot)->akku; b.akku
245 
246  dotprecision a = x + x.err;
247  dotprecision b = y - y.err;
248 
249  int leftsign = sign(a);
250  int rightsign = sign(b);
251  if (leftsign != rightsign)
252  res = (leftsign < rightsign);
253  else if (leftsign == 0)
254  res = true;
255  else
256  {
257  int leftstart = (a_intg)a.akku[A_BEGIN];
258  int leftend = (a_intg)a.akku[A_END];
259  int rightstart = (a_intg)b.akku[A_BEGIN];
260  int rightend = (a_intg)b.akku[A_END];
261 
262  if (leftend < rightstart)
263  res = (leftsign == -1);
264  else if (rightend < leftstart)
265  res = (leftsign != -1);
266  else
267  {
268  cont = true;
269 
270  // ---------------------------------------------------------
271  // hat a Mantissenelemente vor b ==> false
272 
273  while (cont && leftstart < rightstart)
274  {
275  cont = (a.akku[leftstart++] == ZERO);
276  if (!cont)
277  res = false;
278  }
279  // ---------------------------------------------------------
280  // hat b Mantissenelemente vor a ==> true
281 
282  while (cont && rightstart < leftstart)
283  {
284  cont = (b.akku[rightstart++] == ZERO);
285  if (!cont)
286  res = true;
287  }
288 
289  // ---------------------------------------------------------
290  // ein gemeinsames Mantissenelement a <= b ?
291 
292  while (cont && leftstart <= leftend && leftstart <= rightend)
293  {
294  cont = (a.akku[leftstart] == b.akku[rightstart]);
295  if (!cont)
296  {
297  res = (a.akku[leftstart] <= b.akku[rightstart]);
298  }
299  leftstart++,rightstart++;
300  }
301 
302  // ---------------------------------------------------------
303  // hat a weitere Mantissenelemente ==> false
304 
305  while (cont && leftstart <= leftend)
306  {
307  cont = (a.akku[leftstart++] == ZERO);
308  if (!cont)
309  res = false;
310  }
311  // ---------------------------------------------------------
312  // hat b weitere Mantissenelemente ==> true
313 
314  while (cont && rightstart <= rightend)
315  {
316  cont = (b.akku[rightstart++] == ZERO);
317  if (!cont)
318  res = true;
319  }
320 
321  if (cont) // --------------------------------
322  res = true; // Mantissen waren gleich
323  else // --------------------------------
324  if (leftsign == -1)// Mantissen waren unterschiedlich
325  res = !res; // negiere Vergleichsresultat falls
326  // Akku's negativ waren
327  }
328  }
329 
330  return res;
331 }
332 
333 bool operator !=(const dotprecision &a,const dotprecision &b) throw() { return !(a==b); }
334 bool operator <(const dotprecision &a,const dotprecision &b) throw() { return !(b<=a); }
335 bool operator >(const dotprecision &a,const dotprecision &b) throw() { return !(a<=b); }
336 bool operator >=(const dotprecision &a,const dotprecision &b) throw() { return (b<=a); }
337 
338 bool operator ==(const real &r,const dotprecision &d) throw() { return dotprecision(r)==d; }
339 bool operator !=(const real &r,const dotprecision &d) throw() { return dotprecision(r)!=d; }
340 bool operator <(const real &r,const dotprecision &d) throw() { return dotprecision(r)<d; }
341 bool operator >(const real &r,const dotprecision &d) throw() { return dotprecision(r)>d; }
342 bool operator <=(const real &r,const dotprecision &d) throw() { return dotprecision(r)<=d; }
343 bool operator >=(const real &r,const dotprecision &d) throw() { return dotprecision(r)>=d; }
344 
345 bool operator ==(const dotprecision &d,const real &r) throw() { return d==dotprecision(r); }
346 bool operator !=(const dotprecision &d,const real &r) throw() { return d!=dotprecision(r); }
347 bool operator <(const dotprecision &d,const real &r) throw() { return d<dotprecision(r); }
348 bool operator >(const dotprecision &d,const real &r) throw() { return d>dotprecision(r); }
349 bool operator <=(const dotprecision &d,const real &r) throw() { return d<=dotprecision(r); }
350 bool operator >=(const dotprecision &d,const real &r) throw() { return d>=dotprecision(r); }
351 
352 bool operator !(const dotprecision &a) throw() { return sign(a)==0; }
353 
354 
355 void rnd (const dotprecision& d, real& r, rndtype type) throw()
356 {
357  if(type==RND_NEXT) {
358  *((a_real *)(&r))=d_stan((Dotprecision)d.akku);
359  } else if(type==RND_UP) {
360  *((a_real *)(&r))=d_stau((Dotprecision)d.akku);
361  r = addu(r,d.err);
362  } else {
363  *((a_real *)(&r))=d_stad((Dotprecision)d.akku);
364  r = subd(r,d.err);
365  }
366 }
367 
368 void rnd (const dotprecision& d, real& l, real& u) throw()
369 {
370  *((a_real *)(&l))=d_stad(d.akku);
371  *((a_real *)(&u))=d_stau(d.akku);
372 
373  l = subd(l,d.err);
374  u = addu(u,d.err);
375 }
376 
377 void rnd (const dotprecision& d, interval& x) throw()
378 {
379  real a,b;
380  rnd(d,a,b);
381  x = interval(a,b);
382 }
383 
384 real rnd (const dotprecision& d, rndtype type) throw()
385 {
386  real r;
387  if(type==RND_NEXT) {
388  *((a_real *)(&r))=d_stan((Dotprecision)d.akku);
389  } else if(type==RND_UP) {
390  *((a_real *)(&r))=d_stau((Dotprecision)d.akku);
391  r = addu(r,d.err);
392  } else {
393  *((a_real *)(&r))=d_stad((Dotprecision)d.akku);
394  r = subd(r,d.err);
395  }
396  return r;
397 }
398 
399 /*std::string & operator << (std::string & s, const dotprecision & d) throw()
400 {
401  //std::ostringstream o(s);
402 
403  if (ioflags.isset(IOFlags::realformat))
404  {
405 
406  real rl,ru;
407  rnd (d, rl,ru);
408  s += "dot(";
409  //s << SaveOpt << RndDown;
410  // o << rl;
411  //sprintf (&sh[strlen(sh)] << rl, ", "); sh << RndUp;
412  //
413  s+=",";
414  // sprintf (&sh[strlen(sh)] << ru, ")"); sh << RestoreOpt;
415  // o << ru;
416 
417  s+=")";
418  // strcpy (s, sh);
419  } else
420  {
421  //real r=rnd(d);
422 
423  //o << r;
424  s+="<zahl>";
425 
426  / *
427  rndtype rnd;
428 
429  unsigned long length, formatflag, addblanks;
430  unsigned long FracDigits = dotdigits;
431 
432  if (ioflags.isset(ioflags::rndup)) rnd = RND_UP;
433  else if (ioflags.isset(ioflags::rnddown)) rnd = RND_DOWN;
434  else rnd = RND_NEXT;
435 
436  if (ioflags.isset(IOFlags::variable))
437  formatflag = dotwidth;
438  else if (ioflags.isset(IOFlags::varfixwidth))
439  formatflag = dotwidth, FracDigits = -FracDigits;
440  else
441  formatflag = (ioflags.isset(IOFlags::fixed)) ? 0 : -1;
442 
443  // d_outp (str = dm, this->akku, formatflag, digits, rnd, &length);
444  {
445  long dexpo,bdp,len;
446  long expo,i,digits,IntDigits,DecPlaces,vz;
447  long HoldWidth = (FracDigits < 0);
448 
449  if (HoldWidth) FracDigits = -FracDigits;
450 
451  // Kehre noetigenfalls Rundungsrichtung um
452  if (d.sign() < 0)
453  {
454  rnd = -rnd;
455  }
456 
457  bdp = 2+A_I_DIGITS;
458  len = bdp+1;
459  if (c[A_END] > A_D_P) len += B_LENGTH * ((a_intg)c[A_END] - A_D_P);
460 
461 
462  }
463  * /
464 
465  }
466  //s=o.str();
467  return s;
468 }
469 
470 
471 std::ostream & operator << (std::ostream & o, const dotprecision & d) throw()
472 {
473  std::string buffer;
474  buffer << d;
475  o << d;
476  return o;
477 
478 / *
479  // nur real-Ausgabe
480  real rl,ru;
481  *((a_real *)(&rl))=d_stad(d.akku);
482  *((a_real *)(&ru))=d_stau(d.akku);
483 
484  o << "dot(" << rl << "," << ru << ") " << ru-rl;
485 * /
486 }
487 
488 std::istream & operator >> (std::istream & i,dotprecision & d) throw()
489 {
490  // Dies ist noch _schlecht_!
491  double b;
492  i >> b;
493  d=b;
494  return i;
495 }
496 */
497 dotprecision & dotprecision::negdot(void) throw()
498 {
499  this->akku[A_SIGN] = 1 - this->akku[A_SIGN];
500  return *this;
501 }
502 
503 int sign(const dotprecision &a) throw()
504 {
505  if(a.akku[A_BEGIN]==ZERO)
506  return 0;
507  return (a.akku[A_SIGN] ? -1 : 1);
508 }
509 
510 dotprecision abs(const dotprecision &a) throw()
511 {
512  if(sign(a)>=0)
513  return a;
514  dotprecision tmp(a);
515  tmp.negdot();
516  return tmp;
517 }
518 
519 dotprecision & accumulate (dotprecision& dot, const real& a, const real& b) throw()
520 {
521  // nur dann accumulieren, wenn noetig
522  if (!!a && !!b)
523  d_padd (&dot.akku,*(a_real*)&a,*(a_real*)&b);
524  return dot;
525 }
526 
527 } // namespace cxsc
528 
cxsc::dotprecision::operator=
dotprecision & operator=(const dotprecision &)
Implementation of standard assigning operator.
Definition: dot.cpp:84
cxsc::interval
The Scalar Type interval.
Definition: interval.hpp:54
cxsc::abs
ivector abs(const cimatrix_subv &mv)
Returns the absolute value of the matrix.
Definition: cimatrix.inl:737
cxsc::dotprecision
The Data Type dotprecision.
Definition: dot.hpp:111
cxsc
The namespace cxsc, providing all functionality of the class library C-XSC.
Definition: cdot.cpp:29
cxsc::operator+=
cdotprecision & operator+=(cdotprecision &cd, const l_complex &lc)
Implementation of standard algebraic addition and allocation operation.
Definition: cdot.inl:251
cxsc::dotprecision::dotprecision
dotprecision(void)
Constructor of class dotprecision.
Definition: dot.cpp:62
cxsc::real
The Scalar Type real.
Definition: real.hpp:113