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