25#include "cinterval.hpp"
50 return complex(a.re+b.re,a.im+b.im);
55 return complex(a.re-b.re,a.im-b.im);
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; }
131{
return complex(addd(a.re,b.re), addd(a.im,b.im)); }
134{
return complex(addu(a.re,b.re), addu(a.im,b.im)); }
137{
return complex(addd(a.re,b), a.im); }
140{
return complex(addu(a.re,b), a.im); }
143{
return complex(addd(a,b.re), b.im); }
146{
return complex(addu(a,b.re), b.im); }
150{
return complex(subd(a.re,b.re), subd(a.im,b.im)); }
153{
return complex(subu(a.re,b.re), subu(a.im,b.im)); }
156{
return complex(subd(a.re,b), a.im); }
159{
return complex(subu(a.re,b), a.im); }
162{
return complex(subd(a,b.re), -b.im); }
165{
return complex(subu(a,b.re), -b.im); }
170{
return complex( muld(a.re,b), muld(a.im,b) ); }
173{
return complex( mulu(a.re,b), mulu(a.im,b) ); }
176{
return complex( muld(a,b.re), muld(a,b.im) ); }
179{
return complex( mulu(a,b.re), mulu(a,b.im) ); }
184{
return complex( divd(a.re,b), divd(a.im,b) ); }
187{
return complex( divu(a.re,b), divu(a.im,b) ); }
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));
203 accumulate (dot, a.re, b.re);
204 accumulate (dot, -a.im, b.im);
205 rnd (dot, tmp.re, RND_NEXT);
208 accumulate (dot, a.re, b.im);
209 accumulate (dot, a.im, b.re);
210 rnd (dot, tmp.im, RND_NEXT);
221 accumulate (dot, a.re, b.re);
222 accumulate (dot, -a.im, b.im);
223 rnd (dot, tmp.re, RND_DOWN);
226 accumulate (dot, a.re, b.im);
227 accumulate (dot, a.im, b.re);
228 rnd (dot, tmp.im, RND_DOWN);
238 accumulate (dot, a.re, b.re);
239 accumulate (dot, -a.im, b.im);
240 rnd (dot, tmp.re, RND_UP);
243 accumulate (dot, a.re, b.im);
244 accumulate (dot, a.im, b.re);
245 rnd (dot, tmp.im, RND_UP);
251static const int Min_Exp_ = 1074, minexpom = -914,
252 maxexpo1 = 1022, MANT_W = 52;
265 int exa, exb, exc, exd;
278 if ( sign(a) == 0 || sign(b) == 0 )
279 if ( sign(c) == 0 || sign(d) == 0 )
285 if (exc+exd > maxexpo1)
288 if ( exc > exd ) c = comp( mant(c), exc-Min_Exp_ );
289 else d = comp( mant(d), exd-Min_Exp_ );
292 if ( exc+exd < minexpom )
295 if (exc < exd) c = comp( mant(c), exc+Min_Exp_ );
296 else d = comp( mant(d), exd+Min_Exp_ );
301 if ( sign(c) == 0 || sign(d) == 0 )
304 if (exa+exb > maxexpo1)
307 if ( exa > exb ) a = comp( mant(a), exa-Min_Exp_ );
308 else b = comp( mant(b), exb-Min_Exp_ );
311 if (exa+exb < minexpom)
314 if (exa < exb) a = comp( mant(a), exa+Min_Exp_ );
315 else b = comp( mant(b), exb+Min_Exp_ );
322 if (exa+exb > maxexpo1)
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_ );
336 }
else if (exc+exd > maxexpo1)
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_ );
352 if ( exa+exb < minexpom && exc+exd < minexpom )
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_ );
361 accumulate(dot, a, b);
362 accumulate(dot, c, d);
370 p2 = interval( pred(Inf(p2)), succ(Sup(p2)) );
374inline real quotient (real z1, interval z2, real n1,
375 interval n2,
int round,
int zoverfl,
int noverfl)
387 int vorz, anz_scale, ex = 0;
391 if ( zoverfl == -1 && noverfl == 1 )
412 }
else if ( zoverfl==1 && noverfl==-1 )
420 nh = interval(addd(n1, Inf(n2)),
425 accumulate(
id, -q1, n1);
427 accumulate(
id, -q1, n2);
433 q1 = adddown(q1, divd(Inf(q2), Sup(nh)));
436 q1 = q1 + (Inf(q2)+Sup(q2))*0.5/n1;
439 q1 = addup(q1, divu(Sup(q2), Inf(nh)));
451 anz_scale = zoverfl - noverfl;
454 scale = comp(0.5, +1024);
457 else if (anz_scale < 0)
459 scale = comp(0.5, -1022);
469 q1 = multdown(q1, scale);
475 q1 = multup(q1, scale);
482inline complex _c_division(complex a, complex b,
int round)
484 if (0.0 == (
sqr(Re(b))+
sqr(Im(b)))) {
485 cxscthrow(DIV_BY_ZERO(
"complex operator / (const complex&,const complex&)"));
488 int zoverflow, noverflow;
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));
503 return _c_division(a, b, RND_NEXT);
508 return _c_division(a, b, RND_DOWN);
513 return _c_division(a, b, RND_UP);
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);
529 accumulate(dot,a.re,a.re);
530 accumulate(dot,a.im,a.im);
536#ifdef CXSC_FAST_COMPLEX_OPERATIONS
537 return sqrt(Re(z)*Re(z)+Im(z)*Im(z));
complex & operator=(const real &r) noexcept
Implementation of standard assigning operator.
The Data Type dotprecision.
The Scalar Type interval.
The namespace cxsc, providing all functionality of the class library C-XSC.
complex _complex(const real &a) noexcept
Deprecated typecast, which only exist for the reason of compatibility with older versions of C-XSC.
cdotprecision & operator+=(cdotprecision &cd, const l_complex &lc) noexcept
Implementation of standard algebraic addition and allocation operation.
const real MaxReal
Greatest representable floating-point number.
civector operator/(const cimatrix_subv &rv, const cinterval &s) noexcept
Implementation of division operation.
const real minreal
Smallest positive denormalized representable floating-point number.
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.
cinterval sqrt(const cinterval &z) noexcept
Calculates .
void times2pown(cinterval &x, int n) noexcept
Fast multiplication of reference parameter [z] with .
cinterval sqr(const cinterval &z) noexcept
Calculates .
civector operator*(const cimatrix_subv &rv, const cinterval &s) noexcept
Implementation of multiplication operation.
interval sqrtx2y2(const interval &x, const interval &y) noexcept
Calculates .
cimatrix & operator/=(cimatrix &m, const cinterval &c) noexcept
Implementation of division and allocation operation.