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