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
38namespace 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
50unsigned 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
62dotprecision::dotprecision(void) noexcept : 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
105dotprecision::~dotprecision()
106{
107 delete [] akku;
108}
109
110dotprecision operator +(const dotprecision &a,const dotprecision &b) noexcept
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
119dotprecision operator -(const dotprecision &a,const dotprecision &b) noexcept
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
129dotprecision operator +(const dotprecision &d,const real &r) noexcept
130{
131 dotprecision erg(d);
132 d_radd(&(erg.akku),*((a_real *)&r));
133 return erg;
134}
135
136dotprecision operator +(const real &r,const dotprecision &d) noexcept
137{
138 dotprecision erg(d);
139 d_radd(&(erg.akku),*((a_real *)&r));
140 return erg;
141}
142
143dotprecision operator -(const dotprecision &d,const real &r) noexcept
144{
145 dotprecision erg(d);
146 d_radd(&(erg.akku),-*((a_real *)&r));
147 return erg;
148}
149
150dotprecision operator -(const real &r,const dotprecision &d) noexcept
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
171dotprecision & operator -=(dotprecision &d,const real &r) noexcept
172{
173 d_radd(&d.akku,-*((a_real *)&r));
174 return d;
175}
176
177dotprecision & operator -=(dotprecision &d,const dotprecision &d2) noexcept
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
186dotprecision operator -(const dotprecision &a) noexcept
187{
188 dotprecision tmp(a);
189 tmp.negdot();
190 return tmp;
191}
192
193dotprecision operator +(const dotprecision &a) noexcept
194{
195 return a;
196}
197
198// ---------------------------------------------------------------------------
199// Alle Vergleichsoperatoren werden auf die folgenden zwei (==, <=) zurueckgefuehrt!
200
201bool operator ==(const dotprecision &a,const dotprecision &b) noexcept // 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
241bool operator <=(const dotprecision &x,const dotprecision &y) noexcept // 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
333bool operator !=(const dotprecision &a,const dotprecision &b) noexcept { return !(a==b); }
334bool operator <(const dotprecision &a,const dotprecision &b) noexcept { return !(b<=a); }
335bool operator >(const dotprecision &a,const dotprecision &b) noexcept { return !(a<=b); }
336bool operator >=(const dotprecision &a,const dotprecision &b) noexcept { return (b<=a); }
337
338bool operator ==(const real &r,const dotprecision &d) noexcept { return dotprecision(r)==d; }
339bool operator !=(const real &r,const dotprecision &d) noexcept { return dotprecision(r)!=d; }
340bool operator <(const real &r,const dotprecision &d) noexcept { return dotprecision(r)<d; }
341bool operator >(const real &r,const dotprecision &d) noexcept { return dotprecision(r)>d; }
342bool operator <=(const real &r,const dotprecision &d) noexcept { return dotprecision(r)<=d; }
343bool operator >=(const real &r,const dotprecision &d) noexcept { return dotprecision(r)>=d; }
344
345bool operator ==(const dotprecision &d,const real &r) noexcept { return d==dotprecision(r); }
346bool operator !=(const dotprecision &d,const real &r) noexcept { return d!=dotprecision(r); }
347bool operator <(const dotprecision &d,const real &r) noexcept { return d<dotprecision(r); }
348bool operator >(const dotprecision &d,const real &r) noexcept { return d>dotprecision(r); }
349bool operator <=(const dotprecision &d,const real &r) noexcept { return d<=dotprecision(r); }
350bool operator >=(const dotprecision &d,const real &r) noexcept { return d>=dotprecision(r); }
351
352bool operator !(const dotprecision &a) noexcept { return sign(a)==0; }
353
354
355void rnd (const dotprecision& d, real& r, rndtype type) noexcept
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
368void rnd (const dotprecision& d, real& l, real& u) noexcept
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
377void rnd (const dotprecision& d, interval& x) noexcept
378{
379 real a,b;
380 rnd(d,a,b);
381 x = interval(a,b);
382}
383
384real rnd (const dotprecision& d, rndtype type) noexcept
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) noexcept
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
471std::ostream & operator << (std::ostream & o, const dotprecision & d) noexcept
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
488std::istream & operator >> (std::istream & i,dotprecision & d) noexcept
489{
490 // Dies ist noch _schlecht_!
491 double b;
492 i >> b;
493 d=b;
494 return i;
495}
496*/
497dotprecision & dotprecision::negdot(void) noexcept
498{
499 this->akku[A_SIGN] = 1 - this->akku[A_SIGN];
500 return *this;
501}
502
503int sign(const dotprecision &a) noexcept
504{
505 if(a.akku[A_BEGIN]==ZERO)
506 return 0;
507 return (a.akku[A_SIGN] ? -1 : 1);
508}
509
510dotprecision abs(const dotprecision &a) noexcept
511{
512 if(sign(a)>=0)
513 return a;
514 dotprecision tmp(a);
515 tmp.negdot();
516 return tmp;
517}
518
519dotprecision & accumulate (dotprecision& dot, const real& a, const real& b) noexcept
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
The Data Type dotprecision.
Definition dot.hpp:112
dotprecision(void) noexcept
Constructor of class dotprecision.
Definition dot.cpp:62
dotprecision & operator=(const dotprecision &) noexcept
Implementation of standard assigning operator.
Definition dot.cpp:84
The Scalar Type interval.
Definition interval.hpp:55
The Scalar Type real.
Definition real.hpp:114
The namespace cxsc, providing all functionality of the class library C-XSC.
Definition cdot.cpp:29
cdotprecision & operator+=(cdotprecision &cd, const l_complex &lc) noexcept
Implementation of standard algebraic addition and allocation operation.
Definition cdot.inl:251
ivector abs(const cimatrix_subv &mv) noexcept
Returns the absolute value of the matrix.
Definition cimatrix.inl:737