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
28namespace cxsc {
29
30// ---------------- Unary Operators -----------------------------------------
31l_complex operator-(const l_complex& x)
32{
33 return l_complex(-x.re, -x.im);
34}
35
36l_complex operator+(const l_complex& x)
37{
38 return x;
39}
40
41// ----------------- cdotprecision +/- l_complex ----------------------------
42
43cdotprecision operator+(const l_complex& lc, const cdotprecision& cd) noexcept
44{
45 return _cdotprecision(lc) + cd;
46}
47
48cdotprecision operator+(const cdotprecision& cd, const l_complex& lc) noexcept
49{
50 return _cdotprecision(lc) + cd;
51}
52
53cdotprecision operator-(const l_complex& lc, const cdotprecision& cd) noexcept
54{
55 return _cdotprecision(lc) - cd;
56}
57
58cdotprecision operator-(const cdotprecision& cd, const l_complex& lc) noexcept
59{
60 return cd - _cdotprecision(lc);
61}
62
63// ------------------ dotprecision +/- l_complex ----------------------------
64
65cdotprecision operator+(const l_complex& lc, const dotprecision& cd) noexcept
66{
67 return _cdotprecision(lc) + cd;
68}
69
70cdotprecision operator+(const dotprecision& cd, const l_complex& lc) noexcept
71{
72 return _cdotprecision(lc) + cd;
73}
74
75cdotprecision operator-(const l_complex& lc, const dotprecision& cd) noexcept
76{
77 return _cdotprecision(lc) - cd;
78}
79
80cdotprecision 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
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
143static const int maxexpo = 1020;
144
145void 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
180void 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
209void 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
237void 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
346void 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
395l_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
434l_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
448l_complex divn (const l_complex & a, const l_complex & b)
449{
450 return _c_division(a,b,RND_NEXT);
451}
452
453l_complex divd (const l_complex & a, const l_complex & b)
454{
455 return _c_division(a,b,RND_DOWN);
456}
457
458l_complex divu (const l_complex & a, const l_complex & b)
459{
460 return _c_division(a,b,RND_UP);
461}
462
463l_complex operator / (const l_complex &a, const l_complex &b) noexcept
464{
465 return divn(a,b);
466}
467
468int StagPrec(const l_complex& lc) noexcept
469{
470 return StagPrec(lc.re);
471}
472
473void 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
482void 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
491void 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
498void 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
505l_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
513l_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
523l_complex & SetIm(l_complex & a,const l_real & b)
524 { a.im=b; return a; } // See SetRe(...);
525
526l_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
529l_real & Re(l_complex& a) { return a.re; }
530l_real Re(const l_complex& a) { return a.re; }
531l_real & Im(l_complex& a) { return a.im; }
532l_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
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
cdotprecision _cdotprecision(const l_complex &)
Definition cdot.inl:158