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