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