C-XSC - A C++ Class Library for Extended Scientific Computing  2.5.4
l_complex.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 /* CVS $Id: l_complex.cpp,v 1.21 2014/01/30 17:23:46 cxsc Exp $ */
24 
25 #include "l_complex.hpp" // corresponding header file
26 #include "l_interval.hpp"
27 
28 namespace cxsc {
29 
30 // ---------------- Unary Operators -----------------------------------------
31 l_complex operator-(const l_complex& x)
32 {
33  return l_complex(-x.re, -x.im);
34 }
35 
36 l_complex operator+(const l_complex& x)
37 {
38  return x;
39 }
40 
41 // ----------------- cdotprecision +/- l_complex ----------------------------
42 
43 cdotprecision operator+(const l_complex& lc, const cdotprecision& cd) noexcept
44 {
45  return _cdotprecision(lc) + cd;
46 }
47 
48 cdotprecision operator+(const cdotprecision& cd, const l_complex& lc) noexcept
49 {
50  return _cdotprecision(lc) + cd;
51 }
52 
53 cdotprecision operator-(const l_complex& lc, const cdotprecision& cd) noexcept
54 {
55  return _cdotprecision(lc) - cd;
56 }
57 
58 cdotprecision operator-(const cdotprecision& cd, const l_complex& lc) noexcept
59 {
60  return cd - _cdotprecision(lc);
61 }
62 
63 // ------------------ dotprecision +/- l_complex ----------------------------
64 
65 cdotprecision operator+(const l_complex& lc, const dotprecision& cd) noexcept
66 {
67  return _cdotprecision(lc) + cd;
68 }
69 
70 cdotprecision operator+(const dotprecision& cd, const l_complex& lc) noexcept
71 {
72  return _cdotprecision(lc) + cd;
73 }
74 
75 cdotprecision operator-(const l_complex& lc, const dotprecision& cd) noexcept
76 {
77  return _cdotprecision(lc) - cd;
78 }
79 
80 cdotprecision operator-(const dotprecision& cd, const l_complex& lc) noexcept
81 {
82  return cd - _cdotprecision(lc);
83 }
84 
85 // ---------------- l_complex * l_complex -----------------------------------
87  noexcept
88 {
89  l_real r, i;
90  dotprecision dot(0.0);
91 
92  accumulate(dot,a.re,b.re);
93  accumulate(dot,-a.im,b.im);
94  r = dot;
95 
96  dot = 0.0;
97  accumulate(dot,a.im,b.re);
98  accumulate(dot,a.re,b.im);
99  i = dot;
100 
101  return( l_complex(r,i) );
102 }
103 
104 // ------------------ l_complex * complex -----------------------------------
106  noexcept
107 {
108  l_real r, i;
109  dotprecision dot(0.0);
110 
111  accumulate(dot,a.re,Re(b));
112  accumulate(dot,-a.im,Im(b));
113  r = dot; // r is real-part
114 
115  dot = 0.0;
116  accumulate(dot,a.im,Re(b));
117  accumulate(dot,a.re,Im(b));
118  i = dot; // i is imaginary part
119 
120  return( l_complex(r,i) );
121 }
122 
123 l_complex operator * ( const complex& b, const l_complex& a )
124  noexcept
125 {
126  l_real r, i;
127  dotprecision dot(0.0);
128 
129  accumulate(dot,a.re,Re(b));
130  accumulate(dot,-a.im,Im(b));
131  r = dot; // r is real part
132 
133  dot = 0.0;
134  accumulate(dot,a.im,Re(b));
135  accumulate(dot,a.re,Im(b));
136  i = dot; // i is imaginary part
137 
138  return( l_complex(r,i) );
139 }
140 
141 // ---------------- l_complex / l_complex -----------------------------------
142 
143 static const int maxexpo = 1020;
144 
145 void skale_down_exp(int ex1, int ex2, int D, int& d1, int& d2)
146 // l_interval x1,x2 sind Punktintervalle mit den Zweierexponenten
147 // ex1,ex2 > -maxint, d.h. x1,x2 <> 0.0;
148 // x1*x2 ist mit 2^D, D<0, so zu skalieren, dass fuer das Produkt
149 // x1*x2 moeglichst wenig Stellen verloren gehen! x1 ist dazu mit
150 // 2^d1 und x2 ist mit 2^d2 zu multiplizieren. Es muss also gelten:
151 // d1+d2 = D<0;
152 // Diese Funktion berechnet zu den gegebenen Zweierexponenten
153 // ex1 und ex2 die Exponenten d1 und d2.
154 // ex1 berechnet sich z.B. durch: ex1 = expo(x1).
155 // Blomquist, 25.10.2006;
156 
157 {
158  bool change(false);
159  int Diff, D1, c;
160  d1 = 0; d2 = 0;
161 
162  if (D<0)
163  {
164  if (ex2>ex1)
165  { c = ex1; ex1 = ex2; ex2 = c; change = true; }
166  // ab jetzt gilt: ex1 >= ex2;
167  c = ex1 + D;
168  if (c<ex2)
169  {
170  Diff = ex2 - c; D1 = Diff/2;
171  d1 = D + D1; d2 = D - d1; // d1+d2 == D<0;
172  }
173  else d1 = D; // d2 = 0, s.o.
174  if (change)
175  { c = d1; d1 = d2; d2 = c; }
176  }
177 
178 } // skale_down_exp(...)
179 
180 void skale_up_exp1(int ex1, int ex2, int& fillin, int& d1, int& d2)
181 // l_interval x1,x2 sind Punktintervalle mit den Zweierexponenten
182 // ex1,ex2 > -maxint, d.h. x1,x2 <> 0.0;
183 // Das Maximum (ex1+ex2) der beiden gegebenen Exponentensummen
184 // ist <= 1022; Die Punktintervalle x1,x2 sind so mit 2^d1 bzw. mit
185 // 2^d2 zu multiplizieren, dass gilt.
186 // I. 2*|(x1*x2)| < MaxReal;
187 // II. Bei dem Produkt x1*x2 duerfen moeglichst wenig Stellen
188 // verloren gehen.
189 // Blomquist, 25.10.2006;
190 {
191  bool change(false);
192  int c, pot2;
193  d1 = 0; d2 = 0;
194  fillin = maxexpo - (ex1+ex2);
195  if (fillin>0)
196  {
197  if (ex2>ex1)
198  { c = ex1; ex1 = ex2; ex2 = c; change = true; }
199  // ab jetzt gilt: ex1 >= ex2;
200  pot2 = maxexpo - ex2;
201  // Um pot2 kann x2 ohne Overflow hochskaliert werden.
202  if (fillin <= pot2) d2 = fillin;
203  else { d2 = pot2; d1 = fillin - pot2; }
204  if (change)
205  { c = d1; d1 = d2; d2 = c; }
206  }
207 } // skale_up_exp1(...)
208 
209 void skale_up_exp2(int ex1, int ex2, int fillin, int& d1, int& d2)
210 // l_interval x1,x2 sind Punktintervalle mit den Zweierexponenten
211 // ex1,ex2 > -maxint, d.h. x1,x2 <> 0.0;
212 // Das Minimum (ex1+ex2) der beiden gegebenen Exponentensummen
213 // ist <= 1022; Die Punktintervalle x1,x2 sind so mit 2^d1 bzw. mit
214 // 2^d2 zu multiplizieren, dass gilt.
215 // I. 2*|(x1*x2)| < MaxReal;
216 // II. Bei dem Produkt x1*x2 duerfen moeglichst wenig Stellen
217 // verloren gehen.
218 // Blomquist, 25.10.2006;
219 {
220  bool change(false);
221  int c, pot2;
222  d1 = 0; d2 = 0;
223  if (fillin>0)
224  {
225  if (ex2>ex1)
226  { c = ex1; ex1 = ex2; ex2 = c; change = true; }
227  // ab jetzt gilt: ex1 >= ex2;
228  pot2 = 1022 - ex2; // 1022 ??
229  // Um pot2 kann x2 ohne Overflow hochskaliert werden.
230  if (fillin <= pot2) d2 = fillin;
231  else { d2 = pot2; d1 = fillin - pot2; }
232  if (change)
233  { c = d1; d1 = d2; d2 = c; }
234  }
235 } // skale_up_exp2(...)
236 
237 void product(const l_real& a, const l_real& b, const l_real& c,
238  const l_real& d, int& ex, l_interval& res)
239 // Calulation of an inclusion of: a*b + c*d
240 {
241  l_real a1(a), b1(b), c1(c), d1(d);
242 
243  a1 += 0.0; b1 += 0.0; c1 += 0.0; d1 += 0.0;
244  int ex_a1( expo(a1[1]) ), ex_b1( expo(b1[1]) ),
245  ex_c1( expo(c1[1]) ), ex_d1( expo(d1[1]) ), m,p,D1,D2;
246  l_interval a_i(a1),b_i(b1),c_i(c1),d_i(d1),li; // point intervals!
247  idotprecision Akku(0.0);
248  ex = expo(0.0); res = 0.0; // Initialization for a*b + c*d == 0;
249 
250  if ( ex_a1 == ex || ex_b1 == ex ) // a*b == 0;
251  if ( ex_c1 == ex || ex_d1 == ex ) ; // a*b == c*d == 0.0;
252  else
253  {// a*b == 0; c*d != 0;
254  m = ex_c1 + ex_d1;
255  if (m > maxexpo)
256  {
257  p = maxexpo - m;
258  skale_down_exp(ex_c1, ex_d1, p, D1, D2);
259  Times2pown(c_i,D1);
260  Times2pown(d_i,D2);
261  Akku = 0.0;
262  accumulate(Akku,c_i,d_i);
263  res = Akku;
264  ex = -p;
265  }
266  else
267  {
268  skale_up_exp1(ex_c1, ex_d1, p, D1, D2);
269  Times2pown(c_i,D1);
270  Times2pown(d_i,D2);
271  Akku = 0.0;
272  accumulate(Akku,c_i,d_i);
273  res = Akku;
274  ex = -p;
275  }
276  }
277  else // a*b != 0.0;
278  if ( ex_c1 == ex || ex_d1 == ex )
279  {// a*b != 0.0; c*d == 0.0;
280  m = ex_a1 + ex_b1;
281  if (m > maxexpo)
282  {
283  p = maxexpo - m;
284  skale_down_exp(ex_a1, ex_b1, p, D1, D2);
285  Times2pown(a_i,D1);
286  Times2pown(b_i,D2);
287  Akku = 0.0;
288  accumulate(Akku,a_i,b_i);
289  res = Akku;
290  ex = -p;
291  }
292  else
293  {
294  skale_up_exp1(ex_a1, ex_b1, p, D1, D2);
295  Times2pown(a_i,D1);
296  Times2pown(b_i,D2);
297  Akku = 0.0;
298  accumulate(Akku,a_i,b_i);
299  res = Akku;
300  ex = -p;
301  }
302  }
303  else // a*b != 0.0 and c*d != 0.0;
304  {
305  if ( (ex_c1+ex_d1) > (ex_a1+ex_b1) )
306  {
307  li = a_i; a_i = c_i; c_i = li; // a_i <--> c_i
308  li = b_i; b_i = d_i; d_i = li; // b_i <--> d_i
309  m = ex_a1; ex_a1 = ex_c1; ex_c1 = m;
310  m = ex_b1; ex_b1 = ex_d1; ex_d1 = m;
311  }
312  m = ex_a1 + ex_b1;
313  // ab jetzt gilt: m = (ex_a1+ex_b1) >= (ex_c1+ex_d1);
314  if (m > maxexpo)
315  {
316  p = maxexpo - m;
317  skale_down_exp(ex_a1,ex_b1,p,D1,D2);
318  Times2pown(a_i,D1);
319  Times2pown(b_i,D2);
320  skale_down_exp(ex_c1,ex_d1,p,D1,D2);
321  Times2pown(c_i,D1);
322  Times2pown(d_i,D2);
323  Akku = 0.0;
324  accumulate(Akku,a_i,b_i);
325  accumulate(Akku,c_i,d_i);
326  res = Akku;
327  ex = -p;
328  }
329  else
330  {
331  skale_up_exp1(ex_a1, ex_b1, p, D1, D2);
332  Times2pown(a_i,D1);
333  Times2pown(b_i,D2);
334  skale_up_exp2(ex_c1, ex_d1, p, D1, D2);
335  Times2pown(c_i,D1);
336  Times2pown(d_i,D2);
337  Akku = 0.0;
338  accumulate(Akku,a_i,b_i);
339  accumulate(Akku,c_i,d_i);
340  res = Akku;
341  ex = -p;
342  }
343  }
344 } // product(...)
345 
346 void product(const l_real& c, const l_real& d, int& ex, l_interval& res)
347 // Calulation of an inclusion of: c*c + d*d
348 {
349  l_real c1(c), d1(d);
350 
351  c1 += 0.0; d1 += 0.0;
352  int ex_c1( expo(c1[1]) ), ex_d1( expo(d1[1]) ), m;
353 
354  l_interval c_i(c1),d_i(d1),li; // point intervals!
355  idotprecision Akku(0.0);
356  ex = expo(0.0); res = 0.0; // Initialization for c*c + d*d == 0;
357 
358  if (ex_c1 == ex) // c*c == 0;
359  if (ex_d1 == ex) ; // c*c == d*d == 0;
360  else // c*c == 0; d*d != 0;
361  {
362  times2pown(d_i,-ex_d1); // d_i about 1.0;
363  Akku = 0.0;
364  accumulate(Akku,d_i,d_i);
365  res = Akku;
366  ex = 2*ex_d1;
367  }
368  else // c*c != 0;
369  if (ex_d1 == ex)
370  { // c*c != 0; d*d == 0;
371  times2pown(c_i,-ex_c1); // c_i about 1.0;
372  Akku = 0.0;
373  accumulate(Akku,c_i,c_i);
374  res = Akku;
375  ex = 2*ex_c1;
376  }
377  else // c*c != 0; d*d != 0;
378  {
379  if (ex_d1 > ex_c1)
380  {
381  li = c_i; c_i = d_i; d_i = li; // c_i <--> d_i
382  m = ex_c1; ex_c1 = ex_d1; ex_d1 = m;
383  } // c*c >= d*d:
384  times2pown(c_i,-ex_c1); // c_i about 1.0;
385  times2pown(d_i,-ex_c1);
386  Akku = 0.0;
387  accumulate(Akku,c_i,c_i);
388  accumulate(Akku,d_i,d_i);
389  res = Akku;
390  ex = 2*ex_c1;
391  }
392 
393 } // product(...)
394 
395 l_real quotient(const l_interval& z, const l_interval& n, int round,
396  int ex_z, int ex_n)
397 // z is an inclusion of a numerator.
398 // n is an inclusion of a denominator.
399 // quotient(...) calculates with q1 an approximation of z/n
400 // using staggered arithmetic.
401 // Rounding with round (-1,0,+1) is considered.
402 
403 {
404  l_real q1; // return value
405  int ex_diff;
406  l_interval res;
407 
408  if (0.0<=n) // 0 in denominator n leads to error message:
409  {
410  std::cerr << "quotient1(const l_interval& z, const l_interval& n, int round, int ex_z, int ex_n): Division by zero" << std::endl;
411  exit(1);
412  }
413 
414  if ( zero_(z) ) { q1=0; return q1; };
415 
416  ex_diff = ex_z - ex_n;
417  res = z/n;
418  Times2pown(res,ex_diff);
419  switch(round)
420  {
421  case RND_DOWN:
422  q1 = Inf(res);
423  break;
424  case RND_NEXT:
425  q1 = mid(res);
426  break;
427  case RND_UP:
428  q1 = Sup(res);
429  break;
430  } // switch
431  return q1;
432 } // quotient
433 
434 l_complex _c_division(l_complex a, l_complex b, int round)
435 {
436  int ex1, ex2;
437  l_interval z, n;
438  l_complex tmp;
439 
440  product(Re(b),Im(b),ex2,n);
441  product(Re(a),Re(b),Im(a),Im(b),ex1,z);
442  SetRe(tmp, quotient(z,n,round,ex1,ex2));
443  product(Im(a),Re(b),-Re(a),Im(b),ex1,z);
444  SetIm(tmp, quotient(z,n,round,ex1,ex2));
445  return tmp;
446 } // _c_division
447 
448 l_complex divn (const l_complex & a, const l_complex & b)
449 {
450  return _c_division(a,b,RND_NEXT);
451 }
452 
453 l_complex divd (const l_complex & a, const l_complex & b)
454 {
455  return _c_division(a,b,RND_DOWN);
456 }
457 
458 l_complex divu (const l_complex & a, const l_complex & b)
459 {
460  return _c_division(a,b,RND_UP);
461 }
462 
463 l_complex operator / (const l_complex &a, const l_complex &b) noexcept
464 {
465  return divn(a,b);
466 }
467 
468 int StagPrec(const l_complex& lc) noexcept
469 {
470  return StagPrec(lc.re);
471 }
472 
473 void accumulate(cdotprecision& cd, const l_complex& lc1,
474  const l_complex& lc2) noexcept
475 {
476  accumulate(Re(cd),lc1.re,lc2.re);
477  accumulate(Re(cd),-lc1.im,lc2.im);
478  accumulate(Im(cd),lc1.im,lc2.re);
479  accumulate(Im(cd),lc1.re,lc2.im);
480 }
481 
482 void accumulate(cdotprecision& cd, const l_complex& lc,
483  const complex& c) noexcept
484 {
485  accumulate(Re(cd),lc.re,Re(c));
486  accumulate(Re(cd),-lc.im,Im(c));
487  accumulate(Im(cd),lc.im,Re(c));
488  accumulate(Im(cd),lc.re,Im(c));
489 }
490 
491 void accumulate(cdotprecision& cd, const l_complex& lc,
492  const real& r) noexcept
493 {
494  accumulate(Re(cd),lc.re,r);
495  accumulate(Im(cd),lc.im,r);
496 }
497 
498 void accumulate(cdotprecision& cd, const l_complex& lc,
499  const l_real& lr) noexcept
500 {
501  accumulate(Re(cd),lc.re,lr);
502  accumulate(Im(cd),lc.im,lr);
503 }
504 
505 l_real abs2(const l_complex &a) noexcept
506 {
507  dotprecision dot(0.0);
508  accumulate(dot,a.re,a.re);
509  accumulate(dot,a.im,a.im);
510  return l_real(dot);
511 }
512 
513 l_real abs(const l_complex& z) noexcept
514 // Calculation of an approximation of |z|.
515 // In general the maximum precision is stagprec=19, predifined by the used
516 // sqrt-function declared in l_rmath.hpp.
517 // If the difference |exa-exb| of the exponents (base 2) is sufficient high,
518 // precision and accuracy can be choosen greater 19.
519 {
520  return sqrtx2y2(Re(z),Im(z));
521 } // abs(...)
522 
523 l_complex & SetIm(l_complex & a,const l_real & b)
524  { a.im=b; return a; } // See SetRe(...);
525 
526 l_complex & SetRe(l_complex & a,const l_real & b)
527  { a.re=b; return a; } // The real part of a is substituted by b.
528 
529 l_real & Re(l_complex& a) { return a.re; }
530 l_real Re(const l_complex& a) { return a.re; }
531 l_real & Im(l_complex& a) { return a.im; }
532 l_real Im(const l_complex& a) { return a.im; }
533 
535 {
536  real x(Re(a)), y(Im(a));
537  return *this = complex(x,y);
538 }
539 
540 } // end namespace cxsc
541 
542 
543 
544 
545 
546 
The Data Type cdotprecision.
Definition: cdot.hpp:61
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 Multiple-Precision Data Type l_complex.
Definition: l_complex.hpp:46
The Multiple-Precision Data Type l_real.
Definition: l_real.hpp:78
The Scalar Type real.
Definition: real.hpp:114
The namespace cxsc, providing all functionality of the class library C-XSC.
Definition: cdot.cpp:29
civector operator/(const cimatrix_subv &rv, const cinterval &s) noexcept
Implementation of division operation.
Definition: cimatrix.inl:730
cdotprecision _cdotprecision(const l_complex &)
Definition: cdot.inl:158
ivector abs(const cimatrix_subv &mv) noexcept
Returns the absolute value of the matrix.
Definition: cimatrix.inl:737
void times2pown(cinterval &x, int n) noexcept
Fast multiplication of reference parameter [z] with .
Definition: cimath.cpp:2059
cvector mid(const cimatrix_subv &mv) noexcept
Returns the middle of the matrix.
Definition: cimatrix.inl:739
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