C-XSC - A C++ Class Library for Extended Scientific Computing  2.5.4
lx_cinterval.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 
24 /* CVS $Id: lx_cinterval.cpp,v 1.12 2014/01/30 17:23:47 cxsc Exp $ */
25 
26 #include "lx_cinterval.hpp"
27 
28 namespace cxsc {
29 
30 // -----------------------------------------------------------------------
31 // --------- Functions and operators related to type lx_cinterval --------
32 // -----------------------------------------------------------------------
33 
34 l_cinterval & l_cinterval::operator = (const lx_cinterval &a) noexcept
35 // This operator is declared in l_cinterval.hpp !!
36 {
37  l_interval lre,lim;
38  lre = Re(a);
39  lim = Im(a);
40  l_cinterval lc(lre,lim);
41 
42  return (*this) = lc;
43 }
44 
45 cinterval & cinterval::operator = (const lx_cinterval &a) noexcept
46 // This operator is declared in cinterval.hpp !!
47 {
48  l_cinterval lc;
49  cinterval z;
50 
51  lc = a;
52  z = lc;
53 
54  return (*this) = z;
55 }
56 
57 // ------------------------ Input --------------------------------------
58 /*
59 std::string & operator >> (std::string &s, lx_cinterval &a) noexcept
60 // Funktionierende Version mit Zweier-Exponenten
61 {
62  lx_interval Lar, Lai;
63 
64  s = skipwhitespacessinglechar (s, '(');
65  s = s >> SaveOpt >> Lar;
66  s = skipwhitespacessinglechar (s, ',');
67  s = s >> Lai >> RestoreOpt;
68  s = skipwhitespacessinglechar (s, ')');
69  s = skipwhitespaces (s); // erasing whitespaces
70  a = lx_cinterval(Lar,Lai);
71 
72  return s;
73 }
74 */
75 std::string & operator >> (std::string &s, lx_cinterval &a) noexcept
76 // Funktionierende Version mit Zehner-Exponenten, z.B.
77 // s = ( {+1, [4,4.01]} , {-2, [7,7]} )
78 {
79  lx_interval Lar, Lai;
80  string str(s);
81  int i;
82 
83  str = skipwhitespacessinglechar (str, '(');
84  i = str.find("}");
85  str.erase(i+1);
86  str >> SaveOpt >> Lar;
87 
88  i = s.find("}");
89  s.erase(0,i+1);
90  s = skipwhitespacessinglechar (s, ',');
91  s >> Lai >> RestoreOpt;
92  s = "";
93 
94  a = lx_cinterval(Lar,Lai);
95 
96  return s;
97 }
98 
99 void operator >> (const std::string &s, lx_cinterval &a) noexcept
100 {
101 // Writes strings s to variable a of type lx_cinterval;
102  std::string r(s);
103  r >> a;
104 }
105 
106 void operator >> (const char *s, lx_cinterval &a) noexcept
107 {
108  std::string r(s);
109  r >> a;
110 }
111 
112 std::istream & operator >> (std::istream &s, lx_cinterval &a) noexcept
113 {
114  lx_interval Lar, Lai;
115  char c;
116 
117  std::cerr << "Real part: {Exponent to base 10, [a,b]} = ?"
118  << std::endl;
119  s >> Lar;
120  std::cerr << "Img. part: {Exponent to base 10, [a,b]} = ?"
121  << std::endl;
122  s >> Lai >> RestoreOpt;
123  a = lx_cinterval(Lar,Lai);
124 
125  if (!waseolnflag)
126  {
127  skipeolnflag = false; inpdotflag = true;
128  c = skipwhitespaces (s);
129  if (inpdotflag && c != ')')
130  s.putback(c);
131  }
132  return s;
133 }
134 
135 
136 // -----------------------------------------------------------------------------
137 // ---------- Elementary functions related to type lx_cinterval ----------------
138 // -----------------------------------------------------------------------------
139 
140 // --------------------------------------------------------------------------
141 // ------------------------- Help functions ---------------------------------
142 // --------------------------------------------------------------------------
143 
144  int max(int a, int b)
145  {
146  int res(a);
147  if (b>a) res = b;
148  return res;
149  }
150 
151 
152 // --------------------------------------------------------------------------
153 // --------------------------- sqr function ---------------------------------
154 // --------------------------------------------------------------------------
155 
156 lx_cinterval sqr(const lx_cinterval &z) noexcept
157 {
158  lx_interval rez(Re(z)), reza(abs(rez)),
159  imz(Im(z)), imza(abs(imz));
160  lx_real
161  irez = Inf(reza),
162  srez = Sup(reza),
163  iimz = Inf(imza),
164  simz = Sup(imza);
165  lx_interval
166  hxl(irez), hxu(srez), hyl(iimz), hyu(simz);
167  lx_real resxl, resxu;
168  resxl = Inf( (hxl-hyu)*(hxl+hyu) );
169  resxu = Sup( (hxu-hyl)*(hxu+hyl) );
170 
171  hxl = rez * imz;
172  times2pown(hxl,1);
173 
174  return lx_cinterval( lx_interval(resxl,resxu), hxl );
175 } // sqr(...)
176 
177 // -------------------------------------------------------------------------
178 // ---------------------- Begin square root --------------------------------
179 // -------------------------------------------------------------------------
180 
181  lx_interval Sqrt_zpx_m2( const lx_interval &x, const lx_interval &y )
182  // z = x + i*y;
183  // Calculating res = sqrt( 2*(|z| + |x|) );
184  {
185  lx_interval res;
186  res = sqrtx2y2(x,y) + abs(x); // New: 22.05.2008;
187  times2pown(res,1);
188  return sqrt(res);
189  }
190 
191  lx_interval Sqrt_zpx_d2( const lx_interval &x, const lx_interval &y )
192  // z = x + i*y;
193  // Calculating res = sqrt( (|z| + |x|)/2 );
194  {
195  lx_interval res;
196  res = sqrtx2y2(x,y) + abs(x); // New: 22.05.2008;
197  times2pown(res,-1);
198  return sqrt(res);
199  }
200 
201 lx_interval Re_Sqrt_Point(const lx_interval &rez, const lx_interval &imz)
202  // Real part of analytic square root of a POINT INTERVAL only.
203  // Do not use this as a general function
204  // - it is only a subroutine for sqrt(...)
205  // Blomquist, 03.05.2007;
206 {
207  lx_interval res;
208  lx_real irez = Inf( rez ), iimz = Inf( imz );
209 
210  if (eq_zero(iimz))
211  res = (ge_zero(irez))? sqrt(rez) : lx_interval(0);
212  else // iimz <> 0
213  res = (ge_zero(irez))? Sqrt_zpx_d2(rez,imz) :
214  abs(iimz)/Sqrt_zpx_m2(rez,imz);
215  return res;
216 }
217 
218 lx_interval Im_Sqrt_Point(const lx_interval &rez, const lx_interval &imz)
219  // Imaginary part of analytic square root of a POINT INTERVAL only.
220  // Do not use this as a general function
221  // - it is only a subroutine for sqrt(...)
222  // Blomquist, 03.05.2007;
223 {
224  lx_interval res;
225  lx_real irez = Inf( rez ), iimz = Inf( imz );
226 
227  if (eq_zero(iimz))
228  res = (ge_zero(irez))? lx_interval(0) : sqrt(-rez);
229  else
230  if (ge_zero(irez)) res = iimz/Sqrt_zpx_m2(rez,imz);
231  else
232  {
233  res = Sqrt_zpx_d2(rez,imz);
234  if (se_zero(iimz)) res = -res;
235  }
236  return res;
237 }
238 
239 lx_cinterval sqrt(const lx_cinterval& z) noexcept
240 {
241  lx_cinterval y;
242  lx_real
243  irez = Inf( Re(z) ),
244  srez = Sup( Re(z) ),
245  iimz = Inf( Im(z) ),
246  simz = Sup( Im(z) );
247  lx_interval
248  hxl( irez ), hxu( srez ), hyl( iimz ), hyu( simz );
249  lx_real
250  resxl, resxu, resyl, resyu;
251 
252  if ( sm_zero(irez) && sm_zero(iimz) && ge_zero(simz) )
253  cxscthrow(STD_FKT_OUT_OF_DEF(
254  "lx_cinterval sqrt(const lx_cinterval& z); z not in principal branch."));
255 
256  if (ge_zero(iimz))
257  {
258  resxl = Inf( Re_Sqrt_Point( hxl,hyl ) );
259  resxu = Sup( Re_Sqrt_Point( hxu,hyu ) );
260 
261  resyl = Inf( Im_Sqrt_Point( hxu,hyl ) );
262  resyu = Sup( Im_Sqrt_Point( hxl,hyu ) );
263  }
264  else
265  if (se_zero(simz))
266  {
267  resxl = Inf( Re_Sqrt_Point( hxl,hyu ) );
268  resxu = Sup( Re_Sqrt_Point( hxu,hyl ) );
269 
270  resyl = Inf( Im_Sqrt_Point( hxl,hyl ) );
271  resyu = Sup( Im_Sqrt_Point( hxu,hyu ) );
272  }
273  else
274  {
275  resxl = Inf( sqrt(hxl) );
276  resxu = ( -iimz>simz ? Sup( Re_Sqrt_Point( hxu,hyl ) )
277  : Sup( Re_Sqrt_Point( hxu,hyu ) ) );
278 
279  resyl = Inf( Im_Sqrt_Point( hxl,hyl ) );
280  resyu = Sup( Im_Sqrt_Point( hxl,hyu ) );
281  }
282  y = lx_cinterval( lx_interval(resxl,resxu), lx_interval(resyl,resyu) );
283 
284  return y;
285 } // sqrt(...)
286 
287 // sqrt_all(z) computes a list of 2 intervals containing all square roots of z
288 //
289 std::list<lx_cinterval> sqrt_all( const lx_cinterval& z ) noexcept
290 {
291  lx_real
292  irez = Inf( Re(z) ),
293  srez = Sup( Re(z) ),
294  iimz = Inf( Im(z) ),
295  simz = Sup( Im(z) );
296  lx_interval hxl( irez ), hxu( srez ), hyl( iimz ), hyu( simz );
297  lx_real resxl, resxu, resyl, resyu;
298  lx_cinterval w;
299 
300  if( irez < 0.0 && iimz <= 0.0 && simz >= 0.0 )
301  // z contains negative real values
302  {
303  if( iimz == 0.0 )
304  // z in upper half plane
305  // principal values can be used
306  {
307  // min( Re ( sqrt( z ) ) ) in lower left corner
308  // max( Re ( sqrt( z ) ) ) in upper right corner
309  resxl = Inf( Re_Sqrt_Point( hxl, hyl ) );
310  resxu = Sup( Re_Sqrt_Point( hxu, hyu ) );
311  // min( Im ( sqrt( z ) ) ) in lower right corner
312  // max( Im ( sqrt( z ) ) ) in upper left corner
313  resyl = Inf( Im_Sqrt_Point( hxu, hyl ) );
314  resyu = Sup( Im_Sqrt_Point( hxl, hyu ) );
315  }
316  else
317  {
318  if( simz == 0.0 )
319  // z in lower half plane
320  // principal values can be used in lower corners
321  {
322  // z in lower half plane
323  // min( Re ( sqrt( z ) ) ) in upper left corner
324  // max( Re ( sqrt( z ) ) ) in lower right corner
325  resxl = 0;
326  resxu = Sup( Re_Sqrt_Point( hxu, hyl ) );
327  // min( Im ( sqrt( z ) ) ) in lower left corner
328  // max( Im ( sqrt( z ) ) ) in upper right corner
329  resyl = Inf( Im_Sqrt_Point( hxl, hyl ) );
330  resyu = ( srez > 0.0 ? lx_real(0.0,l_real(0)) : -Inf(sqrt(-hxu) ) );
331  }
332  else
333  // 0 is interior point of Im( z )
334  {
335  if( srez <= 0.0 )
336  {
337  // 0 is no interior point of Re (z )
338  // in quadrant III, Re( z ) = Im_Sqrt_Point(|x|,y),
339  // Im( z ) = Re_Sqrt_Point(|x|,y)
340  // min( Re ( sqrt( z ) ) ) in lower right corner = quadrant II/III
341  // max( Re ( sqrt( z ) ) ) in upper right corner = quadrant II
342  resxl = Inf( Im_Sqrt_Point(-hxu, hyl ) );
343  resxu = Sup( Re_Sqrt_Point( hxu, hyu ) );
344  // min( Im ( sqrt( z ) ) ) on real line
345  // max( Im ( sqrt( z ) ) ) in lower or upper left corner
346  resyl = Inf( sqrt( -hxu ) );
347  resyu = ( - iimz > simz ? Sup( Re_Sqrt_Point( -hxl, hyl ) ) :
348  Sup( Im_Sqrt_Point( hxl, hyu ) ) );
349  }
350  else
351  // 0 is interior point of z
352  // here, the principal values apply in all corner points
353  {
354  resxl = 0;
355  // max( Re ( sqrt( z ) ) ) in lower or upper right corner
356  resxu = ( - iimz > simz ? Sup( Re_Sqrt_Point( hxu, hyl ) ) :
357  Sup( Re_Sqrt_Point( hxu, hyu ) ) );
358  // min( Im ( sqrt( z ) ) ) in lower left corner
359  // max( Im ( sqrt( z ) ) ) in upper left corner
360  resyl = Inf( Im_Sqrt_Point( hxl, hyl ) );
361  resyu = Sup( Im_Sqrt_Point( hxl, hyu ) );
362  }
363  }
364  }
365  w = lx_cinterval( lx_interval(resxl,resxu), lx_interval(resyl,resyu ) );
366  }
367  else
368  // sqrt( z ) is well-defined
369  w = sqrt( z );
370 
371  std::list<lx_cinterval> res;
372  res.push_back( w );
373  res.push_back( -w );
374 
375  return res;
376 }
377 //-- end sqrt_all -------------------------------------------------------------
378 // ---------------------- End square root -----------------------------------
379 
380 // ***************************************************************************
381 // ***************************************************************************
382 // *** Single-valued functions ****
383 // ***************************************************************************
384 // ***************************************************************************
385 
386 
387 // ***************************************************************************
388 // *** Power operator pow is not listed here, since it relies on the ****
389 // *** (multi-valued) logarithm ****
390 // ***************************************************************************
391 
392 
393 // ***************************************************************************
394 // *** The hyperbolic functions exp, sin, cos, sinh, cosh are separable: ****
395 // *** Their real and imaginary parts are products of real functions ****
396 // ***************************************************************************
397 // *** With Re(z)=x, Im(z)=y : ****
398 // *** ****
399 // *** exp : Re(exp(z)) = exp(x) * cos(y) ****
400 // *** Im(exp(z)) = exp(x) * sin(y) ****
401 // *** ****
402 // *** sin : Re(sin(z)) = sin(x) * cosh(y) ****
403 // *** Im(sin(x)) = cos(x) * sinh(y) ****
404 // *** ****
405 // *** cos : Re(cos(z)) = cos(x) * cosh(y) ****
406 // *** Im(sin(x)) = -sin(x) * sinh(y) ****
407 // *** ****
408 // *** sinh : Re(sinh(z)) = sinh(x) * cos(y) ****
409 // *** Im(sinh(z)) = cosh(x) * sin(y) ****
410 // *** ****
411 // *** cosh : Re(cosh(z)) = cosh(x) * cos(y) ****
412 // *** Im(cosh(z)) = sinh(x) * sin(y) ****
413 // *** ****
414 // ***************************************************************************
415 
416 // -------------------------- exp(...) --------------------------------------
417 
418  lx_cinterval exp(const lx_cinterval& z) noexcept
419  {
420  int stagsave = stagprec,
421  stagmax = 39;
422  if (stagprec > stagmax)
423  stagprec = stagmax;
424 
425  lx_interval lreal(Re(z)), limg(Im(z));
426  lx_cinterval y;
427  lx_interval
428  A( exp(lreal) ),
429  B( limg );
430  y = lx_cinterval( A*cos( B ) , A*sin( B ) );
431 
432  stagprec = stagsave;
433  y = adjust(y);
434 
435  return y;
436  }
437 
438  lx_cinterval exp2(const lx_cinterval& z) noexcept
439  {
440  return exp(z*Ln2_lx_interval());
441  }
442 
443  lx_cinterval exp10(const lx_cinterval& z) noexcept
444  {
445  return exp(z*Ln10_lx_interval());
446  }
447 
448 // -------------------------- cos(...) --------------------------------------
449 
450  lx_cinterval cos(const lx_cinterval& z) noexcept
451  {
452  int stagsave = stagprec,
453  stagmax = 39;
454  if (stagprec > stagmax)
455  stagprec = stagmax;
456 
457  lx_interval lreal(Re(z)),limg(Im(z));
458  lx_cinterval y;
459  lx_interval
460  A( lreal ),
461  B( limg );
462  y = lx_cinterval( cos( A )*cosh( B ) , -sin( A )*sinh( B ) );
463 
464  stagprec = stagsave;
465  y = adjust(y);
466 
467  return y;
468  }
469 
470 // -------------------------- sin(...) --------------------------------------
471 
472  lx_cinterval sin(const lx_cinterval& z) noexcept
473  {
474  int stagsave = stagprec,
475  stagmax = 39;
476  if (stagprec > stagmax)
477  stagprec = stagmax;
478 
479  lx_interval lreal(Re(z)),limg(Im(z));
480  lx_cinterval y;
481  lx_interval
482  A( lreal ),
483  B( limg );
484  y = lx_cinterval( sin( A )*cosh( B ) , cos( A )*sinh( B ) );
485 
486  stagprec = stagsave;
487  y = adjust(y);
488 
489  return y;
490  }
491 
492 // -------------------------- cosh(...) -------------------------------------
493 
494  lx_cinterval cosh(const lx_cinterval& z) noexcept
495  {
496  int stagsave = stagprec,
497  stagmax = 39;
498  if (stagprec > stagmax)
499  stagprec = stagmax;
500 
501  lx_interval lreal(Re(z)),limg(Im(z));
502  lx_cinterval y;
503  lx_interval
504  A( lreal ),
505  B( limg );
506  y = lx_cinterval( cos( B )*cosh( A ) , sin( B )*sinh( A ) );
507 
508  stagprec = stagsave;
509  y = adjust(y);
510 
511  return y;
512  }
513 
514 // -------------------------- sinh(...) -------------------------------------
515 
516  lx_cinterval sinh(const lx_cinterval& z) noexcept
517  {
518  int stagsave = stagprec,
519  stagmax = 39;
520  if (stagprec > stagmax)
521  stagprec = stagmax;
522 
523  lx_interval lreal(Re(z)),limg(Im(z));
524  lx_cinterval y;
525  lx_interval
526  A( lreal ),
527  B( limg );
528  y = lx_cinterval( cos( B )*sinh( A ) , sin( B )*cosh( A ) );
529 
530  stagprec = stagsave;
531  y = adjust(y);
532 
533  return y;
534  }
535 
536 //-----------------------------------------------------------------------------
537  //
538 // Part I: Multi-valued functions
539  //
540 //-----------------------------------------------------------------------------
541  //
542  //
543 
544 lx_interval Atan(const lx_interval& y,const lx_interval& x)
545 // Calculating an inclusion of atan(y/x) with x<>[0,0];
546 // This help function must only be used for POINT intervals y,x !!
547 // This function avoids internal overflow by calculating y/x.
548 // Blomquist, 25.05.2008;
549 {
550  lx_interval res(0),u;
551  l_interval yl(li_part(y));
552  int ex_yl(expo_gr(yl)),
553  signy(sign(Inf(y))), signx(sign(Inf(x)));
554  real ex_y(expo(Inf(y))), ex_x(expo(Inf(x)));
555  bool neg;
556 
557  if (ex_yl > -1000000)
558  { // y != 0;
559  neg = signy * signx == -1;
560  if (ex_y > ex_x+4197)
561  { // We now assume y>0 and x>0:
562  res = Pi_lx_interval();
563  times2pown(res,-1);
564  if (neg) res = -res;
565  }
566  else
567  {
568  if (ex_x >= 9007199254738946.0)
569  // Now y/x can generate an error!
570  {
571  if (ex_y<-5217)
572  {
573  res = lx_interval( lx_real(0,l_real(0)),
574  lx_real(-Max_Int_R,l_real(minreal)) );
575  if (neg) res = -res;
576  }
577  else // ex_y >= -5217
578  {
579  res = x;
580  times2pown(res,-2045);
581  u = y;
582  times2pown(u,-2045);
583  res = atan(u/res); // Division now without an error!
584  }
585  }
586  else
587  if (ex_x <= -9007199254738891.0)
588  {
589  res = x;
590  times2pown(res,2101);
591  u = y;
592  times2pown(u,2101);
593  res = atan(u/res);
594  }
595  else
596  res = atan(y/x);
597  }
598  }
599 
600  return res;
601 }
602 
603 lx_interval Atan(const lx_interval& y, const lx_real& x)
604 {
605  return Atan(y,lx_interval(x));
606 }
607 
608 //-----------------------------------------------------------------------------
609 //
610 // Section 1: Argument functions
611 //
612 //-----------------------------------------------------------------------------
613 
614 //-- Arg: analytic argument function -------------------------------- 040916 --
615 //
616 // (i) Arg(Z) subset (-pi,pi).
617 // (ii) Arg([0,0]) = 0.
618 // (iii) Arg(Z) is undefined if Z contains negative real numbers.
619 // (iv) Otherwise, Arg(Z) is the interval hull of { Arg(z) | z in Z, z<>0 }.
620 //
621 // atan is the inverse function of tan(t), t in (-pi/2,pi/2).
622 //
623 lx_interval Arg(const lx_cinterval& z) noexcept
624 {
625  lx_real
626  srez = Sup( Re(z) ),
627  irez = Inf( Re(z) ),
628  simz = Sup( Im(z) ),
629  iimz = Inf( Im(z) );
630 
631  lx_interval
632  hxl(irez), hxu(srez), hyl(iimz), hyu(simz), Pid2; // Point intervals
633 
634  lx_real resl, resu;
635 
636  Pid2 = Pid2_lx_interval();
637 
638  if( iimz > 0.0 )
639  // case I: Im(z) > 0
640  {
641  resl = ( srez > 0.0 ? Inf( Atan( hyl,hxu ) ) :
642  ( srez < 0.0 ? Inf( Atan( hyu,hxu ) + Pi_lx_interval() ) :
643  Inf( Pid2 ) ) );
644  resu = ( irez > 0.0 ? Sup( Atan( hyu,hxl ) ) :
645  ( irez < 0.0 ? Sup( Atan( hyl,hxl ) + Pi_lx_interval() ) :
646  Sup( Pid2 ) ) );
647  return lx_interval( resl, resu );
648  }
649  else
650  {
651  if( simz < 0.0 )
652  // case II: Im(z) < 0
653  {
654  resl = ( irez < 0.0 ? Inf( Atan( hyu,hxl ) - Pi_lx_interval() ) :
655  ( irez > 0.0 ? Inf( Atan( hyl,hxl ) ) : -Sup( Pid2 ) ) );
656  resu = ( srez < 0.0 ? Sup( Atan( hyl,hxu ) - Pi_lx_interval() ) :
657  ( srez > 0.0 ? Sup( Atan( hyu,hxu ) ) : -Inf(Pid2) ) );
658  return lx_interval( resl, resu );
659  }
660  else
661  // 0 in Im(z)
662  {
663  if( irez > 0.0 )
664  // case III: Re(z) > 0
665  // z contains positive real values
666  {
667  resl = iimz < 0.0 ? Inf( Atan( hyl,hxl ) ) : lx_real(0.0);
668  return lx_interval( resl, Sup( Atan( hyu,hxl ) ) );
669  }
670  else
671  // z contains nonpositive real numbers
672  {
673  if( irez < 0.0 )
674  {
675  // case IV: z contains negative real numbers
676  cxscthrow (STD_FKT_OUT_OF_DEF("lx_interval Arg(const lx_cinterval& z); z contains negative real numbers"));
677  return lx_interval(0.0);
678  }
679  else
680  // case V: 0 in z, but z doesn't contain negative real numbers
681  {
682  if( srez > 0.0 )
683  // diam( Re(z) > 0.0 )
684  {
685  resl = iimz < 0.0 ? -Sup(Pid2) : lx_real(0.0);
686  resu = simz > 0.0 ? Sup(Pid2) : lx_real(0.0);
687  return lx_interval( resl, resu );
688  }
689  else
690  // Re(z) == 0.0
691  {
692  if( eq_zero(iimz) && eq_zero(simz) )
693  // Z == 0
694  return lx_interval(0.0);
695  else
696  {
697  resl = ( iimz < 0.0 ? - Sup(Pid2) : Inf(Pid2) );
698  resu = ( simz > 0.0 ? Sup(Pid2) : -Inf(Pid2) );
699  return lx_interval( resl, resu );
700  }
701  }
702  }
703  }
704  }
705  }
706 }
707  //
708 //-- end Arg ------------------------------------------------------------------
709 
710 //-- arg: non-analytic argument function --------------------------------------
711  //
712 // (i) arg(Z) is defined for all Z in IC.
713 // (ii) arg([0,0]) = 0.
714 // (iii) arg(Z) subset [-pi,3*pi/2).
715 // (iv) arg(Z) == Arg(Z) if Arg(Z) is well-defined.
716  //
717 
718 lx_interval arg( const lx_cinterval& z ) noexcept
719 {
720  lx_real
721  srez = Sup( Re(z) ),
722  irez = Inf( Re(z) ),
723  simz = Sup( Im(z) ),
724  iimz = Inf( Im(z) );
725 
726  lx_real resl, resu;
727  lx_interval Pid2;
728  Pid2 = Pid2_lx_interval();
729 
730  if( sm_zero(irez) && se_zero(iimz) && ge_zero(simz) )
731  // z contains negative real values
732  {
733  if( gr_zero(srez) )
734  // 0 in z and 0 interior point of Re(z)
735  {
736  resl = ( sm_zero(iimz)? - Sup( Pi_lx_interval() ) : lx_real(0.0) );
737  resu = ( ( sm_zero(iimz) && eq_zero(simz) ) ? lx_real(0.0) :
738  Sup( Pi_lx_interval() ) );
739  return lx_interval( resl, resu );
740  }
741  else
742  { // srez <= 0.0
743  if( iimz == simz )
744  // z is real interval containing no positive values
745  return Pi_lx_interval();
746  else
747  // sup( Re(z) ) <= 0, diam( Im(z) ) > 0
748  {
749  if( eq_zero(srez) )
750  {
751  resl = ( gr_zero(simz)? Inf(Pid2) :
752  -Sup( Pi_lx_interval() ) );
753  resu = ( sm_zero(iimz)?
754  ( gr_zero(simz)? Sup(3*Pid2) : -Inf(Pid2) ) :
755  Sup( Pi_lx_interval() ) );
756  return lx_interval( resl, resu );
757  }
758  else
759  // sup( Re(z) ) < 0, diam( Im(z) ) > 0
760  {
761  lx_interval hyl(iimz), hyu(simz);
762  resl = ( simz > 0.0 ? Inf( Atan( hyu,srez ) + Pi_lx_interval() ) :
763  -Sup( Pi_lx_interval() ) );
764  resu = ( iimz < 0.0 ? ( simz > 0.0 ? Sup( Atan( hyl,srez ) +
765  Pi_lx_interval() ) : Sup( Atan( hyl,srez ) -
766  Pi_lx_interval() ) ) : Sup( Pi_lx_interval() ) );
767  return lx_interval( resl, resu );
768  }
769  }
770  }
771  }
772  else
773  // Arg(z) is well-defined
774  return Arg( z );
775 }
776 //
777 //-- end arg ------------------------------------------------------------------
778 
779 //-------------------- sqrt(z,n): analytic n-th root --------------------------
780  //
781  lx_interval Re_Sqrt_point(const lx_interval& rez, const lx_interval& imz,
782  int n )
783  // before: unsigned int n ---------- 040624 --
784  //
785 // Real part of analytic n-th root.
786 // Do not use this as a general function
787 // - it's only a subroutine for sqrt(z,n) and sqrt_all(z,n).
788 // The calculation is validated but largely overestimating
789 // if (rez+I*imz) is not a complex number.
790 // Blomquist, 25.05.2008;
791  {
792  lx_interval a = sqrtx2y2(rez,imz);
793  if( eq_zero(Sup(a)) )
794  // a == 0
795  return lx_interval(0);
796  else
797  return sqrt(a, n ) *
798  cos( Arg( lx_cinterval( rez, imz ) ) / n );
799  }
800 
801  lx_interval Im_Sqrt_point(const lx_interval& rez, const lx_interval& imz,
802  int n ) // before: unsigned int n --- 040624 --
803  //
804 // Imaginary part of analytic n-th root.
805 // Do not use this as a general function
806 // - it's only a subroutine for sqrt(z,n) and sqrt_all(z,n).
807 // The calculation is validated but largely overestimating
808 // if (rez+I*imz) is not a complex number.
809 // Blomquist, 25.05.2008;
810  {
811  lx_interval a = sqrtx2y2(rez,imz);
812  if( eq_zero(Sup(a)) )
813  // a == 0
814  return lx_interval(0);
815  else
816  return sqrt(a, n ) *
817  sin( Arg( lx_cinterval( rez, imz ) ) / n );
818  }
819 
820 lx_cinterval sqrt(const lx_cinterval& z, int n ) noexcept
821  //
822 // Analytic n-th root function
823 // sqrt(z,n) is undefined if z contains negative real numbers.
824  //
825 {
826  if( n == 0 ) return lx_cinterval(lx_interval(1));
827  if( n == 1 ) return z;
828  if( n == 2 ) return sqrt(z);
829  else
830  {
831  lx_real
832  irez = Inf( Re(z) ),
833  srez = Sup( Re(z) ),
834  iimz = Inf( Im(z) ),
835  simz = Sup( Im(z) );
836  lx_interval hxl( irez ), hxu( srez ), hyl( iimz ), hyu( simz );
837  lx_real resxl, resxu, resyl, resyu;
838 
839  if( sm_zero(irez) && se_zero(iimz) && ge_zero(simz) )
840  {
841  // z contains negative real values
842  cxscthrow(STD_FKT_OUT_OF_DEF("lx_cinterval sqrt(const lx_cinterval& z, int n ); z contains negative real values."));
843  return z;
844  }
845  else
846  {
847  if( sm_zero(simz) )
848  {
849  // z in lower half plane
850  lx_cinterval hres = sqrt( lx_cinterval( Re(z), -Im(z) ), n );
851  return lx_cinterval( Re(hres), -Im(hres) );
852  }
853  else
854  {
855  if( iimz > 0.0 )
856  {
857  // z in upper half plane
858  lx_interval tangle = tan( ( Pi_lx_interval() * n )
859  / ( 2 * ( n-1 ) ) );
860  lx_real tanglel = Inf( tangle ), tangleu = Sup( tangle );
861  //
862  // min( Re( Root( z ) ) )
863  //
864  if ( ge_zero(irez) || Sup( hyl / irez ) <= tanglel )
865  // lower boundary right of phi = n*Pi/(2n-2)
866  // min( Re( Root( z ) ) ) in lower left corner
867  resxl = Inf( Re_Sqrt_point( hxl, hyl, n ) );
868  else
869  {
870  if( sm_zero(srez) && Inf( hyl / srez ) >= tangleu )
871  // lower boundary left of phi = n*Pi/(2n-2)
872  // min( Re( Root( z ) ) ) in lower right corner
873  resxl = Inf( Re_Sqrt_point( hxu, hyl, n ) );
874  else
875  // lower boundary intersects phi = n*Pi/(2n-2)
876  // min( Re( Root( z ) ) ) on intersection
877  resxl = Inf( Re_Sqrt_point( iimz / tangle , hyl, n ) );
878  }
879  //
880  // max( Re( Root( z ) ) )
881  //
882  if ( ge_zero(irez) || Sup( hyu / irez ) <= tanglel )
883  // upper boundary right of phi = n*Pi/(2n-2)
884  // max( Re( Root( z ) ) ) in upper right corner
885  resxu = Sup( Re_Sqrt_point( lx_interval(srez),
886  lx_interval(simz), n ) );
887  else
888  {
889  if ( sm_zero(srez) && Inf( hyu / srez ) >= tangleu )
890  // upper boundary left of phi = n*Pi/(2n-2)
891  // max( Re( Root( z ) ) ) in upper left corner
892  resxu = Sup( Re_Sqrt_point( hxl, hyu, n ) );
893  else
894  // upper boundary intersects phi = n*Pi/(2n-2)
895  // max( Re(Root( z )) ) on upper left or right corner
896  resxu = max( Sup( Re_Sqrt_point( hxl, hyu, n ) ),
897  Sup( Re_Sqrt_point( hxu, hyu, n ) ) );
898  }
899  //
900  // min( Im( Root( z ) ) )
901  //
902  if ( ge_zero(srez) || Sup( hyl / srez ) <= tanglel )
903  // right boundary right of phi = n*Pi/(2n-2)
904  // min( Im( Root( z ) ) ) in lower right corner
905  resyl = Inf( Im_Sqrt_point( hxu, hyl, n ) );
906  else
907  {
908  if( Inf( hyu / srez ) >= tangleu )
909  // right boundary left of phi = n*Pi/(2n-2)
910  // min( Im( Root( z ) ) ) in upper right corner
911  resyl = Inf( Im_Sqrt_point( hxu, hyu, n ) );
912  else
913  // right boundary intersects phi = n*Pi/(2n-2)
914  // min( Im( Root( z ) ) ) on intersection
915  resyl = Inf(Im_Sqrt_point( hxu, srez * tangle, n ));
916  }
917  //
918  // max( Im( Root( z ) ) )
919  //
920  if( ge_zero(irez) || Sup( hyl / irez ) <= tanglel )
921  // left boundary right of phi = n*Pi/(2n-2)
922  // max( Im( Root( z ) ) ) in upper left corner
923  resyu = Sup( Im_Sqrt_point( hxl, hyu, n ) );
924  else
925  {
926  if( Inf( hyu / irez ) >= tangleu )
927  // left boundary left of phi = n*Pi/(2n-2)
928  // max( Im( Root( z ) ) ) in lower left corner
929  resyu = Sup( Im_Sqrt_point( hxl, hyl, n ) );
930  else
931  // left boundary intersects phi = n*Pi/(2n-2)
932  // max( Im(Root( z )) ) on lower or upper left corner
933  resyu = max( Sup( Im_Sqrt_point( hxl, hyl, n ) ),
934  Sup( Im_Sqrt_point( hxl, hyu, n ) ) );
935  }
936  }
937  else
938  {
939  // z intersects positive real axis
940  // min( Re( Root( z ) ) ) on real line
941  // max( Re( Root( z ) ) ) in lower or upper right corner
942  resxl = ( eq_zero(irez)? lx_real(0.0) : Inf( sqrt( hxl, n ) ) );
943  resxu = ( - iimz > simz ? Sup( Re_Sqrt_point( hxu, hyl, n ) ) :
944  Sup( Re_Sqrt_point( hxu, hyu, n ) ) );
945  // min( Im ( sqrt( z ) ) ) in lower left corner
946  // max( Im ( sqrt( z ) ) ) in upper left corner
947  resyl = Inf( Im_Sqrt_point( hxl, hyl, n ) );
948  resyu = Sup( Im_Sqrt_point( hxl, hyu, n ) );
949  }
950  return lx_cinterval( lx_interval( resxl, resxu ),
951  lx_interval( resyl, resyu ) );
952  }
953  }
954  }
955 }
956 //
957 //-- end sqrt(z,n) -----------------------------------------------------------
958 
959 //-- sqrt_all ------------------------------------------------------- 040628 --
960  //
961  std::list<lx_cinterval> sqrt_all( const lx_cinterval& z, int n ) noexcept
962  //
963 // sqrt_all(z,n) computes a list of n intervals containing all n-th roots of z
964  //
965 // For n >=3, computing the optimal interval bounds is very expensive
966 // and thus not considered cost-effective.
967  //
968 // Hence, the polar form is used to calculate sqrt_all(z,n)
969 // (observed overestimation of Re() and Im() in test cases: 5-15%).
970  //
971 // z is enclosed into an interval in polar coordinates
972 // (i.e. a sector of an annulus), from which the roots
973 // (also polar intervals) are computed error-free
974 // (save roundoff, which is enclosed).
975  //
976 // The the final inclusion of the roots into a rectangular complex interval
977 // involves some additional overestimation.
978  //
979  {
980  std::list<lx_cinterval> res;
981 
982  if( n == 0 )
983  {
984  res.push_back( lx_cinterval( lx_interval(0,l_interval(1)),
985  lx_interval(0,l_interval(0)) ) );
986  return res;
987  }
988  else if( n == 1 )
989  {
990  res.push_back(z);
991  return res;
992  }
993  else if( n == 2 ) return sqrt_all( z );
994  else
995  {
996  lx_interval
997  arg_z = arg( z ), root_abs_z = sqrt( abs(z), n );
998 
999  for(int k = 0; k < n; k++)
1000  {
1001  lx_interval arg_k = ( arg_z + 2 * k * Pi_l_interval() ) / n;
1002 
1003  res.push_back( lx_cinterval( root_abs_z * cos( arg_k ),
1004  root_abs_z * sin( arg_k ) ) );
1005  }
1006  return res;
1007  }
1008  }
1009  //
1010 //-- end sqrt_all -------------------------------------------------------------
1011 
1012 
1013 //-----------------------------------------------------------------------------
1014  //
1015 // Section 2: Logarithms
1016  //
1017 //-----------------------------------------------------------------------------
1018 
1019 // ------------------- Ln: analytic natural logarithm -------------------------
1020 //
1021 // Ln(z) is undefined if z contains zero; z must not touch the negative real
1022 // axis from below;
1023 //
1024 lx_cinterval Ln(const lx_cinterval& z) noexcept
1025 { // Blomquist, 24.09.2007;
1026  int stagsave = stagprec,
1027  stagmax = 30;
1028  if (stagprec > stagmax) stagprec = stagmax;
1029 
1030  lx_cinterval y;
1031  lx_real
1032  srez = Sup( Re(z) ),
1033  simz = Sup( Im(z) ),
1034  iimz = Inf( Im(z) );
1035  lx_interval a1( abs(Re(z)) ), a2( abs(Im(z)) );
1036  if ( eq_zero(Inf(a1)) && eq_zero(Inf(a2)) )
1037  // z contains 0
1038  cxscthrow(STD_FKT_OUT_OF_DEF(
1039  "lx_cinterval Ln(const lx_cinterval& z); z contains 0"));
1040  if ( sm_zero(srez) && sm_zero(iimz) && ge_zero(simz) )
1041  cxscthrow(STD_FKT_OUT_OF_DEF(
1042  "lx_cinterval Ln(const lx_cinterval& z); z not allowed"));
1043 
1044  y = lx_cinterval( ln_sqrtx2y2(Re(z),Im(z)), arg(z) );
1045 
1046  stagprec = stagsave;
1047  y = adjust(y);
1048 
1049  return y;
1050 }
1051 //
1052 //-- end Ln -------------------------------------------------------------------
1053 
1054 //-- ln: non-analytic natural logarithm ----------------------------- 040923 --
1055 //
1056 // ln(z) is undefined if z contains zero.
1057 //
1058 lx_cinterval ln(const lx_cinterval& z) noexcept
1059 {
1060  int stagsave = stagprec,
1061  stagmax = 30;
1062  if (stagprec > stagmax) stagprec = stagmax;
1063 
1064  lx_cinterval y;
1065  lx_interval a1( abs(Re(z)) ), a2( abs(Im(z)) );
1066  if ( eq_zero(Inf(a1)) && eq_zero(Inf(a2)) )
1067  // z contains 0
1068  cxscthrow(STD_FKT_OUT_OF_DEF
1069  ("lx_cinterval ln(const lx_cinterval& z); z contains 0"));
1070  y = lx_cinterval( ln_sqrtx2y2(Re(z),Im(z)), arg(z) );
1071 
1072  stagprec = stagsave;
1073  y = adjust(y);
1074 
1075  return y;
1076 }
1077 //
1078 //-- end ln -------------------------------------------------------------------
1079 
1080 // ---------------------- log2, log10: analytic functions ---------------------
1081 //
1082 // log2(z),log10 are undefined if z contains zero; z must not touch the
1083 // negative real axis from below;
1084 //
1085 lx_cinterval log2(const lx_cinterval& z) noexcept
1086 { // Blomquist, 30.11.2008;
1087  return Ln(z) / Ln2_lx_interval();
1088 }
1089 
1090 lx_cinterval log10(const lx_cinterval& z) noexcept
1091 { // Blomquist, 30.11.2008;
1092  return Ln(z) / Ln10_lx_interval();
1093 }
1094 
1095 //-----------------------------------------------------------------------------
1096 //
1097 // Section 3: Power functions
1098 //
1099 //-----------------------------------------------------------------------------
1100 
1101 //-- power_fast ---------------------------------------------------------------
1102 //
1103 // Fast, validated power function for integer powers, based on Moivre's formula.
1104 // If n is not an integer value, an error message is generated.
1105 // Medium amount of overestimation.
1106 //
1107 lx_cinterval power_fast(const lx_cinterval& z, const real& n) noexcept
1108 {
1109  if( n == 0 ) return lx_cinterval(lx_interval(1));
1110  else if( n == 1 ) return z;
1111  else if( n == -1 ) return 1 / z;
1112  else if( n == 2 ) return sqr(z);
1113  else
1114  // n >= 3 or n <= -2
1115  // If z is a large interval, z^n is a distorted set, for which the
1116  // inclusion into a complex rectangle contains large overestimation.
1117  // For example, if n * Arg( z ) intersects the cartesian axes
1118  // more than once then 0 is contained in the rectangular enclosure
1119  // of z^n.
1120  // For n <= -2, also inversion of z or z^n is required;
1121  // both inversions lead to large overestimation of the resulting interval.
1122  //
1123  // In power_fast(z,n), z is enclosed into a sector S of
1124  // an annulus. S^n is again some sector S' of a (different) annulus.
1125  // S' is calculated exactly (apart from roundoff), and then enclosed
1126  // into a rectangle. There is a certain amount of overestimation
1127  // but the calculations are quit cheap.
1128  // The implementation is base on Moivre's formula.
1129  {
1130  lx_interval abs_z = abs(z);
1131 
1132  if( ((n < 0) && eq_zero(Inf(abs_z))) || !(Is_Integer(n)))
1133  // z contains 0
1134  cxscthrow (STD_FKT_OUT_OF_DEF("lx_cinterval power_fast(const lx_cinterval& z, const real& n); z contains 0 or n is not integer."));
1135 
1136  if( eq_zero(Sup(abs_z)) ) return lx_cinterval( lx_interval(0));
1137  else
1138  {
1139  lx_interval arg_z = arg(z);
1140  lx_interval abs_z_n = power(abs_z,n); // Legendre algorithm
1141 
1142  return lx_cinterval( abs_z_n * cos( n * arg_z ),
1143  abs_z_n * sin( n * arg_z ) );
1144  }
1145  }
1146 }
1147 //
1148 //-- end power_fast -----------------------------------------------------------
1149 
1150 //------------------------------------ power ----------------------------------
1151 //
1152 // power function based on the Legendre algotizkm with an integer valued
1153 // exponent n of type real.
1154 //
1155 
1156  lx_cinterval power( const lx_cinterval& x, const real& n ) noexcept
1157  {
1158  if( !(Is_Integer(n)) )
1159  // n is not an integer
1160  cxscthrow(STD_FKT_OUT_OF_DEF("lx_cinterval power(const lx_cinterval& z, const real& n); n is not integer."));
1161 
1162  real one(1.0), zhi(2.0), N(n), r;
1163  double dbl;
1164  lx_cinterval y,neu,X(x);
1165 
1166  if (x == one) y = x;
1167  else
1168  if (N == 0.0) y = one;
1169  else
1170  {
1171  if (N == 1) y = x;
1172  else
1173  if (N == 2) y = sqr(x);
1174  else
1175  {
1176  if (N < 0)
1177  {
1178  X = 1.0/X;
1179  N = -N;
1180  }
1181  // Initialisierung
1182  if ( !Is_Integer(N/2) )
1183  y = X;
1184  else y = one;
1185  neu = sqr(X); // neu = X*X;
1186  do {
1187  dbl = _double(N/zhi);
1188  dbl = floor(dbl);
1189  r = (real) dbl;
1190  if ( !Is_Integer( r/2 ) )
1191  y *= neu;
1192  zhi += zhi;
1193  if (zhi <= N)
1194  neu = sqr(neu); // neu = neu * neu;
1195  } while (zhi <= N);
1196  }
1197  }
1198  return y;
1199  }
1200  //
1201 //-- end power ----------------------------------------------------------------
1202 
1203 //------------------------------------ pow ------------------------------------
1204  //
1205 // Analytic power function for a real interval exponent, based on Ln.
1206  //
1207 
1208  lx_cinterval pow( const lx_cinterval& z, const lx_interval& p ) noexcept
1209  {
1210  return exp( p*Ln(z) );
1211  }
1212 
1213  //
1214 //-- end pow ------------------------------------------------------------------
1215 
1216 //------------------------------------ pow ------------------------------------
1217  //
1218 // Analytic power function for a complex interval exponent, based on Ln.
1219  //
1220 
1221  lx_cinterval pow(const lx_cinterval& z, const lx_cinterval& p) noexcept
1222  {
1223  return exp( p*Ln(z) );
1224  }
1225 
1226 //
1227 //-- end pow ------------------------------------------------------------------
1228 
1229 //-----------------------------------------------------------------------------
1230 //
1231 // Section 2: tan, cot, tanh, coth
1232 //
1233 // The implementation of cot, tanh, and coth is based on tan:
1234 //
1235 // cot( z ) = tan( pi/2 - z )
1236 // tanh( z ) = transp( i * tan( transp( i * z ) )
1237 // coth( z ) = i * cot( i * z ) = i * tan( pi/2 - i * z )
1238 //
1239 //-----------------------------------------------------------------------------
1240 
1241 //---------------------------------- tan --------------------------------------
1242 //
1243 // Complex tangent function
1244 //
1245 
1246 void horizontal_check( const lx_interval& hy, const lx_interval& cosh_2y,
1247  lx_real irez, lx_real srez, const lx_interval& hxl,
1248  const lx_interval& hxu, lx_real& resxl, lx_real& resxu )
1249 //
1250 // Subroutine of tangent function.
1251 // Check intersections with extremal curves of tan on a horizontal boundary.
1252 // This subroutine is only called if an intersection occurs.
1253 // In this case, sinh( 2 * hy ) <> 0.0 (poles are handled before).
1254 //
1255 // There may be 1 or 2 intersections.
1256 // If intersections lie next to a boundary of rez, then it is impossible to
1257 // decide if there are 1 or 2 intersections.
1258 // In this case, 2 intersections are assumed
1259 // (valid enclosure, at the expense of a potential slight overestimation).
1260 //
1261 {
1262  bool both = false, left = false, right = false;
1263  lx_interval hlp, hlp1;
1264 
1265  hlp = Pid2_lx_interval();
1266  hlp1 = lx_interval(srez) - lx_interval(irez);
1267  if (Inf(hlp1) > Inf(hlp))
1268  // 2 intersections
1269  both = true;
1270  else
1271  {
1272  lx_interval
1273  res_l = cos( 2 * hxl ) - cosh_2y,
1274  res_u = cos( 2 * hxu ) - cosh_2y;
1275 
1276  if( Inf( res_l * res_u ) > 0.0 )
1277  // 2 intersections
1278  both = true;
1279  else
1280  {
1281  if( Sup( res_l * res_u ) < 0.0 )
1282  {
1283  // 1 intersection (3 intersections are PI() apart)
1284  // neither of the intervals res_l and res_u contains zero
1285  if( Inf( res_l ) > 0.0 )
1286  // "left" intersection
1287  left = true;
1288  else
1289  // "right" intersection
1290  right = true;
1291  }
1292  else
1293  //
1294  // 1 (or both) intersections lie next to a boundary point
1295  // here, the task is to decide if 2 intersections occurs
1296  // if only 1 intersection takes place, then which one?
1297  //
1298  {
1299  lx_interval
1300  sin_2xl = sin( 2 * hxl ),
1301  sin_2xu = sin( 2 * hxu );
1302 
1303  if( !Disjoint( lx_interval(0), res_l ) )
1304  // intersection on the left boundary
1305  {
1306  if( ge_zero(Inf(sin_2xl)) )
1307  // "left" intersection
1308  {
1309  left = true;
1310  // remove the intersection by changing res_l!
1311  res_l = -lx_interval(1);
1312  }
1313  else
1314  {
1315  if( se_zero(Sup(sin_2xl)) )
1316  // "right" intersection
1317  {
1318  right = true;
1319  // remove the intersection by changing res_l!
1320  res_l = lx_interval(1);
1321  }
1322  else
1323  // zero is interior point of sin_2xl
1324  // if the real sine function has optimal precision,
1325  // this case should never happen
1326  both = true;
1327  }
1328  }
1329 
1330  if( !Disjoint( lx_interval(0), res_u ) )
1331  // intersection on the right boundary
1332  {
1333  if( ge_zero(Inf(sin_2xu)) )
1334  // "left" intersection
1335  {
1336  left = true;
1337  // remove the intersection by changing res_u!
1338  res_u = lx_interval(1);
1339  }
1340  else
1341  {
1342  if( se_zero(Sup(sin_2xu)) )
1343  // "right" intersection
1344  {
1345  right = true;
1346  // remove the intersection by changing res_u!
1347  res_u = -lx_interval(1);
1348  }
1349  else
1350  // zero is interior point of sin_2xu
1351  // if the real sine function has optimal precision,
1352  // this case should never happen
1353  both = true;
1354  }
1355  }
1356  //
1357  // finally, check if there is a second intersection
1358  //
1359  if( Inf( res_l * res_u ) < 0.0 )
1360  both = true;
1361  }
1362  }
1363  }
1364  //
1365  // Calculate extremal values
1366  //
1367  lx_interval re_tan = 1 / sinh( 2 * abs( hy ) );
1368 
1369  // "left" intersection, sin( 2x ) > 0
1370  if( left || both )
1371  {
1372  resxl = min( resxl, Inf( re_tan ) );
1373  resxu = max( resxu, Sup( re_tan ) );
1374  }
1375 
1376  // "right" intersection, sin( 2x ) < 0
1377  if( right || both )
1378  {
1379  resxl = min( resxl, -Sup( re_tan ) );
1380  resxu = max( resxu, -Inf( re_tan ) );
1381  }
1382 } // end horizontal_check
1383 
1384 
1385 lx_cinterval Tan(const lx_cinterval& z) noexcept
1386 {
1387  lx_cinterval y;
1388  lx_interval
1389  rez = Re(z), // rez = z.re(),
1390  imz = Im(z); // imz = z.im();
1391 
1392  lx_real
1393  irez = Inf(rez),
1394  srez = Sup(rez),
1395  iimz = Inf(imz),
1396  simz = Sup(imz);
1397 
1398  lx_interval hxl(irez), hxu(srez), hyl(iimz), hyu(simz);
1399 
1400  lx_real resxl, resxu, resyl, resyu;
1401  //
1402  // 1st: check for poles
1403  //
1404  if( ( !Disjoint( lx_interval(0), imz ) ) &&
1405  ( !Disjoint( lx_interval(0), cos( rez ) ) ) )
1406  cxscthrow (STD_FKT_OUT_OF_DEF(
1407  "lx_cinterval tan(const lx_cinterval& z); Pole(s) in z"));
1408  //
1409  // 2nd: real part
1410  //
1411  // evaluate tan on vertical boundaries
1412  //
1413  lx_interval
1414  cos_2rez = cos( 2 * rez ), sinh_imz_2 = sqr( sinh( imz ) );
1415 
1416  lx_interval
1417  re_tan_l=sin( 2*hxl ) / ( 2*( sqr( cos( hxl ) ) + sinh_imz_2 ) ),
1418  re_tan_u=sin( 2*hxu ) / ( 2*( sqr( cos( hxu ) ) + sinh_imz_2 ) );
1419 
1420  resxl = min( Inf( re_tan_l ), Inf( re_tan_u ) );
1421  resxu = max( Sup( re_tan_l ), Sup( re_tan_u ) );
1422 
1423  //
1424  // look for extremal values on horizontal boundaries
1425  // if one of the horizontal boundaries is the x-axes,
1426  // then the complex tangent is the real tangent on this
1427  // boundary, and due to monotonicity, its range
1428  // is already included in the ranges of the vertical boundaries
1429  //
1430  if( irez < srez )
1431  {
1432  lx_interval
1433  cosh_2yl = - 1 / cosh( 2 * hyl ),
1434  cosh_2yu = - 1 / cosh( 2 * hyu );
1435 
1436  if( !Disjoint( cos_2rez, cosh_2yl ) && iimz != 0.0 )
1437  //extremal curve intersects lower boundary
1438  horizontal_check(hyl,cosh_2yl,irez,srez,hxl,hxu,resxl,resxu);
1439 
1440  if( !Disjoint( cos_2rez, cosh_2yu ) && simz != 0.0 )
1441  //extremal curve intersects upper boundary
1442  horizontal_check(hyu,cosh_2yu,irez,srez,hxl,hxu,resxl,resxu);
1443  }
1444  //
1445  // 3rd: imaginary part
1446  //
1447  // evaluate tan on horizontal boundaries
1448  //
1449  lx_interval cos_rez_2 = sqr( cos( rez ) );
1450 
1451  lx_interval
1452  im_tan_l = sinh( 2*hyl ) / (2*(cos_rez_2 + sqr( sinh( hyl ) ))),
1453  im_tan_u = sinh( 2*hyu ) / (2*(cos_rez_2 + sqr( sinh( hyu ) )));
1454 
1455  resyl = min( Inf( im_tan_l ), Inf( im_tan_u ) );
1456  resyu = max( Sup( im_tan_l ), Sup( im_tan_u ) );
1457 
1458  //
1459  // look for extremal values on vertical boundaries
1460  // here, the situation is simpler than for the real part
1461  // if 2 intersections with extremal curves appear ,
1462  // one lies in the lower half plane, the other in the upper half plane
1463  //
1464  lx_interval
1465  cos_2xl = cos( 2 * hxl ), cos_2xu = cos( 2 * hxu );
1466  lx_interval im_tan;
1467 
1468  if( iimz < 0.0 )
1469  // intersection in lower half plane?
1470  {
1471  lx_interval
1472  imz_h = lx_interval( iimz, min( simz, lx_real(0.0) ) ),
1473  cosh_2imz = - 1 / cosh( 2 * imz_h );
1474 
1475  if( ( !Disjoint( cosh_2imz, cos_2xl ) ) )
1476  //extremal curve intersects left boundary
1477  //in this case, sin( 2 * xl ) <> 0.0 (no poles here!)
1478  {
1479  im_tan = - 1 / abs( sin( 2 * hxl ) );
1480  resyl = min( resyl, Inf( im_tan ) );
1481  resyu = max( resyu, Sup( im_tan ) );
1482  }
1483  if( ( !Disjoint( cosh_2imz, cos_2xu ) ) )
1484  //extremal curve intersects right boundary
1485  //in this case, sin( 2 * xu ) <> 0.0 (no poles here!)
1486  {
1487  im_tan = - 1 / abs( sin( 2 * hxu ) );
1488  resyl = min( resyl, Inf( im_tan ) );
1489  resyu = max( resyu, Sup( im_tan ) );
1490  }
1491  }
1492  if( simz > 0.0 )
1493  // intersection in upper half plane?
1494  {
1495  lx_interval
1496  imz_h = lx_interval( max( iimz, lx_real(0.0) ), simz ),
1497  cosh_2imz = - 1 / cosh( 2 * imz_h );
1498 
1499  if( ( !Disjoint( cosh_2imz, cos_2xl ) ) )
1500  //extremal curve intersects left boundary
1501  //in this case, sin( 2 * xl ) <> 0.0 (no poles here!)
1502  {
1503  im_tan = + 1 / abs( sin( 2 * hxl ) );
1504  resyl = min( resyl, Inf( im_tan ) );
1505  resyu = max( resyu, Sup( im_tan ) );
1506  }
1507  if( ( !Disjoint( cosh_2imz, cos_2xu ) ) )
1508  //extremal curve intersects right boundary
1509  //in this case, sin( 2 * xu ) <> 0.0 (no poles here!)
1510  {
1511  im_tan = + 1 / abs( sin( 2 * hxu ) );
1512  resyl = min( resyl, Inf( im_tan ) );
1513  resyu = max( resyu, Sup( im_tan ) );
1514  }
1515  }
1516 
1517  y = lx_cinterval( lx_interval(resxl,resxu ),lx_interval(resyl,resyu ) );
1518 
1519  return y;
1520 } // Tan
1521 
1522 lx_cinterval tan(const lx_cinterval& z) noexcept
1523 {
1524 // tan(z) has the poles z_k = pi*(1+2k)/2; |k| in {0,1,2,3,...}.
1525 // z = z_k + eps = pi*(1+2k)/2 + eps; With |eps|<<1 k can be calculated
1526 // by: k = Re(z)/pi - 0.5; With this k the complex value eps is given
1527 // by: eps = z - pi*(1+2k)/2; pi = 3.1415926... ;
1528 // It holds:
1529 // tan(z) = tan(z_k+eps) = tan[pi*(1+2k)/2 + eps]
1530 // = tan[pi/2 + pi*k + eps] = tan[pi/2 + eps] = -1 / tan(eps);
1531 // Definitions: u = Re(eps); u = abs(u); v = Im(eps); v = abs(v);
1532 // if (Sup(u)<S && Sup(v)<S) tan(z) = -1 / Tan(eps);
1533 // else tan(z) = Tan(z); S = 1e-15;
1534 // Thus, near the poles tan(z) is calculated in higher accuracy with
1535 // -1 / Tan(eps);
1536 // Blomquist, 27.09.2007;
1537  const real S = 1e-15;
1538  int stagsave = stagprec,
1539  stagmax = 39;
1540  if (stagprec > stagmax)
1541  stagprec = stagmax;
1542 
1543  lx_cinterval y,eps;
1544  l_interval rezl,ul,vl;
1545 
1546  lx_interval
1547  rez = Re(z), Pih; // rez = z.re(),
1548  rezl = rez;
1549  interval re,u_,v_;
1550  re = rezl;
1551  real x(mid(re)),k, pi(3.14159265358979323);
1552  lx_interval u,v;
1553 
1554  int s;
1555  x = x/pi - 0.5;
1556  s = sign(x);
1557  k = (s>=0)? cutint(x+0.5) : cutint(x-0.5);
1558  if (k == 9007199254740992.0)
1559  cxscthrow (STD_FKT_OUT_OF_DEF(
1560  "lx_cinterval tan(const lx_cinterval& z); z out of range"));
1561  Pih = Pi_lx_interval();
1562  rez = k*Pih;
1563  times2pown(Pih,-1); // Pih = pi/2
1564  rez = rez + Pih; // rez: inclusion of k*pi + pi/2
1565  eps = z - rez;
1566 
1567  u = Re(eps); u = abs(u);
1568  v = Im(eps); v = abs(v);
1569  ul = u; vl = v;
1570  u_ = ul; v_ = vl;
1571  y = (Sup(u_)<S && Sup(v_)<S)? -lx_cinterval(1) / Tan(eps) : Tan(z);
1572 
1573  stagprec = stagsave;
1574  y = adjust(y);
1575 
1576  return y;
1577  } // tan()
1578 
1579  lx_cinterval cot(const lx_cinterval& z) noexcept
1580  {
1581 // cot(z) has the poles z_k = k*pi; |k| in {0,1,2,3,...}.
1582 // z = z_k + eps = k*pi + eps; With |eps|<<1 k can be calculated
1583 // by: k = Re(z)/pi; With this k the komplex value eps is given
1584 // by: eps = z - k*pi; pi = 3.1415926... ;
1585 // It holds:
1586 // cot(z) = cot(z_k+eps) = cot(k*pi + eps)
1587 // = cot(eps) = 1 / tan(eps);
1588 // Definitions: u = Re(eps); u = abs(u); v = Im(eps); v = abs(v);
1589 // if (Sup(u)<S && Sup(v)<S) cot(z) = 1 / tan(eps);
1590 // else cot(z) = tan(pi/2 - z); S = 1e-15;
1591 // Thus, near the poles cot(z) is calculated in higher accuracy with
1592 // 1 / tan(eps);
1593 // Blomquist, 28.09.2007;
1594  const real S = 1e-15;
1595  int stagsave = stagprec,
1596  stagmax = 39;
1597  if (stagprec > stagmax)
1598  stagprec = stagmax;
1599 
1600  lx_cinterval y,eps;
1601  l_interval rezl,ul,vl;
1602 
1603  lx_interval rez = Re(z);
1604  rezl = rez;
1605  interval re,u_,v_;
1606  re = rezl;
1607  real x(mid(re)), k, pi(3.14159265358979323);
1608  lx_interval u,v;
1609 
1610  int s;
1611  x = x/pi;
1612  s = sign(x);
1613  k = (s>=0)? cutint(x+0.5) : cutint(x-0.5);
1614  if (k == 9007199254740992.0)
1615  cxscthrow (STD_FKT_OUT_OF_DEF(
1616  "lx_cinterval cot(const lx_cinterval& z); z out of range"));
1617  eps = z - Pi_lx_interval()*k;
1618  u = Re(eps); u = abs(u);
1619  v = Im(eps); v = abs(v);
1620  ul = u; vl = v;
1621  u_ = ul; v_ = vl;
1622  if (Sup(u)<S && Sup(v)<S)
1623  y = 1 / Tan(eps);
1624  else
1625  {
1626  u = Pi_lx_interval();
1627  times2pown(u,-1); // u = pi/2
1628  y = Tan(u - z);
1629  }
1630 
1631  stagprec = stagsave;
1632  y = adjust(y);
1633 
1634  return y;
1635  } // cot
1636 
1637 //---------------------------------- tanh -----------------------------------
1638 //
1639 // tanh( z ) = transp( i * tan( transp( i * z ) )
1640 //
1641 lx_cinterval tanh(const lx_cinterval& z) noexcept
1642 {
1643  int stagsave = stagprec,
1644  stagmax = 39;
1645  if (stagprec > stagmax) stagprec = stagmax;
1646 
1647  lx_cinterval res = tan( lx_cinterval( Im(z), Re(z) ) ),y;
1648  y = lx_cinterval( Im(res), Re(res) );
1649 
1650  stagprec = stagsave;
1651  y = adjust(y);
1652 
1653  return y;
1654 }
1655 //
1656 //-- end tanh -------------------------------------------------------
1657 
1658 //-- coth -----------------------------------------------------------
1659 //
1660 // coth( z ) = i * cot( i * z );
1661 //
1662 
1663 lx_cinterval coth(const lx_cinterval& z) noexcept
1664 { // coth( z ) = i * cot( i * z );
1665  int stagsave = stagprec,
1666  stagmax = 39;
1667  if (stagprec > stagmax) stagprec = stagmax;
1668 
1669  lx_cinterval zh = lx_cinterval( -Im(z), Re(z) ); // zh = i*z;
1670  lx_cinterval res = cot(zh);
1671  zh = lx_cinterval( -Im(res), Re(res) );
1672 
1673  stagprec = stagsave;
1674  zh = adjust(zh);
1675 
1676  return zh;
1677 }
1678 //
1679 //-- end coth -----------------------------------------------------------------
1680 
1681 //-----------------------------------------------------------------------------
1682 //-----------------------------------------------------------------------------
1683 //
1684 // Section 5: asin, acos, asinh, acosh
1685 //
1686 // The implementation of acos, asinh, and acosh is based on asin:
1687 //
1688 // acos( z ) = pi / 2 - asin( z )
1689 // asinh( z ) = i * asin( - i * z )
1690 // acosh( z ) = i * acos( z ) = +/- i * ( pi / 2 - asin( z ) )
1691 //
1692 // Only principal values are computed.
1693 //
1694 //-----------------------------------------------------------------------------
1695 
1696 //------------------------------------ asin -----------------------------------
1697 //
1698 // Analytic inverse sine function
1699 //
1700 lx_interval f_aux_asin(const lx_interval& x, const lx_interval& y)
1701 //
1702 // auxiliary function for evaluating the imaginary part of asin(z);
1703 // f_aux_asin(z) = ( |z+1| + |z-1| ) / 2; z = x + i*y;
1704 // = ( sqrt[(x+1)^2 + y^2] + sqrt[(x-1)^2 + y^2] )/2;
1705 // Blomquist, 01.10.2007;
1706  {
1707  lx_interval res;
1708 // interval sqr_y = sqr( y );
1709 // interval res = ( sqrt( sqr( x + 1.0 ) + sqr_y ) +
1710 // sqrt( sqr( x - 1.0 ) + sqr_y ) ) / 2.0;
1711  res = abs(x);
1712  if (y != 0.0 || Inf(res) < 1.0)
1713  {
1714  res = sqrtx2y2(x+1.0,y) + sqrtx2y2(x-1.0,y); // Blomquist;
1715  times2pown(res,-1);
1716  }
1717 
1718  // Now we correct a possible overestimation of the lower bound
1719  // of res.
1720  // It holds:
1721  // (sqrt[(x+1)^2 + y^2] + sqrt[(x-1)^2 + y^2])/2 >= Max(1,|x|);
1722  lx_real hlb = max( lx_real(1.0), abs( Sup(x) ) );
1723  if( Inf(res) < hlb ) // this is an invalid overestimation!
1724  res = lx_interval( hlb, Sup(res) );
1725  return res;
1726  }
1727 
1728  lx_interval ACOSH_f_aux(const lx_interval& x, const lx_interval& y)
1729 // Function for calculating the imaginary part of asin(z),
1730 // z = x + i*y;
1731 // Notations:
1732 // f_aux_asin(x,y) = alpha(x,y);
1733 // alpha(x,y) := ( sqrt[(x+1)^2 + y^2] + sqrt[(x-1)^2 + y^2] )/2;
1734 // return values:
1735 // 1. res = acosh( f_aux_asin(x,y) ), if |x|>2 or |y|>2;
1736 // 2. res = ln[ alpha + sqrt(alpha^2 - 1) ]; d = alpha(x,y) -1
1737 // = ln[1 + sqrt(d)*(sqrt(d) + sqrt(2+d))], if it holds:
1738 // |x|<=2 und |y|<=2; x,y are point intervals only!
1739 // Blomquist, 02.06.2008;
1740 
1741  {
1742  lx_interval res, xa(abs(x)), ya(abs(y)), S1, S2, S3;
1743  lx_real rx((Inf(xa))), ry(Inf(ya));
1744 
1745  if (rx>2.0 || ry>2.0) res = acosh( f_aux_asin(x,y) );
1746  else
1747  { // Now with improvements:
1748  if (rx == 1.0)
1749  {
1750  res = ya/2; S2 = res/2;
1751  S1 = ya*(0.5 + S2/(sqrt1px2(res)+1)); // S1 = delta
1752  res = lnp1(S1 + sqrt(S1*(2+S1)));
1753  }
1754  else
1755  if (rx<1.0) // 0 <= x < +1;
1756  {
1757  S1 = 1+xa;
1758  S1 = 1/(sqrtx2y2(S1,ya) + S1);
1759  S2 = 1-xa;
1760  S2 = 1/(sqrtx2y2(S2,ya) + S2);
1761  S1 = S1 + S2;
1762  times2pown(S1,-1); // S1 = {...}/2
1763  S2 = ya * sqrt(S1); // S2 = sqrt(delta)
1764  res = lnp1( S2*(S2 + sqrt(2+sqr(ya)*S1)) );
1765  }
1766  else // 1 < x <= 2;
1767  {
1768  if (ya == 0)
1769  {
1770  S1 = xa-1;
1771  if (eq_zero(Inf(S1)))
1772  S1 = lx_interval(Inf(lx_interval(-2097,l_interval(1))),
1773  Sup(S1));
1774  }
1775  else // 1 < x <= 2 and 0 < |y| <= 2
1776  {
1777  S1 = xa-1;
1778  if (eq_zero(Inf(S1)))
1779  S1 = lx_interval(Inf(lx_interval(-2097,l_interval(1))),
1780  Sup(S1));
1781  S2 = ya;
1782  times2pown(S2,1048);
1783  if (Inf(S1)>Inf(S2))
1784  S1 *= lx_interval( lx_real(1),Sup(One_p_lx_interval()) );
1785  else
1786  {
1787  res = sqr(ya);
1788  S3 = res / (sqrtx2y2(S1,ya) + S1); // S3 = v;
1789  S2 = 1+xa;
1790  S2 = res / (sqrtx2y2(S2,ya) + S2); // S2 = u;
1791  times2pown(S1,1); // S1 = w=2*(x-1)
1792  S1 = (S3 + S2) + S1; // S1 = u+v+w
1793  times2pown(S1,-1); // S1 = delta
1794  }
1795  }
1796  res = lnp1( S1+sqrt(S1*(2+S1)) );
1797  }
1798  }
1799  return res;
1800  } // ACOSH_f_aux
1801 
1802  lx_interval Asin_beta(const lx_interval& x, const lx_interval& y )
1803 // Calculating the improved real part of asin(z); Blomquist 22.01.2007;
1804 // Re(asin(z)) = asin[ 2x/(sqrt[(x+1)^2+y^2] + sqrt[(x-1)^2+y^2]) ]=asin[beta];
1805 // Improvements for beta --> +1 and beta --> -1 are necessary because of
1806 // nearly vertical tangents of the real asin(t)-function for |t|-->+1;
1807 // This function should only be used for POINT intervals x,y! z = x + i*y;
1808 // Blomquist, 30.05.2008;
1809  {
1810  const real c1 = 0.75;
1811  bool neg_x, fertig(false);
1812  lx_real Infxa(Inf(x)), L_y(Inf(y));
1813  int ex_xl(expo_gr(lr_part(Infxa))),
1814  ex_yl(expo_gr(lr_part(L_y)));
1815  real ex_x(expo(Infxa)), ex_y(expo(L_y));
1816  lx_interval res,beta,abs_beta,delta,w_delta,xa,Ne,Kla;
1817 
1818  if (ex_xl<-1000000)
1819  res = 0;
1820  else
1821  { // x != 0
1822  if ( ex_yl>-1000000 && ex_y>=9007199254739972.0-ex_yl &&
1823  ex_x <= 2097-ex_xl )
1824  { // It holds: |x/y| <= 2^(-9007199254737871) and |y| != 0
1825  // res in [0,eps] or res in [-eps,0], with
1826  // eps = 2^(-9007199254737870);
1827  xa = x;
1828  neg_x = Inf(x)<0;
1829  if (neg_x) xa = -xa; // Inf(xa) >0 :
1830  res = lx_interval( lx_real(0), lx_real(-9007199254737870.0,l_real(1)) );
1831  if (neg_x)
1832  res = -res;
1833  }
1834  else
1835  {
1836  Kla = sqrtx2y2(1-x,y);
1837  Ne = sqrtx2y2(1+x,y) + Kla;
1838  beta = x / ( Ne/2 );
1839  if (Inf(beta)<-1) Inf(beta) = -1;
1840  if (Sup(beta)> 1) Sup(beta) = 1;
1841  abs_beta = abs(beta);
1842  if (Inf(abs_beta) < c1)
1843  res = asin(beta); // Normal calculation
1844  else
1845  { // Inf(abs_beta)>=c1; Calculation now with improvements:
1846  Ne = Ne + 2.0;
1847  // Ne = 2 + sqrt[(1+x)^2 + y^2] + sqrt[(1-x)^2 + y^2];
1848  xa = x;
1849  neg_x = Inf(x)<0;
1850  if (neg_x) xa = -xa; // Inf(xa) >0 :
1851  Infxa = Inf(xa);
1852  if (Infxa > 1.0)
1853  { // x > +1;
1854  if (y == 0.0)
1855  fertig = true; // due to delta == 0;
1856  else
1857  {
1858  beta = xa - 1;
1859  delta = sqrt(beta);
1860  delta = lx_interval(lx_real(-2100,7.4699)) * delta;
1861  // delta = 7.4699*(2^-2100)*sqrt(x-1);
1862  if (Sup(abs(y)) < Inf(delta))
1863  fertig = true;
1864  else
1865  {
1866  delta = sqr(y/beta); // delta = (y/(x-1))^2
1867  delta = beta * sqrtp1m1(delta);
1868  times2pown(Ne,-1);
1869  delta /= Ne;
1870  }
1871  }
1872  }
1873  else
1874  if (Infxa == 1.0)
1875  { // x = 1;
1876  delta = abs(y);
1877  times2pown(delta,1);
1878  delta /= Ne;
1879  }
1880  else
1881  { // 0.75 <= x < 1;
1882  if (y == 0.0)
1883  delta = 1 - xa;
1884  else // 0.75 <= x < 1 and y != 0:
1885  {
1886  beta = 1 - xa; // beta > 0;
1887  delta = Kla + beta;
1888  times2pown(delta,1);
1889  delta /= (2+Ne);
1890  }
1891  }
1892  res = Pi_lx_interval();
1893  times2pown(res,-1); // res = pi/2;
1894  if (!fertig)
1895  res -= asin( sqrt(delta*(2-delta)) );
1896  if (neg_x) res = -res;
1897  }
1898  }
1899  } // x != 0
1900 
1901  return res;
1902  } // Asin_beta(...)
1903 
1904 lx_cinterval asin(const lx_cinterval& z) noexcept
1905 {
1906  int stagsave = stagprec,
1907  stagmax = 30;
1908  if (stagprec>stagmax) stagprec = stagmax;
1909 
1910  lx_cinterval res;
1911  lx_interval rez = Re(z), imz = Im(z);
1912 
1913  lx_real
1914  irez = Inf(rez),
1915  srez = Sup(rez),
1916  iimz = Inf(imz),
1917  simz = Sup(imz);
1918 
1919  lx_interval hxl(irez), hxu(srez), hyl(iimz), hyu(simz);
1920 
1921  lx_real resxl, resxu, resyl, resyu;
1922 
1923  bool bl = (sm_zero(iimz)) && (gr_zero(simz)),
1924  raxis = (eq_zero(iimz)) && (eq_zero(simz));
1925 
1926  //
1927  // 1st: check for singularities
1928  //
1929  if ( (irez<-1 && (bl || (sm_zero(iimz) && eq_zero(simz)))) ||
1930  (srez >1 && (bl || (eq_zero(iimz) && gr_zero(simz)))) )
1931  cxscthrow(STD_FKT_OUT_OF_DEF(
1932  "lx_cinterval asin(const lx_cinterval& z); z contains singularities."));
1933 
1934  //
1935  // 2nd: real part
1936  //
1937  if (sm_zero(iimz) && gr_zero(simz))
1938  // z intersects [-1,1]
1939  {
1940  if (se_zero(irez))
1941  resxl = Inf( asin(hxl) );
1942  else
1943  resxl = Inf(Asin_beta(hxl,lx_interval( max(-iimz,simz) )) );
1944  if (sm_zero(srez))
1945  resxu = Sup(Asin_beta(hxu,lx_interval( max(-iimz,simz) )) );
1946  else
1947  resxu = Sup( asin(hxu) );
1948  }
1949  else
1950  {
1951  if ( (ge_zero(iimz) && ge_zero(irez)) ||
1952  (se_zero(simz) && se_zero(irez)) )
1953  // left boundary in quadrants I or III
1954  // min( Re( z ) ) in upper left corner
1955  resxl = Inf( Asin_beta(hxl,hyu) ); // Blomquist, 19.06.2005;
1956  else
1957  // left boundary in quadrants II or IV
1958  // min( Re( z ) ) in lower left corner
1959  resxl = Inf( Asin_beta(hxl,hyl) ); // Blomquist, 19.06.2005;
1960  if ( (ge_zero(iimz) && ge_zero(srez)) || (se_zero(simz) && se_zero(srez) ) )
1961  // right boundary in quadrants I or III
1962  // max( Re( z ) ) in lower right corner
1963  resxu = Sup( Asin_beta(hxu,hyl) ); // Blomquist, 19.06.2005;
1964  else
1965  // right boundary in quadrants II or IV
1966  // max( Re( z ) ) in upper right corner
1967  resxu = Sup( Asin_beta(hxu,hyu) ); // Blomquist, 19.06.2005;
1968  }
1969 
1970  //
1971  // 3rd: imaginary part
1972  //
1973  if (raxis)
1974  { // Interval argument is now a subset of the real axis.
1975  // Blomquist, 16.06.2005;
1976  if (sm_zero(srez)) resyl = Inf( ACOSH_f_aux( hxu, hyu ));
1977  else resyl = -Sup( ACOSH_f_aux( hxu, hyu ));
1978  if (gr_zero(irez)) resyu = -Inf( ACOSH_f_aux( hxl, hyu ));
1979  else resyu = Sup( ACOSH_f_aux( hxl, hyu ));
1980  }
1981  else
1982  if (se_zero(simz))
1983  // z in lower half plane
1984  // min( Im( z ) ) in point with max |z|
1985  // max( Im( z ) ) in point with min |z|
1986  {
1987  if (irez < -srez)
1988  // most of z in quadrant III
1989  {
1990  resyl = -Sup( ACOSH_f_aux( hxl, hyl ) );
1991  if (sm_zero(srez))
1992  resyu = -Inf( ACOSH_f_aux( hxu, hyu ) );
1993  else
1994  resyu = -Inf( ACOSH_f_aux( lx_interval(0),hyu ) );
1995  }
1996  else
1997  // most of z in quadrant IV
1998  {
1999  resyl = -Sup( ACOSH_f_aux( hxu, hyl ) );
2000  if (gr_zero(irez))
2001  resyu = -Inf( ACOSH_f_aux( hxl, hyu ) );
2002  else
2003  resyu = -Inf( ACOSH_f_aux( lx_interval(0),hyu ) );
2004  }
2005  }
2006  else
2007  if (ge_zero(iimz))
2008  // z in upper half plane
2009  // min( Im( z ) ) in point with min |z|
2010  // max( Im( z ) ) in point with max |z|
2011  {
2012  if( irez < -srez ) // if( irez + srez < 0.0 )
2013  // most of z in quadrant II
2014  {
2015  resyu = Sup( ACOSH_f_aux( hxl, hyu ) );
2016  if (sm_zero(srez))
2017  resyl = Inf( ACOSH_f_aux( hxu, hyl ) );
2018  else
2019  resyl = Inf( ACOSH_f_aux( lx_interval(0), hyl ) );
2020  }
2021  else
2022  // most of z in quadrant I
2023  {
2024  resyu = Sup( ACOSH_f_aux( hxu, hyu ) );
2025  if (gr_zero(irez))
2026  resyl = Inf( ACOSH_f_aux( hxl, hyl ) );
2027  else
2028  resyl = Inf( ACOSH_f_aux( lx_interval(0), hyl ) );
2029  }
2030  }
2031  else
2032  // z intersects imaginary axes
2033  // min( Im( z ) ) in point in lower half plane with max |z|
2034  // max( Im( z ) ) in point in upper half plane with max |z|
2035  {
2036  if( irez < -srez ) // if( irez + srez < 0.0 )
2037  // most of z in quadrants II and IV
2038  {
2039  resyl = - Sup( ACOSH_f_aux( hxl, hyl ) );
2040  resyu = Sup( ACOSH_f_aux( hxl, hyu ) );
2041  }
2042  else
2043  {
2044  resyl = - Sup( ACOSH_f_aux( hxu, hyl ) );
2045  resyu = Sup( ACOSH_f_aux( hxu, hyu ) );
2046  }
2047  }
2048 
2049  res = lx_cinterval( lx_interval(resxl,resxu),lx_interval(resyl,resyu) );
2050  stagprec = stagsave;
2051  res = adjust(res);
2052  return res;
2053 }
2054 //
2055 //-- end asin -----------------------------------------------------------------
2056 
2057 //-----------------------------------------------------------------------------
2058 //----------------------------- acos ------------------------------------------
2059 
2060 lx_interval BETA_xy(const lx_interval& x, const lx_interval& y)
2061 // Calculation of beta = x / ( sqrt[(x+1)^2 + y^2] + sqrt[(x-1)^2 + y^2]/2 );
2062 // Only for the internal use for calculating acos(beta);
2063  {
2064  lx_interval res,xa;
2065  l_interval xl,yl;
2066  real exx,exy;
2067  int exxl,exyl;
2068  interval z;
2069  bool neg;
2070 
2071  xa = x;
2072  xl = li_part(xa); yl = li_part(y);
2073  neg = Inf(xl)<0;
2074  if (neg) xa = -xa;
2075  exx = expo(xa); exxl = expo_gr(Inf(xl));
2076  exy = expo(y); exyl = expo_gr(Inf(yl));
2077 
2078  if (exyl<-100000) // y = 0;
2079  if (Inf(xa)<1) res = xa;
2080  else res = 1.0;
2081  else // y != 0;
2082  if (exxl<-100000) // x = 0
2083  res = 0;
2084  else
2085  {
2086  z = interval(exy) + (exyl-exxl-2103);
2087  if (exx<Inf(z)) res = 0;
2088  else
2089  res = xa / f_aux_asin(xa,y);
2090  }
2091  if (Sup(res)>1) res = 1.0;
2092 
2093  if (neg) res = -res;
2094  return res;
2095  } // BETA_xy(...)
2096 
2097 lx_interval Asin_arg(const lx_interval& x, const lx_interval& y )
2098 // Asin_arg calculats for point intervals x,y with
2099 // beta := 2*|x| / ( sqrt((x+1)^2+y^2) + sqrt((x-1)^2+y^2) )
2100 // and delta := 1 - |beta| an inclusion of:
2101 // arcsin[ sqrt(delta)*sqrt(2-delta) ].
2102 // The point interval x may be negative!
2103 // Blomquist, 06.06.2008;
2104  {
2105  const real c2 = -9007199254739968.0;
2106  lx_interval res,xa,F_xa,S1,S2,xm1,xp1,w_delta,T;
2107  lx_real Infxa;
2108  bool neg_x;
2109  xa = x;
2110  neg_x = Inf(x)<0;
2111  if (neg_x) xa = -xa; // Inf(xa) > 0 :
2112  Infxa = Inf(xa);
2113 
2114  if (Infxa > 1.0)
2115  { // x > +1;
2116  if (y == 0.0)
2117  res = 0.0;
2118  else // y != 0;
2119  {
2120  F_xa = xa;
2121  times2pown(F_xa,c2);
2122 
2123  if (Inf(abs(y)) > Inf(F_xa))
2124  // y > x * 2^(-9007199254739968):
2125  // Now the calculation of w_delta according to (1) is
2126  // possible and additional the calculation of the total
2127  // argument of the asin-function:
2128  {
2129  xm1 = xa-1.0;
2130  if (Inf(xm1)<=0)
2131  xm1 = lx_interval(Inf( lx_interval(-2097,l_interval(1)) ),
2132  Sup(xm1));
2133  xp1 = xa+1.0;
2134  S1 = sqrtx2y2(xm1,y);
2135  S2 = sqrtx2y2(xp1,y);
2136  w_delta = sqrt(2 + S1 + S2);
2137  w_delta = w_delta * sqrt(S1+xm1);
2138  w_delta = Sqrt2_lx_interval()*abs(y) / w_delta;
2139  T = (Sqrt2_lx_interval()-w_delta)*(Sqrt2_lx_interval()+w_delta);
2140  T = w_delta * sqrt(T);
2141  res = asin( T );
2142  }
2143  else
2144  {
2145  // y <= x * 2^(-9007199254739968)
2146  if (Inf(xa) >= 2)
2147  {
2148  real exx,exxl,exy,exyl;
2149  interval z;
2150  double dbl;
2151 
2152  exx = expo(xa); exy = expo(y);
2153  exxl = expo_gr(li_part(xa));
2154  exyl = expo_gr(li_part(y));
2155 
2156  z = interval(exy)-interval(exx) + (exyl-exxl+1079);
2157  // 1079 = 1074 + 5; see documentation
2158  exx = Sup(z);
2159  if (exx<-Max_Int_R)
2160  exyl = -Max_Int_R;
2161  else // Sup(z) >= -9007199254740991
2162  {
2163  if (diam(z)==0) exyl = Sup(z);
2164  else
2165  {
2166  dbl = floor(_double(exx));
2167  exyl = real(dbl) + 1; // exyl is integer!
2168  }
2169  }
2170 
2171  res = lx_interval(lx_real(0.0),
2172  lx_real(exyl,l_real(minreal)));
2173  }
2174  else // 1 < Inf(xa) < 2;
2175  {
2176  T = sqrt(xa*(xa-1.0));
2177  if (Inf(T) <= 0)
2178  {
2179  times2pown(xa,-2097);
2180  xa = sqrt(xa);
2181  times2pown(xa,3000);
2182  res = abs(y);
2183  times2pown(res,3001);
2184  res = res / xa;
2185  }
2186  else
2187  {
2188  res = abs(y);
2189  times2pown(res,1);
2190  res = res/T;
2191  }
2192  res = lx_interval(lx_real(0.0),Sup(res));
2193  }
2194  }
2195  }
2196  }
2197  else
2198  if (Infxa == 1.0)
2199  { // x = 1;
2200  T = abs(y);
2201  res = 2 + T + sqrtx2y2(lx_interval(0,l_interval(2)),T);
2202  times2pown(T,1);
2203  T = T / res;
2204  res = asin(sqrt(T*(2-T)));
2205  }
2206  else
2207  { // 0.75 <= x < 1;
2208  if (y == 0.0)
2209  {
2210  T = sqrt1mx2(xa);
2211  res = asin(T);
2212  }
2213  else
2214  {
2215  T = 1 - xa;
2216  S1 = sqrtx2y2(T,y);
2217  S2 = sqrtx2y2(x+1,y);
2218  res = S1 + T;
2219  times2pown(res,1);
2220  T = 2 + S1 + S2;
2221  T = res / T;
2222  T = sqrt( T*(2-T) );
2223  res = asin(T);
2224  }
2225  }
2226  return res;
2227 
2228  } // Asin_arg
2229 
2230 lx_interval Acos_beta(const lx_interval& x, const lx_interval& y)
2231 // Calculating the improved real part of acos(z); Blomquist 05.06.2005;
2232 // Re(acos(z)) = acos[ 2x / (sqrt[(x+1)^2+y^2] + sqrt[(x-1)^2+y^2]) ]
2233  {
2234  const real c1 = 0.75;
2235  lx_interval res(0),beta;
2236  beta = BETA_xy(x,y);
2237  if (Sup(beta)<c1)
2238  if (Sup(beta)<-c1)
2239  { // Improvement for beta --> -1
2240  res = Pi_lx_interval() - Asin_arg(x,y);
2241  }
2242  else res = acos(beta); // Normal calculation
2243  else // Sup(beta)>=c1
2244  res = Asin_arg(x,y);
2245  return res;
2246 } // Acos_beta(...)
2247 
2248 //
2249 lx_cinterval acos(const lx_cinterval& z) noexcept
2250 {
2251  lx_interval
2252  rez = Re(z),
2253  imz = Im(z);
2254 
2255  lx_real
2256  irez = Inf(rez),
2257  srez = Sup(rez),
2258  iimz = Inf(imz),
2259  simz = Sup(imz);
2260 
2261  lx_interval
2262  hxl(irez), hxu(srez), hyl(iimz), hyu(simz);
2263 
2264  bool bl = iimz< 0.0 && simz>0.0,
2265  raxis = iimz==0.0 && simz==0;
2266  lx_real resxl, resxu, resyl, resyu;
2267  //
2268  // 1st: check for singularities
2269  //
2270  if( (irez<-1 && (bl || (iimz<0 && simz==0))) ||
2271  (srez >1 && (bl || (iimz==0 && simz>0))) )
2272  cxscthrow(STD_FKT_OUT_OF_DEF("lx_cinterval acos(const lx_cinterval& z); z contains singularities."));
2273  //
2274  // 2nd: real part
2275  //
2276  // Blomquist, 05.02.2007;
2277 
2278  if( iimz < 0.0 && simz > 0.0 )
2279  // z intersects [-1,1] on the x-axis
2280  {
2281  if( irez <= 0.0 ) resxu = Sup( acos( hxl ) );
2282  else resxu = Sup( Acos_beta(hxl,lx_interval( max(-iimz,simz) )) );
2283 
2284  if( srez < 0.0 )
2285  resxl = Inf( Acos_beta(hxu,lx_interval(max(-iimz,simz))) );
2286  else resxl = Inf( acos( hxu ) );
2287  }
2288  else
2289  {
2290  if (irez<0 && srez>0)
2291  // z intersects the posizive or negative y-axis
2292  if (ge_zero(iimz))
2293  {
2294  resxl = Inf( Acos_beta(hxu,hyl) );
2295  resxu = Sup( Acos_beta(hxl,hyl) );
2296  }
2297  else
2298  {
2299  resxl = Inf( Acos_beta(hxu,hyu) );
2300  resxu = Sup( Acos_beta(hxl,hyu) );
2301  }
2302  else
2303  {
2304  if ( (ge_zero(iimz) && ge_zero(irez)) ||
2305  (se_zero(simz) && sm_zero(irez)) )
2306  // left boundary in quadrants I or III
2307  // min( Re( z ) ) in lower right corner
2308  resxl = Inf( Acos_beta(hxu,hyl) );
2309  else
2310  // left boundary in quadrants II or IV
2311  // min( Re( z ) ) in upper right corner
2312  resxl = Inf( Acos_beta(hxu,hyu) );
2313 
2314  if ( (ge_zero(iimz) && gr_zero(srez)) ||
2315  (se_zero(simz) && se_zero(srez)) )
2316  // right boundary in quadrants I or III
2317  // max( Re( z ) ) in upper left corner
2318  resxu = Sup( Acos_beta(hxl,hyu) );
2319  else
2320  // right boundary in quadrants II or IV
2321  // max( Re( z ) ) in lower left corner
2322  resxu = Sup( Acos_beta(hxl,hyl) );
2323  }
2324  }
2325 
2326  //
2327  // 3rd: imaginary part
2328  //
2329  if (raxis)
2330  { // Interval argument is now a subset of the real axis.
2331  // Blomquist, 16.06.2005;
2332  if (srez<0.0) resyl = Inf( ACOSH_f_aux( hxu, hyu ));
2333  else resyl = -Sup( ACOSH_f_aux( hxu, hyu ));
2334  if (irez>0.0) resyu = -Inf( ACOSH_f_aux( hxl, hyu ));
2335  else resyu = Sup( ACOSH_f_aux( hxl, hyu ));
2336  }
2337  else
2338  if( simz <= 0.0 )
2339  // z in lower half plane
2340  // min( Im( z ) ) in point with max |z|
2341  // max( Im( z ) ) in point with min |z|
2342  {
2343  if (irez < -srez)
2344  // most of z in quadrant III
2345  {
2346  resyl = -Sup( ACOSH_f_aux( hxl, hyl ) );
2347  if( srez < 0.0 )
2348  resyu = -Inf( ACOSH_f_aux( hxu, hyu ) );
2349  else
2350  resyu = -Inf( ACOSH_f_aux( lx_interval(0), hyu ) );
2351  }
2352  else
2353  // most of z in quadrant IV
2354  {
2355  resyl = -Sup( ACOSH_f_aux( hxu, hyl ) );
2356  if( irez > 0.0 )
2357  resyu = -Inf( ACOSH_f_aux( hxl, hyu ) );
2358  else
2359  resyu = -Inf( ACOSH_f_aux( lx_interval(0), hyu ) );
2360  }
2361  }
2362  else
2363  if( ge_zero(iimz) )
2364  // z in upper half plane
2365  // min( Im( z ) ) in point with min |z|
2366  // max( Im( z ) ) in point with max |z|
2367  {
2368  if( irez < -srez ) // if( irez + srez < 0.0 )
2369  // most of z in quadrant II
2370  {
2371  resyu = Sup( ACOSH_f_aux( hxl, hyu ) );
2372  if( sm_zero(srez) )
2373  resyl = Inf( ACOSH_f_aux( hxu, hyl ) );
2374  else
2375  resyl = Inf( ACOSH_f_aux( lx_interval(0), hyl ) );
2376  }
2377  else
2378  // most of z in quadrant I
2379  {
2380  resyu = Sup( ACOSH_f_aux( hxu, hyu ) );
2381  if( irez > 0.0 )
2382  resyl = Inf( ACOSH_f_aux( hxl, hyl ) );
2383  else
2384  resyl = Inf( ACOSH_f_aux( lx_interval(0), hyl ) );
2385  }
2386  }
2387  else
2388  // z intersects imaginary axes
2389  // min( Im( z ) ) in point in lower half plane with max |z|
2390  // max( Im( z ) ) in point in upper half plane with max |z|
2391  {
2392  if( irez < -srez ) // if( irez + srez < 0.0 )
2393  // most of z in quadrants II and IV
2394  {
2395  resyl = -Sup( ACOSH_f_aux( hxl, hyl ) );
2396  resyu = Sup( ACOSH_f_aux( hxl, hyu ) );
2397  }
2398  else
2399  {
2400  resyl = -Sup( ACOSH_f_aux( hxu, hyl ) );
2401  resyu = Sup( ACOSH_f_aux( hxu, hyu ) );
2402  }
2403  }
2404 
2405  return lx_cinterval( lx_interval( resxl, resxu ),-lx_interval( resyl, resyu ) );
2406 
2407 } // acos(...)
2408 
2409 //-- end acos -----------------------------------------------------------------
2410 
2411 //-- asinh --------------------------------------------------------------------
2412 //
2413 lx_cinterval asinh( const lx_cinterval& z ) noexcept
2414 //
2415 // asinh( z ) = i * asin( -i * z )
2416 //
2417 {
2418  lx_cinterval res = asin( lx_cinterval( Im(z), -Re(z) ) );
2419  return lx_cinterval( -Im(res), Re(res) );
2420 }
2421 //
2422 //-- end asinh ----------------------------------------------------------------
2423 
2424 //-- acosh --------------------------------------------------------------------
2425 //
2426 lx_cinterval acosh( const lx_cinterval& z ) noexcept
2427 //
2428 // acosh( z ) = i * acos( z ) = +/- i * ( pi / 2 - asin( z ) )
2429 //
2430 {
2431  lx_interval
2432  rez = Re(z),
2433  imz = Im(z);
2434 
2435  lx_real
2436  irez = Inf(rez),
2437  srez = Sup(rez),
2438  iimz = Inf(imz),
2439  simz = Sup(imz);
2440 
2441  lx_interval hxl(irez), hxu(srez), hyl(iimz), hyu(simz);
2442 
2443  lx_real resxl, resxu, resyl, resyu;
2444 
2445  //
2446  // 1st: check for singularities
2447  //
2448  if( ( se_zero(iimz) && ge_zero(simz) ) && ( irez < 1.0 ) )
2449  cxscthrow(STD_FKT_OUT_OF_DEF(
2450  "lx_cinterval acosh( const lx_cinterval& z ); z contains singularities."));
2451  // With this restriction the complex interval argument and the real axis
2452  // must not have any common point, if irez < +1;
2453  // So for example the negative real axis must not be touched from above if
2454  // irez<1, although this should be possible if the principal branch is
2455  // considered! So the above restriction is too widely in
2456  // some cases; Blomquist, 21.06.2005;
2457  //
2458  // 2nd: z in upper half plane (or on the real axis)
2459  // acosh( z ) = + i * ( pi / 2 - asin( z ) )
2460  //
2461  if( gr_zero(iimz) )
2462  {
2463  lx_cinterval res = acos(z);
2464  return lx_cinterval( -Im(res),Re(res) );
2465  }
2466  //
2467  // 3rd: z in lower half plane
2468  // acosh( z ) = - i * ( pi / 2 - asin( z ) )
2469  //
2470  if( sm_zero(simz) )
2471  {
2472  // cinterval res = HALFPI() - asin( z );
2473  lx_cinterval res = acos(z); // Blomquist, 14.06.2005
2474  return lx_cinterval( Im(res), -Re(res) );
2475  }
2476  //
2477  // z intersects [1,infinity)
2478  //
2479  // real part
2480  // minimum on the left on real axes, maximum in lower or upper right corner
2481  //
2482  resxl = Inf( acosh( hxl ) );
2483  lx_interval ytilde( max( -iimz, simz ) );
2484  // resxu = Sup( acosh( f_aux_asin( hxu, ytilde ) ) );
2485  resxu = Sup( ACOSH_f_aux(hxu,ytilde) ); // Blomquist, 14.06.2005;
2486  //
2487  // imaginary part
2488  // minimum in lower left corner, maximum in upper left corner
2489  //
2490  // resyl = -Sup( acos( hxl / f_aux_asin( hxl, hyl ) ) );
2491  resyl = -Sup( Acos_beta(hxl,hyl) ); // Blomquist, 14.06.2005;
2492  // resyu = Sup( acos( hxl / f_aux_asin( hxl, hyu ) ) );
2493  resyu = Sup( Acos_beta(hxl,hyu) );
2494 
2495  return lx_cinterval(lx_interval( resxl, resxu ),lx_interval( resyl, resyu ));
2496 }
2497 //
2498 //-- end acosh ----------------------------------------------------------------
2499 
2500 
2501 //-- atan ---------------------------------------------------------------------
2502 //
2503 
2504 bool sign_test(const lx_interval& x, int s_org)
2505 // Only for internal use by function re_atan(...)
2506 {
2507  bool bl1,bl2,alter;
2508  if (diam(li_part(x))>0)
2509  {
2510  bl1 = eq_zero(Sup(x)) && (s_org==1);
2511  bl2 = eq_zero(Inf(x)) && (s_org==-1);
2512  alter = bl1 || bl2;
2513  }
2514  else
2515  alter = sign(Sup(li_part(x))) != s_org;
2516  return alter;
2517 } // sign_test(...)
2518 
2519 void re_atan( const lx_interval& y, lx_interval& x, lx_interval& res )
2520 // Calculating an inclusion res of 1 - y^2 - x^2;
2521 // x is always a point interval and y =[y1,y2], with y1 <= y2;
2522 // Blomquist, 10.06.2008;
2523  {
2524  const real c1 = 4503599627369982.0;
2525  lx_interval ya(abs(y)), One(lx_interval(0,l_interval(1))), xa;
2526  lx_real R,U;
2527  real ex_x,ex_xl,ex_y,ex_yl,s;
2528  int signx;
2529  bool alter;
2530 
2531  ex_x = expo(x);
2532  ex_xl = expo_gr(Inf(li_part(x)));
2533  ex_y = expo(y);
2534  ex_yl = expo_gr(Sup(li_part(ya)));
2535  ya = y;
2536  signx = sign(Inf(li_part(x)));
2537 
2538  if (ex_xl < -1000000)
2539  res = 1.0;
2540  else // x <> 0:
2541  if (ex_yl < -1000000) // y == [0,0] and x <> [0,0];
2542  if(ex_x > c1)
2543  {
2544  res = 1/x - x;
2545  x = -1; // Eigentlich sollte hier x = +1 stehen ???
2546  }
2547  else // y == [0,0] and x^2 without overflow
2548  {
2549  res = (1-x)*(1+x);
2550  if (Sup(res)>1) // to avoid an overestimation
2551  res = lx_interval(Inf(res),lx_real(1.0));
2552  }
2553  else // y <> [0,0] and x <> [0,0];
2554  if (ex_x>c1) // |x| --> infty;
2555  if (ex_y>c1) // |x| --> infty and |y| --> infty;
2556  {
2557  s = c1 - max(ex_x,ex_y); // s < 0; and s is an integer value!
2558  times2pown_neg(One,s); // 2^s
2559  times2pown_neg(x,s); // x*2^s
2560  times2pown(ya,s); // y*2^s
2561  res = (One-x)*(One+x) - sqr(ya);
2562  times2pown_neg(x,s); // x*2^2s
2563  }
2564  else // |x| --> infty and |y| not--> infty;
2565  {
2566  res = sqr(y); xa = abs(x);
2567  if (res == 1.0)
2568  res = xa;
2569  else
2570  {
2571  res = res - 1.0; // res = y^2 - 1;
2572  R = Sup(abs(res));
2573  ex_y = expo(R);
2574  ex_yl = expo_gr(lr_part(R));
2575  ex_xl = ex_xl - 1051;
2576  if (ex_y+ex_yl < 2*(ex_x+ex_xl))
2577  res = xa *
2578  lx_interval(Inf(One_m_lx_interval()),
2579  Sup(One_p_lx_interval()));
2580  else
2581  res = res/xa + xa;
2582  }
2583  x = (Sup(x)>0)? -1.0 : 1.0;
2584  }
2585  else // |x| not to infty
2586  if (ex_y>c1) // |y| --> infty;
2587  {
2588  s = c1 - ex_y; //
2589  res = (x-1)*(x+1);
2590  R = Sup( abs(res) );
2591  times2pown_neg(x,s); // x*2^s
2592  times2pown(ya,s); // y*2^s
2593  if (eq_zero(R))
2594  res = sqr(ya);
2595  else
2596  { // R > 0;
2597  ya = sqr(ya); // ya = (y*2^s)^2;
2598  times2pown_neg(res,s);
2599  times2pown_neg(res,s);
2600  // res = {(x-1)*(x+1)}*2^2s;
2601  res = res + ya;
2602  }
2603  times2pown_neg(x,s); // (x*2^s)*2^s = x*2^(2s)
2604  x = -x;
2605  }
2606  else // |x|,|y| not to +infty
2607  if (ex_y<-c1) // |y| ---> 0, y<>[0,0];
2608  if (abs(x)==1)
2609  {
2610  times2pown(x,9007199254738894.0);
2611  x = -x;
2612  times2pown(ya,4503599627369447.0);
2613  res = sqr(ya);
2614  }
2615  else
2616  {
2617  res = (x-1)*(x+1) + sqr(ya);
2618  x = -x;
2619  }
2620  else // |x|,|y| not to +infty and |y| not to 0
2621  if (ex_x < -c1) // now: |x| --> 0
2622  {
2623  res = (y-1)*(y+1);
2624  if (res == 0.0) // alpha = -1/x;
2625  {
2626  x = -1;
2627  res = x;
2628  }
2629  else
2630  {
2631  res += sqr(x);
2632  x = -x;
2633  }
2634  }
2635  else // x^2 and y2 can now be evaluated without any
2636  // integer overflow!
2637  res = (1-x)*(1+x) - sqr(y);
2638  alter = sign_test(x,signx);
2639  if (alter)
2640  {
2641  x = -x;
2642  res = -res;
2643  }
2644 } // re_atan
2645 
2646 void re_vert( const lx_real& x, const lx_interval& hx,
2647  const lx_real& rew_inf, const lx_real& rew_sup,
2648  lx_real& resxl, lx_real& resxu )
2649 //
2650 // Subroutine of analytic inverse tangent function.
2651 // Evaluate real part on a vertical boundary.
2652 //
2653 {
2654  if( eq_zero(x) )
2655  // singularities have been handled before, hence Re( w ) > 0
2656  {
2657  resxl = 0.0;
2658  resxu = 0.0;
2659  }
2660  else
2661  {
2662  lx_interval hx2(hx),Pid2,Pid4;
2663  Pid4 = Pid4_lx_interval();
2664  times2pown(hx2,1);
2665  if( gr_zero(x) )
2666  // w in quadrants I and/or II
2667  // atan is the inverse function of tan(t), t in (-pi/2,pi/2).
2668  {
2669  resxl = gr_zero(rew_sup)? Inf( Atan( hx2,rew_sup )/2.0 )
2670  : ( sm_zero(rew_sup)? Inf( (Atan( hx2,rew_sup ) + Pi_lx_interval() )/2.0 )
2671  : Inf( Pid4 ) );
2672 
2673  resxu = gr_zero(rew_inf)? Sup( Atan( hx2,rew_inf )/2.0 )
2674  : ( sm_zero(rew_inf)? Sup( (Atan( hx2,rew_inf ) + Pi_lx_interval())/2.0 )
2675  : Sup( Pid4 ) );
2676  }
2677  else
2678  // w in quadrants III and/or IV
2679  {
2680  resxl = sm_zero(rew_inf)? Inf( (Atan( hx2,rew_inf ) - Pi_lx_interval())/2.0 )
2681  : ( gr_zero(rew_inf)? Inf( Atan( hx2,rew_inf )/2.0 )
2682  : -Sup( Pid4 ) );
2683  resxu = sm_zero(rew_sup)? Sup( (Atan( hx2,rew_sup ) - Pi_lx_interval())/2.0 )
2684  : ( gr_zero(rew_sup)? Sup( Atan( hx2,rew_sup )/2.0 )
2685  : -Inf( Pid4 ) );
2686  }
2687  }
2688 } // re_vert
2689 
2690 lx_interval T_atan(const lx_real& x)
2691 // Calculating an inclusion of
2692 // ln[ 1+2/(sqrt(1+x^2)-1) ].
2693 // x will be handeld as a point interval [x]=[x,x], with x>0.
2694 // Blomquist, 10.06.2008;
2695 {
2696  const real c1 = 4503599627367859.0;
2697  lx_interval res,
2698  ix(x); // ix is point interval with x>0;
2699  real ex_x(expo(ix));
2700 
2701  if (ex_x<-c1)
2702  {
2703  res = sqrt1px2(ix) + 1;
2704  times2pown(res,1);
2705  res += sqr(ix); // res = 2(1+sqrt(1+x^2)) + x^2;
2706  ix = ln(ix);
2707  times2pown(ix,1); // ix = 2*ln(x);
2708  res = ln(res) - ix;
2709  }
2710  else // -c1 <= ex_x
2711  if (ex_x<1150) // -c1 <= ex_x < 1150
2712  res = lnp1(2/sqrtp1m1(sqr(ix)));
2713  else res = lnp1(2/(sqrt1px2(ix)-1));
2714 
2715  return res;
2716 } // T_atan
2717 
2718 lx_interval Q_atan(const lx_interval& x, const lx_interval& y)
2719 {
2720 // x: abs(Re(z)); So x is a real interval x = [x1,x2], with 0<=x1<x2.
2721 // y: Inf(Im(z)); So y is a point interval, with y >= 0.
2722 // Q_atan returns an inclusion of ln[1 + 4y/(x^2+(1-y)^2)]
2723 // Tested in detail; Blomquist, 13.06.2008;
2724 
2725  const real c1 = 4503599627367859.0,
2726  c2 = 9007199254740990.0,
2727  c3 = 4503599627370495.0;
2728  const lx_real S = lx_real(-c2,MinReal);
2729 
2730  double Dbl;
2731  int r;
2732  lx_real R;
2733  lx_interval res(0.0);
2734  real ex_y(expo(y)), ex_yl(expo_gr(li_part(y))),
2735  ex_x(expo(Sup(x))),ex_xl,exx,exxl,dbl,s,up,exy,exyl;
2736  interval dbli,z;
2737 
2738  ex_xl = expo_gr(lr_part(Sup(x)));
2739  if (ex_xl<-100000) ex_x = 0;
2740  if (ex_yl>-100000) // y = [y,y], y>0;
2741  {
2742  if (y==1.0)
2743  { // y = 1;
2744  if (ex_x < -c1)
2745  {
2746  res = ln(x);
2747  times2pown(res,1); // res = 2*ln(x)
2748  res = ln(4+sqr(x)) - res;
2749  }
2750  else
2751  if (ex_x > c1)
2752  {
2753  R = Inf(x);
2754  exx = expo(R);
2755  exxl = expo_gr(lr_part(R));
2756  if (514-exxl < -c3+exx)
2757  res = lx_interval(lx_real(0),S);
2758  else
2759  if (3 - exxl >= -c3 + exx) // dbl >= -c2
2760  {
2761  dbl = -2*(exx+exxl-3); // without any rounding!
2762  res = lx_interval(lx_real(0),lx_real(dbl,l_real(1)));
2763  }
2764  else // dbl < -c2
2765  {
2766  s = c2-exx; s = (s - exx) -2*exxl + 6;
2767  r = (int) _double(s);
2768  res = lx_interval(lx_real(0),lx_real(-c2,comp(0.5,r+1)));
2769  }
2770  }
2771  else // x^2 now without integer overflow
2772  res = lnp1(sqr(2/x));
2773  }
2774  else // Now: y>0 and y != [1,1] and 0<=x1<=x2;
2775  {
2776  if (ex_x>c1 || ex_y>c1 )
2777  {
2778  R = Inf(x);
2779  exx = expo(R); exxl = expo_gr(lr_part(R));
2780  R = Inf(y-1.0); // R != 0:
2781  exy = expo(R); exyl = expo_gr(lr_part(R));
2782 
2783  if (exxl<-100000) // x1 = 0;
2784  dbli = 2*( interval(exy) + (exyl-2) );
2785  else
2786  {
2787  z = 2*( interval(exx) + (exxl-2) );
2788  dbli = 2*( interval(exy) + (exyl-2) );
2789  if (Sup(z) > Sup(dbli)) dbli = z;
2790  }
2791  dbli = (interval(exy) - dbli) + (exyl + 3);
2792  dbl = Sup(dbli);
2793  // Now it holds: Sup{4y/(x^2+(1-y)^2)} < 2^dbl;
2794  // and this upper bound 2^dbl is guaranteed,
2795  // even if dbl is not an integer value!
2796  up = Inf( interval(-c2) - 1022);
2797  if (dbl<up)
2798  res = lx_interval(lx_real(0),S);
2799  else // dbl >= up
2800  if (dbl>=-c2)
2801  {
2802  if (!(Is_Integer(dbl)))
2803  { // rounding to the next greater integer value.
2804  Dbl = floor(_double(dbl)) + 1; // Dbl: integer value
2805  dbl = (real) Dbl;
2806  }
2807  res = lx_interval(lx_real(0),lx_real(dbl,l_real(1)));
2808  }
2809  else
2810  {
2811  s = Sup(interval(c2) + dbl);
2812  if ( !(Is_Integer(s)) )
2813  { // rounding to the next greater integer value.
2814  Dbl = floor(_double(s)) + 1;
2815  r = (int) Dbl;
2816  }
2817  else r = (int) _double(s);
2818  res = lx_interval(lx_real(0),lx_real(-c2,comp(0.5,r+1)));
2819  }
2820  }
2821  else // Now it holds: y>0 and y!=1 and 0<=x1<=x2;
2822  // and: x^2, (1-y)^2 produce no overflow for
2823  // x2, |1-y| --> +infty;
2824  // furthermore |1-y| produces no integer overflow
2825  // for y --> +1.
2826  {
2827  res = y;
2828  times2pown(res,2); // res = 4*y;
2829  res = res/(sqr(x) + sqr(1-y));
2830  res = lnp1(res); // res: inclusion of Q(x,y);
2831  }
2832  }
2833  }
2834  return res;
2835 } // Q_atan
2836 
2837 lx_cinterval atan( const lx_cinterval& z ) noexcept
2838 {
2839  int stagsave = stagprec,
2840  stagmax = 30;
2841  if (stagprec>stagmax) stagprec = stagmax;
2842 
2843  lx_cinterval res;
2844  lx_interval
2845  rez = Re(z),
2846  imz = Im(z);
2847 
2848  lx_real
2849  irez = Inf(rez),
2850  srez = Sup(rez),
2851  iimz = Inf(imz),
2852  simz = Sup(imz);
2853 
2854  lx_interval // all theses intervals are point intervals!
2855  hxl(irez), hxu(srez), hyl(iimz), hyu(simz);
2856 
2857  lx_real
2858  resxl, resxu, resyl, resyu;
2859  //
2860  // 1st: check for singularities
2861  //
2862  if( (se_zero(irez) && ge_zero(srez)) && ( iimz <= -1.0 || simz >= 1.0 ) )
2863  cxscthrow(STD_FKT_OUT_OF_DEF("lx_cinterval atan( const lx_cinterval& z ); points of the branch cuts are not allowed in z."));
2864  //
2865  // 2nd: real part
2866  // Re( atan( z ) ) = Arg( w ) / 2, where w = 1 - x^2 - y^2 + i * 2x )
2867  //
2868  // evaluate atan on vertical boundaries
2869  //
2870  lx_interval
2871  // y_sqr = sqr( imz ),
2872  // rew_l = (1 - y_sqr) - sqr( hxl ),
2873  // Blomquist; before: rew_l = 1 - sqr(hxl) - y_sqr,
2874  // rew_u = (1 - y_sqr) - sqr( hxu );
2875  // Blomquist; before: rew_u = 1 - sqr(hxu) - y_sqr;
2876  rew_l, rew_u;
2877 // ------------------------------ Blomquist ---------------------------------
2878 // ---------- Improvements for Im(z) = [1,1] or Im(z) = [-1,-1]
2879  bool sqrImz_1 = (iimz==simz) && (iimz==1.0 || iimz==-1.0);
2880  // Test for Im(z) = [1,1] or [-1,-1]
2881 
2882  if (sqrImz_1)
2883  {
2884  rew_l = -abs(hxl); hxl = lx_interval(0,l_interval(sign(irez)));
2885  rew_u = -abs(hxu); hxu = lx_interval(0,l_interval(sign(srez)));
2886  }
2887  else
2888  {
2889  // rew_l = (1 - sqr( imz )) - sqr( hxl );
2890  re_atan(imz,hxl,rew_l);
2891  // rew_u = (1 - sqr( imz )) - sqr( hxu );
2892  re_atan(imz,hxu,rew_u);
2893  }
2894 // ------------------------------ Blomquist; 22.02.05; --------------------
2895 
2896  //
2897  // left boundary
2898  //
2899  lx_real rew_inf = Inf( rew_l );
2900  lx_real rew_sup = Sup( rew_l );
2901  re_vert( irez, hxl, rew_inf, rew_sup, resxl, resxu );
2902 
2903  //
2904  // right boundary
2905  //
2906  rew_inf = Inf( rew_u );
2907  rew_sup = Sup( rew_u );
2908  lx_real res_l, res_u;
2909  re_vert( srez, hxu, rew_inf, rew_sup, res_l, res_u );
2910 
2911  if (res_l<resxl) resxl = res_l;
2912  if (res_u>resxu) resxu = res_u;
2913 
2914  //
2915  // look for extremal values on horizontal boundaries
2916  // since atan( x+iy ) = atan( x-iy ),
2917  // intersections can be considered in the upper half plane
2918  //
2919  lx_real abs_y_min = Inf( abs( imz ) );
2920  if( abs_y_min > 1.0 )
2921  {
2922  lx_interval
2923  abs_hyl = lx_interval( abs_y_min ),
2924  // abs_hxl = sqrt( sqr( abs_hyl ) - 1.0 );
2925  abs_hxl = sqrtx2m1(abs_hyl); // Blomquist;
2926 
2927  if( Sup( abs_hxl ) > irez && Inf( abs_hxl ) < srez )
2928  // extremal curve intersects lower boundary of x+i|y| in quadrant I
2929  // intersection in Q I or Q IV: update minimum
2930  // resxl = inf( atan( abs_y_min / abs_hxl ) ) / 2.0;
2931  resxl = Inf( (Pi_lx_interval() - atan( 1.0 / abs_hxl ))/2.0 );
2932  else
2933  if ( -Inf( abs_hxl ) > irez && -Sup( abs_hxl ) < srez )
2934  // extremal curve intersects lower boundary of x+i|y| in quadrant II
2935  // intersection in Q II or Q III: update maximum
2936  resxu = Sup( (atan( 1.0 / abs_hxl ) - Pi_lx_interval())/2.0 );
2937  }
2938  // 3rd: imaginary part
2939  // Im( atan( z ) ) = +/- Ln( 1 +/- 4y/( x^2 + (1 -/+ y)^2 ) ) / 4
2940  //
2941  // evaluate atan on horizontal boundaries
2942  lx_interval
2943  abs_rez = abs(rez),
2944  im_atan_l, im_atan_u;
2945 
2946  if ( sm_zero(iimz) )
2947 // im_atan_l = -ln( 1 - 4 * hyl / ( x_sqr + sqr( 1 + hyl ) ) ) / 4.0;
2948 // im_atan_l = -lnp1(-4 * hyl / ( x_sqr + sqr( 1 + hyl ) )) / 4.0; // Blomquist
2949  im_atan_l = -Q_atan(abs_rez,-hyl); // Blomquist
2950  else
2951 // im_atan_l = ln( 1 + 4 * hyl / ( x_sqr + sqr( 1 - hyl ) ) ) / 4.0;
2952  im_atan_l = Q_atan(abs_rez,hyl); // Blomquist
2953  times2pown(im_atan_l,-2); // Division by 4
2954  if ( sm_zero(simz) )
2955 // im_atan_u = -ln( 1 - 4 * hyu / ( x_sqr + sqr( 1 + hyu ) ) ) / 4.0;
2956 // im_atan_u = -lnp1(-4 * hyu / ( x_sqr + sqr( 1 + hyu ) ) ) / 4.0; // Blomquist
2957  im_atan_u = -Q_atan(abs_rez,-hyu); // Blomquist
2958  else
2959 // im_atan_u = ln( 1 + 4 * hyu / ( x_sqr + sqr( 1 - hyu ) ) ) / 4.0;
2960  im_atan_u = Q_atan(abs_rez,hyu); // Blomquist
2961  times2pown(im_atan_u,-2); // Division by 4
2962  resyl = min( Inf( im_atan_l ), Inf( im_atan_u ) );
2963  resyu = max( Sup( im_atan_l ), Sup( im_atan_u ) );
2964  //
2965  // look for extremal values on vertical boundaries,
2966  // if vertical boundaries intersect extremal curves
2967  //
2968  lx_real abs_x_min = Inf( abs( rez ) );
2969  lx_interval
2970  x_extr = lx_interval( abs_x_min ),
2971 // y_extr = sqrt( 1.0 + sqr( x_extr ) );
2972  y_extr = sqrt1px2(x_extr); // Blomquist;
2973 
2974  if( Inf( y_extr ) < simz && Sup( y_extr ) > iimz )
2975  // extremal curve intersects left boundary of |x|+iy in quadrant I
2976  // update maximum
2977  // resyu = Sup( ln( 1 + 4 * y_extr / ( sqr( x_extr ) + sqr( 1 - y_extr ) ) ) ) / 4.0;
2978  // resyu = Sup( T_atan(abs_x_min)/4.0 ); // Blomquist
2979  {
2980  rez = T_atan(abs_x_min);
2981  times2pown(rez,-2);
2982  resyu = Sup(rez);
2983  }
2984 
2985  if( -Sup(y_extr) < simz && -Inf(y_extr) > iimz )
2986  // extremal curve intersects left boundary of |x|+iy in quadrant IV
2987  // update minimum
2988  // resyl = -Sup( ln( 1 + 4 * y_extr / ( sqr( x_extr ) + sqr( 1 - y_extr ) ) ) ) / 4.0;
2989  // resyl = -Sup( T_atan(abs_x_min)/4.0 ); // Blomquist
2990  {
2991  rez = T_atan(abs_x_min);
2992  times2pown(rez,-2);
2993  resyl = -Sup(rez);
2994  }
2995 
2996  res = lx_cinterval( lx_interval(resxl,resxu),lx_interval(resyl,resyu) );
2997  stagprec = stagsave;
2998  res = adjust(res);
2999  return res;
3000 }
3001 //
3002 //-- end atan -----------------------------------------------------------------
3003 
3004 lx_cinterval acot( const lx_cinterval& z ) noexcept
3005 {
3006  int stagsave = stagprec,
3007  stagmax = 30;
3008  if (stagprec>stagmax) stagprec = stagmax;
3009 
3010  lx_cinterval res;
3011  lx_interval
3012  rez = Re(z),
3013  imz = Im(z);
3014 
3015  lx_real
3016  irez = Inf(rez),
3017  srez = Sup(rez),
3018  iimz = Inf(imz),
3019  simz = Sup(imz);
3020 
3021  lx_interval // all theses intervals are point intervals!
3022  hxl(irez), hxu(srez), hyl(iimz), hyu(simz);
3023 
3024  lx_real
3025  resxl, resxu, resyl, resyu;
3026  //
3027  // 1st: check for singularities
3028  //
3029  if( (se_zero(irez) && ge_zero(srez)) && ( iimz <= 1.0 || simz >= -1.0 ) )
3030  cxscthrow(STD_FKT_OUT_OF_DEF("lx_cinterval acot( const lx_cinterval& z ); points of the branch cuts are not allowed in z."));
3031  //
3032  // 2nd: real part
3033  // Re( atan( z ) ) = Arg( w ) / 2, where w = 1 - x^2 - y^2 + i * 2x )
3034  //
3035  // evaluate atan on vertical boundaries
3036  //
3037  lx_interval
3038 // y_sqr = sqr( imz ),
3039 // rew_l = (1 - y_sqr) - sqr( hxl ), // Blomquist; before: rew_l = 1 - sqr(hxl) - y_sqr,
3040 // rew_u = (1 - y_sqr) - sqr( hxu ); // Blomquist; before: rew_u = 1 - sqr(hxu) - y_sqr;
3041  rew_l, rew_u;
3042 
3043 // ------------------------------ Blomquist ---------------------------------
3044 // ---------- Improvements for Im(z) = [1,1] or Im(z) = [-1,-1]
3045  bool sqrImz_1 = (iimz==simz) && (iimz==1.0 || iimz==-1.0);
3046  // Test for Im(z) = [1,1] or [-1,-1]
3047 
3048  if (sqrImz_1)
3049  {
3050  rew_l = abs(hxl); hxl = lx_interval(0,l_interval(sign(irez)));
3051  rew_u = abs(hxu); hxu = lx_interval(0,l_interval(sign(srez)));
3052  }
3053  else
3054  {
3055  // rew_l = (1 - sqr( imz )) - sqr( hxl );
3056  re_atan(imz,hxl,rew_l);
3057  rew_l = -rew_l;
3058  // rew_u = (1 - sqr( imz )) - sqr( hxu );
3059  re_atan(imz,hxu,rew_u);
3060  rew_u = -rew_u;
3061  }
3062 // ------------------------------ Blomquist; 22.02.05; --------------------
3063 
3064  //
3065  // left boundary
3066  //
3067  lx_real rew_inf = Inf( rew_l );
3068  lx_real rew_sup = Sup( rew_l );
3069  re_vert( irez, hxl, rew_inf, rew_sup, resxl, resxu );
3070 
3071  //
3072  // right boundary
3073  //
3074  rew_inf = Inf( rew_u );
3075  rew_sup = Sup( rew_u );
3076  lx_real res_l, res_u;
3077  re_vert( srez, hxu, rew_inf, rew_sup, res_l, res_u );
3078 
3079  if (res_l<resxl) resxl = res_l;
3080  if (res_u>resxu) resxu = res_u;
3081 
3082  //
3083  // look for extremal values on horizontal boundaries
3084  // since atan( x+iy ) = atan( x-iy ),
3085  // intersections can be considered in the upper half plane
3086  //
3087  lx_real abs_y_min = Inf( abs( imz ) );
3088 
3089  if( abs_y_min > 1.0 )
3090  {
3091  lx_interval
3092  abs_hyl = lx_interval( abs_y_min ),
3093 // abs_hxl = sqrt( sqr( abs_hyl ) - 1.0 );
3094  abs_hxl = sqrtx2m1(abs_hyl); // Blomquist;
3095 
3096  if( Sup( abs_hxl ) > irez && Inf( abs_hxl ) < srez )
3097  // extremal curve intersects lower boundary of x+i|y| in quadrant I
3098  // intersection in Q I or Q IV: update maximum
3099  // resxl = inf( atan( abs_y_min / abs_hxl ) ) / 2.0;
3100  resxu = Sup( atan( 1.0 / abs_hxl ) / 2.0 );
3101  if ( -Inf( abs_hxl ) > irez && -Sup( abs_hxl ) < srez )
3102  // extremal curve intersects lower boundary of x+i|y| in quadrant II
3103  // intersection in Q II or Q III: update minimum
3104  resxl = -Sup( atan( 1.0 / abs_hxl ) /2.0 );
3105  }
3106 
3107  // 3rd: imaginary part
3108  // Im( atan( z ) ) = +/- Ln( 1 +/- 4y/( x^2 + (1 -/+ y)^2 ) ) / 4
3109  //
3110  // evaluate atan on horizontal boundaries
3111  lx_interval
3112  abs_rez = abs(rez),
3113  im_atan_l, im_atan_u;
3114 
3115  if ( sm_zero(iimz) )
3116 // im_atan_l = -ln( 1 - 4 * hyl / ( x_sqr + sqr( 1 + hyl ) ) ) / 4.0;
3117 // im_atan_l = -lnp1(-4 * hyl / ( x_sqr + sqr( 1 + hyl ) )) / 4.0; // Blomquist
3118  im_atan_l = -Q_atan(abs_rez,-hyl); // Blomquist
3119  else
3120 // im_atan_l = ln( 1 + 4 * hyl / ( x_sqr + sqr( 1 - hyl ) ) ) / 4.0;
3121  im_atan_l = Q_atan(abs_rez,hyl); // Blomquist
3122  times2pown(im_atan_l,-2); // Division by 4
3123 
3124  if ( sm_zero(simz) )
3125 // im_atan_u = -ln( 1 - 4 * hyu / ( x_sqr + sqr( 1 + hyu ) ) ) / 4.0;
3126 // im_atan_u = -lnp1(-4 * hyu / ( x_sqr + sqr( 1 + hyu ) ) ) / 4.0; // Blomquist
3127  im_atan_u = -Q_atan(abs_rez,-hyu); // Blomquist
3128  else
3129 // im_atan_u = ln( 1 + 4 * hyu / ( x_sqr + sqr( 1 - hyu ) ) ) / 4.0;
3130  im_atan_u = Q_atan(abs_rez,hyu); // Blomquist
3131  times2pown(im_atan_u,-2); // Division by 4
3132 
3133  resyl = min( Inf( im_atan_l ), Inf( im_atan_u ) );
3134  resyu = max( Sup( im_atan_l ), Sup( im_atan_u ) );
3135  //
3136  // look for extremal values on vertical boundaries,
3137  // if vertical boundaries intersect extremal curves
3138  //
3139  lx_real abs_x_min = Inf( abs( rez ) );
3140  lx_interval
3141  x_extr = lx_interval( abs_x_min ),
3142 // y_extr = sqrt( 1.0 + sqr( x_extr ) );
3143  y_extr = sqrt1px2(x_extr); // Blomquist;
3144 
3145  if( Inf( y_extr ) < simz && Sup( y_extr ) > iimz )
3146  // extremal curve intersects left boundary of |x|+iy in quadrant I
3147  // update maximum
3148  // resyu = Sup( ln( 1 + 4 * y_extr / ( sqr( x_extr ) + sqr( 1 - y_extr ) ) ) ) / 4.0;
3149  // resyu = Sup( T_atan(abs_x_min)/4.0 ); // Blomquist
3150  {
3151  rez = T_atan(abs_x_min);
3152  times2pown(rez,-2);
3153  resyu = Sup(rez);
3154  }
3155 
3156  if( -Sup(y_extr) < simz && -Inf(y_extr) > iimz )
3157  // extremal curve intersects left boundary of |x|+iy in quadrant IV
3158  // update minimum
3159  // resyl = -Sup( ln( 1 + 4 * y_extr / ( sqr( x_extr ) + sqr( 1 - y_extr ) ) ) ) / 4.0;
3160  // resyl = -Sup( T_atan(abs_x_min)/4.0 ); // Blomquist
3161  {
3162  rez = T_atan(abs_x_min);
3163  times2pown(rez,-2);
3164  resyl = -Sup(rez);
3165  }
3166 
3167  res = lx_cinterval( lx_interval(resxl,resxu),lx_interval(-resyu,-resyl) );
3168  stagprec = stagsave;
3169  res = adjust(res);
3170  return res;
3171 }
3172 //
3173 //-- end acot -----------------------------------------------------------------
3174 
3175 //-- atanh --------------------------------------------------------------------
3176 //
3177 lx_cinterval atanh( const lx_cinterval& z ) noexcept
3178 //
3179 // atanh( z ) = - i * atan( i * z )
3180 //
3181 {
3182  lx_cinterval res = atan( lx_cinterval( -Im(z), Re(z) ) );
3183  return lx_cinterval( Im(res), -Re(res) );
3184 }
3185 //
3186 //-- end atanh ----------------------------------------------------------------
3187 
3188 //-- acoth --------------------------------------------------------------------
3189 //
3190 lx_cinterval acoth( const lx_cinterval& z ) noexcept
3191 //
3192 // acoth( z ) = i * acot( i * z )
3193 //
3194 {
3195  lx_cinterval res = acot( lx_cinterval( -Im(z), Re(z) ) );
3196  return lx_cinterval( -Im(res), Re(res) );
3197 }
3198 //
3199 //-- end acoth ----------------------------------------------------------------
3200 
3201 // ---- sqrt1px2 --------------------------------------------------------------
3202 //
3203 lx_cinterval sqrt1px2(const lx_cinterval& z) noexcept
3204 // sqrt(1+z^2);
3205 // Blomquist, 03.07.2008;
3206 {
3207  const lx_real c = lx_real(1600,l_real(1));
3208  int stagsave = stagprec,
3209  stagmax = 30; // l_imath.cpp: sqrt() uses stagmax = 30;
3210  if (stagprec>stagmax) stagprec = stagmax;
3211 
3212  lx_cinterval res;
3213  lx_interval absz(abs(z));
3214  lx_real Inf_absz(Inf(absz));
3215 
3216  if (Inf_absz > c)
3217  {
3218  absz = 1 / lx_interval(Inf_absz);
3219  Inf_absz = Sup(absz);
3220  res = lx_cinterval(lx_interval(-Inf_absz,Inf_absz),
3221  lx_interval(-Inf_absz,Inf_absz));
3222  // res is the correcture interval, i.e.
3223  // z + res or -z + res is the inclusionof sqrt{1+z^2}
3224  res = (Inf(Re(z))>=0)? z + res : -z + res;
3225  }
3226  else
3227  {
3228  res = lx_cinterval(lx_interval(0,l_interval(0)),
3229  lx_interval(0,l_interval(1))); // res = i
3230  if ( Sup(abs(z-res))<0.5 || Sup(abs(z+res))<0.5 )
3231  {
3232  res = lx_cinterval(-Im(z),Re(z)); // Res = i*z;
3233  // (1 - i*z)*(1 + i*z) = 1+z^2;
3234  res = sqrt( (1-res)*(1+res) );
3235  }
3236  else
3237  res = sqrt(1+sqr(z));
3238  }
3239  if (sm_zero(Inf(Re(res))))
3240  res = lx_cinterval(lx_interval(lx_real(0),Sup(Re(res))),Im(res));
3241  stagprec = stagsave;
3242  res = adjust(res);
3243  return res;
3244 }
3245 // -- end sqrt1px2 ------------------------------------------------------------
3246 
3247 lx_cinterval sqrt1mx2(const lx_cinterval& z) noexcept
3248 // sqrt(1-z^2);
3249 // Blomquist, 03.07.2008;
3250 {
3251  const lx_real c = lx_real(1600,l_real(1));
3252  int stagsave = stagprec,
3253  stagmax = 30; // l_imath.cpp: sqrt() uses stagmax = 30;
3254  if (stagprec>stagmax) stagprec = stagmax;
3255 
3256  lx_cinterval res,u;
3257  lx_interval absz(abs(z));
3258  lx_real Inf_absz(Inf(absz));
3259 
3260  if (Inf_absz > c)
3261  {
3262  absz = 1 / lx_interval(Inf_absz);
3263  Inf_absz = Sup(absz);
3264  res = lx_cinterval(lx_interval(-Inf_absz,Inf_absz),
3265  lx_interval(-Inf_absz,Inf_absz)); // res = Delta
3266  u = lx_cinterval(-Im(z),Re(z)); // u = i*z;
3267  // res is the correcture interval, i.e.
3268  // i*z + res or -i*z + res is the inclusion of sqrt{1-z^2}
3269  res = (Inf(Im(z))>=0)? -u + res : u + res;
3270  }
3271  else
3272  {
3273  res = 1-z; u = 1+z;
3274  res = (Sup(abs(res))<0.5 || Sup(abs(u))<0.5)? sqrt(res*u) : sqrt(1-sqr(z));
3275  }
3276  if (sm_zero(Inf(Re(res))))
3277  res = lx_cinterval(lx_interval(lx_real(0),Sup(Re(res))),Im(res));
3278  stagprec = stagsave;
3279  res = adjust(res);
3280  return res;
3281 }
3282 
3283 lx_cinterval sqrtx2m1(const lx_cinterval& z) noexcept
3284 // sqrt(z^2-1);
3285 // Blomquist, 03.07.2008;
3286 {
3287  const lx_real c = lx_real(1600,l_real(1));
3288  int stagsave = stagprec,
3289  stagmax = 30; // l_imath.cpp: sqrt() uses stagmax = 30;
3290  if (stagprec>stagmax) stagprec = stagmax;
3291 
3292  lx_cinterval res,u;
3293  lx_interval absz(abs(z));
3294  lx_real Inf_absz(Inf(absz));
3295 
3296  if (Inf_absz > c)
3297  {
3298  absz = 1 / lx_interval(Inf_absz);
3299  Inf_absz = Sup(absz);
3300  res = lx_cinterval(lx_interval(-Inf_absz,Inf_absz),
3301  lx_interval(-Inf_absz,Inf_absz)); // res = Delta
3302  // res is the correcture interval, i.e.
3303  res = (Inf(Re(z))>=0)? z + res : -z + res;
3304  }
3305  else
3306  {
3307  res = z-1; u = z+1;
3308  res = (Sup(abs(res))<0.5 || Sup(abs(u))<0.5)? sqrt(res*u) : sqrt(sqr(z)-1);
3309  }
3310 
3311  if (sm_zero(Inf(Re(res))))
3312  res = lx_cinterval(lx_interval(lx_real(0),Sup(Re(res))),Im(res));
3313 
3314  stagprec = stagsave;
3315  res = adjust(res);
3316  return res;
3317 }
3318 
3319 lx_cinterval sqrtp1m1(const lx_cinterval& z) noexcept
3320 // sqrt(1+z)-1;
3321 // Blomquist, 08.07.2008;
3322 {
3323  const real c = 0.125;
3324  int stagsave = stagprec,
3325  stagmax = 30; // l_imath.cpp: sqrt() uses stagmax = 30;
3326  if (stagprec>stagmax) stagprec = stagmax;
3327 
3328  lx_cinterval res;
3329  lx_interval absz(abs(z));
3330  lx_real Sup_absz(Sup(absz));
3331 
3332  if (Sup_absz < c)
3333  res = z / (sqrt(z+1) + 1);
3334  else
3335  res = sqrt(z+1) - 1;
3336 
3337  stagprec = stagsave;
3338  res = adjust(res);
3339  return res;
3340 }
3341 
3342 lx_cinterval expm1(const lx_cinterval& z) noexcept
3343 // exp(z) - 1;
3344 // Blomquist, 09.08.2008;
3345 {
3346  int stagsave = stagprec,
3347  stagmax = 30; // l_imath.cpp: sqrt() uses stagmax = 30;
3348  if (stagprec>stagmax) stagprec = stagmax;
3349 
3350  const lx_interval cancl_test = lx_interval(0,l_interval(0.995,1.005));
3351  lx_interval rez(Re(z)), imz(Im(z));
3352  lx_interval exp_x, sin_y, cos_y, h, Xt;
3353  lx_cinterval res;
3354 
3355  exp_x = exp(rez);
3356  sin_y = sin(imz);
3357  cos_y = cos(imz);
3358 
3359  h = exp_x*cos_y;
3360  if (h < cancl_test && cos_y < cancl_test)
3361  {
3362  h = lnp1(-sqr(sin_y));
3363  times2pown(h,-1);
3364  // h = 0.5 * ln(1-sqr( sin(y) ));
3365  h = expm1(rez+h); // Cancellation also possible here!
3366  }
3367  else h = h - 1; // Cancellation possible here (real part)!
3368  // h: Real part;
3369  imz = exp_x * sin_y;
3370  res = lx_cinterval(h,imz);
3371 
3372  stagprec = stagsave;
3373  res = adjust(res);
3374 
3375  return res;
3376 }
3377 
3378 lx_cinterval lnp1(const lx_cinterval& z) noexcept
3379 {// ln(1+z);
3380  // Blomquist, 11.08.2008;
3381  int stagsave = stagprec,
3382  stagmax = 30;
3383  if (stagprec > stagmax) stagprec = stagmax;
3384  const real c1 = 1e-128;
3385  lx_cinterval y;
3386  lx_interval abs_z(abs(z));
3387  lx_real
3388  srez = Sup( Re(z) ),
3389  simz = Sup( Im(z) ),
3390  iimz = Inf( Im(z) );
3391 
3392  if (-1 <= z) // z contains -1
3393  cxscthrow(STD_FKT_OUT_OF_DEF(
3394  "lx_cinterval lnp1(const lx_cinterval& z); z contains -1"));
3395  if ( srez<-1 && iimz<0 && simz>=0 )
3396  cxscthrow(STD_FKT_OUT_OF_DEF(
3397  "lx_cinterval lnp1(const lx_cinterval& z); z not allowed"));
3398  if (Sup(abs_z) < c1)
3399  {
3400  abs_z = Re(z);
3401  abs_z = lnp1( abs_z*(2+abs_z) + sqr(Im(z)) );
3402  times2pown(abs_z,-1);
3403  y = lx_cinterval( abs_z, arg(1+z) );
3404  }
3405  else
3406  y = Ln(1+z);
3407  stagprec = stagsave;
3408  y = adjust(y);
3409 
3410  return y;
3411 }
3412 
3413 //-- pow_all ------
3414 //
3415 // Non-analytic power function for real l_interval exponent.
3416 //
3417 // If 0 \not\in z, then compute four rectangular intervals that comprehend
3418 // an annulus which contains all values zeta^phi, zeta in z, phi in p.
3419 //
3420 // If 0 in z and negative reals in p, then abort execution
3421 // (potential modification: return the entire complex plane).
3422 //
3423 std::list<lx_cinterval> pow_all( const lx_cinterval& z,
3424  const lx_interval& p ) noexcept
3425 {
3426  lx_interval abs_z = abs(z);
3427 
3428  if( 0.0 < Inf( abs_z ) )
3429  {
3430  lx_interval abs_z_p = exp( p * ln( abs_z ) );
3431 
3432  // Inner and outer radii of the annulus are inf/sup( abs_z_n )
3433  // Inscribed square has side length sqrt( 2 ) * rad_1
3434  lx_interval rad_1 = Sqrt2r_l_interval() * lx_interval(Inf( abs_z_p ));
3435  lx_interval rad_2 = lx_interval(Sup( abs_z_p ));
3436 
3437  std::list<lx_cinterval> res;
3438 
3439  // 4 intervals covering the annulus, counter-clockwise
3440  res.push_back(lx_cinterval(lx_interval( Inf( rad_1 ), Sup( rad_2 ) ),
3441  lx_interval( -Sup( rad_1 ), Sup( rad_2 ) )));
3442  res.push_back(lx_cinterval(lx_interval( -Sup( rad_2 ), Sup( rad_1 ) ),
3443  lx_interval( Inf( rad_1 ), Sup( rad_2 ) ) ));
3444  res.push_back(lx_cinterval(lx_interval( -Sup( rad_2 ), -Inf( rad_1 ) ),
3445  lx_interval( -Sup( rad_2 ), Sup(rad_1) ) ) );
3446  res.push_back(lx_cinterval(lx_interval( -Sup( rad_1 ), Sup( rad_2 ) ),
3447  lx_interval( -Sup( rad_2 ), -Inf(rad_1) ) ) );
3448 
3449  return res;
3450  }
3451  else
3452  // 0 in z
3453  {
3454  if( Inf( p ) > 0.0 )
3455  //
3456  // All values zeta^phi, zeta in z, phi in p lie in a disk
3457  //
3458  {
3459  lx_interval abs_z_p = exp( p * ln( lx_interval( Sup( abs_z ) ) ) );
3460  lx_real rad_p = Sup( abs_z_p );
3461 
3462  std::list<lx_cinterval> res;
3463 
3464  res.push_back( lx_cinterval( lx_interval( -rad_p, rad_p ),
3465  lx_interval( -rad_p, rad_p ) ) );
3466  return res;
3467  }
3468  else
3469  {
3470  //
3471  // The set zeta^phi, zeta in z, phi in p is unbounded
3472  // if inf( p ) < 0. 0^p is undefined for p <= 0.
3473  //
3474  cxscthrow(STD_FKT_OUT_OF_DEF(
3475  "pow_all(lx_cinterval, lx_interval); 0^p is undefined for p <= 0."));
3476  std::list<lx_cinterval> res;
3477  return res;
3478  }
3479  }
3480 }
3481 //
3482 //-- end pow_all --------------------------------------------------------------
3483 
3484 } // namespace cxsc
cxsc::Sqrt2_lx_interval
lx_interval Sqrt2_lx_interval() noexcept
Enclosure-Interval for .
Definition: lx_interval.cpp:4720
cxsc::power
cinterval power(const cinterval &z, int n) noexcept
Calculates .
Definition: cimath.cpp:1941
cxsc::mid
cvector mid(const cimatrix_subv &mv) noexcept
Returns the middle of the matrix.
Definition: cimatrix.inl:739
cxsc::MinReal
const real MinReal
Smallest normalized representable floating-point number.
Definition: real.cpp:62
cxsc::Ln
cinterval Ln(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:829
cxsc::lnp1
cinterval lnp1(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:867
cxsc::asinh
cinterval asinh(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:2718
cxsc::sqrt1px2
cinterval sqrt1px2(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:1071
cxsc::One_p_lx_interval
lx_interval One_p_lx_interval() noexcept
Enclosure-Interval for .
Definition: lx_interval.cpp:6297
cxsc::coth
cinterval coth(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:578
cxsc::sin
cinterval sin(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:215
cxsc::sqrt_all
std::list< cinterval > sqrt_all(const cinterval &z)
Calculates and returns all possible solutions.
Definition: cimath.cpp:1176
cxsc::minreal
const real minreal
Smallest positive denormalized representable floating-point number.
Definition: real.cpp:63
cxsc::cot
cinterval cot(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:538
cxsc::exp10
cinterval exp10(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:172
cxsc::Is_Integer
bool Is_Integer(const real &x)
Returns 1 if x is an integer value and if .
Definition: lx_real.inl:63
cxsc::log10
cinterval log10(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:903
cxsc::arg
interval arg(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:741
cxsc::tan
cinterval tan(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:393
cxsc::interval
The Scalar Type interval.
Definition: interval.hpp:55
cxsc::Ln10_lx_interval
lx_interval Ln10_lx_interval() noexcept
Enclosure-Interval for .
Definition: lx_interval.cpp:4279
cxsc::Pi_lx_interval
lx_interval Pi_lx_interval() noexcept
Enclosure-Interval for .
Definition: lx_interval.cpp:4058
cxsc::expo_gr
int expo_gr(const l_interval &x)
Definition: l_interval.inl:522
cxsc::tanh
cinterval tanh(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:565
cxsc::l_cinterval::operator=
l_cinterval & operator=(const real &) noexcept
Implementation of standard assigning operator.
Definition: l_cinterval.inl:82
cxsc::sqrt
cinterval sqrt(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:1007
cxsc::ln
cinterval ln(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:851
cxsc::log2
cinterval log2(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:898
cxsc::abs
ivector abs(const cimatrix_subv &mv) noexcept
Returns the absolute value of the matrix.
Definition: cimatrix.inl:737
cxsc::Pid4_lx_interval
lx_interval Pid4_lx_interval() noexcept
Enclosure-Interval for .
Definition: lx_interval.cpp:6446
cxsc::sqrtx2m1
cinterval sqrtx2m1(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:1109
cxsc::Sqrt2r_l_interval
l_interval Sqrt2r_l_interval() noexcept
Enclosure-Interval for .
Definition: l_imath.cpp:2233
cxsc::Arg
interval Arg(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:654
cxsc::exp
cinterval exp(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:159
cxsc::acos
cinterval acos(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:2553
cxsc::One_m_lx_interval
lx_interval One_m_lx_interval() noexcept
Enclosure-Interval for .
Definition: lx_interval.cpp:6260
cxsc::l_interval
The Multiple-Precision Data Type l_interval.
Definition: l_interval.hpp:72
cxsc::cosh
cinterval cosh(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:223
cxsc::Pi_l_interval
l_interval Pi_l_interval() noexcept
Enclosure-Interval for .
Definition: l_imath.cpp:1423
cxsc::sinh
cinterval sinh(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:231
cxsc::ln_sqrtx2y2
interval ln_sqrtx2y2(const interval &x, const interval &y) noexcept
Calculates .
Definition: imath.cpp:581
cxsc::pow
cinterval pow(const cinterval &z, const interval &p) noexcept
Calculates .
Definition: cimath.cpp:2074
cxsc::exp2
cinterval exp2(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:167
cxsc::acoth
cinterval acoth(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:3330
cxsc
The namespace cxsc, providing all functionality of the class library C-XSC.
Definition: cdot.cpp:29
cxsc::l_real
The Multiple-Precision Data Type l_real.
Definition: l_real.hpp:78
cxsc::sqrtp1m1
cinterval sqrtp1m1(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:1054
cxsc::expm1
cinterval expm1(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:177
cxsc::diam
cvector diam(const cimatrix_subv &mv) noexcept
Returns the diameter of the matrix.
Definition: cimatrix.inl:738
cxsc::atanh
cinterval atanh(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:3317
cxsc::pow_all
std::list< cinterval > pow_all(const cinterval &z, const interval &p) noexcept
Calculates and returns all possible solutions.
Definition: cimath.cpp:2107
cxsc::l_cinterval
The Multiple-Precision Data Type l_cinterval.
Definition: l_cinterval.hpp:54
cxsc::times2pown
void times2pown(cinterval &x, int n) noexcept
Fast multiplication of reference parameter [z] with .
Definition: cimath.cpp:2059
cxsc::cos
cinterval cos(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:207
cxsc::atan
cinterval atan(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:2938
cxsc::cinterval
The Scalar Type cinterval.
Definition: cinterval.hpp:55
cxsc::acot
cinterval acot(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:3130
cxsc::Pid2_lx_interval
lx_interval Pid2_lx_interval() noexcept
Enclosure-Interval for .
Definition: lx_interval.cpp:6558
cxsc::cinterval::operator=
cinterval & operator=(const real &) noexcept
Implementation of standard assigning operator.
Definition: cinterval.inl:53
cxsc::Ln2_lx_interval
lx_interval Ln2_lx_interval() noexcept
Enclosure-Interval for .
Definition: lx_interval.cpp:4169
cxsc::Disjoint
int Disjoint(const interval &a, const interval &b)
Checks arguments for disjointness.
Definition: interval.cpp:288
cxsc::power_fast
cinterval power_fast(const cinterval &z, int n) noexcept
Calculates .
Definition: cimath.cpp:1520
cxsc::sqrt1mx2
cinterval sqrt1mx2(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:1140
cxsc::acosh
cinterval acosh(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:2732
cxsc::cutint
real cutint(const real &x) noexcept
Returns the truncated integer part of x.
Definition: lx_real.inl:364
cxsc::sqr
cinterval sqr(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:3342
cxsc::sqrtx2y2
interval sqrtx2y2(const interval &x, const interval &y) noexcept
Calculates .
Definition: imath.cpp:80
cxsc::real
The Scalar Type real.
Definition: real.hpp:114
cxsc::asin
cinterval asin(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:2311