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