C-XSC - A C++ Class Library for Extended Scientific Computing 2.5.4
complex.inl
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: complex.inl,v 1.26 2014/01/30 18:13:52 cxsc Exp $ */
25#include "cinterval.hpp"
26#include "idot.hpp"
27
28namespace cxsc {
29// ---- Constructors ----------------------------------------------
30
31inline complex & complex::operator= (const real & r) noexcept
32{
33 re=r;im=0;
34 return *this;
35}
36
37 // ---- Std.Operators ---------------------------------------
38inline complex operator -(const complex &a) noexcept
39{
40 return complex(-a.re,-a.im);
41}
42
43inline complex operator +(const complex &a) noexcept
44{
45 return a;
46}
47
48inline complex operator +(const complex &a,const complex &b) noexcept
49{
50 return complex(a.re+b.re,a.im+b.im);
51}
52
53inline complex operator -(const complex &a,const complex &b) noexcept
54{
55 return complex(a.re-b.re,a.im-b.im);
56}
57
58inline complex & operator +=(complex &a, const complex &b) noexcept { return a=a+b; }
59inline complex & operator -=(complex &a, const complex &b) noexcept { return a=a-b; }
60inline complex & operator *=(complex &a, const complex &b) noexcept { return a=a*b; }
61inline complex & operator /=(complex &a, const complex &b) noexcept { return a=a/b; }
62
63inline complex operator +(const complex & a,const real & b) noexcept
64{
65 return complex(a.re+b,a.im);
66}
67
68inline complex operator +(const real & a,const complex & b) noexcept
69{
70 return complex(a+b.re,b.im);
71}
72
73inline complex operator -(const complex & a,const real & b) noexcept
74{
75 return complex(a.re-b,a.im);
76}
77
78inline complex operator -(const real & a,const complex & b) noexcept
79{
80 return complex(a-b.re,-b.im);
81}
82
83inline complex operator *(const complex & a,const real & b) noexcept
84{
85// return a*_complex(b);
86 return complex(a.re*b,a.im*b); // Blomquist, 07.11.02;
87}
88
89inline complex operator *(const real & a,const complex & b) noexcept
90{
91// return _complex(a)*b;
92 return complex(a*b.re, a*b.im); // Blomquist, 07.11.02;
93}
94
95inline complex operator /(const complex & a,const real & b) noexcept
96{
97// return a/_complex(b);
98 return complex(a.re/b, a.im/b); // Blomquist, 07.11.02;
99}
100
101inline complex operator /(const real & a,const complex & b) noexcept
102{
103 return _complex(a)/b;
104}
105
106inline complex & operator +=(complex & a, const real & b) noexcept { return a=a+b; }
107inline complex & operator -=(complex & a, const real & b) noexcept { return a=a-b; }
108inline complex & operator *=(complex & a, const real & b) noexcept { return a=a*b; }
109inline complex & operator /=(complex & a, const real & b) noexcept { return a=a/b; }
110
111 // ---- Comp.Operat. ---------------------------------------
112inline bool operator! (const complex & a) noexcept { return !a.re && !a.im; }
113inline bool operator== (const complex & a, const complex & b) noexcept { return a.re==b.re && a.im==b.im; }
114inline bool operator!= (const complex & a, const complex & b) noexcept { return a.re!=b.re || a.im!=b.im; }
115inline bool operator== (const complex & a, const real & b) noexcept { return !a.im && a.re==b; }
116inline bool operator== (const real & a, const complex & b) noexcept { return !b.im && a==b.re; }
117inline bool operator!= (const complex & a, const real & b) noexcept { return !!a.im || a.re!=b; }
118inline bool operator!= (const real & a, const complex & b) noexcept { return !!b.im || a!=b.re; }
119
120 // ---- Others -------------------------------------------
121
122inline complex conj(const complex & a) noexcept { return complex(a.re,-a.im); }
123
124
125// ----------- Directed Rounding, Blomquist -------------------------------
126// ------------------------------------------------------------------------
127
128 // -------------------- addition --------------------------------
129
130inline complex addd(const complex& a, const complex& b) noexcept
131{ return complex(addd(a.re,b.re), addd(a.im,b.im)); }
132
133inline complex addu(const complex& a, const complex& b) noexcept
134{ return complex(addu(a.re,b.re), addu(a.im,b.im)); }
135
136inline complex addd(const complex& a, const real& b) noexcept
137{ return complex(addd(a.re,b), a.im); }
138
139inline complex addu(const complex& a, const real& b) noexcept
140{ return complex(addu(a.re,b), a.im); }
141
142inline complex addd(const real& a, const complex& b) noexcept
143{ return complex(addd(a,b.re), b.im); }
144
145inline complex addu(const real& a, const complex& b) noexcept
146{ return complex(addu(a,b.re), b.im); }
147 // ----------------- subtraction: ----------------------------
148
149inline complex subd(const complex& a, const complex& b) noexcept
150{ return complex(subd(a.re,b.re), subd(a.im,b.im)); }
151
152inline complex subu(const complex& a, const complex& b) noexcept
153{ return complex(subu(a.re,b.re), subu(a.im,b.im)); }
154
155inline complex subd(const complex& a, const real& b) noexcept
156{ return complex(subd(a.re,b), a.im); }
157
158inline complex subu(const complex& a, const real& b) noexcept
159{ return complex(subu(a.re,b), a.im); }
160
161inline complex subd(const real& a, const complex& b) noexcept
162{ return complex(subd(a,b.re), -b.im); }
163
164inline complex subu(const real& a, const complex& b) noexcept
165{ return complex(subu(a,b.re), -b.im); }
166
167 // --------------- multiplikation ------------------------
168
169inline complex muld(const complex &a, const real &b) noexcept
170{ return complex( muld(a.re,b), muld(a.im,b) ); }
171
172inline complex mulu(const complex &a, const real &b) noexcept
173{ return complex( mulu(a.re,b), mulu(a.im,b) ); }
174
175inline complex muld(const real &a, const complex &b) noexcept
176{ return complex( muld(a,b.re), muld(a,b.im) ); }
177
178inline complex mulu(const real &a, const complex &b) noexcept
179{ return complex( mulu(a,b.re), mulu(a,b.im) ); }
180
181 // -------------- division ---------------------------------
182
183inline complex divd(const complex &a, const real &b) noexcept
184{ return complex( divd(a.re,b), divd(a.im,b) ); }
185
186inline complex divu(const complex &a, const real &b) noexcept
187{ return complex( divu(a.re,b), divu(a.im,b) ); }
188
189inline complex divd(const real &a, const complex &b) noexcept
190{ return divd(_complex(a),b); }
191
192inline complex divu(const real &a, const complex &b) noexcept
193{ return divu(_complex(a),b); }
194
195inline complex operator *(const complex &a,const complex &b) noexcept
196{
197#ifdef CXSC_FAST_COMPLEX_OPERATIONS
198 return complex(Re(a)*Re(b)-Im(a)*Im(b), Re(a)*Im(b)+Im(a)*Re(b));
199#else
200 complex tmp;
201 dotprecision dot(0.0);
202
203 accumulate (dot, a.re, b.re);
204 accumulate (dot, -a.im, b.im);
205 rnd (dot, tmp.re, RND_NEXT);
206
207 dot = 0.0;
208 accumulate (dot, a.re, b.im);
209 accumulate (dot, a.im, b.re);
210 rnd (dot, tmp.im, RND_NEXT);
211
212 return tmp;
213#endif
214}
215
216inline complex muld(const complex &a, const complex &b) noexcept
217{ // Blomquist 07.11.02;
218 complex tmp;
219 dotprecision dot(0.0);
220
221 accumulate (dot, a.re, b.re);
222 accumulate (dot, -a.im, b.im);
223 rnd (dot, tmp.re, RND_DOWN);
224
225 dot = 0.0;
226 accumulate (dot, a.re, b.im);
227 accumulate (dot, a.im, b.re);
228 rnd (dot, tmp.im, RND_DOWN);
229
230 return tmp;
231}
232
233inline complex mulu(const complex &a, const complex &b) noexcept
234{ // Blomquist 07.11.02;
235 complex tmp;
236 dotprecision dot(0.0);
237
238 accumulate (dot, a.re, b.re);
239 accumulate (dot, -a.im, b.im);
240 rnd (dot, tmp.re, RND_UP);
241
242 dot = 0.0;
243 accumulate (dot, a.re, b.im);
244 accumulate (dot, a.im, b.re);
245 rnd (dot, tmp.im, RND_UP);
246
247 return tmp;
248}
249
250
251static const int Min_Exp_ = 1074, minexpom = -914,
252 maxexpo1 = 1022, MANT_W = 52;
253
254inline void product(real a, real b, real c, real d,
255 int& overfl, real& p1, interval& p2)
256// New version of function product(...) from Blomquist, 26.10.02;
257// Input data: a,b,c,d; Output data: overfl, p1, p2;
258// In case of overfl=0 the interval p1+p2 is an inclusion of a*b + c*d;
259// overfl=1 (overflow) signalizes that p1+p2 must be multiplied with 2^1074
260// to be an inclusion of a*b + c*d;
261// overfl=-1 (underflow) signalizes that p1+p2 must be multiplied with
262// 2^-1074 to be an inclusion of a*b + c*d;
263
264{
265 int exa, exb, exc, exd; // Exp. von a-d
266 dotprecision dot;
267 int inexact;
268
269 overfl = 0;
270 inexact = 0; // false
271
272 dot = 0.0;
273 exa = expo(a);
274 exb = expo(b);
275 exc = expo(c);
276 exd = expo(d);
277
278 if ( sign(a) == 0 || sign(b) == 0 ) // a * b == 0.0
279 if ( sign(c) == 0 || sign(d) == 0 ) // a * b == c * d == 0
280 // dot := #(0);
281 ; // No Operation necessary
282 else
283 {
284 // a * b == 0; c * d != 0;
285 if (exc+exd > maxexpo1)
286 {
287 // overflow !
288 if ( exc > exd ) c = comp( mant(c), exc-Min_Exp_ );
289 else d = comp( mant(d), exd-Min_Exp_ );
290 overfl = 1;
291 } else
292 if ( exc+exd < minexpom )
293 { // undeflow;
294 // c = comp( mant(c), exc+Min_Exp_ ); kann Overflow erzeugen!
295 if (exc < exd) c = comp( mant(c), exc+Min_Exp_ );
296 else d = comp( mant(d), exd+Min_Exp_ );
297 overfl = -1;
298 }
299 accumulate(dot,c,d);
300 } else // a,b != 0
301 if ( sign(c) == 0 || sign(d) == 0 )
302 {
303 // a*b != 0, c * d == 0
304 if (exa+exb > maxexpo1)
305 {
306 // overflow !
307 if ( exa > exb ) a = comp( mant(a), exa-Min_Exp_ );
308 else b = comp( mant(b), exb-Min_Exp_ );
309 overfl = 1;
310 } else
311 if (exa+exb < minexpom)
312 { // undeflow;
313 // a = comp( mant(a), exa+Min_Exp_ ); kann Overflow erzeugen!
314 if (exa < exb) a = comp( mant(a), exa+Min_Exp_ );
315 else b = comp( mant(b), exb+Min_Exp_ );
316 overfl = -1;
317 }
318 accumulate(dot,a,b);
319 } else
320 {
321 // a,b,c,d != 0
322 if (exa+exb > maxexpo1)
323 { // overflow bei a*b
324 if ( exa > exb ) a = comp( mant(a), exa-Min_Exp_ );
325 else b = comp( mant(b), exb-Min_Exp_ );
326 if (exc > MANT_W) c = comp( mant(c), exc-Min_Exp_ );
327 else if (exd > MANT_W)
328 d = comp( mant(d), exd-Min_Exp_ );
329 else
330 {
331 // underflow wegen Skalierung bei c*d
332 c = 0.0;
333 inexact = 1; // true
334 }
335 overfl = 1; // Hat vorher gefehlt!! Blomquist, 24.10.02;
336 } else if (exc+exd > maxexpo1)
337 {
338 // overflow bei c*d
339 if ( exc > exd ) c = comp( mant(c), exc-Min_Exp_ );
340 else d = comp( mant(d), exd-Min_Exp_ );
341 if (exa > MANT_W) a = comp( mant(a), exa-Min_Exp_ );
342 else if (exb > MANT_W)
343 b = comp( mant(b), exb-Min_Exp_ );
344 else
345 {
346 // underflow wegen Skalierung bei a*b
347 a = 0.0;
348 inexact = 1; // true
349 }
350 overfl = 1;
351 } else
352 if ( exa+exb < minexpom && exc+exd < minexpom )
353 {
354 // underflow bei a*b und bei c*d
355 if (exa < exb) a = comp( mant(a), exa+Min_Exp_ );
356 else b = comp( mant(b), exb+Min_Exp_ );
357 if (exc < exd) c = comp( mant(c), exc+Min_Exp_ );
358 else d = comp( mant(d), exd+Min_Exp_ );
359 overfl = -1;
360 }
361 accumulate(dot, a, b);
362 accumulate(dot, c, d);
363 }
364
365 p1 = rnd(dot);
366 dot -= p1;
367 rnd(dot,p2); // Blomquist, 07.11.02;
368
369 if (inexact)
370 p2 = interval( pred(Inf(p2)), succ(Sup(p2)) );
371
372} // end product
373
374inline real quotient (real z1, interval z2, real n1,
375 interval n2, int round, int zoverfl, int noverfl)
376// z1+z2 is an inclusion of a numerator.
377// n1+n2 is an inclusion of a denominator.
378// quotient(...) calculates with q1 an approximation of (z1+z2)/(n1+n2)
379// using staggered arithmetic.
380// Rounding with round (-1,0,+1) is considered.
381// zoverfl and noverfl are considered by scaling back with 2^1074 or 2^-1074
382// n1+n2 > 0 is assumed.
383{
384 real q1=0, scale;
385 interval q2, nh;
386 idotprecision id;
387 int vorz, anz_scale, ex = 0;
388
389 vorz = sign(z1); // n1,n2 > 0 is assumed!!
390 // so the sign of q1 ist sign of z1.
391 if ( zoverfl == -1 && noverfl == 1 )
392 {
393 // result in the undeflow range:
394 switch (round)
395 {
396 case RND_DOWN:
397 if (vorz >= 0)
398 q1 = 0.0;
399 else
400 q1 = -minreal; // Blomquist: MinReal --> minreal;
401 break;
402 case RND_NEXT:
403 q1 = 0.0;
404 break;
405 case RND_UP:
406 if (vorz <= 0)
407 q1 = 0.0;
408 else
409 q1 = minreal; // Blomquist: MinReal --> minreal;
410 break;
411 } // switch
412 } else if ( zoverfl==1 && noverfl==-1 )
413 {
414 // result in the overflow range:
415 if (vorz >= 0) q1 = MaxReal+MaxReal; // Overflow
416 else q1 = -MaxReal-MaxReal; // Overflow
417 } else
418 {
419 q1 = divd(z1, n1); // down, to get q2 >= 0
420 nh = interval(addd(n1, Inf(n2)),
421 addu(n1, Sup(n2)));
422
423 // q2:= ##( z1 - q1*n1 + z2 - q1*n2 );
424 id = z1;
425 accumulate(id, -q1, n1);
426 id += z2;
427 accumulate(id, -q1, n2);
428 q2 = rnd(id);
429
430 switch (round) // Considering the rounding before scaling
431 {
432 case RND_DOWN:
433 q1 = adddown(q1, divd(Inf(q2), Sup(nh)));
434 break;
435 case RND_NEXT:
436 q1 = q1 + (Inf(q2)+Sup(q2))*0.5/n1;
437 break;
438 case RND_UP:
439 q1 = addup(q1, divu(Sup(q2), Inf(nh)));
440 break;
441 } // switch
442
443
444 // scaling back, if numerator|denominator - over-|underflow :
445
446 // actuell as follows:
447 // q1:= comp( mant(q1), expo(q1) + (zoverfl-noverfl)*1074 );
448 // The scaling with 2^1074 must be done in two steps with 2^ex
449 // and with the factor scale:
450
451 anz_scale = zoverfl - noverfl; // |anz_scale| <= 1;
452 if (anz_scale > 0)
453 {
454 scale = comp(0.5, +1024);
455 ex = MANT_W-1;
456 }
457 else if (anz_scale < 0)
458 {
459 scale = comp(0.5, -1022);
460 ex = -MANT_W+1;
461 }
462 else scale = 1.0; // ex = 0 is already initialized for this case
463
464
465 if (ex) times2pown(q1,ex); // EXACT part scaling, if ex != 0.
466 switch (round) // correct rounding with the second factor scale:
467 {
468 case RND_DOWN:
469 q1 = multdown(q1, scale);
470 break;
471 case RND_NEXT:
472 q1 = q1 * scale;
473 break;
474 case RND_UP:
475 q1 = multup(q1, scale);
476 break;
477 } // switch
478 }
479 return q1;
480} // end of quotient
481
482inline complex _c_division(complex a, complex b, int round)
483{
484 if (0.0 == (sqr(Re(b))+sqr(Im(b)))) {
485 cxscthrow(DIV_BY_ZERO("complex operator / (const complex&,const complex&)"));
486 }
487
488 int zoverflow, noverflow;
489 real z1, n1;
490 interval z2, n2;
491 complex tmp;
492
493 product (Re(b), Re(b), Im(b), Im(b), noverflow, n1, n2);
494 product (Re(a), Re(b), Im(a), Im(b), zoverflow, z1, z2);
495 SetRe (tmp, quotient (z1, z2, n1, n2, round, zoverflow, noverflow));
496 product (Im(a), Re(b), -Re(a), Im(b), zoverflow, z1, z2);
497 SetIm (tmp, quotient (z1, z2, n1, n2, round, zoverflow, noverflow));
498 return tmp;
499}
500
501inline complex divn (const complex & a, const complex & b)
502{ // Blomquist: vorher c_divd(...), 07.11.02;
503 return _c_division(a, b, RND_NEXT);
504}
505
506inline complex divd (const complex & a, const complex & b)
507{ // Blomquist: vorher c_divd(...), 07.11.02;
508 return _c_division(a, b, RND_DOWN);
509}
510
511inline complex divu (const complex & a, const complex & b)
512{ // Blomquist: vorher c_divu(...), 07.11.02;
513 return _c_division(a, b, RND_UP);
514}
515
516inline complex operator / (const complex &a,const complex &b) noexcept
517{
518#ifdef CXSC_FAST_COMPLEX_OPERATIONS
519 real q = Re(b)*Re(b) + Im(b)*Im(b);
520 return complex((Re(a)*Re(b)+Im(a)*Im(b))/q, (Im(a)*Re(b)-Re(a)*Im(b))/q);
521#else
522 return divn(a,b);
523#endif
524}
525
526inline real abs2(const complex &a) noexcept
527{
528 dotprecision dot(0.0);
529 accumulate(dot,a.re,a.re);
530 accumulate(dot,a.im,a.im);
531 return rnd(dot);
532}
533
534inline real abs (complex z) noexcept
535{ // calculation of |z|; Blomquist 06.12.02;
536#ifdef CXSC_FAST_COMPLEX_OPERATIONS
537 return sqrt(Re(z)*Re(z)+Im(z)*Im(z));
538#else
539 return sqrtx2y2(Re(z),Im(z));
540#endif
541}
542
543
544} // namespace cxsc
The Scalar Type complex.
Definition complex.hpp:50
complex & operator=(const real &r) noexcept
Implementation of standard assigning operator.
Definition complex.inl:31
The Data Type dotprecision.
Definition dot.hpp:112
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
complex _complex(const real &a) noexcept
Deprecated typecast, which only exist for the reason of compatibility with older versions of C-XSC.
Definition complex.hpp:366
cdotprecision & operator+=(cdotprecision &cd, const l_complex &lc) noexcept
Implementation of standard algebraic addition and allocation operation.
Definition cdot.inl:251
const real MaxReal
Greatest representable floating-point number.
Definition real.cpp:65
civector operator/(const cimatrix_subv &rv, const cinterval &s) noexcept
Implementation of division operation.
Definition cimatrix.inl:730
const real minreal
Smallest positive denormalized representable floating-point number.
Definition real.cpp:63
cimatrix & operator*=(cimatrix &m, const cinterval &c) noexcept
Implementation of multiplication and allocation operation.
ivector abs(const cimatrix_subv &mv) noexcept
Returns the absolute value of the matrix.
Definition cimatrix.inl:737
cinterval sqrt(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:1007
void times2pown(cinterval &x, int n) noexcept
Fast multiplication of reference parameter [z] with .
Definition cimath.cpp:2059
cinterval sqr(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:3342
civector operator*(const cimatrix_subv &rv, const cinterval &s) noexcept
Implementation of multiplication operation.
Definition cimatrix.inl:731
interval sqrtx2y2(const interval &x, const interval &y) noexcept
Calculates .
Definition imath.cpp:80
cimatrix & operator/=(cimatrix &m, const cinterval &c) noexcept
Implementation of division and allocation operation.