25#include "l_complex.hpp"
26#include "l_interval.hpp"
92 accumulate(dot,a.re,b.re);
93 accumulate(dot,-a.im,b.im);
97 accumulate(dot,a.im,b.re);
98 accumulate(dot,a.re,b.im);
111 accumulate(dot,a.re,Re(b));
112 accumulate(dot,-a.im,Im(b));
116 accumulate(dot,a.im,Re(b));
117 accumulate(dot,a.re,Im(b));
129 accumulate(dot,a.re,Re(b));
130 accumulate(dot,-a.im,Im(b));
134 accumulate(dot,a.im,Re(b));
135 accumulate(dot,a.re,Im(b));
143static const int maxexpo = 1020;
145void skale_down_exp(
int ex1,
int ex2,
int D,
int& d1,
int& d2)
165 { c = ex1; ex1 = ex2; ex2 = c; change =
true; }
170 Diff = ex2 - c; D1 = Diff/2;
171 d1 = D + D1; d2 = D - d1;
175 { c = d1; d1 = d2; d2 = c; }
180void skale_up_exp1(
int ex1,
int ex2,
int& fillin,
int& d1,
int& d2)
194 fillin = maxexpo - (ex1+ex2);
198 { c = ex1; ex1 = ex2; ex2 = c; change =
true; }
200 pot2 = maxexpo - ex2;
202 if (fillin <= pot2) d2 = fillin;
203 else { d2 = pot2; d1 = fillin - pot2; }
205 { c = d1; d1 = d2; d2 = c; }
209void skale_up_exp2(
int ex1,
int ex2,
int fillin,
int& d1,
int& d2)
226 { c = ex1; ex1 = ex2; ex2 = c; change =
true; }
230 if (fillin <= pot2) d2 = fillin;
231 else { d2 = pot2; d1 = fillin - pot2; }
233 { c = d1; d1 = d2; d2 = c; }
237void product(
const l_real& a,
const l_real& b,
const l_real& c,
238 const l_real& d,
int& ex, l_interval& res)
241 l_real a1(a), b1(b), c1(c), d1(d);
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;
247 idotprecision Akku(0.0);
248 ex = expo(0.0); res = 0.0;
250 if ( ex_a1 == ex || ex_b1 == ex )
251 if ( ex_c1 == ex || ex_d1 == ex ) ;
258 skale_down_exp(ex_c1, ex_d1, p, D1, D2);
262 accumulate(Akku,c_i,d_i);
268 skale_up_exp1(ex_c1, ex_d1, p, D1, D2);
272 accumulate(Akku,c_i,d_i);
278 if ( ex_c1 == ex || ex_d1 == ex )
284 skale_down_exp(ex_a1, ex_b1, p, D1, D2);
288 accumulate(Akku,a_i,b_i);
294 skale_up_exp1(ex_a1, ex_b1, p, D1, D2);
298 accumulate(Akku,a_i,b_i);
305 if ( (ex_c1+ex_d1) > (ex_a1+ex_b1) )
307 li = a_i; a_i = c_i; c_i = li;
308 li = b_i; b_i = d_i; d_i = li;
309 m = ex_a1; ex_a1 = ex_c1; ex_c1 = m;
310 m = ex_b1; ex_b1 = ex_d1; ex_d1 = m;
317 skale_down_exp(ex_a1,ex_b1,p,D1,D2);
320 skale_down_exp(ex_c1,ex_d1,p,D1,D2);
324 accumulate(Akku,a_i,b_i);
325 accumulate(Akku,c_i,d_i);
331 skale_up_exp1(ex_a1, ex_b1, p, D1, D2);
334 skale_up_exp2(ex_c1, ex_d1, p, D1, D2);
338 accumulate(Akku,a_i,b_i);
339 accumulate(Akku,c_i,d_i);
346void product(
const l_real& c,
const l_real& d,
int& ex, l_interval& res)
351 c1 += 0.0; d1 += 0.0;
352 int ex_c1( expo(c1[1]) ), ex_d1( expo(d1[1]) ), m;
354 l_interval c_i(c1),d_i(d1),li;
355 idotprecision Akku(0.0);
356 ex = expo(0.0); res = 0.0;
364 accumulate(Akku,d_i,d_i);
373 accumulate(Akku,c_i,c_i);
381 li = c_i; c_i = d_i; d_i = li;
382 m = ex_c1; ex_c1 = ex_d1; ex_d1 = m;
387 accumulate(Akku,c_i,c_i);
388 accumulate(Akku,d_i,d_i);
395l_real quotient(
const l_interval& z,
const l_interval& n,
int round,
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;
414 if ( zero_(z) ) { q1=0;
return q1; };
416 ex_diff = ex_z - ex_n;
418 Times2pown(res,ex_diff);
434l_complex _c_division(l_complex a, l_complex b,
int round)
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));
450 return _c_division(a,b,RND_NEXT);
455 return _c_division(a,b,RND_DOWN);
460 return _c_division(a,b,RND_UP);
470 return StagPrec(lc.re);
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);
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));
492 const real& r)
noexcept
494 accumulate(Re(cd),lc.re,r);
495 accumulate(Im(cd),lc.im,r);
499 const l_real& lr)
noexcept
501 accumulate(Re(cd),lc.re,lr);
502 accumulate(Im(cd),lc.im,lr);
508 accumulate(dot,a.re,a.re);
509 accumulate(dot,a.im,a.im);
524 { a.im=b;
return a; }
527 { a.re=b;
return a; }
536 real x(Re(a)), y(Im(a));
The Data Type cdotprecision.
complex & operator=(const real &r) noexcept
Implementation of standard assigning operator.
The Data Type dotprecision.
The Multiple-Precision Data Type l_complex.
The Multiple-Precision Data Type l_real.
The namespace cxsc, providing all functionality of the class library C-XSC.
civector operator/(const cimatrix_subv &rv, const cinterval &s) noexcept
Implementation of division operation.
ivector abs(const cimatrix_subv &mv) noexcept
Returns the absolute value of the matrix.
void times2pown(cinterval &x, int n) noexcept
Fast multiplication of reference parameter [z] with .
cvector mid(const cimatrix_subv &mv) noexcept
Returns the middle of the matrix.
civector operator*(const cimatrix_subv &rv, const cinterval &s) noexcept
Implementation of multiplication operation.
interval sqrtx2y2(const interval &x, const interval &y) noexcept
Calculates .
cdotprecision _cdotprecision(const l_complex &)