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 
28 namespace cxsc {
29 // ---- Constructors ----------------------------------------------
30 
31 inline complex & complex::operator= (const real & r) throw()
32 {
33  re=r;im=0;
34  return *this;
35 }
36 
37  // ---- Std.Operators ---------------------------------------
38 inline complex operator -(const complex &a) throw ()
39 {
40  return complex(-a.re,-a.im);
41 }
42 
43 inline complex operator +(const complex &a) throw ()
44 {
45  return a;
46 }
47 
48 inline complex operator +(const complex &a,const complex &b) throw()
49 {
50  return complex(a.re+b.re,a.im+b.im);
51 }
52 
53 inline complex operator -(const complex &a,const complex &b) throw()
54 {
55  return complex(a.re-b.re,a.im-b.im);
56 }
57 
58 inline complex & operator +=(complex &a, const complex &b) throw() { return a=a+b; }
59 inline complex & operator -=(complex &a, const complex &b) throw() { return a=a-b; }
60 inline complex & operator *=(complex &a, const complex &b) throw() { return a=a*b; }
61 inline complex & operator /=(complex &a, const complex &b) throw() { return a=a/b; }
62 
63 inline complex operator +(const complex & a,const real & b) throw()
64 {
65  return complex(a.re+b,a.im);
66 }
67 
68 inline complex operator +(const real & a,const complex & b) throw()
69 {
70  return complex(a+b.re,b.im);
71 }
72 
73 inline complex operator -(const complex & a,const real & b) throw()
74 {
75  return complex(a.re-b,a.im);
76 }
77 
78 inline complex operator -(const real & a,const complex & b) throw()
79 {
80  return complex(a-b.re,-b.im);
81 }
82 
83 inline complex operator *(const complex & a,const real & b) throw()
84 {
85 // return a*_complex(b);
86  return complex(a.re*b,a.im*b); // Blomquist, 07.11.02;
87 }
88 
89 inline complex operator *(const real & a,const complex & b) throw()
90 {
91 // return _complex(a)*b;
92  return complex(a*b.re, a*b.im); // Blomquist, 07.11.02;
93 }
94 
95 inline complex operator /(const complex & a,const real & b) throw()
96 {
97 // return a/_complex(b);
98  return complex(a.re/b, a.im/b); // Blomquist, 07.11.02;
99 }
100 
101 inline complex operator /(const real & a,const complex & b) throw()
102 {
103  return _complex(a)/b;
104 }
105 
106 inline complex & operator +=(complex & a, const real & b) throw() { return a=a+b; }
107 inline complex & operator -=(complex & a, const real & b) throw() { return a=a-b; }
108 inline complex & operator *=(complex & a, const real & b) throw() { return a=a*b; }
109 inline complex & operator /=(complex & a, const real & b) throw() { return a=a/b; }
110 
111  // ---- Comp.Operat. ---------------------------------------
112 inline bool operator! (const complex & a) throw() { return !a.re && !a.im; }
113 inline bool operator== (const complex & a, const complex & b) throw() { return a.re==b.re && a.im==b.im; }
114 inline bool operator!= (const complex & a, const complex & b) throw() { return a.re!=b.re || a.im!=b.im; }
115 inline bool operator== (const complex & a, const real & b) throw() { return !a.im && a.re==b; }
116 inline bool operator== (const real & a, const complex & b) throw() { return !b.im && a==b.re; }
117 inline bool operator!= (const complex & a, const real & b) throw() { return !!a.im || a.re!=b; }
118 inline bool operator!= (const real & a, const complex & b) throw() { return !!b.im || a!=b.re; }
119 
120  // ---- Others -------------------------------------------
121 
122 inline complex conj(const complex & a) throw() { return complex(a.re,-a.im); }
123 
124 
125 // ----------- Directed Rounding, Blomquist -------------------------------
126 // ------------------------------------------------------------------------
127 
128  // -------------------- addition --------------------------------
129 
130 inline complex addd(const complex& a, const complex& b) throw()
131 { return complex(addd(a.re,b.re), addd(a.im,b.im)); }
132 
133 inline complex addu(const complex& a, const complex& b) throw()
134 { return complex(addu(a.re,b.re), addu(a.im,b.im)); }
135 
136 inline complex addd(const complex& a, const real& b) throw()
137 { return complex(addd(a.re,b), a.im); }
138 
139 inline complex addu(const complex& a, const real& b) throw()
140 { return complex(addu(a.re,b), a.im); }
141 
142 inline complex addd(const real& a, const complex& b) throw()
143 { return complex(addd(a,b.re), b.im); }
144 
145 inline complex addu(const real& a, const complex& b) throw()
146 { return complex(addu(a,b.re), b.im); }
147  // ----------------- subtraction: ----------------------------
148 
149 inline complex subd(const complex& a, const complex& b) throw()
150 { return complex(subd(a.re,b.re), subd(a.im,b.im)); }
151 
152 inline complex subu(const complex& a, const complex& b) throw()
153 { return complex(subu(a.re,b.re), subu(a.im,b.im)); }
154 
155 inline complex subd(const complex& a, const real& b) throw()
156 { return complex(subd(a.re,b), a.im); }
157 
158 inline complex subu(const complex& a, const real& b) throw()
159 { return complex(subu(a.re,b), a.im); }
160 
161 inline complex subd(const real& a, const complex& b) throw()
162 { return complex(subd(a,b.re), -b.im); }
163 
164 inline complex subu(const real& a, const complex& b) throw()
165 { return complex(subu(a,b.re), -b.im); }
166 
167  // --------------- multiplikation ------------------------
168 
169 inline complex muld(const complex &a, const real &b) throw()
170 { return complex( muld(a.re,b), muld(a.im,b) ); }
171 
172 inline complex mulu(const complex &a, const real &b) throw()
173 { return complex( mulu(a.re,b), mulu(a.im,b) ); }
174 
175 inline complex muld(const real &a, const complex &b) throw()
176 { return complex( muld(a,b.re), muld(a,b.im) ); }
177 
178 inline complex mulu(const real &a, const complex &b) throw()
179 { return complex( mulu(a,b.re), mulu(a,b.im) ); }
180 
181  // -------------- division ---------------------------------
182 
183 inline complex divd(const complex &a, const real &b) throw()
184 { return complex( divd(a.re,b), divd(a.im,b) ); }
185 
186 inline complex divu(const complex &a, const real &b) throw()
187 { return complex( divu(a.re,b), divu(a.im,b) ); }
188 
189 inline complex divd(const real &a, const complex &b) throw()
190 { return divd(_complex(a),b); }
191 
192 inline complex divu(const real &a, const complex &b) throw()
193 { return divu(_complex(a),b); }
194 
195 inline complex operator *(const complex &a,const complex &b) throw()
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 
216 inline complex muld(const complex &a, const complex &b) throw()
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 
233 inline complex mulu(const complex &a, const complex &b) throw()
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 
251 static const int Min_Exp_ = 1074, minexpom = -914,
252  maxexpo1 = 1022, MANT_W = 52;
253 
254 inline 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 
374 inline 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 
482 inline complex _c_division(complex a, complex b, int round) throw(DIV_BY_ZERO)
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 
501 inline 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 
506 inline 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 
511 inline 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 
516 inline complex operator / (const complex &a,const complex &b) throw()
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 
526 inline real abs2(const complex &a) throw()
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 
534 inline real abs (complex z) throw()
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
cxsc::minreal
const real minreal
Smallest positive denormalized representable floating-point number.
Definition: real.cpp:63
cxsc::operator*=
cimatrix & operator*=(cimatrix &m, const cinterval &c)
Implementation of multiplication and allocation operation.
Definition: cimatrix.inl:1605
cxsc::interval
The Scalar Type interval.
Definition: interval.hpp:54
cxsc::_complex
complex _complex(const real &a)
Deprecated typecast, which only exist for the reason of compatibility with older versions of C-XSC.
Definition: complex.hpp:366
cxsc::sqrt
cinterval sqrt(const cinterval &z)
Calculates .
Definition: cimath.cpp:1007
cxsc::abs
ivector abs(const cimatrix_subv &mv)
Returns the absolute value of the matrix.
Definition: cimatrix.inl:737
cxsc::operator*
civector operator*(const cimatrix_subv &rv, const cinterval &s)
Implementation of multiplication operation.
Definition: cimatrix.inl:731
cxsc::dotprecision
The Data Type dotprecision.
Definition: dot.hpp:111
cxsc::operator/=
cimatrix & operator/=(cimatrix &m, const cinterval &c)
Implementation of division and allocation operation.
Definition: cimatrix.inl:1623
cxsc::times2pown
void times2pown(cinterval &x, int n)
Fast multiplication of reference parameter [z] with .
Definition: cimath.cpp:2059
cxsc::sqrtx2y2
interval sqrtx2y2(const interval &x, const interval &y)
Calculates .
Definition: imath.cpp:80
cxsc
The namespace cxsc, providing all functionality of the class library C-XSC.
Definition: cdot.cpp:29
cxsc::MaxReal
const real MaxReal
Greatest representable floating-point number.
Definition: real.cpp:65
cxsc::operator+=
cdotprecision & operator+=(cdotprecision &cd, const l_complex &lc)
Implementation of standard algebraic addition and allocation operation.
Definition: cdot.inl:251
cxsc::sqr
cinterval sqr(const cinterval &z)
Calculates .
Definition: cimath.cpp:3342
cxsc::operator/
civector operator/(const cimatrix_subv &rv, const cinterval &s)
Implementation of division operation.
Definition: cimatrix.inl:730
cxsc::complex
The Scalar Type complex.
Definition: complex.hpp:49
cxsc::complex::operator=
complex & operator=(const real &r)
Implementation of standard assigning operator.
Definition: complex.inl:31
cxsc::real
The Scalar Type real.
Definition: real.hpp:113