C-XSC - A C++ Class Library for Extended Scientific Computing 2.5.4
l_imath.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_imath.cpp,v 1.38 2014/01/30 17:23:46 cxsc Exp $ */
25
26#include <math.h>
27#include "l_imath.hpp"
28#include "imath.hpp"
29#include "rmath.hpp"
30
31namespace cxsc {
32
33#define CXSC_One 1.
34#define CXSC_Zero 0.
35#define CXSC_MinusOne -1.
36
37l_interval pow(const l_interval & x, const l_interval & e)
38{
39 int stagsave = stagprec,
40 stagmax = 19,
41 intexp;
42
43 bool fertig;
44
45 l_interval y;
46 interval dx = interval(x),
47 de = interval(e),
48 einfachgenau;
49 real supabsde = Sup(abs(de));
50
51 einfachgenau = pow(dx,de);
52
53 fertig = false;
54 if (Inf(de) == Sup(de)) // &&
55 if (supabsde < 32768.0)
56 {
57 intexp = int(_double(real(Sup(e))));
58 if (real(intexp) == Sup(e))
59 {
60 y = power(x,intexp); // Integer-Potenz wesentlich schneller
61 fertig = true;
62 }
63 }
64
65 if (!fertig)
66 {
67 if (Inf(x) < 0.0)
68 cxscthrow(ERROR_LINTERVAL_STD_FKT_OUT_OF_DEF("l_interval pow(const l_interval & x, const l_interval & e)"));
69 else if (Inf(dx) == Sup(dx) && Sup(dx) == CXSC_One)
70 y = x;
71 else if (Inf(de) == Sup(de) && Sup(de) == CXSC_One)
72 y = x;
73 else if (Inf(de) == Sup(de) && Sup(de) == CXSC_Zero)
74 y = 1.0;
75 else
76 {
77 if (stagprec < stagmax)
78 stagprec++;
79 else
80 stagprec = stagmax;
81 y = exp(e*ln(x));
82 stagprec = stagsave;
83 y = adjust(y);
84 y = y & einfachgenau;
85 }
86 }
87
88 return y;
89}
90
91l_interval power(const l_interval &x, int n) // Power(x,n)
92{
93 int stagsave = stagprec,
94 stagmax = 19;
95
96 bool neg = false;
97
98 long int zhi = 2;
99 interval dx = interval(x),
100 einfachgenau;
101 l_interval y, neu;
102
103 einfachgenau = Power(dx,n);
104
105 if (Inf(dx) == Sup(dx) && Sup(dx) == CXSC_One)
106 y = x;
107 else if (n == 0)
108 y = adjust(l_interval(1.0));
109 else
110 {
111 if (stagprec < stagmax)
112 stagprec++;
113 else
114 stagprec = stagmax;
115
116 if (n == 1)
117 y = x;
118 else if (n == 2)
119 y = sqr(x);
120 else
121 {
122 if (n < 0)
123 {
124 neg = true;
125 n = -n;
126 }
127 // Initialisierung
128 if (n%2)
129 y = x;
130 else
131 y = l_interval(1.0); // Praezision wird bei 1 Mult. auf
132 // aktuellen Wert gesetzt;
133 // Berechnugn durch binaere Darstellung der n
134 neu = sqr(x); // neu = x*x;
135 do {
136 if ((n/zhi)%2) y *= neu;
137 zhi += zhi;
138 if (zhi <= n) // letzte Mult. entfaellt --> schneller
139 neu *= neu;
140 } while (zhi <= n);
141
142 if (neg)
143 y = 1.0/(y);
144 }
145 stagprec = stagsave;
146 y = adjust(y);
147 y = y & einfachgenau;
148 }
149
150 return y;
151}
152
153l_interval sqr(const l_interval & x) // Sqr(x)
154{
155 l_interval y;
156
157 if (Inf(x) >= 0.0) /* result = [x.inf*x.inf, x.sup*x.sup] */
158 y = x * x;
159 else if (Sup(x)<= 0.0)
160 {
161 /* result = [x.sup*x.sup, x.inf*x.inf] */
162 y = l_interval (-Sup(x) , -Inf(x));
163 y = (y) * (y);
164 } else
165 {
166 /* result = [0.0, max(x.sup*x.sup,x.inf*x.inf)] */
167 if (abs(Inf(x)) >= abs(Sup(x)))
168 y = l_interval(0.0, abs(Inf(x)));
169 else
170 y = l_interval(0.0, abs(Sup(x)));
171 y = (y) * (y);
172 }
173 return y;
174}
175
177 // Sqrt(x)
178{ // Blomquist: scaling with 2^ex is necessary if expo(Sup(dx)) is too small!
179 int stagsave = stagprec,
180 stagmax = 30,
181 stagsave2;
182 interval dx = interval(x),
183 einfachgenau;
184 l_interval a1,y,t,mt;
185 bool Inf_Zero;
186
187 einfachgenau = sqrt(dx);
188
189 if (Inf(x) < 0.0)
190 cxscthrow(ERROR_LINTERVAL_STD_FKT_OUT_OF_DEF("l_interval sqrt(const l_interval & x)"));
191 else if (Inf(dx) == Sup(dx) &&
192 (Sup(dx) == CXSC_Zero || Sup(dx) == CXSC_One))
193 y = x;
194 else
195 { // scaling necessary if exponent ex < 0
196 l_interval x1 = x;
197 Inf_Zero = (Inf(dx)==0);
198 if (Inf_Zero) x1 = Sup(x1);
199 int ex = expo(Sup(dx));
200 if (ex>0) ex = 0; // scaling unnecessary if ex>0
201 else
202 {
203 ex = -ex; // ex >= 0
204 if (ex > 1023) ex = 1023; // ex==1023 is sufficient
205 if (ex%2) ex--; // ex>=0 is even
206 }
207 if (ex) times2pown(x1,ex); // ex >= 0 ---> exact scaling!
208 // Interval-Newton-methode: y = m(y)-f(m(y))/f'(y)
209 t = sqrt(interval(x1));
210 if (stagprec < stagmax)
211 stagsave2 = stagprec+1;
212 else
213 stagsave2 = stagmax;
214 stagprec = 1;
215 while (stagprec < stagsave2)
216 {
217 stagprec += stagprec;
218 if (stagprec > stagmax)
219 stagprec = stagmax;
220 mt = mid(t); times2pown(t,1);
221 t = mt-((mt*mt-x1)/t);
222 }
223 if (ex) times2pown(t,-ex/2); // ex!=0 --> backscaling with 2^(-ex/2)
224 stagprec = stagsave; // restore the previous precision
225 y = adjust(t); // matching to previous stagprec.
226 if (Inf_Zero) SetInf(y,0.0);
227 y = y & einfachgenau; // seeking optimal inclusion with intersection
228 }
229 return y;
230}
231
232l_interval sqrt(const l_interval &x, int n)
233 // Sqrt(x,n)
234// -1073741822 <= n <= +1073741823, sonst autom. Programm-Abbruch
235// sqrt(x,n) jetzt mit Skalierung --> Hohe Genauigkeit in allen Bereichen!
236// Blomquist, 28.12.03;
237{
238 int stagsave = stagprec,
239 stagshort, staglong, i, ex, N,
240 stagmax = 19,
241 max = 2;
242 l_interval my, fy, corr, y, xx;
243 interval dx = interval(x), einfachgenau;
244
245 einfachgenau = sqrt(dx,n);
246
247 if (Inf(x) < 0.0)
248 cxscthrow(ERROR_LINTERVAL_STD_FKT_OUT_OF_DEF
249 ("l_interval sqrt(const l_interval &x, int n)"));
250 else if (stagprec == 1)
251 y = pow(dx, interval(1.0)/n);
252 else if (Inf(dx) == Sup(dx) && (Sup(dx) == CXSC_Zero ||
253 Sup(dx) == CXSC_One))
254 y = x;
255 else
256 {
257 if (stagprec < stagmax)
258 stagprec++;
259 else stagprec = stagmax;
260
261 while (max < stagprec)
262 max += max; // quadratische Konvergenz kann
263 // nicht eingehalten werden,
264 max += max; // deshalb eine Schleife mehr
265
266 xx = x;
267 ex = expo(Sup(dx));
268 N = -ex;
269 // Skalierung mit 2^N, so dass dx,xx etwa = 1; ----------------
270 times2pown(dx,N);
271 if (N>1023) {
272 times2pown(xx,1023);
273 times2pown(xx,N-1023);
274 }
275 else times2pown(xx,N);
276 // Skalierung beendet, Blomquist 28.12.03 ---------------------
277
278 y = pow(dx, interval(1.0/n));
279 stagprec = 1;
280 // Intervall-Newton-Verfahren
281 for (i = 2; i <= max; i += i)
282 {
283 // Verdoppelung der Genauigkeit:
284 stagshort = stagprec;
285 stagprec += stagprec; // Verdoppelung hier
286 if (stagprec > stagmax)
287 stagprec = stagmax;
288 staglong = stagprec;
289 my = l_interval(mid(y));
290 fy = power(my, n)-xx;
291 // Fehlerauswertung nur in halber Genauigkeit notwendig!
292 stagprec = stagshort;
293 corr = fy/(real(n)*power(y, n-1));
294 stagprec = staglong;
295 // Fehlerkorrektur in normaler Genauigkeit:
296 y = my-corr;
297 }
298 // Rueckskalierung mit dem Faktor 2^(-N/n):
299 fy = l_interval(-N)/n;
300 y *= exp(fy*li_ln2()); // li_ln2() = ln(2)
301
302 stagprec = stagsave;
303 y = adjust(y); // Anpassung an Ausgangs-stagprec = stagsave.
304 y = y & einfachgenau; // Falls y breiter ist als einfachgenau
305 }
306
307 return y;
308} // sqrt(x,n)
309
311// Calculation of an optimal inclusion of sqrt(1+x^2); Blomquist, 13.12.02;
312// With stagmax=19 we get about 16*19=304 exact decimal digits.
313{ // First step: Calculation of sqrt(1+x*x) in simple type interval:
314 interval einfach = sqrt1px2( interval(x) ); // Only interval types
315 int stagsave, stagmax=19;
316 stagsave = stagprec;
317 if (stagprec > stagmax) stagprec = stagmax;
318 // Second step: Inclusion of sqrt(1+x^2) with stagprec <= 19
319 const int exmax=512;
320 l_interval y = abs(x);
321 l_real lr = Sup( 1 + l_interval(Sup(y)) );
322 interval z = interval(Sup(x));
323 int ex = expo(Sup(z));
324 if (ex > exmax)
325 { // scaling to avoid overflow by sqr(y):
326 ex = exmax - ex; // ex = 512 - ex;
327 times2pown(y,ex); // scaling to avoid overflow by sqr(y)
328 y = sqrt( comp(0.5,2*ex+1) + sqr(y) ); // sqr(y) without overflow!
329 times2pown(y,-ex); // backscaling of the result with 2^(-ex)
330 } else
331 { // no scaling necessary:
332 y = sqrt(1+sqr(x));
333 }
334 if (Inf(y)<1.0) SetInf(y,1.0); // improvement, Blomquist, 25.02.07;
335 if (Sup(y)>lr) SetSup(y,lr); // improvement, Blomquist, 26.02.07;
336 stagprec = stagsave; // restore the old stagprec value
337 y = adjust(y); // y gets the previous precision
338 y = einfach & y; // This intersection delivers for intervals x with large
339 // diameters the optimal result inclusion
340 return y;
341} // sqrt1px2(...)
342
343l_interval sqrtx2y2(const l_interval& x, const l_interval& y) noexcept
344// Inclusion of sqrt(x^2+y^2); Blomquist, 14.12.02;
345{
346 interval ia = abs(interval(x)), ib = abs(interval(y)), einfach;
347 einfach = sqrtx2y2(ia,ib); // Inclusion only with type interval.
348 if (!einfach) return l_interval(0.0);
349
350 int stagsave=stagprec, stagmax=19;
351 if (stagprec>stagmax) stagprec = stagmax;
352
353 l_interval a=abs(x), b=abs(y), r;
354 int exa=expo(Sup(ia)), exb=expo(Sup(ib)), ex;
355 if (exb > exa)
356 { // Permutation of a,b:
357 r = a; a = b; b = r;
358 ex = exa; exa = exb; exb = ex;
359 }
360 ex = 511 - exa; exa = 0;
361 if (ex>1022)
362 { exa = ex - 1022; ex = 1022; } // ex > 1022 --> scaling in two steps
363 times2pown(a,ex); // is necessary!
364 if (exa) times2pown(a,exa);
365 times2pown(b,ex); // First step: scaling b with 2^ex;
366 if (exa) times2pown(b,exa); // Second step: scaling b with 2^exa;
367 r = sqrt(a*a + b*b);
368 times2pown(r,-ex); // Backscaling, first step
369 if (exa) times2pown(r,-exa); // Backscaling, second step
370 stagprec = stagsave;
371 r = adjust(r);
372 r = einfach & r;
373 return r;
374} // sqrtx2y2
375
377// sqrtp1m1(x) calculates an inclusion of sqrt(x+1)-1;
378// Blomquist, 05.08.03;
379{
380 int stagsave=stagprec, stagmax=19;
381 stagprec++;
382 if (stagprec>stagmax) stagprec = stagmax;
383 l_interval y,tmp;
384 interval z = interval(x); // z is an inclusion of x
385 real r = Inf(z);
386 if (r < -1)
387 cxscthrow(STD_FKT_OUT_OF_DEF("l_interval sqrtp1m1(const l_interval&)"));
388 const real c = 1e-10;
389 tmp = x+1;
390 y = x<=interval(-c,c) ? x / (sqrt(tmp)+1) : sqrt(tmp)-1;
391 stagprec = stagsave;
392 y = adjust(y);
393 return y;
394} // sqrtp1m1
395
396
397l_interval sin(const l_interval & x) // Sin(x)
398{
399 int stagsave = stagprec,
400 stagmax = 19;
401 l_interval pihalbe,
402 y;
403 interval dx = interval(x),
404 einfachgenau;
405
406 einfachgenau = sin(dx);
407
408 if (stagprec == 1)
409 y = sin(dx);
410 else if (Inf(dx) == Sup(dx) && Sup(dx) == CXSC_Zero)
411 y = adjust(l_interval(0.0));
412 else
413 {
414 if (stagprec < stagmax)
415 stagprec++;
416 else
417 stagprec = stagmax;
418
419// pihalbe = 2.0*atan(l_interval(1.0));
420 pihalbe = li_pi4();
421 times2pown(pihalbe,1); // Blomquist, 05.12.03;
422 y = x-pihalbe;
423
424 try {
425 y = cos(y);
426 }
427 catch(const ERROR_LINTERVAL_FAK_OVERFLOW &) // Damit Funktionsname stimmt
428 {
429 cxscthrow(ERROR_LINTERVAL_FAK_OVERFLOW("l_interval sin(const l_interval & x)"));
430 }
431
432
433 stagprec = stagsave;
434 y = adjust(y);
435 y = y & einfachgenau;
436 }
437
438 return y;
439}
440
441l_interval cos(const l_interval & x) // Cos(x)
442{
443 long int mm = 6;
444 int stagsave = stagprec,
445 stagmax = 19,
446 n = 0,
447 digits = 53,
448 sign = 0,
449 degree, k,
450 m2;
451
452 bool xinf = false,
453 xsup = false,
454 fertig = false;
455
456 real abst, m, eps,
457// bas, tn, t4,
458 lneps, lnt,
459 zk,
460 zhn = 1.0,
461 lnb = 0.69314718,
462 fak = 720.0; // 6!
463 interval dx = interval(x),
464 einfachgenau,
465 dt, extr, error;
466 l_interval zwopi, t, t2, p, ph, test,
467 y;
468
469 einfachgenau = cos(dx);
470 if (stagprec == 1)
471 y = cos(dx);
472 else if (Inf(dx) == Sup(dx) && Sup(dx) == CXSC_Zero)
473 y = adjust(l_interval(1.0));
474 else
475 {
476 if (stagprec < stagmax)
477 stagprec++;
478 else
479 stagprec = stagmax;
480
481 // stagprec++;
482 // zwopi = 8.0*atan(l_interval(1.0));
483 zwopi = li_pi4();
484 times2pown(zwopi,3); // Blomquist, 05.12.03;
485 // stagprec--;
486
487 // erste Argumentreduktion ( cos(x) = cos(x+2k*pi) )
488 if (Sup(zwopi) < Sup(abs(x)))
489 {
490 m = floor(_double(real(Sup((x/zwopi))))); // floor rundet zur naechsten
491 t = x-m*zwopi; // ganzen Zahl kleiner x
492 if (Sup(zwopi) < Sup(abs(t)))
493 { // das ganze nochmal, falls m uebergelaufen ist!
494 // ohne Sup wird Inclusion geprueft!
495 m = floor(_double(real(Sup((t/zwopi)))));// rundet zur naechsten Zahl < x
496 t = t-m*zwopi;
497 }
498 } else
499 t = x;
500
501 // ueberpruefen, ob Maximum oder Minimum im Inneren des Intervalls
502 extr = interval(2.0/zwopi*t);
503 m2 = int(double(floor(_double(Sup(extr)))));
504 if (interval(real(m2)) <= extr)
505 {
506 if (interval(real(m2-1),real(m2)) <= extr)
507 {
508 y = l_interval(-1.0,1.0);
509 fertig = true;
510 } else
511 if (m2%2)
512 xinf = TRUE;
513 else
514 xsup = TRUE;
515 }
516
517 if (!(fertig))
518 {
519 // zweite Argumentreduktion
520 dt = interval(t);
521 // eps = 0.01/(stagprec*stagprec*stagprec);
522 eps = 0.01;
523 while (Sup(abs(dt))/zhn >= eps)
524 {
525 n++;
526 zhn += zhn;
527 }
528 t /= zhn;
529
530 // Abschaetzung der Genauigkeit
531
532 t2 = t*t;
533 abst = real(Sup(abs(t)));
534 if (abst < MinReal)
535 abst = MinReal; // nur zur Sicherheit
536 lnt = ln(abst);
537 lneps = (1.0-digits*stagprec)*lnb;
538 while (lneps-(real(mm)*lnt+ln(2.0/fak)) <= 0.0)
539 {
540 mm += 4;
541 if (mm > 170)
542 {
543 // 170! = 7.26*10^306
544 cxscthrow(ERROR_LINTERVAL_FAK_OVERFLOW("l_interval cos(const l_interval & x)"));
545 mm = 170;
546 break;
547 }
548 fak *= mm*(mm-1.0)*(mm-2.0)*(mm-3.0);
549 }
550 /*
551 bas = real(Inf(power(l_interval(2.0),(1-digits*stagprec))));
552 tn = abs(real(Sup(power(t,6))));
553 t4 = real(mid(power(t,4)));
554 while (bas-2.0*tn/fak <= 0.0)
555 {
556 mm += 4;
557 if (mm > 170)
558 { // 170! = 7.26*10^306
559 errmon (ERROR_LINTERVAL(FAKOVERFLOW));
560 mm = 170;
561 break;
562 }
563 tn *= t4;
564 fak *= mm*(mm-1.0)*(mm-2.0)*(mm-3.0);
565 }
566 */
567
568 degree = mm-2; // Achtung mm := 2n+2 !
569
570 // Polynomauswertung
571
572 sign = (degree/2)%2;
573 zk = real(degree)*real(degree-1);
574 if (sign)
575 p = -t2/zk;
576 else
577 p = t2/zk;
578 for (k = degree-2; k >= 2; k -= 2)
579 {
580 sign = 1-sign;
581 if (sign)
582 p -= 1.0;
583 else
584 p += 1.0;
585 zk = real(k)*real(k-1);
586 p *= t2/zk;
587 }
588
589 error = pow(interval(2.0), interval(1.0-digits*stagprec))
590 * interval(-0.5,0.5);
591 p = l_interval(1.0)+error+p;
592
593 // Rueckgaengigmachung der zweiten Argumentreduktion
594 for (int i = 0; i < n; i++)
595 p = 2.0*p*p-1.0;
596
597 stagprec = stagsave;
598 y = adjust(p);
599 if (Inf(y) < -1.0)
600 SetInf(y,-1.0);
601 if (Sup(y) > 1.0)
602 SetSup(y,1.0);
603 if (xinf)
604 SetInf(y,-1.0);
605 if (xsup)
606 SetSup(y,1.0);
607 }
608 y = y & einfachgenau;
609 }
610 return y;
611}
612
613l_interval tan(const l_interval & x) // Tan(x)
614{
615 interval dx = interval(x),
616 einfachgenau;
617 l_interval s, c, y;
618
619 einfachgenau = tan(dx);
620
621 if (stagprec == 1)
622 y = tan(dx);
623 else if (Inf(dx) == Sup(dx) && Sup(dx) == CXSC_Zero)
624 y = adjust(l_interval(0.0));
625 else
626 {
627 try
628 {
629 c = cos(x);
630 }
631 catch(const ERROR_LINTERVAL_FAK_OVERFLOW &) // Damit Funktionsname stimmt
632 {
633 cxscthrow(ERROR_LINTERVAL_FAK_OVERFLOW("l_interval tan(const l_interval &x)"));
634 }
635
636 if (interval(0.0) <= c)
637 {
638 cxscthrow(ERROR_LINTERVAL_STD_FKT_OUT_OF_DEF("l_interval tan(const l_interval &x)"));
639 }
640 try
641 {
642 s = sin(x);
643 }
644 catch(const ERROR_LINTERVAL_FAK_OVERFLOW &)
645 {
646 cxscthrow(ERROR_LINTERVAL_FAK_OVERFLOW("l_interval tan(const l_interval &x)"));
647 }
648 stagprec++;
649 y = s/c;
650 stagprec--;
651 y = adjust(y);
652 y = y & einfachgenau;
653 }
654 return y;
655}
656
657l_interval cot(const l_interval & x) // Cot(x)
658{
659 interval dx = interval(x),
660 einfachgenau;
661 l_interval s, c, y;
662
663 einfachgenau = cot(dx);
664
665 if (stagprec == 1)
666 y = tan(dx);
667 else if (Inf(dx) == Sup(dx) && Sup(dx) == CXSC_Zero)
668 y = adjust(l_interval(0.0));
669 else
670 {
671 try
672 {
673 s = sin(x);
674 }
675 catch(const ERROR_LINTERVAL_FAK_OVERFLOW &) // Damit Funktionsname stimmt
676 {
677 cxscthrow(ERROR_LINTERVAL_FAK_OVERFLOW("l_interval cot(const l_interval &x)"));
678 }
679
680 if (interval(0.0) <= s)
681 {
682 cxscthrow(ERROR_LINTERVAL_STD_FKT_OUT_OF_DEF("l_interval cot(const l_interval &x)"));
683 }
684
685 try
686 {
687 c = cos(x);
688 }
689 catch(const ERROR_LINTERVAL_FAK_OVERFLOW &) // Damit Funktionsname stimmt
690 {
691 cxscthrow(ERROR_LINTERVAL_FAK_OVERFLOW("l_interval cot(const l_interval &x)"));
692 }
693
694 c = cos(x);
695 stagprec++;
696 y = c/s;
697 stagprec--;
698 y = adjust(y);
699 y = y & einfachgenau;
700 }
701
702 return y;
703}
704
705l_interval asin(const l_interval & x) // ASin(x)
706{
707 l_interval t, ta, u, pihalbe,
708 y;
709 interval dx = interval(x),
710 einfachgenau;
711 real supabsdx = Sup(abs(dx)),
712 infdx = Inf(dx),
713 supdx = Sup(dx);
714
715 einfachgenau = asin(dx);
716
717 stagprec++;
718// pihalbe = 2.0*atan(l_interval(1.0));
719 pihalbe = li_pi4();
720 times2pown(pihalbe,1); // Blomquist, 05.12.03;
721 stagprec--;
722
723 if (Inf(x) < CXSC_MinusOne || Sup(x) > CXSC_One)
724 cxscthrow(ERROR_LINTERVAL_STD_FKT_OUT_OF_DEF("l_interval asin(const l_interval & x)"));
725 else if (stagprec == 1)
726 y = asin(dx);
727 else if (infdx == supdx && supdx == CXSC_Zero)
728 y = adjust(l_interval(0.0));
729 else if (infdx == supdx && supabsdx == CXSC_One)
730 {
731 if (supdx == 1.0)
732 y = pihalbe;
733 else
734 y = -pihalbe;
735 } else
736 {
737 stagprec++;
738 try
739 {
740 if (supabsdx <= 0.75)
741 u = x;
742 else
743 u = 2.0*x*sqrt((1.0-x)*(1.0+x));
744 t = u/sqrt((1.0-u)*(1.0+u));
745 stagprec--;
746 ta = atan(t);
747 }
748 catch(const ERROR_LINTERVAL_STD_FKT_OUT_OF_DEF &)
749 {
750 cxscthrow(ERROR_LINTERVAL_STD_FKT_OUT_OF_DEF("l_interval asin(const l_interval & x)"));
751 }
752 stagprec++;
753 if (supabsdx <= 0.75)
754 y = ta;
755 else if (Sup(t) >= 0.0)
756 y = pihalbe-0.5*ta;
757 else
758 y = -pihalbe-0.5*ta;
759 stagprec--;
760 y = adjust(y);
761 y = y & einfachgenau;
762 }
763
764 return y;
765}
766
767l_interval acos(const l_interval & x) // ACos(x)
768{
769 bool neg=false;
770 l_interval pi, s, y;
771 interval dx = interval(x),
772 einfachgenau;
773 real supabsdx = Sup(abs(dx)),
774 infdx = Inf(dx),
775 supdx = Sup(dx);
776
777 try {
778
779 einfachgenau = acos(dx);
780
781 stagprec++;
782// pi = 4.0*atan(l_interval(1.0));
783 pi = li_pi4();
784 times2pown(pi,2); // Blomquist, 05.12.03;
785 stagprec--;
786 if (Inf(x) < CXSC_MinusOne || Sup(x) > CXSC_One)
787 cxscthrow(ERROR_LINTERVAL_STD_FKT_OUT_OF_DEF());
788 else if (stagprec == 1)
789 y = acos(dx);
790 else if (infdx == supdx && supabsdx == CXSC_One)
791 {
792 if (supdx == 1.0)
793 y = adjust(l_interval(0.0));
794 else
795 y = adjust(pi);
796 } else
797 {
798 if (supdx < 0.0)
799 {
800 y = -x;
801 neg = true;
802 } else
803 {
804 y = x;
805 neg = false;
806 }
807 stagprec++;
808 if (supabsdx > 0.75)
809 {
810 y = (1.0-(y))*(1.0+(y));
811 // stagprec--;
812 y = asin(sqrt(y));
813 // stagprec++;
814 if (neg)
815 y = pi-(y);
816 } else
817 {
818 // stagprec--;
819 s = asin(y);
820 // stagprec++;
821 if (neg)
822 y = 0.5*pi+s;
823 else
824 y = 0.5*pi-s;
825 }
826
827 // Fehler in der Systemumgebung, deshalb die Aufblaehung des Intervalls
828 real err1, err2;
829 interval error;
830 err1 = 5.0*real(diam(y));
831 err2 = 5.0*power(10.0,-16*(stagprec-1)-1);
832 if (err1 > err2)
833 error = interval(-err1, err1);
834 else
835 error = interval(-err2, err2);
836 y += error;
837
838 stagprec--;
839 y = adjust(y);
840 y = y & einfachgenau;
841 }
842 }
843 catch(const ERROR_LINTERVAL_STD_FKT_OUT_OF_DEF &)
844 {
845 cxscthrow(ERROR_LINTERVAL_STD_FKT_OUT_OF_DEF("l_interval acos(const l_interval & x)"));
846 }
847 return y;
848}
849
850static real CXSC_ln2[21]; // CXSC_ln2[0], ... CXSC_ln2[20]
851static bool CXSC_ln2_initialized = false;
852
853// l_interval li_ln2() noexcept
855// Inclusion of ln(2), Blomquist, 04.12.03;
856{
857 l_interval y;
858 int stagsave = stagprec,
859 stagmax = 20;
860
861 if (!CXSC_ln2_initialized)
862 {
863 std::string str;
864 std::cout << SaveOpt;
865 std::cout << Hex;
866 str = "+162E42FEFA39EFe3FE";
867 str >> CXSC_ln2[0];
868 str = "+1ABC9E3B39803Fe3C7";
869 str >> CXSC_ln2[1];
870 str = "+17B57A079A1934e390";
871 str >> CXSC_ln2[2];
872 str = "-1ACE93A4EBE5D1e35A";
873 str >> CXSC_ln2[3];
874 str = "-123A2A82EA0C24e324";
875 str >> CXSC_ln2[4];
876 str = "+1D881B7AEB2615e2ED";
877 str >> CXSC_ln2[5];
878 str = "+19552FB4AFA1B1e2B7";
879 str >> CXSC_ln2[6];
880 str = "+1DA5D5C6B82704e27E";
881 str >> CXSC_ln2[7];
882 str = "+14427573B29117e247";
883 str >> CXSC_ln2[8];
884 str = "-191F6B05A4D7A7e211";
885 str >> CXSC_ln2[9];
886 str = "-1DB5173AE53426e1DB";
887 str >> CXSC_ln2[10];
888 str = "+11317C387EB9EBe1A3";
889 str >> CXSC_ln2[11];
890 str = "-190F13B267F137e16D";
891 str >> CXSC_ln2[12];
892 str = "+16FA0EC7657F75e137";
893 str >> CXSC_ln2[13];
894 str = "-1234C5E1398A6Be101";
895 str >> CXSC_ln2[14];
896 str = "+1195EBBF4D7A70e0CA";
897 str >> CXSC_ln2[15];
898 str = "+18192432AFD0C4e094";
899 str >> CXSC_ln2[16];
900 str = "-1A1BE38BA4BA4De05E";
901 str >> CXSC_ln2[17];
902 str = "-1D7860151CFC06e024";
903 str >> CXSC_ln2[18];
904 str = "+1000032847ED6Fe000";
905 str >> CXSC_ln2[19];
906 str = "+1000032847ED70e000";
907 str >> CXSC_ln2[20];
908
909 CXSC_ln2_initialized = true;
910 std::cout << RestoreOpt;
911 }
912 stagprec = stagmax;
913 y = adjust(l_interval(0.0));
914 for (int i=0; i <= stagmax; i++)
915 y.data[i] = CXSC_ln2[i];
916 stagprec = stagsave;
917 y = adjust(y);
918 return y;
919}
920
921static real CXSC_ln10[21]; // CXSC_ln10[0], ... CXSC_ln10[20]
922static bool CXSC_ln10_initialized = false;
923
924// l_interval li_ln10() noexcept
926{
927 l_interval y;
928 int stagsave = stagprec,
929 stagmax = 20;
930
931 if (!CXSC_ln10_initialized)
932 {
933 std::string str;
934 std::cout << SaveOpt;
935 std::cout << Hex;
936 str = "+126BB1BBB55516e400";
937 str >> CXSC_ln10[0];
938 str = "-1F48AD494EA3E9e3CA";
939 str >> CXSC_ln10[1];
940 str = "-19EBAE3AE0260Ce394";
941 str >> CXSC_ln10[2];
942 str = "-12D10378BE1CF1e35E";
943 str >> CXSC_ln10[3];
944 str = "+10403E05AE52C6e328";
945 str >> CXSC_ln10[4];
946 str = "-1FA509CAFDF466e2F0";
947 str >> CXSC_ln10[5];
948 str = "-1C79A1FE9D0795e2BA";
949 str >> CXSC_ln10[6];
950 str = "+1058C448308218e284";
951 str >> CXSC_ln10[7];
952 str = "-1D250470877BFDe24D";
953 str >> CXSC_ln10[8];
954 str = "-1AE92987D3075De215";
955 str >> CXSC_ln10[9];
956 str = "-1D5CDBB8626956e1DF";
957 str >> CXSC_ln10[10];
958 str = "-13C4F27CE0410Ae1A9";
959 str >> CXSC_ln10[11];
960 str = "+1B3AC12ACF1BE9e173";
961 str >> CXSC_ln10[12];
962 str = "+1161BB49D219C8e13D";
963 str >> CXSC_ln10[13];
964 str = "-110D6613293728e107";
965 str >> CXSC_ln10[14];
966 str = "+142163A4CDA351e0CF";
967 str >> CXSC_ln10[15];
968 str = "+1E2713D6C22C16e097";
969 str >> CXSC_ln10[16];
970 str = "-15090EF85CB0ADe05E";
971 str >> CXSC_ln10[17];
972 str = "-1C5B3E859F876Ee027";
973 str >> CXSC_ln10[18];
974 str = "-10000703552C52e000";
975 str >> CXSC_ln10[19];
976 str = "-10000703552C51e000";
977 str >> CXSC_ln10[20];
978
979 CXSC_ln10_initialized = true;
980 std::cout << RestoreOpt;
981 }
982 stagprec = stagmax;
983 y = adjust(l_interval(0.0));
984 for (int i=0; i <= stagmax; i++)
985 y.data[i] = CXSC_ln10[i];
986 stagprec = stagsave;
987 y = adjust(y);
988 return y;
989}
990
991static real CXSC_Rln10[21]; // CXSC_Rln10[0], ... CXSC_Rln10[20]
992static bool CXSC_Rln10_initialized = false;
993
994// l_interval li_Rln10() noexcept
996{
997 l_interval y;
998 int stagsave = stagprec,
999 stagmax = 20;
1000
1001 if (!CXSC_Rln10_initialized)
1002 {
1003 std::string str;
1004 std::cout << SaveOpt;
1005 std::cout << Hex;
1006 str = "+1BCB7B1526E50Ee3FD";
1007 str >> CXSC_Rln10[0];
1008 str = "+195355BAAAFAD3e3C6";
1009 str >> CXSC_Rln10[1];
1010 str = "+1EE191F71A3012e38F";
1011 str >> CXSC_Rln10[2];
1012 str = "+17268808E8FCB5e358";
1013 str >> CXSC_Rln10[3];
1014 str = "+13DE3A94F1D509e320";
1015 str >> CXSC_Rln10[4];
1016 str = "+1DF42805E7E524e2E9";
1017 str >> CXSC_Rln10[5];
1018 str = "+11AAC96323250Be2B3";
1019 str >> CXSC_Rln10[6];
1020 str = "-1CE63884C058E4e27D";
1021 str >> CXSC_Rln10[7];
1022 str = "-1A1C82EA3969BAe247";
1023 str >> CXSC_Rln10[8];
1024 str = "+1B4F6686AD7A33e211";
1025 str >> CXSC_Rln10[9];
1026 str = "-1B97C8035FFC70e1DB";
1027 str >> CXSC_Rln10[10];
1028 str = "+1630771369962Ee1A0";
1029 str >> CXSC_Rln10[11];
1030 str = "-1E15BD37B295AFe16A";
1031 str >> CXSC_Rln10[12];
1032 str = "-132484B432318Be134";
1033 str >> CXSC_Rln10[13];
1034 str = "+15430212AE68C0e0FE";
1035 str >> CXSC_Rln10[14];
1036 str = "+1351923B322731e0C8";
1037 str >> CXSC_Rln10[15];
1038 str = "+11F934D794D64Fe092";
1039 str >> CXSC_Rln10[16];
1040 str = "+13E4B475D9FF20e05B";
1041 str >> CXSC_Rln10[17];
1042 str = "+185D9B63ED9A24e025";
1043 str >> CXSC_Rln10[18];
1044 str = "+1000035B8CA18Ce000";
1045 str >> CXSC_Rln10[19];
1046 str = "+1000035B8CA18De000";
1047 str >> CXSC_Rln10[20];
1048
1049 CXSC_Rln10_initialized = true;
1050 std::cout << RestoreOpt;
1051 }
1052 stagprec = stagmax;
1053 y = adjust(l_interval(0.0));
1054 for (int i=0; i <= stagmax; i++)
1055 y.data[i] = CXSC_Rln10[i];
1056 stagprec = stagsave;
1057 y = adjust(y);
1058 return y;
1059}
1060
1061static real CXSC_pi4[21];
1062static bool CXSC_pi4_initialized = false;
1063
1064// l_interval li_pi4() noexcept
1066{
1067 l_interval y;
1068 int stagsave = stagprec,
1069 stagmax = 20;
1070
1071 if (!CXSC_pi4_initialized)
1072 {
1073 std::string str;
1074 std::cout << SaveOpt;
1075 std::cout << Hex;
1076 str = "+1921FB54442D18e3FE";
1077 str >> CXSC_pi4[0];
1078 str = "+11A62633145C06e3C8";
1079 str >> CXSC_pi4[1];
1080 str = "+1C1CD129024E08e393";
1081 str >> CXSC_pi4[2];
1082 str = "+114CF98E804178e35E";
1083 str >> CXSC_pi4[3];
1084 str = "-159C4EC64DDAECe327";
1085 str >> CXSC_pi4[4];
1086 str = "+1410F31C6809BCe2F2";
1087 str >> CXSC_pi4[5];
1088 str = "-106AE64C32C5BCe2BB";
1089 str >> CXSC_pi4[6];
1090 str = "-1C99FA9EB241B4e286";
1091 str >> CXSC_pi4[7];
1092 str = "-1D791603D95252e24E";
1093 str >> CXSC_pi4[8];
1094 str = "-1571EDD0DBD254e218";
1095 str >> CXSC_pi4[9];
1096 str = "-133B4302721768e1E2";
1097 str >> CXSC_pi4[10];
1098 str = "+10BA698DFB5AC2e1AD";
1099 str >> CXSC_pi4[11];
1100 str = "+1FFAE5B7A035C0e178";
1101 str >> CXSC_pi4[12];
1102 str = "-1211C79404A576e143";
1103 str >> CXSC_pi4[13];
1104 str = "-1816945836FBA0e10D";
1105 str >> CXSC_pi4[14];
1106 str = "-1DA700CDB6BCCEe0D8";
1107 str >> CXSC_pi4[15];
1108 str = "+11ECE45B3DC200e0A3";
1109 str >> CXSC_pi4[16];
1110 str = "+1F2E2858EFC166e06D";
1111 str >> CXSC_pi4[17];
1112 str = "+1B4906C38ABA73e036";
1113 str >> CXSC_pi4[18];
1114 str = "+19A458FEA3F493e000";
1115 str >> CXSC_pi4[19];
1116 str = "+19A458FEA3F494e000";
1117 str >> CXSC_pi4[20];
1118
1119 CXSC_pi4_initialized = true;
1120 std::cout << RestoreOpt;
1121 }
1122 stagprec = stagmax;
1123 y = adjust(l_interval(0.0));
1124 for (int i=0; i <= stagmax; i++)
1125 y.data[i] = CXSC_pi4[i];
1126 stagprec = stagsave;
1127 y = adjust(y);
1128 return y;
1129}
1130
1131static real CXSC_sqrt2[21];
1132static bool CXSC_sqrt2_initialized = false;
1133
1134// l_interval li_sqrt2() noexcept
1136{
1137 l_interval y;
1138 int stagsave = stagprec,
1139 stagmax = 20;
1140
1141 if (!CXSC_sqrt2_initialized)
1142 {
1143 std::string str;
1144 std::cout << SaveOpt;
1145 std::cout << Hex;
1146 str = "+16A09E667F3BCDe3FF";
1147 str >> CXSC_sqrt2[0];
1148 str = "-1BDD3413B26456e3C9";
1149 str >> CXSC_sqrt2[1];
1150 str = "+157D3E3ADEC175e393";
1151 str >> CXSC_sqrt2[2];
1152 str = "+12775099DA2F59e35B";
1153 str >> CXSC_sqrt2[3];
1154 str = "+160CCE64552BF2e322";
1155 str >> CXSC_sqrt2[4];
1156 str = "+1821D5C5161D46e2E9";
1157 str >> CXSC_sqrt2[5];
1158 str = "-1C032046F8498Ee2B3";
1159 str >> CXSC_sqrt2[6];
1160 str = "+1EE950BC8738F7e27B";
1161 str >> CXSC_sqrt2[7];
1162 str = "-1AC3FDBC64E103e245";
1163 str >> CXSC_sqrt2[8];
1164 str = "+13B469101743A1e20D";
1165 str >> CXSC_sqrt2[9];
1166 str = "+15E3E9CA60B38Ce1D7";
1167 str >> CXSC_sqrt2[10];
1168 str = "+11BC337BCAB1BDe19C";
1169 str >> CXSC_sqrt2[11];
1170 str = "-1BBA5DEE9D6E7De166";
1171 str >> CXSC_sqrt2[12];
1172 str = "-1438DD083B1CC4e130";
1173 str >> CXSC_sqrt2[13];
1174 str = "+1B56A28E2EDFA7e0FA";
1175 str >> CXSC_sqrt2[14];
1176 str = "+1CCB2A634331F4e0C4";
1177 str >> CXSC_sqrt2[15];
1178 str = "-1BD9056876F83Ee08D";
1179 str >> CXSC_sqrt2[16];
1180 str = "-1234FA22AB6BEFe057";
1181 str >> CXSC_sqrt2[17];
1182 str = "+19040CA4A81395e020";
1183 str >> CXSC_sqrt2[18];
1184 str = "-1000002A493818e000";
1185 str >> CXSC_sqrt2[19];
1186 str = "-1000002A493817e000";
1187 str >> CXSC_sqrt2[20];
1188
1189 CXSC_sqrt2_initialized = true;
1190 std::cout << RestoreOpt;
1191 }
1192 stagprec = stagmax;
1193 y = adjust(l_interval(0.0));
1194 for (int i=0; i <= stagmax; i++)
1195 y.data[i] = CXSC_sqrt2[i];
1196 stagprec = stagsave;
1197 y = adjust(y);
1198 return y;
1199}
1200
1201static real CXSC_sqrt5[21]; // CXSC_sqrt5[0], ... CXSC_sqrt5[20]
1202static bool CXSC_sqrt5_initialized = false;
1203
1205// Inclusion of sqrt(5), Blomquist, 30.11.2008;
1206{
1207 l_interval y;
1208 int stagsave = stagprec,
1209 stagmax = 20;
1210
1211 if (!CXSC_sqrt5_initialized)
1212 {
1213 std::string str;
1214 std::cout << SaveOpt;
1215 std::cout << Hex;
1216 str = "+11E3779B97F4A8e400";
1217 str >> CXSC_sqrt5[0];
1218 str = "-1F506319FCFD19e3C9";
1219 str >> CXSC_sqrt5[1];
1220 str = "+1B906821044ED8e393";
1221 str >> CXSC_sqrt5[2];
1222 str = "-18BB1B5C0F272Ce35B";
1223 str >> CXSC_sqrt5[3];
1224 str = "+11D0C18E952768e324";
1225 str >> CXSC_sqrt5[4];
1226 str = "-1E9D585B0901F9e2EB";
1227 str >> CXSC_sqrt5[5];
1228 str = "-1C7DD252073EC0e2B5";
1229 str >> CXSC_sqrt5[6];
1230 str = "-1FCEF21EDAF7FAe27F";
1231 str >> CXSC_sqrt5[7];
1232 str = "+160EB25D20799Be241";
1233 str >> CXSC_sqrt5[8];
1234 str = "-1C90F95285168Fe208";
1235 str >> CXSC_sqrt5[9];
1236 str = "+1E1DFA160E75BCe1D2";
1237 str >> CXSC_sqrt5[10];
1238 str = "-10A08E66CB368Ce196";
1239 str >> CXSC_sqrt5[11];
1240 str = "+1C5371682CADD1e160";
1241 str >> CXSC_sqrt5[12];
1242 str = "-1998100220F4EDe129";
1243 str >> CXSC_sqrt5[13];
1244 str = "+1C6771A0968663e0F3";
1245 str >> CXSC_sqrt5[14];
1246 str = "+1DFB9E3C86CA7Ce0BD";
1247 str >> CXSC_sqrt5[15];
1248 str = "-18AE38ED5304B1e086";
1249 str >> CXSC_sqrt5[16];
1250 str = "+182A5FEC507706e050";
1251 str >> CXSC_sqrt5[17];
1252 str = "-1B5191A18C5647e018";
1253 str >> CXSC_sqrt5[18];
1254 str = "+100000000F9D52e000";
1255 str >> CXSC_sqrt5[19];
1256 str = "+100000000F9D53e000";
1257 str >> CXSC_sqrt5[20];
1258
1259 CXSC_sqrt5_initialized = true;
1260 std::cout << RestoreOpt;
1261 }
1262 stagprec = stagmax;
1263 y = adjust(l_interval(0.0));
1264 for (int i=0; i <= stagmax; i++)
1265 y.data[i] = CXSC_sqrt5[i];
1266 stagprec = stagsave;
1267 y = adjust(y);
1268 return y;
1269}
1270
1271// ************************************************************************
1272
1273static real CXSC_sqrt7[21]; // CXSC_sqrt7[0], ... CXSC_sqrt7[20]
1274static bool CXSC_sqrt7_initialized = false;
1275
1277// Inclusion of sqrt(7), Blomquist, 30.11.2008;
1278{
1279 l_interval y;
1280 int stagsave = stagprec,
1281 stagmax = 20;
1282
1283 if (!CXSC_sqrt7_initialized)
1284 {
1285 std::string str;
1286 std::cout << SaveOpt;
1287 std::cout << Hex;
1288 str = "+152A7FA9D2F8EAe400";
1289 str >> CXSC_sqrt7[0];
1290 str = "-121C62B033C079e3CA";
1291 str >> CXSC_sqrt7[1];
1292 str = "-177CAAD6200612e391";
1293 str >> CXSC_sqrt7[2];
1294 str = "-1EFA880DC72D64e359";
1295 str >> CXSC_sqrt7[3];
1296 str = "-171D206D5B1A4Ce31F";
1297 str >> CXSC_sqrt7[4];
1298 str = "+119392FA9B0494e2E6";
1299 str >> CXSC_sqrt7[5];
1300 str = "+17BB8A64890057e2AD";
1301 str >> CXSC_sqrt7[6];
1302 str = "-17E89300383DDEe277";
1303 str >> CXSC_sqrt7[7];
1304 str = "+130FB7AF68A6FBe241";
1305 str >> CXSC_sqrt7[8];
1306 str = "+1322281D303D36e209";
1307 str >> CXSC_sqrt7[9];
1308 str = "+1996109A16D3B1e1D3";
1309 str >> CXSC_sqrt7[10];
1310 str = "+1F239C301DFBB4e19C";
1311 str >> CXSC_sqrt7[11];
1312 str = "-1B5CA40AB771A2e163";
1313 str >> CXSC_sqrt7[12];
1314 str = "-1675711487FEAAe12A";
1315 str >> CXSC_sqrt7[13];
1316 str = "+122CB7FA26ABA5e0F4";
1317 str >> CXSC_sqrt7[14];
1318 str = "+1059211B7D5398e0BD";
1319 str >> CXSC_sqrt7[15];
1320 str = "-10F15BFA46EB7Fe087";
1321 str >> CXSC_sqrt7[16];
1322 str = "+15AB71566CE72Be051";
1323 str >> CXSC_sqrt7[17];
1324 str = "-1386BDCA3845C7e01A";
1325 str >> CXSC_sqrt7[18];
1326 str = "+10000000AC4BC7e000";
1327 str >> CXSC_sqrt7[19];
1328 str = "+10000000AC4BC8e000";
1329 str >> CXSC_sqrt7[20];
1330
1331 CXSC_sqrt7_initialized = true;
1332 std::cout << RestoreOpt;
1333 }
1334 stagprec = stagmax;
1335 y = adjust(l_interval(0.0));
1336 for (int i=0; i <= stagmax; i++)
1337 y.data[i] = CXSC_sqrt7[i];
1338 stagprec = stagsave;
1339 y = adjust(y);
1340 return y;
1341}
1342
1343
1344// ************************************************************************
1345
1346static real CXSC_ln2r[21]; // CXSC_ln2r[0], ... CXSC_ln2r[20]
1347static bool CXSC_ln2r_initialized = false;
1348
1350// Inclusion of 1/ln(2), Blomquist, 12.09.06;
1351{
1352 l_interval y;
1353 int stagsave = stagprec,
1354 stagmax = 20;
1355
1356 if (!CXSC_ln2r_initialized)
1357 {
1358 std::string str;
1359 std::cout << SaveOpt;
1360 std::cout << Hex;
1361 str = "+171547652B82FEe3FF";
1362 str >> CXSC_ln2r[0];
1363 str = "+1777D0FFDA0D24e3C7";
1364 str >> CXSC_ln2r[1];
1365 str = "-160BB8A5442AB9e391";
1366 str >> CXSC_ln2r[2];
1367 str = "-14B52D3BA6D74De359";
1368 str >> CXSC_ln2r[3];
1369 str = "+19A342648FBC39e323";
1370 str >> CXSC_ln2r[4];
1371 str = "-1E0455744994EEe2ED";
1372 str >> CXSC_ln2r[5];
1373 str = "+1B25EEB82D7C16e2B7";
1374 str >> CXSC_ln2r[6];
1375 str = "+1F5485CF306255e281";
1376 str >> CXSC_ln2r[7];
1377 str = "-1EC07680A1F958e24B";
1378 str >> CXSC_ln2r[8];
1379 str = "-106326680EB5B6e215";
1380 str >> CXSC_ln2r[9];
1381 str = "-1B3D04C549BC98e1DF";
1382 str >> CXSC_ln2r[10];
1383 str = "+1EABCEAD10305Be1A9";
1384 str >> CXSC_ln2r[11];
1385 str = "-14440C57D7AB97e170";
1386 str >> CXSC_ln2r[12];
1387 str = "-17185D42A4E6D6e139";
1388 str >> CXSC_ln2r[13];
1389 str = "-1F332B5BE48526e101";
1390 str >> CXSC_ln2r[14];
1391 str = "+12CE4F199E108De0CB";
1392 str >> CXSC_ln2r[15];
1393 str = "-18DAFCC6077F2Ae092";
1394 str >> CXSC_ln2r[16];
1395 str = "+19ABB71EC25E12e05B";
1396 str >> CXSC_ln2r[17];
1397 str = "-11473D7A3366BDe022";
1398 str >> CXSC_ln2r[18];
1399 str = "-1000004977D38Be000";
1400 str >> CXSC_ln2r[19];
1401 str = "-1000004977D38Ae000";
1402 str >> CXSC_ln2r[20];
1403
1404 CXSC_ln2r_initialized = true;
1405 std::cout << RestoreOpt;
1406 }
1407 stagprec = stagmax;
1408 y = adjust(l_interval(0.0));
1409 for (int i=0; i <= stagmax; i++)
1410 y.data[i] = CXSC_ln2r[i];
1411// y[i+1] = CXSC_ln2r[i];
1412 stagprec = stagsave;
1413 y = adjust(y);
1414 return y;
1415}
1416
1417
1418// **********************************************************************
1419
1420static real CXSC_Pi[21]; // CXSC_Pi[0], ... CXSC_Pi[20]
1421static bool CXSC_Pi_initialized = false;
1422
1424// Inclusion of Pi, Blomquist, 12.09.06;
1425{
1426 l_interval y;
1427 int stagsave = stagprec,
1428 stagmax = 20;
1429
1430 if (!CXSC_Pi_initialized)
1431 {
1432 std::string str;
1433 std::cout << SaveOpt;
1434 std::cout << Hex;
1435 str = "+1921FB54442D18e400";
1436 str >> CXSC_Pi[0];
1437 str = "+11A62633145C07e3CA";
1438 str >> CXSC_Pi[1];
1439 str = "-1F1976B7ED8FBCe392";
1440 str >> CXSC_Pi[2];
1441 str = "+14CF98E804177De35C";
1442 str >> CXSC_Pi[3];
1443 str = "+131D89CD9128A5e326";
1444 str >> CXSC_Pi[4];
1445 str = "+10F31C6809BBDFe2EC";
1446 str >> CXSC_Pi[5];
1447 str = "+1519B3CD3A431Be2B5";
1448 str >> CXSC_Pi[6];
1449 str = "+18158536F92F8Ae27E";
1450 str >> CXSC_Pi[7];
1451 str = "+1BA7F09AB6B6A9e246";
1452 str >> CXSC_Pi[8];
1453 str = "-1EDD0DBD2544CFe20E";
1454 str >> CXSC_Pi[9];
1455 str = "+179FB1BD1310BAe1D7";
1456 str >> CXSC_Pi[10];
1457 str = "+1A637ED6B0BFF6e1A1";
1458 str >> CXSC_Pi[11];
1459 str = "-1A485FCA40908Ee16A";
1460 str >> CXSC_Pi[12];
1461 str = "-1E501295D98169e133";
1462 str >> CXSC_Pi[13];
1463 str = "-1160DBEE83B4E0e0FD";
1464 str >> CXSC_Pi[14];
1465 str = "-19B6D799AE131Ce0C5";
1466 str >> CXSC_Pi[15];
1467 str = "+16CF70801F2E28e08F";
1468 str >> CXSC_Pi[16];
1469 str = "+163BF0598DA483e059";
1470 str >> CXSC_Pi[17];
1471 str = "+1871574E69A459e023";
1472 str >> CXSC_Pi[18];
1473 str = "-10000005702DB4e000";
1474 str >> CXSC_Pi[19];
1475 str = "-10000005702DB3e000";
1476 str >> CXSC_Pi[20];
1477
1478 CXSC_Pi_initialized = true;
1479 std::cout << RestoreOpt;
1480 }
1481 stagprec = stagmax;
1482 y = adjust(l_interval(0.0));
1483 for (int i=0; i <= stagmax; i++)
1484 y.data[i] = CXSC_Pi[i];
1485// y[i+1] = CXSC_Pi[i];
1486 stagprec = stagsave;
1487 y = adjust(y);
1488 return y;
1489}
1490
1491// **********************************************************************
1492
1493static real CXSC_Pid2[21]; // CXSC_Pid2[0], ... CXSC_Pid2[20]
1494static bool CXSC_Pid2_initialized = false;
1495
1497// Inclusion of Pi/2, Blomquist, 12.09.06;
1498{
1499 l_interval y;
1500 int stagsave = stagprec,
1501 stagmax = 20;
1502
1503 if (!CXSC_Pid2_initialized)
1504 {
1505 std::string str;
1506 std::cout << SaveOpt;
1507 std::cout << Hex;
1508 str = "+1921FB54442D18e3FF";
1509 str >> CXSC_Pid2[0];
1510 str = "+11A62633145C07e3C9";
1511 str >> CXSC_Pid2[1];
1512 str = "-1F1976B7ED8FBCe391";
1513 str >> CXSC_Pid2[2];
1514 str = "+14CF98E804177De35B";
1515 str >> CXSC_Pid2[3];
1516 str = "+131D89CD9128A5e325";
1517 str >> CXSC_Pid2[4];
1518 str = "+10F31C6809BBDFe2EB";
1519 str >> CXSC_Pid2[5];
1520 str = "+1519B3CD3A431Be2B4";
1521 str >> CXSC_Pid2[6];
1522 str = "+18158536F92F8Ae27D";
1523 str >> CXSC_Pid2[7];
1524 str = "+1BA7F09AB6B6A9e245";
1525 str >> CXSC_Pid2[8];
1526 str = "-1EDD0DBD2544CFe20D";
1527 str >> CXSC_Pid2[9];
1528 str = "+179FB1BD1310BAe1D6";
1529 str >> CXSC_Pid2[10];
1530 str = "+1A637ED6B0BFF6e1A0";
1531 str >> CXSC_Pid2[11];
1532 str = "-1A485FCA40908Ee169";
1533 str >> CXSC_Pid2[12];
1534 str = "-1E501295D98169e132";
1535 str >> CXSC_Pid2[13];
1536 str = "-1160DBEE83B4E0e0FC";
1537 str >> CXSC_Pid2[14];
1538 str = "-19B6D799AE131Ce0C4";
1539 str >> CXSC_Pid2[15];
1540 str = "+16CF70801F2E28e08E";
1541 str >> CXSC_Pid2[16];
1542 str = "+163BF0598DA483e058";
1543 str >> CXSC_Pid2[17];
1544 str = "+1871574E69A459e022";
1545 str >> CXSC_Pid2[18];
1546 str = "-10000002B816DAe000";
1547 str >> CXSC_Pid2[19];
1548 str = "-10000002B816D0e000";
1549 str >> CXSC_Pid2[20];
1550
1551 CXSC_Pid2_initialized = true;
1552 std::cout << RestoreOpt;
1553 }
1554 stagprec = stagmax;
1555 y = adjust(l_interval(0.0));
1556 for (int i=0; i <= stagmax; i++)
1557 y.data[i] = CXSC_Pid2[i];
1558// y[i+1] = CXSC_Pid2[i];
1559 stagprec = stagsave;
1560 y = adjust(y);
1561 return y;
1562}
1563
1564// **********************************************************************
1565
1566static real CXSC_Pi2[21]; // CXSC_Pi2[0], ... CXSC_Pi2[20]
1567static bool CXSC_Pi2_initialized = false;
1568
1570// Inclusion of 2*Pi, Blomquist, 12.09.06;
1571{
1572 l_interval y;
1573 int stagsave = stagprec,
1574 stagmax = 20;
1575
1576 if (!CXSC_Pi2_initialized)
1577 {
1578 std::string str;
1579 std::cout << SaveOpt;
1580 std::cout << Hex;
1581 str = "+1921FB54442D18e401";
1582 str >> CXSC_Pi2[0];
1583 str = "+11A62633145C07e3CB";
1584 str >> CXSC_Pi2[1];
1585 str = "-1F1976B7ED8FBCe393";
1586 str >> CXSC_Pi2[2];
1587 str = "+14CF98E804177De35D";
1588 str >> CXSC_Pi2[3];
1589 str = "+131D89CD9128A5e327";
1590 str >> CXSC_Pi2[4];
1591 str = "+10F31C6809BBDFe2ED";
1592 str >> CXSC_Pi2[5];
1593 str = "+1519B3CD3A431Be2B6";
1594 str >> CXSC_Pi2[6];
1595 str = "+18158536F92F8Ae27F";
1596 str >> CXSC_Pi2[7];
1597 str = "+1BA7F09AB6B6A9e247";
1598 str >> CXSC_Pi2[8];
1599 str = "-1EDD0DBD2544CFe20F";
1600 str >> CXSC_Pi2[9];
1601 str = "+179FB1BD1310BAe1D8";
1602 str >> CXSC_Pi2[10];
1603 str = "+1A637ED6B0BFF6e1A2";
1604 str >> CXSC_Pi2[11];
1605 str = "-1A485FCA40908Ee16B";
1606 str >> CXSC_Pi2[12];
1607 str = "-1E501295D98169e134";
1608 str >> CXSC_Pi2[13];
1609 str = "-1160DBEE83B4E0e0FE";
1610 str >> CXSC_Pi2[14];
1611 str = "-19B6D799AE131Ce0C6";
1612 str >> CXSC_Pi2[15];
1613 str = "+16CF70801F2E28e090";
1614 str >> CXSC_Pi2[16];
1615 str = "+163BF0598DA483e05A";
1616 str >> CXSC_Pi2[17];
1617 str = "+1871574E69A459e024";
1618 str >> CXSC_Pi2[18];
1619 str = "-1000000AE05B67e000";
1620 str >> CXSC_Pi2[19];
1621 str = "-1000000AE05B66e000";
1622 str >> CXSC_Pi2[20];
1623
1624 CXSC_Pi2_initialized = true;
1625 std::cout << RestoreOpt;
1626 }
1627 stagprec = stagmax;
1628 y = adjust(l_interval(0.0));
1629 for (int i=0; i <= stagmax; i++)
1630 y.data[i] = CXSC_Pi2[i];
1631// y[i+1] = CXSC_Pi2[i];
1632 stagprec = stagsave;
1633 y = adjust(y);
1634 return y;
1635}
1636
1637// **********************************************************************
1638
1639static real CXSC_Pid3[21]; // CXSC_Pid3[0], ... CXSC_Pid3[20]
1640static bool CXSC_Pid3_initialized = false;
1641
1643// Inclusion of Pi/3, Blomquist, 12.09.06;
1644{
1645 l_interval y;
1646 int stagsave = stagprec,
1647 stagmax = 20;
1648
1649 if (!CXSC_Pid3_initialized)
1650 {
1651 std::string str;
1652 std::cout << SaveOpt;
1653 std::cout << Hex;
1654 str = "+10C152382D7366e3FF";
1655 str >> CXSC_Pid3[0];
1656 str = "-1EE6913347C2A6e3C9";
1657 str >> CXSC_Pid3[1];
1658 str = "-14BBA47A9E5FD2e391";
1659 str >> CXSC_Pid3[2];
1660 str = "-1CCAEF65529B02e35B";
1661 str >> CXSC_Pid3[3];
1662 str = "+197CB7BCC18B87e324";
1663 str >> CXSC_Pid3[4];
1664 str = "-13EBBDA1FF3058e2EE";
1665 str >> CXSC_Pid3[5];
1666 str = "-11D10CB320F4D1e2B6";
1667 str >> CXSC_Pid3[6];
1668 str = "+1958EB892987ECe27F";
1669 str >> CXSC_Pid3[7];
1670 str = "+167C54B11CF247e249";
1671 str >> CXSC_Pid3[8];
1672 str = "+12C2E985923A44e210";
1673 str >> CXSC_Pid3[9];
1674 str = "+1945484A2DD81Fe1D8";
1675 str >> CXSC_Pid3[10];
1676 str = "+1197A9E475D54Fe1A0";
1677 str >> CXSC_Pid3[11];
1678 str = "-1E181FEE158585e16A";
1679 str >> CXSC_Pid3[12];
1680 str = "+1047FCE7066A6Ee134";
1681 str >> CXSC_Pid3[13];
1682 str = "+1D1A8602EA0C85e0FE";
1683 str >> CXSC_Pid3[14];
1684 str = "+14430C5998BF34e0C8";
1685 str >> CXSC_Pid3[15];
1686 str = "+173BF40AAD43D9e091";
1687 str >> CXSC_Pid3[16];
1688 str = "-137B014DDEDCF5e05B";
1689 str >> CXSC_Pid3[17];
1690 str = "-1A5F1B210EE7C5e022";
1691 str >> CXSC_Pid3[18];
1692 str = "+100000A8DA9B6Ee000";
1693 str >> CXSC_Pid3[19];
1694 str = "+100000A8DA9B6Fe000";
1695 str >> CXSC_Pid3[20];
1696
1697 CXSC_Pid3_initialized = true;
1698 std::cout << RestoreOpt;
1699 }
1700 stagprec = stagmax;
1701 y = adjust(l_interval(0.0));
1702 for (int i=0; i <= stagmax; i++)
1703 y.data[i] = CXSC_Pid3[i];
1704// y[i+1] = CXSC_Pid3[i];
1705 stagprec = stagsave;
1706 y = adjust(y);
1707 return y;
1708}
1709
1710
1711// **********************************************************************
1712
1713static real CXSC_Pir[21]; // CXSC_Pir[0], ... CXSC_Pir[20]
1714static bool CXSC_Pir_initialized = false;
1715
1717// Inclusion of 1/Pi, Blomquist, 12.09.06;
1718{
1719 l_interval y;
1720 int stagsave = stagprec,
1721 stagmax = 20;
1722
1723 if (!CXSC_Pir_initialized)
1724 {
1725 std::string str;
1726 std::cout << SaveOpt;
1727 std::cout << Hex;
1728 str = "+145F306DC9C883e3FD";
1729 str >> CXSC_Pir[0];
1730 str = "-16B01EC5417056e3C7";
1731 str >> CXSC_Pir[1];
1732 str = "-16447E493AD4CEe391";
1733 str >> CXSC_Pir[2];
1734 str = "+1E21C820FF28B2e35B";
1735 str >> CXSC_Pir[3];
1736 str = "-1508510EA79237e324";
1737 str >> CXSC_Pir[4];
1738 str = "+1B8E909374B802e2EC";
1739 str >> CXSC_Pir[5];
1740 str = "-1B6D115F62E6DEe2B6";
1741 str >> CXSC_Pir[6];
1742 str = "-180F10A71A76B3e27F";
1743 str >> CXSC_Pir[7];
1744 str = "+1CFBA208D7D4BBe248";
1745 str >> CXSC_Pir[8];
1746 str = "-12EDEC598E3F65e210";
1747 str >> CXSC_Pir[9];
1748 str = "-1741037D8CDC54e1D9";
1749 str >> CXSC_Pir[10];
1750 str = "+1CC1A99CFA4E42e1A3";
1751 str >> CXSC_Pir[11];
1752 str = "+17E2EF7E4A0EC8e16C";
1753 str >> CXSC_Pir[12];
1754 str = "-1DA00087E99FC0e130";
1755 str >> CXSC_Pir[13];
1756 str = "-10D0EE74A5F593e0FA";
1757 str >> CXSC_Pir[14];
1758 str = "+1F6D367ECF27CBe0C2";
1759 str >> CXSC_Pir[15];
1760 str = "+136E9E8C7ECD3De089";
1761 str >> CXSC_Pir[16];
1762 str = "-100AE9456C229Ce053";
1763 str >> CXSC_Pir[17];
1764 str = "-141A0E84C2F8C6e01A";
1765 str >> CXSC_Pir[18];
1766 str = "-1000000010EB5Be000";
1767 str >> CXSC_Pir[19];
1768 str = "-1000000010EB5Ae000";
1769 str >> CXSC_Pir[20];
1770
1771 CXSC_Pir_initialized = true;
1772 std::cout << RestoreOpt;
1773 }
1774 stagprec = stagmax;
1775 y = adjust(l_interval(0.0));
1776 for (int i=0; i <= stagmax; i++)
1777 y.data[i] = CXSC_Pir[i];
1778// y[i+1] = CXSC_Pir[i];
1779 stagprec = stagsave;
1780 y = adjust(y);
1781 return y;
1782}
1783
1784// **********************************************************************
1785
1786static real CXSC_Pi2r[21]; // CXSC_Pi2r[0], ... CXSC_Pi2r[20]
1787static bool CXSC_Pi2r_initialized = false;
1788
1790// Inclusion of 1/(2*Pi), Blomquist, 12.09.06;
1791{
1792 l_interval y;
1793 int stagsave = stagprec,
1794 stagmax = 20;
1795
1796 if (!CXSC_Pi2r_initialized)
1797 {
1798 std::string str;
1799 std::cout << SaveOpt;
1800 std::cout << Hex;
1801 str = "+145F306DC9C883e3FC";
1802 str >> CXSC_Pi2r[0];
1803 str = "-16B01EC5417056e3C6";
1804 str >> CXSC_Pi2r[1];
1805 str = "-16447E493AD4CEe390";
1806 str >> CXSC_Pi2r[2];
1807 str = "+1E21C820FF28B2e35A";
1808 str >> CXSC_Pi2r[3];
1809 str = "-1508510EA79237e323";
1810 str >> CXSC_Pi2r[4];
1811 str = "+1B8E909374B802e2EB";
1812 str >> CXSC_Pi2r[5];
1813 str = "-1B6D115F62E6DEe2B5";
1814 str >> CXSC_Pi2r[6];
1815 str = "-180F10A71A76B3e27E";
1816 str >> CXSC_Pi2r[7];
1817 str = "+1CFBA208D7D4BBe247";
1818 str >> CXSC_Pi2r[8];
1819 str = "-12EDEC598E3F65e20F";
1820 str >> CXSC_Pi2r[9];
1821 str = "-1741037D8CDC54e1D8";
1822 str >> CXSC_Pi2r[10];
1823 str = "+1CC1A99CFA4E42e1A2";
1824 str >> CXSC_Pi2r[11];
1825 str = "+17E2EF7E4A0EC8e16B";
1826 str >> CXSC_Pi2r[12];
1827 str = "-1DA00087E99FC0e12F";
1828 str >> CXSC_Pi2r[13];
1829 str = "-10D0EE74A5F593e0F9";
1830 str >> CXSC_Pi2r[14];
1831 str = "+1F6D367ECF27CBe0C1";
1832 str >> CXSC_Pi2r[15];
1833 str = "+136E9E8C7ECD3De088";
1834 str >> CXSC_Pi2r[16];
1835 str = "-100AE9456C229Ce052";
1836 str >> CXSC_Pi2r[17];
1837 str = "-141A0E84C2F8C6e019";
1838 str >> CXSC_Pi2r[18];
1839 str = "-100000000875AEe000";
1840 str >> CXSC_Pi2r[19];
1841 str = "-100000000875ADe000";
1842 str >> CXSC_Pi2r[20];
1843
1844 CXSC_Pi2r_initialized = true;
1845 std::cout << RestoreOpt;
1846 }
1847 stagprec = stagmax;
1848 y = adjust(l_interval(0.0));
1849 for (int i=0; i <= stagmax; i++)
1850 y.data[i] = CXSC_Pi2r[i];
1851// y[i+1] = CXSC_Pi2r[i];
1852 stagprec = stagsave;
1853 y = adjust(y);
1854 return y;
1855}
1856
1857
1858// **********************************************************************
1859
1860static real CXSC_SqrtPi[21]; // CXSC_SqrtPi[0], ... CXSC_SqrtPi[20]
1861static bool CXSC_SqrtPi_initialized = false;
1862
1864// Inclusion of Sqrt(Pi), Blomquist, 12.09.06;
1865{
1866 l_interval y;
1867 int stagsave = stagprec,
1868 stagmax = 20;
1869
1870 if (!CXSC_SqrtPi_initialized)
1871 {
1872 std::string str;
1873 std::cout << SaveOpt;
1874 std::cout << Hex;
1875 str = "+1C5BF891B4EF6Be3FF";
1876 str >> CXSC_SqrtPi[0];
1877 str = "-1618F13EB7CA89e3C9";
1878 str >> CXSC_SqrtPi[1];
1879 str = "-1B1F0071B7AAE4e391";
1880 str >> CXSC_SqrtPi[2];
1881 str = "-1389B5A46BDFE8e35A";
1882 str >> CXSC_SqrtPi[3];
1883 str = "-160AF5C5C89448e324";
1884 str >> CXSC_SqrtPi[4];
1885 str = "-14835F07122994e2E8";
1886 str >> CXSC_SqrtPi[5];
1887 str = "+1CEC283C18EE8Fe2B2";
1888 str >> CXSC_SqrtPi[6];
1889 str = "-13ADEBB9223CA8e27B";
1890 str >> CXSC_SqrtPi[7];
1891 str = "+1454912430D291e245";
1892 str >> CXSC_SqrtPi[8];
1893 str = "-1E8B2345020EF6e20F";
1894 str >> CXSC_SqrtPi[9];
1895 str = "-17262982556291e1D8";
1896 str >> CXSC_SqrtPi[10];
1897 str = "+1196FA9B140CABe1A1";
1898 str >> CXSC_SqrtPi[11];
1899 str = "-175EEE59D91D39e16B";
1900 str >> CXSC_SqrtPi[12];
1901 str = "+1789268B7D9D48e130";
1902 str >> CXSC_SqrtPi[13];
1903 str = "+17162E2F06B89Ce0FA";
1904 str >> CXSC_SqrtPi[14];
1905 str = "+1EC9C08F40A3DBe0C3";
1906 str >> CXSC_SqrtPi[15];
1907 str = "+1B6048DD0729E2e08D";
1908 str >> CXSC_SqrtPi[16];
1909 str = "+1471CF4C33FF6Be056";
1910 str >> CXSC_SqrtPi[17];
1911 str = "+1D75FBD8B36F94e020";
1912 str >> CXSC_SqrtPi[18];
1913 str = "+1000002D74B3A2e000";
1914 str >> CXSC_SqrtPi[19];
1915 str = "+1000002D74B3A3e000";
1916 str >> CXSC_SqrtPi[20];
1917
1918 CXSC_SqrtPi_initialized = true;
1919 std::cout << RestoreOpt;
1920 }
1921 stagprec = stagmax;
1922 y = adjust(l_interval(0.0));
1923 for (int i=0; i <= stagmax; i++)
1924 y.data[i] = CXSC_SqrtPi[i];
1925// y[i+1] = CXSC_SqrtPi[i];
1926 stagprec = stagsave;
1927 y = adjust(y);
1928 return y;
1929}
1930
1931
1932// **********************************************************************
1933
1934static real CXSC_Sqrt2Pi[21]; // CXSC_Sqrt2Pi[0], ... CXSC_Sqrt2Pi[20]
1935static bool CXSC_Sqrt2Pi_initialized = false;
1936
1938// Inclusion of sqrt(2*Pi), Blomquist, 12.09.06;
1939{
1940 l_interval y;
1941 int stagsave = stagprec,
1942 stagmax = 20;
1943
1944 if (!CXSC_Sqrt2Pi_initialized)
1945 {
1946 std::string str;
1947 std::cout << SaveOpt;
1948 std::cout << Hex;
1949 str = "+140D931FF62706e400";
1950 str >> CXSC_Sqrt2Pi[0];
1951 str = "-1A6A0D6F814637e3CA";
1952 str >> CXSC_Sqrt2Pi[1];
1953 str = "-1311D073060ACEe394";
1954 str >> CXSC_Sqrt2Pi[2];
1955 str = "+16000B50DC2F41e35B";
1956 str >> CXSC_Sqrt2Pi[3];
1957 str = "+16EF75CA45A834e324";
1958 str >> CXSC_Sqrt2Pi[4];
1959 str = "+19BDB2B4C39342e2EC";
1960 str >> CXSC_Sqrt2Pi[5];
1961 str = "+1F5582E2063EE6e2B5";
1962 str >> CXSC_Sqrt2Pi[6];
1963 str = "+183F879BEA150Ce27C";
1964 str >> CXSC_Sqrt2Pi[7];
1965 str = "-1F1EA3CA289B00e244";
1966 str >> CXSC_Sqrt2Pi[8];
1967 str = "-1699CDA77736F9e20D";
1968 str >> CXSC_Sqrt2Pi[9];
1969 str = "-11A379D298B55Ee1D4";
1970 str >> CXSC_Sqrt2Pi[10];
1971 str = "-1A6DDB0152BA94e19E";
1972 str >> CXSC_Sqrt2Pi[11];
1973 str = "-1957E2E58A02FEe167";
1974 str >> CXSC_Sqrt2Pi[12];
1975 str = "-1D6160F18E604De131";
1976 str >> CXSC_Sqrt2Pi[13];
1977 str = "+1311860CDF7215e0F8";
1978 str >> CXSC_Sqrt2Pi[14];
1979 str = "+12271F44C50274e0C1";
1980 str >> CXSC_Sqrt2Pi[15];
1981 str = "-100BF5C5497A21e08A";
1982 str >> CXSC_Sqrt2Pi[16];
1983 str = "+1E94B6E6AD51E2e052";
1984 str >> CXSC_Sqrt2Pi[17];
1985 str = "-1C910B5F3D27CEe019";
1986 str >> CXSC_Sqrt2Pi[18];
1987 str = "+100000007C99B0e000";
1988 str >> CXSC_Sqrt2Pi[19];
1989 str = "+100000007C99B1e000";
1990 str >> CXSC_Sqrt2Pi[20];
1991
1992 CXSC_Sqrt2Pi_initialized = true;
1993 std::cout << RestoreOpt;
1994 }
1995 stagprec = stagmax;
1996 y = adjust(l_interval(0.0));
1997 for (int i=0; i <= stagmax; i++)
1998 y.data[i] = CXSC_Sqrt2Pi[i];
1999// y[i+1] = CXSC_Sqrt2Pi[i];
2000 stagprec = stagsave;
2001 y = adjust(y);
2002 return y;
2003}
2004
2005
2006// **********************************************************************
2007
2008static real CXSC_SqrtPir[21]; // CXSC_SqrtPir[0], ... CXSC_SqrtPir[20]
2009static bool CXSC_SqrtPir_initialized = false;
2010
2012// Inclusion of 1/sqrt(Pi), Blomquist, 12.09.06;
2013{
2014 l_interval y;
2015 int stagsave = stagprec,
2016 stagmax = 20;
2017
2018 if (!CXSC_SqrtPir_initialized)
2019 {
2020 std::string str;
2021 std::cout << SaveOpt;
2022 std::cout << Hex;
2023 str = "+120DD750429B6De3FE";
2024 str >> CXSC_SqrtPir[0];
2025 str = "+11AE3A914FED80e3C6";
2026 str >> CXSC_SqrtPir[1];
2027 str = "-13CBBEBF65F145e38F";
2028 str >> CXSC_SqrtPir[2];
2029 str = "-1E0C574632F53Ee358";
2030 str >> CXSC_SqrtPir[3];
2031 str = "-1E6633BE9E7F15e322";
2032 str >> CXSC_SqrtPir[4];
2033 str = "+1CF859270F1141e2EB";
2034 str >> CXSC_SqrtPir[5];
2035 str = "-1FE4FB499C328Ae2B4";
2036 str >> CXSC_SqrtPir[6];
2037 str = "-10B82C446DC78De27D";
2038 str >> CXSC_SqrtPir[7];
2039 str = "-1878B089078800e247";
2040 str >> CXSC_SqrtPir[8];
2041 str = "-13DAEADA9E233Ee20F";
2042 str >> CXSC_SqrtPir[9];
2043 str = "+1137197A708BD2e1D9";
2044 str >> CXSC_SqrtPir[10];
2045 str = "-109009506D5BA2e19E";
2046 str >> CXSC_SqrtPir[11];
2047 str = "+17C9F0B5951E94e168";
2048 str >> CXSC_SqrtPir[12];
2049 str = "-1735F4949633A4e131";
2050 str >> CXSC_SqrtPir[13];
2051 str = "-146014DBC90D0Ee0FB";
2052 str >> CXSC_SqrtPir[14];
2053 str = "+1CAB0B222EEEA0e0C5";
2054 str >> CXSC_SqrtPir[15];
2055 str = "+1B1C750754B40Ae08F";
2056 str >> CXSC_SqrtPir[16];
2057 str = "-16B2CD2E72C16Ee057";
2058 str >> CXSC_SqrtPir[17];
2059 str = "-148C024FF194B2e021";
2060 str >> CXSC_SqrtPir[18];
2061 str = "+10000073E19B74e000";
2062 str >> CXSC_SqrtPir[19];
2063 str = "+10000073E19B75e000";
2064 str >> CXSC_SqrtPir[20];
2065
2066 CXSC_SqrtPir_initialized = true;
2067 std::cout << RestoreOpt;
2068 }
2069 stagprec = stagmax;
2070 y = adjust(l_interval(0.0));
2071 for (int i=0; i <= stagmax; i++)
2072 y.data[i] = CXSC_SqrtPir[i];
2073// y[i+1] = CXSC_SqrtPir[i];
2074 stagprec = stagsave;
2075 y = adjust(y);
2076 return y;
2077}
2078
2079
2080// **********************************************************************
2081
2082static real CXSC_Sqrt2Pir[21]; // CXSC_Sqrt2Pir[0], ... CXSC_Sqrt2Pir[20]
2083static bool CXSC_Sqrt2Pir_initialized = false;
2084
2086// Inclusion of 1/sqrt(2*Pi), Blomquist, 12.09.06;
2087{
2088 l_interval y;
2089 int stagsave = stagprec,
2090 stagmax = 20;
2091
2092 if (!CXSC_Sqrt2Pir_initialized)
2093 {
2094 std::string str;
2095 std::cout << SaveOpt;
2096 std::cout << Hex;
2097 str = "+19884533D43651e3FD";
2098 str >> CXSC_Sqrt2Pir[0];
2099 str = "-1CBC0D30EBFD15e3C7";
2100 str >> CXSC_Sqrt2Pir[1];
2101 str = "-1C7402C7D60CFBe38F";
2102 str >> CXSC_Sqrt2Pir[2];
2103 str = "+12706D8C0471B5e357";
2104 str >> CXSC_Sqrt2Pir[3];
2105 str = "-1FF6718B45881De321";
2106 str >> CXSC_Sqrt2Pir[4];
2107 str = "-13AABB82C248DCe2EB";
2108 str >> CXSC_Sqrt2Pir[5];
2109 str = "-1458A899162EE4e2B2";
2110 str >> CXSC_Sqrt2Pir[6];
2111 str = "-14EBD8868F41EBe27B";
2112 str >> CXSC_Sqrt2Pir[7];
2113 str = "+13278E993445F1e243";
2114 str >> CXSC_Sqrt2Pir[8];
2115 str = "-1CC019F5F4780Ae20D";
2116 str >> CXSC_Sqrt2Pir[9];
2117 str = "+147CE4B4ECDBD7e1D7";
2118 str >> CXSC_Sqrt2Pir[10];
2119 str = "-19A3DCC6A3534Be19F";
2120 str >> CXSC_Sqrt2Pir[11];
2121 str = "+11379A7BA8CB0Ae169";
2122 str >> CXSC_Sqrt2Pir[12];
2123 str = "-12D909C875312Ee132";
2124 str >> CXSC_Sqrt2Pir[13];
2125 str = "+1C1CEC4882C77Be0FB";
2126 str >> CXSC_Sqrt2Pir[14];
2127 str = "-14C4078263DF36e0C5";
2128 str >> CXSC_Sqrt2Pir[15];
2129 str = "+1AB3FC8D2AB243e08F";
2130 str >> CXSC_Sqrt2Pir[16];
2131 str = "+17B9172454310Ae059";
2132 str >> CXSC_Sqrt2Pir[17];
2133 str = "-1444B6B781B7F2e023";
2134 str >> CXSC_Sqrt2Pir[18];
2135 str = "-100001DB5C6774e000";
2136 str >> CXSC_Sqrt2Pir[19];
2137 str = "-100001DB5C6773e000";
2138 str >> CXSC_Sqrt2Pir[20];
2139
2140 CXSC_Sqrt2Pir_initialized = true;
2141 std::cout << RestoreOpt;
2142 }
2143 stagprec = stagmax;
2144 y = adjust(l_interval(0.0));
2145 for (int i=0; i <= stagmax; i++)
2146 y.data[i] = CXSC_Sqrt2Pir[i];
2147// y[i+1] = CXSC_Sqrt2Pir[i];
2148 stagprec = stagsave;
2149 y = adjust(y);
2150 return y;
2151}
2152
2153
2154// **********************************************************************
2155
2156static real CXSC_Pip2[21]; // CXSC_Pip2[0], ... CXSC_Pip2[20]
2157static bool CXSC_Pip2_initialized = false;
2158
2160// Inclusion of Pi^2, Blomquist, 12.09.06;
2161{
2162 l_interval y;
2163 int stagsave = stagprec,
2164 stagmax = 20;
2165
2166 if (!CXSC_Pip2_initialized)
2167 {
2168 std::string str;
2169 std::cout << SaveOpt;
2170 std::cout << Hex;
2171 str = "+13BD3CC9BE45DEe402";
2172 str >> CXSC_Pip2[0];
2173 str = "+1692B71366CC04e3CC";
2174 str >> CXSC_Pip2[1];
2175 str = "+18358E10ACD480e396";
2176 str >> CXSC_Pip2[2];
2177 str = "-1F2F5DD7997DDFe35F";
2178 str >> CXSC_Pip2[3];
2179 str = "+129E39B47B884Ee324";
2180 str >> CXSC_Pip2[4];
2181 str = "-12CF7459DD5DAFe2EE";
2182 str >> CXSC_Pip2[5];
2183 str = "-11842F87B5FE0Fe2B8";
2184 str >> CXSC_Pip2[6];
2185 str = "+1FFD8A79616A21e282";
2186 str >> CXSC_Pip2[7];
2187 str = "+12492A6663E899e24C";
2188 str >> CXSC_Pip2[8];
2189 str = "-1A15F4352CC511e215";
2190 str >> CXSC_Pip2[9];
2191 str = "-1301AA1792FF3Ce1DE";
2192 str >> CXSC_Pip2[10];
2193 str = "+122B6F31626EFEe1A8";
2194 str >> CXSC_Pip2[11];
2195 str = "+1B317FA13BDD8Fe172";
2196 str >> CXSC_Pip2[12];
2197 str = "+16F83B49040075e13C";
2198 str >> CXSC_Pip2[13];
2199 str = "-1B1890A945FE17e106";
2200 str >> CXSC_Pip2[14];
2201 str = "+12DCD389B96CDBe0D0";
2202 str >> CXSC_Pip2[15];
2203 str = "-1743F5DDE2F157e097";
2204 str >> CXSC_Pip2[16];
2205 str = "-153F96FFD4AEB5e060";
2206 str >> CXSC_Pip2[17];
2207 str = "+13CD6F5847D569e028";
2208 str >> CXSC_Pip2[18];
2209 str = "+10001471E79A7Be000";
2210 str >> CXSC_Pip2[19];
2211 str = "+10001471E79A8Be000";
2212 str >> CXSC_Pip2[20];
2213
2214 CXSC_Pip2_initialized = true;
2215 std::cout << RestoreOpt;
2216 }
2217 stagprec = stagmax;
2218 y = adjust(l_interval(0.0));
2219 for (int i=0; i <= stagmax; i++)
2220 y.data[i] = CXSC_Pip2[i];
2221// y[i+1] = CXSC_Pip2[i];
2222 stagprec = stagsave;
2223 y = adjust(y);
2224 return y;
2225}
2226
2227
2228// **********************************************************************
2229
2230static real CXSC_Sqrt2r[21]; // CXSC_Sqrt2r[0], ... CXSC_Sqrt2r[20]
2231static bool CXSC_Sqrt2r_initialized = false;
2232
2234// Inclusion of 1/sqrt(2), Blomquist, 12.09.06;
2235{
2236 l_interval y;
2237 int stagsave = stagprec,
2238 stagmax = 20;
2239
2240 if (!CXSC_Sqrt2r_initialized)
2241 {
2242 std::string str;
2243 std::cout << SaveOpt;
2244 std::cout << Hex;
2245 str = "+16A09E667F3BCDe3FE";
2246 str >> CXSC_Sqrt2r[0];
2247 str = "-1BDD3413B26456e3C8";
2248 str >> CXSC_Sqrt2r[1];
2249 str = "+157D3E3ADEC175e392";
2250 str >> CXSC_Sqrt2r[2];
2251 str = "+12775099DA2F59e35A";
2252 str >> CXSC_Sqrt2r[3];
2253 str = "+160CCE64552BF2e321";
2254 str >> CXSC_Sqrt2r[4];
2255 str = "+1821D5C5161D46e2E8";
2256 str >> CXSC_Sqrt2r[5];
2257 str = "-1C032046F8498Ee2B2";
2258 str >> CXSC_Sqrt2r[6];
2259 str = "+1EE950BC8738F7e27A";
2260 str >> CXSC_Sqrt2r[7];
2261 str = "-1AC3FDBC64E103e244";
2262 str >> CXSC_Sqrt2r[8];
2263 str = "+13B469101743A1e20C";
2264 str >> CXSC_Sqrt2r[9];
2265 str = "+15E3E9CA60B38Ce1D6";
2266 str >> CXSC_Sqrt2r[10];
2267 str = "+11BC337BCAB1BDe19B";
2268 str >> CXSC_Sqrt2r[11];
2269 str = "-1BBA5DEE9D6E7De165";
2270 str >> CXSC_Sqrt2r[12];
2271 str = "-1438DD083B1CC4e12F";
2272 str >> CXSC_Sqrt2r[13];
2273 str = "+1B56A28E2EDFA7e0F9";
2274 str >> CXSC_Sqrt2r[14];
2275 str = "+1CCB2A634331F4e0C3";
2276 str >> CXSC_Sqrt2r[15];
2277 str = "-1BD9056876F83Ee08C";
2278 str >> CXSC_Sqrt2r[16];
2279 str = "-1234FA22AB6BEFe056";
2280 str >> CXSC_Sqrt2r[17];
2281 str = "+19040CA4A81395e01F";
2282 str >> CXSC_Sqrt2r[18];
2283 str = "-10000015249C0Ce000";
2284 str >> CXSC_Sqrt2r[19];
2285 str = "-10000015249C0Be000";
2286 str >> CXSC_Sqrt2r[20];
2287
2288 CXSC_Sqrt2r_initialized = true;
2289 std::cout << RestoreOpt;
2290 }
2291 stagprec = stagmax;
2292 y = adjust(l_interval(0.0));
2293 for (int i=0; i <= stagmax; i++)
2294 y.data[i] = CXSC_Sqrt2r[i];
2295// y[i+1] = CXSC_Sqrt2r[i];
2296 stagprec = stagsave;
2297 y = adjust(y);
2298 return y;
2299}
2300
2301
2302// **********************************************************************
2303
2304static real CXSC_Sqrt3[21]; // CXSC_Sqrt3[0], ... CXSC_Sqrt3[20]
2305static bool CXSC_Sqrt3_initialized = false;
2306
2308// Inclusion of sqrt(3), Blomquist, 12.09.06;
2309{
2310 l_interval y;
2311 int stagsave = stagprec,
2312 stagmax = 20;
2313
2314 if (!CXSC_Sqrt3_initialized)
2315 {
2316 std::string str;
2317 std::cout << SaveOpt;
2318 std::cout << Hex;
2319 str = "+1BB67AE8584CAAe3FF";
2320 str >> CXSC_Sqrt3[0];
2321 str = "+1CEC95D0B5C1E3e3C9";
2322 str >> CXSC_Sqrt3[1];
2323 str = "-1F11DB689F2CCFe391";
2324 str >> CXSC_Sqrt3[2];
2325 str = "+13DA4798C720A6e35B";
2326 str >> CXSC_Sqrt3[3];
2327 str = "+121B9169B89243e325";
2328 str >> CXSC_Sqrt3[4];
2329 str = "-1813508751212Be2EC";
2330 str >> CXSC_Sqrt3[5];
2331 str = "-1B3D547B775C1Ee2B5";
2332 str >> CXSC_Sqrt3[6];
2333 str = "-19D986D92E2F0Ae27C";
2334 str >> CXSC_Sqrt3[7];
2335 str = "+1A34334CE806B6e245";
2336 str >> CXSC_Sqrt3[8];
2337 str = "+1A383B9E122E61e20F";
2338 str >> CXSC_Sqrt3[9];
2339 str = "+1C61D736F2F6F2e1D8";
2340 str >> CXSC_Sqrt3[10];
2341 str = "-10AF49233F9250e1A1";
2342 str >> CXSC_Sqrt3[11];
2343 str = "-1558A109EC0523e16A";
2344 str >> CXSC_Sqrt3[12];
2345 str = "+1F799D4D4FF2BCe134";
2346 str >> CXSC_Sqrt3[13];
2347 str = "-1AD7B219E34EDBe0FE";
2348 str >> CXSC_Sqrt3[14];
2349 str = "+15AB940B6677E3e0C8";
2350 str >> CXSC_Sqrt3[15];
2351 str = "-1D9B2A8203B8F0e091";
2352 str >> CXSC_Sqrt3[16];
2353 str = "-1DB0C8975A3834e05B";
2354 str >> CXSC_Sqrt3[17];
2355 str = "-1BCAAB3F6BE884e025";
2356 str >> CXSC_Sqrt3[18];
2357 str = "+100000531C2B6Ce000";
2358 str >> CXSC_Sqrt3[19];
2359 str = "+100000531C2B6De000";
2360 str >> CXSC_Sqrt3[20];
2361
2362 CXSC_Sqrt3_initialized = true;
2363 std::cout << RestoreOpt;
2364 }
2365 stagprec = stagmax;
2366 y = adjust(l_interval(0.0));
2367 for (int i=0; i <= stagmax; i++)
2368 y.data[i] = CXSC_Sqrt3[i];
2369// y[i+1] = CXSC_Sqrt3[i];
2370 stagprec = stagsave;
2371 y = adjust(y);
2372 return y;
2373}
2374
2375
2376// **********************************************************************
2377
2378static real CXSC_Sqrt3d2[21]; // CXSC_Sqrt3d2[0], ... CXSC_Sqrt3d2[20]
2379static bool CXSC_Sqrt3d2_initialized = false;
2380
2382// Inclusion of Sqrt(3)/2, Blomquist, 12.09.06;
2383{
2384 l_interval y;
2385 int stagsave = stagprec,
2386 stagmax = 20;
2387
2388 if (!CXSC_Sqrt3d2_initialized)
2389 {
2390 std::string str;
2391 std::cout << SaveOpt;
2392 std::cout << Hex;
2393 str = "+1BB67AE8584CAAe3FE";
2394 str >> CXSC_Sqrt3d2[0];
2395 str = "+1CEC95D0B5C1E3e3C8";
2396 str >> CXSC_Sqrt3d2[1];
2397 str = "-1F11DB689F2CCFe390";
2398 str >> CXSC_Sqrt3d2[2];
2399 str = " +13DA4798C720A6e35A";
2400 str >> CXSC_Sqrt3d2[3];
2401 str = "+121B9169B89243e324";
2402 str >> CXSC_Sqrt3d2[4];
2403 str = " -1813508751212Be2EB";
2404 str >> CXSC_Sqrt3d2[5];
2405 str = "-1B3D547B775C1Ee2B4";
2406 str >> CXSC_Sqrt3d2[6];
2407 str = "-19D986D92E2F0Ae27B";
2408 str >> CXSC_Sqrt3d2[7];
2409 str = "+1A34334CE806B6e244";
2410 str >> CXSC_Sqrt3d2[8];
2411 str = "+1A383B9E122E61e20E";
2412 str >> CXSC_Sqrt3d2[9];
2413 str = "+1C61D736F2F6F2e1D7";
2414 str >> CXSC_Sqrt3d2[10];
2415 str = "-10AF49233F9250e1A0";
2416 str >> CXSC_Sqrt3d2[11];
2417 str = " -1558A109EC0523e169";
2418 str >> CXSC_Sqrt3d2[12];
2419 str = "+1F799D4D4FF2BCe133";
2420 str >> CXSC_Sqrt3d2[13];
2421 str = "-1AD7B219E34EDBe0FD";
2422 str >> CXSC_Sqrt3d2[14];
2423 str = "+15AB940B6677E3e0C7";
2424 str >> CXSC_Sqrt3d2[15];
2425 str = "-1D9B2A8203B8F0e090";
2426 str >> CXSC_Sqrt3d2[16];
2427 str = "-1DB0C8975A3834e05A";
2428 str >> CXSC_Sqrt3d2[17];
2429 str = "-1BCAAB3F6BE884e024";
2430 str >> CXSC_Sqrt3d2[18];
2431 str = "+100000298E15B6e000";
2432 str >> CXSC_Sqrt3d2[19];
2433 str = "+100000298E15B7e000";
2434 str >> CXSC_Sqrt3d2[20];
2435
2436 CXSC_Sqrt3d2_initialized = true;
2437 std::cout << RestoreOpt;
2438 }
2439 stagprec = stagmax;
2440 y = adjust(l_interval(0.0));
2441 for (int i=0; i <= stagmax; i++)
2442 y.data[i] = CXSC_Sqrt3d2[i];
2443// y[i+1] = CXSC_Sqrt3d2[i];
2444 stagprec = stagsave;
2445 y = adjust(y);
2446 return y;
2447}
2448
2449
2450// **********************************************************************
2451
2452static real CXSC_Sqrt3r[21]; // CXSC_Sqrt3r[0], ... CXSC_Sqrt3r[20]
2453static bool CXSC_Sqrt3r_initialized = false;
2454
2456// Inclusion of 1/Sqrt(3), Blomquist, 12.09.06;
2457{
2458 l_interval y;
2459 int stagsave = stagprec,
2460 stagmax = 20;
2461
2462 if (!CXSC_Sqrt3r_initialized)
2463 {
2464 std::string str;
2465 std::cout << SaveOpt;
2466 std::cout << Hex;
2467 str = "+1279A74590331Ce3FE";
2468 str >> CXSC_Sqrt3r[0];
2469 str = "+134863E0792BEDe3C8";
2470 str >> CXSC_Sqrt3r[1];
2471 str = "-1A82F9E6C53222e392";
2472 str >> CXSC_Sqrt3r[2];
2473 str = "-1CB0F41134253Ae35C";
2474 str >> CXSC_Sqrt3r[3];
2475 str = "+1859ED919EC30Be326";
2476 str >> CXSC_Sqrt3r[4];
2477 str = "+1454874FB1F3F4e2EF";
2478 str >> CXSC_Sqrt3r[5];
2479 str = "-1DE69C6D3D2741e2B9";
2480 str >> CXSC_Sqrt3r[6];
2481 str = "+17EEC450C48BE1e283";
2482 str >> CXSC_Sqrt3r[7];
2483 str = "-16F743EEE65D53e24D";
2484 str >> CXSC_Sqrt3r[8];
2485 str = "-1887B505D7E7C2e215";
2486 str >> CXSC_Sqrt3r[9];
2487 str = "-1484D2E10C1161e1DE";
2488 str >> CXSC_Sqrt3r[10];
2489 str = "-1A0B1F86177FB7e1A8";
2490 str >> CXSC_Sqrt3r[11];
2491 str = "+1FE389D3F2C54Ee170";
2492 str >> CXSC_Sqrt3r[12];
2493 str = "+1F29F77C671544e13A";
2494 str >> CXSC_Sqrt3r[13];
2495 str = "-16CE74ED77D9BEe104";
2496 str >> CXSC_Sqrt3r[14];
2497 str = "-1E38708FF0CCB5e0CE";
2498 str >> CXSC_Sqrt3r[15];
2499 str = "-1F13BCC70157D1e098";
2500 str >> CXSC_Sqrt3r[16];
2501 str = "+17EC34CF9B1930e062";
2502 str >> CXSC_Sqrt3r[17];
2503 str = "-117A638EFF3A8Be02B";
2504 str >> CXSC_Sqrt3r[18];
2505 str = "-10016A8EF69C32e000";
2506 str >> CXSC_Sqrt3r[19];
2507 str = "-10016A8EF69C31e000";
2508 str >> CXSC_Sqrt3r[20];
2509
2510 CXSC_Sqrt3r_initialized = true;
2511 std::cout << RestoreOpt;
2512 }
2513 stagprec = stagmax;
2514 y = adjust(l_interval(0.0));
2515 for (int i=0; i <= stagmax; i++)
2516 y.data[i] = CXSC_Sqrt3r[i];
2517// y[i+1] = CXSC_Sqrt3r[i];
2518 stagprec = stagsave;
2519 y = adjust(y);
2520 return y;
2521}
2522
2523
2524// **********************************************************************
2525
2526static real CXSC_LnPi[21]; // CXSC_LnPi[0], ... CXSC_LnPi[20]
2527static bool CXSC_LnPi_initialized = false;
2528
2530// Inclusion of ln(Pi), Blomquist, 12.09.06;
2531{
2532 l_interval y;
2533 int stagsave = stagprec,
2534 stagmax = 20;
2535
2536 if (!CXSC_LnPi_initialized)
2537 {
2538 std::string str;
2539 std::cout << SaveOpt;
2540 std::cout << Hex;
2541 str = "+1250D048E7A1BDe3FF";
2542 str >> CXSC_LnPi[0];
2543 str = "+17ABF2AD8D5088e3C6";
2544 str >> CXSC_LnPi[1];
2545 str = "-16CCF43244818Ae38E";
2546 str >> CXSC_LnPi[2];
2547 str = "+1F9303719C0176e358";
2548 str >> CXSC_LnPi[3];
2549 str = "+15DF52611CB54Ee322";
2550 str >> CXSC_LnPi[4];
2551 str = "-1D9056E74F8C97e2EC";
2552 str >> CXSC_LnPi[5];
2553 str = "+100B095B6C2E1Ae2B5";
2554 str >> CXSC_LnPi[6];
2555 str = "-18C7557878A9E7e27F";
2556 str >> CXSC_LnPi[7];
2557 str = "+1B9BBBB4F4CEE7e248";
2558 str >> CXSC_LnPi[8];
2559 str = "+1B477FCC702F86e212";
2560 str >> CXSC_LnPi[9];
2561 str = "+141F1344A31799e1DC";
2562 str >> CXSC_LnPi[10];
2563 str = "+1B6740BE95CD58e1A6";
2564 str >> CXSC_LnPi[11];
2565 str = "-1F2C63904D27DBe16E";
2566 str >> CXSC_LnPi[12];
2567 str = "+1426F00B933976e136";
2568 str >> CXSC_LnPi[13];
2569 str = "+125703BE5FAA20e100";
2570 str >> CXSC_LnPi[14];
2571 str = "-1DADAE5397F95Be0C9";
2572 str >> CXSC_LnPi[15];
2573 str = "+17C9D110381543e091";
2574 str >> CXSC_LnPi[16];
2575 str = "-1259230E627FCAe05B";
2576 str >> CXSC_LnPi[17];
2577 str = "+191CEAB6B13A33e024";
2578 str >> CXSC_LnPi[18];
2579 str = "+10000109D49A14e000";
2580 str >> CXSC_LnPi[19];
2581 str = "+10000109D49A15e000";
2582 str >> CXSC_LnPi[20];
2583
2584 CXSC_LnPi_initialized = true;
2585 std::cout << RestoreOpt;
2586 }
2587 stagprec = stagmax;
2588 y = adjust(l_interval(0.0));
2589 for (int i=0; i <= stagmax; i++)
2590 y.data[i] = CXSC_LnPi[i];
2591// y[i+1] = CXSC_LnPi[i];
2592 stagprec = stagsave;
2593 y = adjust(y);
2594 return y;
2595}
2596
2597
2598// **********************************************************************
2599
2600static real CXSC_Ln2Pi[21]; // CXSC_Ln2Pi[0], ... CXSC_Ln2Pi[20]
2601static bool CXSC_Ln2Pi_initialized = false;
2602
2604// Inclusion of ln(2*Pi), Blomquist, 12.09.06;
2605{
2606 l_interval y;
2607 int stagsave = stagprec,
2608 stagmax = 20;
2609
2610 if (!CXSC_Ln2Pi_initialized)
2611 {
2612 std::string str;
2613 std::cout << SaveOpt;
2614 std::cout << Hex;
2615 str = "+1D67F1C864BEB5e3FF";
2616 str >> CXSC_Ln2Pi[0];
2617 str = "-165B5A1B7FF5DFe3C9";
2618 str >> CXSC_Ln2Pi[1];
2619 str = "-1B7F70C13DC1CCe392";
2620 str >> CXSC_Ln2Pi[2];
2621 str = "+13458B4DDEC6A3e35C";
2622 str >> CXSC_Ln2Pi[3];
2623 str = "+133DAA155D2130e324";
2624 str >> CXSC_Ln2Pi[4];
2625 str = "-18A007FC5E501Be2EE";
2626 str >> CXSC_Ln2Pi[5];
2627 str = "-15406FA3AA9644e2B4";
2628 str >> CXSC_Ln2Pi[6];
2629 str = "-13E8D52A392CC9e27E";
2630 str >> CXSC_Ln2Pi[7];
2631 str = "-1A43099131E88De248";
2632 str >> CXSC_Ln2Pi[8];
2633 str = "-114835B6623C4De212";
2634 str >> CXSC_Ln2Pi[9];
2635 str = "-1ABB7858CF827Ae1DC";
2636 str >> CXSC_Ln2Pi[10];
2637 str = "+1D8D7045A5A495e1A6";
2638 str >> CXSC_Ln2Pi[11];
2639 str = "+1A26094B3F6FC5e16F";
2640 str >> CXSC_Ln2Pi[12];
2641 str = "-1EF27932D0E3D0e137";
2642 str >> CXSC_Ln2Pi[13];
2643 str = "-12128804136AB6e100";
2644 str >> CXSC_Ln2Pi[14];
2645 str = "+15F8A4AC0BEE17e0C7";
2646 str >> CXSC_Ln2Pi[15];
2647 str = "+1892F2A5B69B5Fe091";
2648 str >> CXSC_Ln2Pi[16];
2649 str = "+1CC7C09477ADCEe05B";
2650 str >> CXSC_Ln2Pi[17];
2651 str = "-116DD579AF074Ae022";
2652 str >> CXSC_Ln2Pi[18];
2653 str = "+100000321C8783e000";
2654 str >> CXSC_Ln2Pi[19];
2655 str = "+100000321C8784e000";
2656 str >> CXSC_Ln2Pi[20];
2657
2658 CXSC_Ln2Pi_initialized = true;
2659 std::cout << RestoreOpt;
2660 }
2661 stagprec = stagmax;
2662 y = adjust(l_interval(0.0));
2663 for (int i=0; i <= stagmax; i++)
2664 y.data[i] = CXSC_Ln2Pi[i];
2665// y[i+1] = CXSC_Ln2Pi[i];
2666 stagprec = stagsave;
2667 y = adjust(y);
2668 return y;
2669}
2670
2671
2672// **********************************************************************
2673
2674static real CXSC_E[21]; // CXSC_E[0], ... CXSC_E[20]
2675static bool CXSC_E_initialized = false;
2676
2678// Inclusion of e=exp(1), Blomquist, 12.09.06;
2679{
2680 l_interval y;
2681 int stagsave = stagprec,
2682 stagmax = 20;
2683
2684 if (!CXSC_E_initialized)
2685 {
2686 std::string str;
2687 std::cout << SaveOpt;
2688 std::cout << Hex;
2689 str = "+15BF0A8B145769e400";
2690 str >> CXSC_E[0];
2691 str = "+14D57EE2B1013Ae3CA";
2692 str >> CXSC_E[1];
2693 str = "-1618713A31D3E2e392";
2694 str >> CXSC_E[2];
2695 str = "+1C5A6D2B53C26De35C";
2696 str >> CXSC_E[3];
2697 str = "-1F75CDE60219B6e326";
2698 str >> CXSC_E[4];
2699 str = "-188C76D93041A1e2EF";
2700 str >> CXSC_E[5];
2701 str = "+12FE363630C75Ee2B9";
2702 str >> CXSC_E[6];
2703 str = "-1C25F937F544EEe283";
2704 str >> CXSC_E[7];
2705 str = "-1E852C20E12A2Ae24D";
2706 str >> CXSC_E[8];
2707 str = "-14D4F6DE605705e212";
2708 str >> CXSC_E[9];
2709 str = "-1F3225EF539355e1D8";
2710 str >> CXSC_E[10];
2711 str = "-16109728625547e1A2";
2712 str >> CXSC_E[11];
2713 str = "-194301506D94CFe16C";
2714 str >> CXSC_E[12];
2715 str = "-1879C78F8CBA44e136";
2716 str >> CXSC_E[13];
2717 str = "-1D5976250C1018e0FD";
2718 str >> CXSC_E[14];
2719 str = "+1C877C56284DABe0C7";
2720 str >> CXSC_E[15];
2721 str = "+1E73530ACCA4F5e091";
2722 str >> CXSC_E[16];
2723 str = "-1F161A150FD53Ae05B";
2724 str >> CXSC_E[17];
2725 str = "+159927DB0E8845e022";
2726 str >> CXSC_E[18];
2727 str = "+10000094BB2C8Ee000";
2728 str >> CXSC_E[19];
2729 str = "+10000094BB2C8Fe000";
2730 str >> CXSC_E[20];
2731
2732 CXSC_E_initialized = true;
2733 std::cout << RestoreOpt;
2734 }
2735 stagprec = stagmax;
2736 y = adjust(l_interval(0.0));
2737 for (int i=0; i <= stagmax; i++)
2738 y.data[i] = CXSC_E[i];
2739// y[i+1] = CXSC_E[i];
2740 stagprec = stagsave;
2741 y = adjust(y);
2742 return y;
2743}
2744
2745
2746// **********************************************************************
2747
2748static real CXSC_Er[21]; // CXSC_Er[0], ... CXSC_Er[20]
2749static bool CXSC_Er_initialized = false;
2750
2752// Inclusion of 1/e, Blomquist, 12.09.06;
2753{
2754 l_interval y;
2755 int stagsave = stagprec,
2756 stagmax = 20;
2757
2758 if (!CXSC_Er_initialized)
2759 {
2760 std::string str;
2761 std::cout << SaveOpt;
2762 std::cout << Hex;
2763 str = "+178B56362CEF38e3FD";
2764 str >> CXSC_Er[0];
2765 str = "-1CA8A4270FADF5e3C6";
2766 str >> CXSC_Er[1];
2767 str = "-1837912B3FD2AAe390";
2768 str >> CXSC_Er[2];
2769 str = "-152711999FB68Ce35A";
2770 str >> CXSC_Er[3];
2771 str = "-17AD7C1289274Ee324";
2772 str >> CXSC_Er[4];
2773 str = "+17E8E56842B705e2E6";
2774 str >> CXSC_Er[5];
2775 str = "-1D24CB13796C2De2B0";
2776 str >> CXSC_Er[6];
2777 str = "-1456AABDA5C8F2e279";
2778 str >> CXSC_Er[7];
2779 str = "+1229F03C6276DDe243";
2780 str >> CXSC_Er[8];
2781 str = "-1569CFC4F53109e20D";
2782 str >> CXSC_Er[9];
2783 str = "-155B63C9B68091e1D5";
2784 str >> CXSC_Er[10];
2785 str = "+1580CF14DC087Ce19F";
2786 str >> CXSC_Er[11];
2787 str = "+1F9FF222313669e168";
2788 str >> CXSC_Er[12];
2789 str = "+15BC9CB1A22487e132";
2790 str >> CXSC_Er[13];
2791 str = "-1857E415C89B13e0FB";
2792 str >> CXSC_Er[14];
2793 str = "+13DF75706E3643e0C5";
2794 str >> CXSC_Er[15];
2795 str = "+13BDF5B7646234e08D";
2796 str >> CXSC_Er[16];
2797 str = "+1C956A5A3BE55De057";
2798 str >> CXSC_Er[17];
2799 str = "-167243FE9CD95Ee020";
2800 str >> CXSC_Er[18];
2801 str = "+1000002F30CCDBe000";
2802 str >> CXSC_Er[19];
2803 str = "+1000002F30CCDCe000";
2804 str >> CXSC_Er[20];
2805
2806 CXSC_Er_initialized = true;
2807 std::cout << RestoreOpt;
2808 }
2809 stagprec = stagmax;
2810 y = adjust(l_interval(0.0));
2811 for (int i=0; i <= stagmax; i++)
2812 y.data[i] = CXSC_Er[i];
2813// y[i+1] = CXSC_Er[i];
2814 stagprec = stagsave;
2815 y = adjust(y);
2816 return y;
2817}
2818
2819
2820// **********************************************************************
2821
2822static real CXSC_Ep2[21]; // CXSC_Ep2[0], ... CXSC_Ep2[20]
2823static bool CXSC_Ep2_initialized = false;
2824
2826// Inclusion of e^2, Blomquist, 12.09.06;
2827{
2828 l_interval y;
2829 int stagsave = stagprec,
2830 stagmax = 20;
2831
2832 if (!CXSC_Ep2_initialized)
2833 {
2834 std::string str;
2835 std::cout << SaveOpt;
2836 std::cout << Hex;
2837 str = "+1D8E64B8D4DDAEe401";
2838 str >> CXSC_Ep2[0];
2839 str = "-19E62E22EFCA4Ce3CA";
2840 str >> CXSC_Ep2[1];
2841 str = "+1577508F5CF5EDe394";
2842 str >> CXSC_Ep2[2];
2843 str = "-186EF0294C2511e35E";
2844 str >> CXSC_Ep2[3];
2845 str = "+177D109F148782e327";
2846 str >> CXSC_Ep2[4];
2847 str = "+166BBC354AB700e2F0";
2848 str >> CXSC_Ep2[5];
2849 str = "-1273AEC0115969e2BA";
2850 str >> CXSC_Ep2[6];
2851 str = "-1C5AE00D3BEEF1e284";
2852 str >> CXSC_Ep2[7];
2853 str = "+15ACA3FDC9595Fe24C";
2854 str >> CXSC_Ep2[8];
2855 str = "-113FCDFE2B1F0Ce215";
2856 str >> CXSC_Ep2[9];
2857 str = "+10EEDFD1AE90C9e1DF";
2858 str >> CXSC_Ep2[10];
2859 str = "+1D2CB8EDC7078Be1A9";
2860 str >> CXSC_Ep2[11];
2861 str = "+11827A19F175F8e173";
2862 str >> CXSC_Ep2[12];
2863 str = "-10267512A9BFB2e13C";
2864 str >> CXSC_Ep2[13];
2865 str = "-19A1E2FC413AE3e105";
2866 str >> CXSC_Ep2[14];
2867 str = "+1170C7A5981ADBe0CF";
2868 str >> CXSC_Ep2[15];
2869 str = "-1FC991480067CFe099";
2870 str >> CXSC_Ep2[16];
2871 str = "-12E9A54CF5CFB5e062";
2872 str >> CXSC_Ep2[17];
2873 str = "-166FA6C468910Ae02A";
2874 str >> CXSC_Ep2[18];
2875 str = "+100043EA6DC142e000";
2876 str >> CXSC_Ep2[19];
2877 str = "+100043EA6DC143e000";
2878 str >> CXSC_Ep2[20];
2879
2880 CXSC_Ep2_initialized = true;
2881 std::cout << RestoreOpt;
2882 }
2883 stagprec = stagmax;
2884 y = adjust(l_interval(0.0));
2885 for (int i=0; i <= stagmax; i++)
2886 y.data[i] = CXSC_Ep2[i];
2887// y[i+1] = CXSC_Ep2[i];
2888 stagprec = stagsave;
2889 y = adjust(y);
2890 return y;
2891}
2892
2893
2894// **********************************************************************
2895
2896static real CXSC_Ep2r[21]; // CXSC_Ep2r[0], ... CXSC_Ep2r[20]
2897static bool CXSC_Ep2r_initialized = false;
2898
2900// Inclusion of 1/e^2, Blomquist, 12.09.06;
2901{
2902 l_interval y;
2903 int stagsave = stagprec,
2904 stagmax = 20;
2905
2906 if (!CXSC_Ep2r_initialized)
2907 {
2908 std::string str;
2909 std::cout << SaveOpt;
2910 std::cout << Hex;
2911 str = "+1152AAA3BF81CCe3FC";
2912 str >> CXSC_Ep2r[0];
2913 str = "-1809224547B4BFe3C6";
2914 str >> CXSC_Ep2r[1];
2915 str = "-16A8E079134F13e390";
2916 str >> CXSC_Ep2r[2];
2917 str = "+14564CACF0994Ee358";
2918 str >> CXSC_Ep2r[3];
2919 str = "+1B796438129AF8e322";
2920 str >> CXSC_Ep2r[4];
2921 str = "-1ACFED57EF2AE5e2EC";
2922 str >> CXSC_Ep2r[5];
2923 str = "-1A968CBDBB5D9De2B5";
2924 str >> CXSC_Ep2r[6];
2925 str = "+1A7238CBD97B71e27C";
2926 str >> CXSC_Ep2r[7];
2927 str = "-146C53DB77BB01e245";
2928 str >> CXSC_Ep2r[8];
2929 str = "-1EEC161C3EBBD7e20C";
2930 str >> CXSC_Ep2r[9];
2931 str = "-12D084DC157ACEe1D5";
2932 str >> CXSC_Ep2r[10];
2933 str = "+12A61F46883347e19F";
2934 str >> CXSC_Ep2r[11];
2935 str = "+1993BAF10CAE0Be164";
2936 str >> CXSC_Ep2r[12];
2937 str = "+1F9224351178FFe12E";
2938 str >> CXSC_Ep2r[13];
2939 str = "-1C366D1C7BA64Ae0F7";
2940 str >> CXSC_Ep2r[14];
2941 str = "-17D9938EFA4657e0C0";
2942 str >> CXSC_Ep2r[15];
2943 str = "+1B6668DF0C1286e08A";
2944 str >> CXSC_Ep2r[16];
2945 str = "+1F7A4FFC9B48C6e050";
2946 str >> CXSC_Ep2r[17];
2947 str = "+1F3E3AF6F17591e01A";
2948 str >> CXSC_Ep2r[18];
2949 str = "+100000006C7831e000";
2950 str >> CXSC_Ep2r[19];
2951 str = "+100000006C7832e000";
2952 str >> CXSC_Ep2r[20];
2953
2954 CXSC_Ep2r_initialized = true;
2955 std::cout << RestoreOpt;
2956 }
2957 stagprec = stagmax;
2958 y = adjust(l_interval(0.0));
2959 for (int i=0; i <= stagmax; i++)
2960 y.data[i] = CXSC_Ep2r[i];
2961// y[i+1] = CXSC_Ep2r[i];
2962 stagprec = stagsave;
2963 y = adjust(y);
2964 return y;
2965}
2966
2967
2968// **********************************************************************
2969
2970static real CXSC_EpPi[21]; // CXSC_EpPi[0], ... CXSC_EpPi[20]
2971static bool CXSC_EpPi_initialized = false;
2972
2974// Inclusion of e^(Pi), Blomquist, 12.09.06;
2975{
2976 l_interval y;
2977 int stagsave = stagprec,
2978 stagmax = 20;
2979
2980 if (!CXSC_EpPi_initialized)
2981 {
2982 std::string str;
2983 std::cout << SaveOpt;
2984 std::cout << Hex;
2985 str = "+1724046EB0933Ae403";
2986 str >> CXSC_EpPi[0];
2987 str = "-184C962DD81952e3CD";
2988 str >> CXSC_EpPi[1];
2989 str = "-12D659C0BCD22Ee396";
2990 str >> CXSC_EpPi[2];
2991 str = "+117496B8A92F91e360";
2992 str >> CXSC_EpPi[3];
2993 str = "+16A8C4203E5FCDe32A";
2994 str >> CXSC_EpPi[4];
2995 str = "-166B11F99A663Be2F4";
2996 str >> CXSC_EpPi[5];
2997 str = "-118EC2076DABB1e2BE";
2998 str >> CXSC_EpPi[6];
2999 str = "+19776E5BEB18A5e288";
3000 str >> CXSC_EpPi[7];
3001 str = "+1AD4091E84B051e252";
3002 str >> CXSC_EpPi[8];
3003 str = "+1E89AA12909B40e21C";
3004 str >> CXSC_EpPi[9];
3005 str = "+1ACE3C0DDBB994e1E3";
3006 str >> CXSC_EpPi[10];
3007 str = "+141EC9379CBBFEe1AC";
3008 str >> CXSC_EpPi[11];
3009 str = "+1FC4E78D00A016e173";
3010 str >> CXSC_EpPi[12];
3011 str = "+1608BE35B9A409e13D";
3012 str >> CXSC_EpPi[13];
3013 str = "-1A0D8AA90EB6B9e103";
3014 str >> CXSC_EpPi[14];
3015 str = "+106FE8AFD21ACFe0CD";
3016 str >> CXSC_EpPi[15];
3017 str = "+1C072FEA1BFCAFe095";
3018 str >> CXSC_EpPi[16];
3019 str = "+1915B9F352EC68e05B";
3020 str >> CXSC_EpPi[17];
3021 str = "-13FA07C37897E9e024";
3022 str >> CXSC_EpPi[18];
3023 str = "-100003D8039138e000";
3024 str >> CXSC_EpPi[19];
3025 str = "-100003D8039137e000";
3026 str >> CXSC_EpPi[20];
3027
3028 CXSC_EpPi_initialized = true;
3029 std::cout << RestoreOpt;
3030 }
3031 stagprec = stagmax;
3032 y = adjust(l_interval(0.0));
3033 for (int i=0; i <= stagmax; i++)
3034 y.data[i] = CXSC_EpPi[i];
3035// y[i+1] = CXSC_EpPi[i];
3036 stagprec = stagsave;
3037 y = adjust(y);
3038 return y;
3039}
3040
3041
3042// **********************************************************************
3043
3044static real CXSC_Ep2Pi[21]; // CXSC_Ep2Pi[0], ... CXSC_Ep2Pi[20]
3045static bool CXSC_Ep2Pi_initialized = false;
3046
3048// Inclusion of e^(2*Pi), Blomquist, 12.09.06;
3049{
3050 l_interval y;
3051 int stagsave = stagprec,
3052 stagmax = 20;
3053
3054 if (!CXSC_Ep2Pi_initialized)
3055 {
3056 std::string str;
3057 std::cout << SaveOpt;
3058 std::cout << Hex;
3059 str = "+10BBEEE9177E19e408";
3060 str >> CXSC_Ep2Pi[0];
3061 str = "+1C7DD9272526B1e3D0";
3062 str >> CXSC_Ep2Pi[1];
3063 str = "+15200F57AB89EDe39A";
3064 str >> CXSC_Ep2Pi[2];
3065 str = "-1FCCB6EDBE9C36e363";
3066 str >> CXSC_Ep2Pi[3];
3067 str = "+1BEA0BF179A589e32D";
3068 str >> CXSC_Ep2Pi[4];
3069 str = "-1F3AD5A6B77F9Ee2F7";
3070 str >> CXSC_Ep2Pi[5];
3071 str = "-1622F702B57637e2C0";
3072 str >> CXSC_Ep2Pi[6];
3073 str = "-100C09AE818734e287";
3074 str >> CXSC_Ep2Pi[7];
3075 str = "+10DA7ADA79EFE6e24D";
3076 str >> CXSC_Ep2Pi[8];
3077 str = "+1FF9BF48B72959e216";
3078 str >> CXSC_Ep2Pi[9];
3079 str = "-17AD7A3F6D2A14e1E0";
3080 str >> CXSC_Ep2Pi[10];
3081 str = "+1FCD4B0FA971E4e1A9";
3082 str >> CXSC_Ep2Pi[11];
3083 str = "+193A2CDC04526Be172";
3084 str >> CXSC_Ep2Pi[12];
3085 str = "-18CBE5FDFAF25Fe13C";
3086 str >> CXSC_Ep2Pi[13];
3087 str = "+1D47EEE171DA93e105";
3088 str >> CXSC_Ep2Pi[14];
3089 str = "-15B0F8DA29DB32e0CF";
3090 str >> CXSC_Ep2Pi[15];
3091 str = "-19207AD7E637D8e097";
3092 str >> CXSC_Ep2Pi[16];
3093 str = "+191CA743F265A6e061";
3094 str >> CXSC_Ep2Pi[17];
3095 str = "+1A15069182EF28e02A";
3096 str >> CXSC_Ep2Pi[18];
3097 str = "-1000EAC5FC05A9e000";
3098 str >> CXSC_Ep2Pi[19];
3099 str = "-1000EAC5FC05A8e000";
3100 str >> CXSC_Ep2Pi[20];
3101
3102 CXSC_Ep2Pi_initialized = true;
3103 std::cout << RestoreOpt;
3104 }
3105 stagprec = stagmax;
3106 y = adjust(l_interval(0.0));
3107 for (int i=0; i <= stagmax; i++)
3108 y.data[i] = CXSC_Ep2Pi[i];
3109// y[i+1] = CXSC_Ep2Pi[i];
3110 stagprec = stagsave;
3111 y = adjust(y);
3112 return y;
3113}
3114
3115
3116// **********************************************************************
3117
3118static real CXSC_EpPid2[21]; // CXSC_EpPid2[0], ... CXSC_EpPid2[20]
3119static bool CXSC_EpPid2_initialized = false;
3120
3122// Inclusion of e^(Pi/2), Blomquist, 12.09.06;
3123{
3124 l_interval y;
3125 int stagsave = stagprec,
3126 stagmax = 20;
3127
3128 if (!CXSC_EpPid2_initialized)
3129 {
3130 std::string str;
3131 std::cout << SaveOpt;
3132 std::cout << Hex;
3133 str = "+133DEDC855935Fe401";
3134 str >> CXSC_EpPid2[0];
3135 str = "+13E45A768FB73Ce3CB";
3136 str >> CXSC_EpPid2[1];
3137 str = "-1FB31CF300FF3Ce395";
3138 str >> CXSC_EpPid2[2];
3139 str = "-1E80D8BEB83F79e35F";
3140 str >> CXSC_EpPid2[3];
3141 str = "-14A3DE039142DDe326";
3142 str >> CXSC_EpPid2[4];
3143 str = "-18792D7A37282Be2EB";
3144 str >> CXSC_EpPid2[5];
3145 str = "-19DF43A5980C28e2B5";
3146 str >> CXSC_EpPid2[6];
3147 str = "-1C6F0F641C0D67e27F";
3148 str >> CXSC_EpPid2[7];
3149 str = "-1779C86C2DB5ACe249";
3150 str >> CXSC_EpPid2[8];
3151 str = "+168521EE91B16Fe211";
3152 str >> CXSC_EpPid2[9];
3153 str = "+12530F905D97BDe1DB";
3154 str >> CXSC_EpPid2[10];
3155 str = "+13498112CB7585e1A5";
3156 str >> CXSC_EpPid2[11];
3157 str = "+1BA4546B13A434e16F";
3158 str >> CXSC_EpPid2[12];
3159 str = "+14FF791C56421Ce138";
3160 str >> CXSC_EpPid2[13];
3161 str = "-1F375C223A2152e102";
3162 str >> CXSC_EpPid2[14];
3163 str = "-126AB0C8C77412e0CC";
3164 str >> CXSC_EpPid2[15];
3165 str = "-1B39C9C0B8C54Ae094";
3166 str >> CXSC_EpPid2[16];
3167 str = "-167741414E31E3e05D";
3168 str >> CXSC_EpPid2[17];
3169 str = "+1DEFB4462546C1e025";
3170 str >> CXSC_EpPid2[18];
3171 str = "-1000010F7B89CDe000";
3172 str >> CXSC_EpPid2[19];
3173 str = "-1000010F7B89CCe000";
3174 str >> CXSC_EpPid2[20];
3175
3176 CXSC_EpPid2_initialized = true;
3177 std::cout << RestoreOpt;
3178 }
3179 stagprec = stagmax;
3180 y = adjust(l_interval(0.0));
3181 for (int i=0; i <= stagmax; i++)
3182 y.data[i] = CXSC_EpPid2[i];
3183// y[i+1] = CXSC_EpPid2[i];
3184 stagprec = stagsave;
3185 y = adjust(y);
3186 return y;
3187}
3188
3189
3190// **********************************************************************
3191
3192static real CXSC_EpPid4[21]; // CXSC_EpPid4[0], ... CXSC_EpPid4[20]
3193static bool CXSC_EpPid4_initialized = false;
3194
3196// Inclusion of e^(Pi/4), Blomquist, 12.09.06;
3197{
3198 l_interval y;
3199 int stagsave = stagprec,
3200 stagmax = 20;
3201
3202 if (!CXSC_EpPid4_initialized)
3203 {
3204 std::string str;
3205 std::cout << SaveOpt;
3206 std::cout << Hex;
3207 str = "+118BD669471CAAe400";
3208 str >> CXSC_EpPid4[0];
3209 str = "+1F0ED609715756e3CA";
3210 str >> CXSC_EpPid4[1];
3211 str = "-1B9C7B871FE1DBe394";
3212 str >> CXSC_EpPid4[2];
3213 str = "+15C0FECE98F209e35D";
3214 str >> CXSC_EpPid4[3];
3215 str = "+18C9FACC5DF3CEe327";
3216 str >> CXSC_EpPid4[4];
3217 str = "+15EDE838B4A399e2EF";
3218 str >> CXSC_EpPid4[5];
3219 str = "-1C7EFACA363051e2B9";
3220 str >> CXSC_EpPid4[6];
3221 str = "-1A1EBEA1646411e283";
3222 str >> CXSC_EpPid4[7];
3223 str = "+1AEF54E68CE03Be24C";
3224 str >> CXSC_EpPid4[8];
3225 str = "-11250CB97FDDBFe212";
3226 str >> CXSC_EpPid4[9];
3227 str = "-169ADC0E65B8A7e1DB";
3228 str >> CXSC_EpPid4[10];
3229 str = "+198A501DB90EDDe1A5";
3230 str >> CXSC_EpPid4[11];
3231 str = "-1586909A3F6365e16E";
3232 str >> CXSC_EpPid4[12];
3233 str = "+1BE542410F8CE7e138";
3234 str >> CXSC_EpPid4[13];
3235 str = "+1E7EEC51889EECe102";
3236 str >> CXSC_EpPid4[14];
3237 str = "-1913C9FC19333Ce0CC";
3238 str >> CXSC_EpPid4[15];
3239 str = "+1112C71EA1E6F0e095";
3240 str >> CXSC_EpPid4[16];
3241 str = "-1C4CCF0F5D1E14e05E";
3242 str >> CXSC_EpPid4[17];
3243 str = "+1AC4A72310FA27e028";
3244 str >> CXSC_EpPid4[18];
3245 str = "-100013EC6A07AEe000";
3246 str >> CXSC_EpPid4[19];
3247 str = "-100013EC6A07ADe000";
3248 str >> CXSC_EpPid4[20];
3249
3250 CXSC_EpPid4_initialized = true;
3251 std::cout << RestoreOpt;
3252 }
3253 stagprec = stagmax;
3254 y = adjust(l_interval(0.0));
3255 for (int i=0; i <= stagmax; i++)
3256 y.data[i] = CXSC_EpPid4[i];
3257// y[i+1] = CXSC_EpPid4[i];
3258 stagprec = stagsave;
3259 y = adjust(y);
3260 return y;
3261}
3262
3263
3264// **********************************************************************
3265
3266static real CXSC_EulerGa[21]; // CXSC_EulerGa[0], ... CXSC_EulerGa[20]
3267static bool CXSC_EulerGa_initialized = false;
3268
3270// Inclusion of EulerGamma, Blomquist, 12.09.06;
3271{
3272 l_interval y;
3273 int stagsave = stagprec,
3274 stagmax = 20;
3275
3276 if (!CXSC_EulerGa_initialized)
3277 {
3278 std::string str;
3279 std::cout << SaveOpt;
3280 std::cout << Hex;
3281 str = "+12788CFC6FB619e3FE";
3282 str >> CXSC_EulerGa[0];
3283 str = "-16CB90701FBFABe3C5";
3284 str >> CXSC_EulerGa[1];
3285 str = "-134A95E3133C51e38F";
3286 str >> CXSC_EulerGa[2];
3287 str = "+19730064300F7De359";
3288 str >> CXSC_EulerGa[3];
3289 str = "-171ECA0084E369e322";
3290 str >> CXSC_EulerGa[4];
3291 str = "-1302FE2B078898e2EC";
3292 str >> CXSC_EulerGa[5];
3293 str = "+192732D88415F4e2B5";
3294 str >> CXSC_EulerGa[6];
3295 str = "+11056AE9132136e27F";
3296 str >> CXSC_EulerGa[7];
3297 str = "-17DC6F12E630A3e249";
3298 str >> CXSC_EulerGa[8];
3299 str = "+175FD4B1BD70F2e212";
3300 str >> CXSC_EulerGa[9];
3301 str = "-19BC9466120C20e1DC";
3302 str >> CXSC_EulerGa[10];
3303 str = "-18FD5699260EADe1A6";
3304 str >> CXSC_EulerGa[11];
3305 str = "-12EA987665551Fe16F";
3306 str >> CXSC_EulerGa[12];
3307 str = "-1FB159BA4A423De138";
3308 str >> CXSC_EulerGa[13];
3309 str = "+1FA543D43BCC60e102";
3310 str >> CXSC_EulerGa[14];
3311 str = "-1E6F04E0F639F6e0C9";
3312 str >> CXSC_EulerGa[15];
3313 str = "-1A23768654F43De091";
3314 str >> CXSC_EulerGa[16];
3315 str = "-14F1C5CB4F55EBe058";
3316 str >> CXSC_EulerGa[17];
3317 str = "+1E71DF52EDAA7Fe020";
3318 str >> CXSC_EulerGa[18];
3319 str = "+1000001C398F9Be000";
3320 str >> CXSC_EulerGa[19];
3321 str = "+1000001C398F9Ce000";
3322 str >> CXSC_EulerGa[20];
3323
3324 CXSC_EulerGa_initialized = true;
3325 std::cout << RestoreOpt;
3326 }
3327 stagprec = stagmax;
3328 y = adjust(l_interval(0.0));
3329 for (int i=0; i <= stagmax; i++)
3330 y.data[i] = CXSC_EulerGa[i];
3331// y[i+1] = CXSC_EulerGa[i];
3332 stagprec = stagsave;
3333 y = adjust(y);
3334 return y;
3335}
3336
3337
3338// **********************************************************************
3339
3340static real CXSC_Catalan[21]; // CXSC_Catalan[0], ... CXSC_Catalan[20]
3341static bool CXSC_Catalan_initialized = false;
3342
3344// Inclusion of Catalan-Constant, Blomquist, 12.09.06;
3345{
3346 l_interval y;
3347 int stagsave = stagprec,
3348 stagmax = 20;
3349
3350 if (!CXSC_Catalan_initialized)
3351 {
3352 std::string str;
3353 std::cout << SaveOpt;
3354 std::cout << Hex;
3355 str = "+1D4F9713E8135De3FE";
3356 str >> CXSC_Catalan[0];
3357 str = "+11485608B8DF4De3C5";
3358 str >> CXSC_Catalan[1];
3359 str = "-12F39C13BC1EC8e38F";
3360 str >> CXSC_Catalan[2];
3361 str = "+1C2FF8094A263Ee357";
3362 str >> CXSC_Catalan[3];
3363 str = "+168F335DBE5370e321";
3364 str >> CXSC_Catalan[4];
3365 str = "+16291BBB16163Ee2E9";
3366 str >> CXSC_Catalan[5];
3367 str = "+124D663F739C43e2B3";
3368 str >> CXSC_Catalan[6];
3369 str = "-136A0725ED0E94e27B";
3370 str >> CXSC_Catalan[7];
3371 str = "-1D3A26F9C06FCEe240";
3372 str >> CXSC_Catalan[8];
3373 str = "-164E42486BFCD2e209";
3374 str >> CXSC_Catalan[9];
3375 str = "14F358CFDEC843e1D3";
3376 str >> CXSC_Catalan[10];
3377 str = "-11EB82210976ABe19D";
3378 str >> CXSC_Catalan[11];
3379 str = "-17D31F6DF5E801e167";
3380 str >> CXSC_Catalan[12];
3381 str = "+13FD19CE3E396Ae131";
3382 str >> CXSC_Catalan[13];
3383 str = "-1C8CBB3852FF3Fe0F8";
3384 str >> CXSC_Catalan[14];
3385 str = "+1A86EB34EAD01Ae0C2";
3386 str >> CXSC_Catalan[15];
3387 str = "+1C68C37800513Be087";
3388 str >> CXSC_Catalan[16];
3389 str = "+1D46EBB334D7C9e050";
3390 str >> CXSC_Catalan[17];
3391 str = "-1944C5E2711625e019";
3392 str >> CXSC_Catalan[18];
3393 str = "-100000005E2172e000";
3394 str >> CXSC_Catalan[19];
3395 str = "-100000005E2171e000";
3396 str >> CXSC_Catalan[20];
3397
3398 CXSC_Catalan_initialized = true;
3399 std::cout << RestoreOpt;
3400 }
3401 stagprec = stagmax;
3402 y = adjust(l_interval(0.0));
3403 for (int i=0; i <= stagmax; i++)
3404 y.data[i] = CXSC_Catalan[i];
3405// y[i+1] = CXSC_Catalan[i];
3406 stagprec = stagsave;
3407 y = adjust(y);
3408 return y;
3409}
3410
3411
3412l_interval atan(const l_interval & x) noexcept // aTan(x)
3413{
3414 int stagsave = stagprec,
3415 stagmax = 19,
3416 digits = 53,
3417 m = 1,
3418 stagsave2,
3419 degree, sign, zk;
3420 real lnt2, lneps, ln3m, eps,
3421 // bas, tn, t3,
3422 zhn = 1.0,
3423 lnb = 0.69314718, // ln(2.0)
3424 ln3 = 1.098612289; // ln(3.0), Blomquist 27.12.03
3425 l_interval t, t2, p,
3426 y;
3427 interval dx = interval(x),
3428 einfachgenau,
3429 err;
3430
3431 einfachgenau = atan(dx);
3432
3433 if (stagprec == 1)
3434 y = atan(dx);
3435 else if (Inf(dx) == Sup(dx) && Sup(dx) == CXSC_Zero)
3436 y = adjust(l_interval(0.0));
3437 else
3438 {
3439 if (dx == interval(1.0)) y = li_pi4(); // y = Pi/4
3440 else
3441 {
3442 if (stagprec < stagmax)
3443 stagprec++;
3444 else
3445 stagprec = stagmax;
3446
3447 // Argumentreduktion
3448 t = x;
3449 eps = 0.01/stagprec; // vorher 0.1, aber zu gross
3450 while (Sup(abs(interval(t))) > eps)
3451 { // t = t/(1.0+sqrt(1.0+t*t));
3452 t = t/(1.0+sqrt1px2(t)); // Blomquist, 13.12.02;
3453 zhn += zhn;
3454 }
3455 t2 = t*t;
3456
3457 // Bestimmung des Approximationsgrades
3458 // Beginn Blomquist, 27.12.03 ------------------------------
3459 err = interval(abs(t2));
3460 if (expo(Sup(err)) < -300) m = 4; // Vermeidet alte Fehlermeldg.
3461 else
3462 {
3463 // lnt2 = Sup(ln(err)); // Erzeugt alten Fehler nicht mehr!! ALT!!
3464 lnt2 = ln(Sup(err)); // Neu, 24.08.12;
3465 ln3m = ln(3.0-real(Sup(t2)));
3466 lneps = (1.0-digits*stagprec)*lnb;
3467 do {
3468 m += 3;
3469 } while (lneps-ln3-real(m+1)*lnt2+ln(2.0*m+3.0)+ln3m <= 0.0);
3470 }
3471 // Ende Blomquist, 27.12.03 ---------------------------------
3472
3473 degree = m;
3474
3475 // Approximation
3476 sign = (degree%2) ? -1 : 1;
3477 zk = sign*(2*degree+1);
3478 p = l_interval(1.0)/real(zk);
3479 for (int k = degree-1; k >= 0; k--)
3480 {
3481 sign = -sign;
3482 zk = sign*(2*k+1);
3483 p = p*t2+l_interval(1.0)/real(zk);
3484 } // letztes t wird beim Fehler mit einmultipliziert
3485
3486 // Fehler bestimmen
3487 stagsave2 = stagprec;
3488 stagprec = 2;
3489 l_interval relerr;
3490 stagprec = stagsave2;
3491 err = pow(interval(2.0), interval(1.0-digits*stagprec))
3492 * interval(-1.0, 1.0); // von 2^(2-d*s) ge\x{201E}ndert !
3493 relerr[1] = 1.0;
3494 relerr[2] = Inf(err);
3495 relerr[3] = Sup(err);
3496
3497 // Argumentreduktion rueckgaengig machen, mit t mult.
3498 // und Fehler einbauen
3499 y = zhn*t*p*relerr;
3500
3501 stagprec = stagsave;
3502 y = adjust(y);
3503 y = y & einfachgenau;
3504 }
3505 }
3506
3507 return y;
3508}
3509
3510l_interval acot(const l_interval &x) noexcept // ACot(x)
3511{
3512 l_interval pihalbe, y;
3513 interval dx = interval(x),
3514 einfachgenau;
3515
3516 einfachgenau = acot(dx);
3517 stagprec++;
3518// pihalbe = 2.0*atan(l_interval(1.0));
3519 pihalbe = li_pi4();
3520 times2pown(pihalbe,1); // Blomquist, 05.12.03;
3521 stagprec--;
3522
3523 if (stagprec == 1)
3524 y = acot(dx);
3525 else if (Inf(dx) == Sup(dx) && Sup(dx) == CXSC_Zero)
3526 y = adjust(pihalbe);
3527 else
3528 {
3529 y = pihalbe-atan(x);
3530 // y = atan(1.0/x); andere M�glichkeit der Berechnung
3531 y = adjust(y);
3532 y = y & einfachgenau;
3533 }
3534
3535 return y;
3536}
3537
3538l_interval exp(const l_interval & x) // exp(x)
3539{
3540 long int n = 2;
3541 int stagsave = stagprec,
3542 stagmax = 19,
3543 degree,
3544 rednum = 0,
3545 digits = 53;
3546 real fak = 2.0,
3547 zhn = 1.0, // kann groesser 32768 werden
3548 tmp, lny, lneps, eps,
3549 lnb = 0.69314718;
3550 l_interval p, t,
3551 y;
3552 interval dx = interval(x),
3553 einfachgenau,
3554 dt, error;
3555 // ilnb =interval(0.6931471, 0.6931572),
3556 // ilny, ilneps;
3557
3558 einfachgenau = exp(dx);
3559
3560 // gueltigen Bereich pruefen
3561 if (stagprec == 1)
3562 y = exp(dx);
3563 else if (Inf(dx) == Sup(dx) && Inf(dx) == CXSC_Zero)
3564 y = adjust(l_interval(1.0));
3565 else if (Inf(dx)<-708.396418532264) y = einfachgenau; // Blomquist,12.07.04
3566 else
3567 {
3568 if (stagprec < stagmax)
3569 stagprec++;
3570 else
3571 stagprec = stagmax;
3572
3573 if (Sup(dx) <= 0.0)
3574 t = -x; // Werte unter -1e308 nicht auswertbar,
3575 else
3576 t = x; // dafuer andere Ergebnisse exakter
3577
3578 dt = interval(t);
3579
3580 // Argumentreduktion
3581 // eps = 0.000001*real(Sup(t))/stagprec;
3582 // eps = 1.0/pow10(stagprec/4);
3583 eps = 0.1;
3584 while(Sup(dt)/zhn > eps)
3585 {
3586 rednum++;
3587 zhn += zhn;
3588 }
3589 t /= l_interval(zhn);
3590
3591 // Taylor Approximation
3592 dt = interval(t);
3593
3594 // Anzahl der Durchlaeufe bei der Approximation
3595 tmp = Sup(abs(dt));
3596 if (MinReal > tmp)
3597 tmp = MinReal;
3598 lny = ln(tmp);
3599 lneps = (1.0-digits*stagprec)*lnb;
3600
3601 while(lneps-lnb+ln(fak)-real(n)*lny <= 0.0)
3602 {
3603 n += 3;
3604 if (n > 170)
3605 { // 170! = 7.26*10306
3606 cxscthrow(ERROR_LINTERVAL_FAK_OVERFLOW("l_interval exp(const l_interval & x)"));
3607 n = 170;
3608 break;
3609 }
3610 fak = fak*real(n)*(n-1.0)*(n-2.0);
3611 }
3612 degree = n;
3613
3614 p = t/real(degree);
3615 int i;
3616 for (i=degree-1; i >= 1; i--)
3617 p = (p+1.0)*t/real(i);
3618
3619 // Fehlerabschaetzung
3620 error = interval(-2.0, 2.0)*pow(dt, interval(real(n)))/fak;
3621
3622 // Fehler einbauen
3623 p += 1.0+l_interval(error);
3624
3625 // Argumentreduktion wird rueckgaengig gemacht
3626 for (i = 1; i <= rednum; i++)
3627 p *= p;
3628
3629 if (Sup(dx) <= 0.0)
3630 p = 1.0/p;
3631
3632 stagprec = stagsave;
3633 y = adjust(p);
3634 y = y & einfachgenau;
3635 }
3636
3637 return y;
3638}
3639
3640l_interval exp2(const l_interval & x) // 2^x
3641{
3642 int stagsave = stagprec,
3643 stagmax = 19;
3644 if (stagprec>stagmax)
3645 stagprec = stagmax;
3646 l_interval y;
3647
3648 y = exp(x*Ln2_l_interval());
3649
3650 stagprec = stagsave;
3651 y = adjust(y);
3652
3653 return y;
3654}
3655
3656l_interval exp10(const l_interval & x) // 10^x
3657{
3658 int stagsave = stagprec,
3659 stagmax = 19;
3660 if (stagprec>stagmax)
3661 stagprec = stagmax;
3662 l_interval y;
3663
3664 y = exp(x*Ln10_l_interval());
3665
3666 stagprec = stagsave;
3667 y = adjust(y);
3668
3669 return y;
3670}
3671
3672l_interval expm1(const l_interval & x) noexcept
3673// exp(x)-1;
3674{
3675 l_interval y(0.0);
3676 real B,S_absdx,Sbru,abserr;
3677 interval dx = interval(x),
3678 einfachgenau,sx,bru;
3679 int stagsave = stagprec,ex,N,
3680 stagmax = 19; // from exponential function
3681 einfachgenau = expm1(dx); // produces all necessary error messages!
3682
3683 if (stagprec == 1) y = einfachgenau;
3684 else
3685 {
3686 if (stagprec>stagmax) stagprec = stagmax;
3687 S_absdx = Sup( abs(dx) );
3688 ex = expo(S_absdx);
3689 if (ex<-49) { // Taylor series
3690 // P_N(x) := x + x^2/2! + ... + x^N/N!
3691 sx = interval(S_absdx); // sx: point interval
3692 N = ex-53*stagprec; // This N is here not the polynomial degree!
3693 if (N<-1022) N = -1022;
3694 if (N > 0) goto Fertig; // because S_absdx = 0;
3695 B = comp(0.5,N); // B: rough bound of the absolute error
3696 // The absolute error should be smaller than B.
3697 N=0; bru = sx; // N is now the polynomial degree!
3698 // Calculation of the polynomia degree N:
3699 // Calculation of the absolute error:
3700 // |absolute error| := |(exp(x)-1) - P_N(x)| <= AE;
3701 // AE = S_absdx^(N+1)/[(N+1)!*(1-S_absdx)]
3702 do
3703 {
3704 N++;
3705 bru = (bru * S_absdx)/(N+1);
3706 Sbru = Sup(bru);
3707 }
3708 while(Sbru > B);
3709 bru = bru * 1.00000001; // for 1/(1-S_absdx)
3710 abserr = Sup(bru);
3711 // |absolute error| <= abserr;
3712 // Caculating an inclusion y of P_N(x):
3713 y = x/N;
3714 for (int i=N-1; i>=1; i--)
3715 y = (y+1)*x/i;
3716 // x + x^2/2! + ... + x^N/N! = P_N(x) included by y
3717 // Conserning the absolute error:
3718 y = y + interval(-abserr,abserr);
3719
3720 } else {
3721 stagprec++;
3722 if (stagprec>stagmax) stagprec = stagmax;
3723 y = exp(x) - 1;
3724 }
3725 Fertig:
3726 stagprec = stagsave; // restore old stagprec value
3727 y = adjust(y); // adjust precision of y to old stagprec value
3728 y = y & einfachgenau;
3729 }
3730 return y;
3731} // expm1
3732
3734// e^(-x^2); Blomquist, 13.04.04;
3735{
3736 int stagsave = stagprec,
3737 stagmax = 19;
3738 l_interval y,z = abs(x);
3739 l_real r1,r2;
3740 interval dx = interval(z),
3741 einfachgenau;
3742 einfachgenau = expmx2(dx);
3743 if (stagprec>stagmax) stagprec = stagmax;
3744
3745 if (stagprec == 1) y = expmx2(dx); // simple precision
3746 else if (Inf(z)>30) y = einfachgenau; // verhindert Overflow bei z*z!
3747 else y = exp(-z*z);
3748
3749 stagprec = stagsave;
3750 y = adjust(y);
3751 y = y & einfachgenau;
3752 return y;
3753} // expmx2
3754
3755
3756
3757
3758int cxsc_zweihoch(int) noexcept;
3759
3760l_interval ln(const l_interval & x) // Ln(x)
3761{
3762 int stagsave = stagprec,
3763 stagmax = 19,
3764 k,
3765 m = 0,
3766 n = 0,
3767 digits = 53, // Anzahl der Mantissenstellen
3768 cmp,
3769 zhn1 = 1,
3770 zhn2 = 1;
3771 bool zwei=false;
3772 real mx, tmp,
3773 // zkp1,
3774 ln2 = 0.69314718,
3775 lnb = 0.69314718,
3776 lny, lneps;
3777 l_interval y, t, t2, p;
3778 interval dx = interval(x),
3779 einfachgenau,
3780 dy, error;
3781
3782 einfachgenau = ln(dx);
3783 // ---------------------------------------------------------------------
3784 if (Sup(dx) > succ(succ(Inf(dx)))) { // Dieser Teil vermeidet eine
3785 y = einfachgenau; y = adjust(y); // Fehlermeldung der pow(..)-Fkt,
3786 goto fertig; // wenn x zu breit ist und die 1
3787 } // enthaelt. Bei zu breiten Inter-
3788 // vallen wird y daher mit einfachgenau berechnet. Blomquist,05.12.03;
3789 // ---------------------------------------------------------------------
3790 if (Inf(x) <= 0.0)
3791 cxscthrow(ERROR_LINTERVAL_STD_FKT_OUT_OF_DEF("l_interval ln(const l_interval & x)"));
3792 else if (stagprec == 1)
3793 y = ln(dx);
3794 else if (Inf(dx) == Sup(dx) && Sup(dx) == CXSC_One)
3795 y = adjust(l_interval(0.0));
3796 else
3797 {
3798 if (stagprec < stagmax)
3799 stagprec++;
3800 else
3801 stagprec = stagmax;
3802
3803 // erste Argumentreduktion
3804 /*
3805 if (Inf(dx) > 1.0)
3806 {
3807 y = 1.0/x;
3808 dx = interval(x);
3809 neg = true;
3810 } else
3811 y = x;
3812 */
3813 y = x;
3814
3815 // zweite Argumentreduktion
3816 mx = Sup(dx);
3817 tmp = 1.0/Sup(dx);
3818 if (tmp > mx)
3819 mx=tmp;
3820 tmp = 1.0/Inf(dx);
3821 if (tmp > mx)
3822 mx = tmp;
3823 if (mx > 1.1)
3824 {
3825 cmp = int(_double(ln(ln(mx)/ln(1.1))/ln2));
3826 if (cmp > 0)
3827 {
3828 if (cmp > 14)
3829 n = 14;
3830 else
3831 n = cmp;
3832 zhn1 = cxsc_zweihoch(n);
3833 }
3834 }
3835// y = sqrt(y, zhn1); // Old slow version
3836 if (zhn1 != 1)
3837 for (int k=1; k<=n; k++)
3838 y = sqrt(y); // Blomquist, 26.06.2007;
3839 mx = abs(real(Sup(y)));
3840 if (mx > 1.1)
3841 {
3842 cmp = int(_double(ln(ln(mx)/ln(1.1))/ln2));
3843 if (cmp > 0)
3844 {
3845 if (cmp > 14)
3846 n = 14;
3847 else
3848 n = cmp;
3849 zhn2 = cxsc_zweihoch(n);
3850 }
3851// y = sqrt(y, zhn2); // Old slow version
3852 if (zhn2 != 1)
3853 for (int k=1; k<=n; k++)
3854 y = sqrt(y); // Blomquist, 26.06.2007;
3855 zwei = TRUE;
3856 }
3857
3858 // dritte Argumentreduktion
3859 t = (y-1.0)/(y+1.0);
3860
3861 // Abschaetzung der Genauigkeit
3862 t2 = t*t;
3863 tmp = Sup(interval(t2));
3864 if (tmp == 0.0)
3865 tmp = MinReal;
3866 dy = tmp;
3867 lny = Sup(ln(dy));
3868 lneps = (1.0-digits*stagprec)*lnb;
3869
3870 do {
3871 m += 2;
3872 } while (lneps-lnb+ln(2.0*m+3.0)-real(m+1)*lny <= 0.0);
3873
3874 // Polynomauswertung
3875 p = 0.0;
3876 for (k = 2*m+1; k >= 1; k -= 2)
3877 p = p*t2+l_interval(2.0)/real(k);
3878 p *= t;
3879
3880 // Fehlereinschliessung
3881 error = interval(-4.0, 4.0)*pow(interval(t2), interval(real(m+1)))/(2.0*m+3.0);
3882
3883 // Fehler einbauen
3884 p = p+error;
3885
3886 // Argumentreduktion 2 wird rueckgaengig gemacht
3887 y = real(zhn1)*p;
3888 if (zwei)
3889 y *= real(zhn2);
3890
3891 // Argumentreduktion 1 wird rueckgaengig gemacht
3892 /* if (neg) y = -(y);*/
3893
3894 stagprec = stagsave;
3895 y = adjust(y);
3896 y = y & einfachgenau;
3897
3898 // Zusaetzliche Fehlerbehandlung noetig, da Underflow-Fehler
3899 // nicht abgefangen werden kann
3900 /*
3901 error = interval(-4.0*pow10(-stagprec*16-1), 4.0*pow10(-stagprec*16-1));
3902 y += error;
3903 */
3904
3905 }
3906 fertig: // Blomquist,05.12.03;
3907 return y;
3908}
3909
3910int cxsc_zweihoch(int n) noexcept
3911{
3912 // liefert 2^n
3913
3914 int res = 1;
3915 switch (n)
3916 {
3917 case 0:
3918 return 1;
3919 case 1:
3920 return 2;
3921 case 2:
3922 return 4;
3923 default:
3924 {
3925 int y = 1,
3926 x = 2,
3927 zhi = 2;
3928
3929 if (n%2)
3930 res = 2;
3931 y = x*x;
3932 do {
3933 if ((n/zhi)%2)
3934 res *= y;
3935 zhi += zhi;
3936 if (zhi <= n)
3937 y *= y;
3938 } while (zhi <= n);
3939 }
3940 }
3941 return res;
3942}
3943
3944l_interval log2(const l_interval & x) // log2(x)
3945{
3946 int stagsave = stagprec,
3947 stagmax = 19;
3948 if (stagprec>stagmax)
3949 stagprec = stagmax;
3950 l_interval y;
3951
3952 y = ln(x)/Ln2_l_interval();
3953
3954 stagprec = stagsave;
3955 y = adjust(y);
3956
3957 return y;
3958}
3959
3960l_interval log10(const l_interval & x) // log10(x)
3961{
3962 int stagsave = stagprec,
3963 stagmax = 19;
3964 if (stagprec>stagmax)
3965 stagprec = stagmax;
3966 l_interval y;
3967
3968 y = ln(x)/Ln10_l_interval();
3969
3970 stagprec = stagsave;
3971 y = adjust(y);
3972
3973 return y;
3974}
3975
3976l_interval sinh(const l_interval & x) // Sinh(x)
3977{
3978 long int n = 1;
3979 int stagsave = stagprec,
3980 stagmax = 19,
3981 sign = 1,
3982 digits = 53,
3983 degree, stagsave2;
3984 real tmp,
3985 lnt, lneps,
3986 fak = 1.0,
3987 lnb = 0.69314718;
3988 l_interval y,
3989 t, t2, p, pot;
3990 interval dy,
3991 dx = interval(x),
3992 einfachgenau,
3993 ibase = interval(2.0),
3994 // ilnb = interval(0.6931471,0.6931472),
3995 // ilnx, ilneps, in,
3996 err;
3997 // ifak = interval(fak);
3998
3999 einfachgenau = sinh(dx);
4000
4001 if (stagprec == 1)
4002 y = sinh(dx);
4003 else if (Inf(dx) == Sup(dx) && Sup(dx) == CXSC_Zero)
4004 y = x;
4005 else
4006 {
4007 if (stagprec < stagmax)
4008 stagprec++;
4009 else
4010 stagprec = stagmax;
4011
4012 // Ausnuztung der Punktsymmetrie
4013 if (Sup(dx) < 0.0)
4014 {
4015 y = -x;
4016 sign = -1;
4017 } else
4018 y = x;
4019
4020 dy = interval(y);
4021
4022 // Bei Werten �ber 0.5 -- Auswertung �ber e-Funktion
4023 if (Sup(dy) > 0.5)
4024 {
4025 try {
4026 t = exp(y);
4027 }
4028 catch(const ERROR_LINTERVAL_FAK_OVERFLOW &)
4029 {
4030 cxscthrow( ERROR_LINTERVAL_FAK_OVERFLOW("l_interval sinh(const l_interval & x)"));
4031 }
4032 y = sign*0.5*(t-1.0/t);
4033 }
4034 // Auswertung �ber Potenzreihenentwicklung
4035 else
4036 {
4037 t = y;
4038 t2 = t*t;
4039
4040 // Absch�tzung der Genauigkeit
4041 tmp = real(Sup(abs(t))); // bei Andreas und Christian "t2"
4042 if (tmp < MinReal)
4043 tmp = MinReal;
4044 lnt = ln(tmp);
4045 lneps = (1.0-digits*stagprec)*lnb;
4046 do {
4047 n += 3;
4048 if (n > 170)
4049 { // 170! = 7.26*10^306
4050 cxscthrow(ERROR_LINTERVAL_FAK_OVERFLOW("l_interval sinh(const l_interval & x)"));
4051 n = 170;
4052 fak *= n*(n-1.0);
4053 break;
4054 }
4055 fak *= n*(n-1.0)*(n-2.0);
4056 } while(lneps+ln(fak)-real(n)*lnt-lnb <= 0.0);
4057 /*
4058 real bas, tn, t2d;
4059 bas = double(real(power(l_real(2.0),1-digits*stagprec)));
4060 tn = abs(double(real(mid(t))));
4061 t2d = double(real(mid(t2)));
4062 do {
4063 n += 2.0;
4064 if (n > 170)
4065 { // 170! = 7.26*10^306
4066 errmon(ERROR_LINTERVAL(FAKOVERFLOW));
4067 n = 170;
4068 fak *= n;
4069 break;
4070 }
4071 tn *= t2d;
4072 fak *= n*(n-1.0);
4073 } while(bas-tn/fak <= 0.0);
4074 */
4075 degree = 2*(n/2+1); // degree ist gerade
4076 // Achtung degree = 2n+2!
4077
4078 // Auswertung mit Horner-Schema
4079 p = 1.0;
4080 for (int i = degree; i >= 2; i -= 2)
4081 {
4082 pot = interval((i+1.0)*i);
4083 p = (p*(t2))/pot+1.0;
4084 }
4085 p *= t;
4086
4087 // Negative Werte zur�ckwandeln
4088 if (sign == -1)
4089 p = -p;
4090
4091 // Fehlerauswertung
4092 stagsave2 = stagprec;
4093 stagprec = 2;
4094 l_interval relerr;
4095 stagprec = stagsave2;
4096 err = pow(ibase , interval(1.0-digits*stagprec)) * interval(-1.0, 1.0);
4097 relerr[1] = 1.0;
4098 relerr[2] = Inf(err);
4099 relerr[3] = Sup(err);
4100
4101 // Fehler einbauen
4102 y = p*relerr;
4103 }
4104 stagprec = stagsave;
4105 y = adjust(y);
4106 y = y & einfachgenau;
4107 }
4108 return y;
4109}
4110
4111l_interval cosh(const l_interval & x) // Cosh(x)
4112{
4113 int stagsave = stagprec,
4114 stagmax = 19;
4115 l_interval y, t;
4116 interval dx = interval(x),
4117 einfachgenau;
4118
4119 einfachgenau = cosh(dx);
4120
4121 if (stagprec == 1)
4122 y = cosh(dx);
4123 else if (Inf(dx) == Sup(dx) && Sup(dx) == CXSC_Zero)
4124 y = adjust(l_interval(1.0));
4125 else
4126 {
4127 if (stagprec < stagmax)
4128 stagprec++;
4129 else
4130 stagprec = stagmax;
4131 if (Sup(dx) <= 0.0)
4132 y = -x;
4133 else
4134 y = x;
4135
4136 try
4137 {
4138 t = exp(y);
4139 }
4140 catch(const ERROR_LINTERVAL_FAK_OVERFLOW &)
4141 {
4142 cxscthrow(ERROR_LINTERVAL_FAK_OVERFLOW("l_interval cosh(const l_interval & x)"));
4143 }
4144
4145 y = 0.5*(t+1.0/t);
4146 if (interval(0.0) <= dx)
4147 SetInf(y,1.0);
4148 stagprec = stagsave;
4149 y = adjust(y);
4150 y = y & einfachgenau;
4151 }
4152
4153 return y;
4154}
4155
4156// Fields for tanh-function
4157int ex_tanh[8] = { 1,-1,-2,-4,-5,-6,-8,-9 };
4158int c_tanhN[8] = { 1,-1,2,-17,62,-1382,21844,-929569 };
4159int c_tanhD[8] = { 1,3,15,315,2835,155925,6081075,638512875 };
4160
4161l_interval tanh(const l_interval & x) noexcept
4162// tanh(x) with Taylor-series and x --> MaxReal; Blomquist 24.05.04;
4163{
4164 l_interval s,t,y;
4165 interval dx = interval(x),
4166 einfachgenau,r,r2;
4167 int stagsave = stagprec,ex,m,N,k,
4168 stagmax = 19;
4169 bool neg;
4170 real Sdx;
4171 einfachgenau = tanh(dx);
4172 if (stagprec>stagmax) stagprec = stagmax;
4173 if (stagprec == 1) y = tanh(dx);
4174 else if (dx==0) y = 0;
4175 else if (0<=dx) y = einfachgenau;
4176 else // 0 not in x and not in dx:
4177 {
4178 t = x;
4179 neg = Sup(dx)<0; Sdx = Sup(dx);
4180 if (neg) {t = -t; dx = -dx; Sdx = Sup(dx);} // Inf(dx) > 0:
4181 ex = expo(Sdx);
4182 if (ex < -70) {
4183 m = ex - 53*stagprec; // m >= -1074 necessary !!
4184 if (m < -1074) m = -1074;
4185 r = interval(Sdx);
4186 r2 = r*r;
4187 N = 0;
4188 do {
4189 N++;
4190 r = r*r2;
4191 k = expo(Sup(r)) + ex_tanh[N];
4192 } while (k>m);
4193 r2 = interval(c_tanhN[N])/c_tanhD[N];
4194 r2 = r * abs(r2); // r2: inclusion of the absolute error
4195 N--; // N: Polynomial degree
4196 y = l_interval(c_tanhN[N])/c_tanhD[N];
4197 s = t*t;
4198 for (k=N-1; k>=0; k--) // Horner-scheme
4199 y = y*s + l_interval(c_tanhN[k])/c_tanhD[k];
4200 y = y * t;
4201 y = y + interval(-Sup(r2),Sup(r2)); // Approxim. error
4202 } else if (ex<-4) {
4203 s = sinh(t); y = cosh(t);
4204 if (stagprec<stagmax) stagprec++;
4205 y = s/y;
4206 } else if (Sdx<352) {
4207 times2pown(t,1);
4208 y = 1 - 2/(exp(t)+1); // --> Sup(y) <= 1 !!
4209 } else {
4210 if (Inf(dx)<352) {
4211 t = Inf(t); times2pown(t,1);
4212 y = 1 - 2/(exp(t)+1);
4213 } else y = l_interval(1) - comp(0.5,-1013);
4214 SetSup(y,1);
4215 }
4216 if (neg) y = -y;
4217 } // else
4218
4219 stagprec = stagsave; // restore old stagprec value
4220 y = adjust(y); // adjust precision of y to old stagprec value
4221 y = y & einfachgenau; // optimal inclusion in case of too
4222 // large relative diameters of y;
4223 return y;
4224} // tanh
4225
4226// ************************************************************************
4227// Fields for coth-function
4228int ex_coth[8] = { -1,-5,-8,-12,-15,-18,-22,-25 };
4229int c_cothN[8] = { 1,-1,2,-1,2,-1382,4,-3617 };
4230real c_cothD[8] = { 1,15,315,1575,31185,212837625,
4231 6081075,54273594375.0 }; // components exactly stored!
4232// ************************************************************************
4233
4234l_interval coth(const l_interval & x) noexcept
4235// coth(x); Blomquist 17.04.04;
4236{
4237 l_interval t,s,c,y;
4238 interval dx = interval(x),
4239 einfachgenau,r,r2;
4240 int stagsave = stagprec,ex,m,N,k,
4241 stagmax = 19;
4242 bool neg;
4243 einfachgenau = coth(dx); // produces all necessary error messages!
4244
4245 if (stagprec == 1) y = coth(dx);
4246 else
4247 {
4248 if (stagprec>stagmax) stagprec = stagmax;
4249 neg = Sup(dx)<0.0;
4250 t = x;
4251 if (neg) { t = -t; dx = -dx; } // Inf(t),Inf(dx) > 0;
4252 ex = expo(Inf(dx));
4253 if (ex<-66) { // Laurent series
4254 m = -ex - 53*stagprec;
4255 r = interval(Sup(dx)); r2 = r*r;
4256 N = 0;
4257 do {
4258 N++;
4259 if (N>7) { y = einfachgenau; goto Fertig; }
4260 r = r*r2;
4261 k = expo(Sup(r)) + ex_coth[N];
4262 } while (k>m);
4263 r2 = interval(c_cothN[N])/c_cothD[N];
4264 r2 = r*abs(r2)/3; // r2: inclusion of the absolute error
4265 N--; // Polynomial degree
4266 y = l_interval(c_cothN[N])/c_cothD[N];
4267 s = t*t;
4268 for (k=N-1; k>=0; k--)
4269 y = y*s + l_interval(c_cothN[k])/c_cothD[k];
4270 y = (y*t)/3 + 1/t;
4271 y = y + interval(-Sup(r2),Sup(r2)); // with approx. error
4272 } else if (ex < 2) {
4273 if (stagprec<stagmax) stagprec++; // stagprec <= 19
4274 c = cosh(t);
4275 s = sinh(t);
4276 y = c/s;
4277 } else if (Inf(dx)<353.0) {
4278 if (stagprec<stagmax) stagprec++; // stagprec <= 19
4279 times2pown(t,1);
4280 y = 1.0 + 2.0/(exp(t)-1.0);
4281 } else { // Inf(dx) >= 353.0
4282 y = l_interval(1.0) + comp(0.5,-1016);
4283 SetInf(y,1.0);
4284 }
4285 if (interval(1.0) <= y) SetInf(y,1.0);
4286 if (neg) y = -y;
4287 Fertig: stagprec = stagsave; // restore old stagprec value
4288 y = adjust(y); // adjust precision of y to old stagprec value
4289 y = y & einfachgenau;
4290 }
4291 return y;
4292} // coth
4293
4294
4295l_interval acosh(const l_interval & x) noexcept
4296// acosh(x); Blomquist 14.04.04;
4297{
4298 int stagsave = stagprec,ex1,ex2,
4299 stagmax = 19;
4300 interval dx = interval(x),
4301 einfachgenau;
4302 l_interval y,t;
4303
4304 einfachgenau = acosh(dx); // false definition range stops program
4305
4306 if (stagprec == 1) y = acosh(dx);
4307 else if (Inf(dx) == Sup(dx) && Sup(dx) == 1.0)
4308 y = adjust(l_interval(0.0));
4309 else
4310 {
4311 if (stagprec < stagmax)
4312 stagprec++;
4313 else stagprec = stagmax;
4314 ex1 = expo(Inf(dx)); ex2 = expo(Sup(dx));
4315 if (ex1>500) { // acosh(x) approx ln(2) + ln(x):
4316 y = li_ln2() + ln(x);
4317 // Absolute approximation error:
4318 t = 1.0/Inf(x);
4319 y += l_interval(-Sup(t*t),0); // approxim. error realized here
4320 } else if (ex2<2) {
4321 t = x-1;
4322 y = lnp1( t+sqrt(t*(2+t)) );
4323 } else y = ln(x+sqrtx2m1(x));
4324
4325 stagprec = stagsave;
4326 y = adjust(y);
4327 y = y & einfachgenau;
4328 }
4329 return y;
4330}
4331
4333
4334// ASinh(x) hier ohne zeitaufwendiges Newton-Verfahren, da die Funktion
4335// sqrt1px2(x) = sqrt(1+x*x) zur Verfuegung steht, die Overflow vermeidet!
4336// Blomquist, 28.12.03;
4337{
4338 int stagsave = stagprec,
4339 stagmax = 19;
4340 l_interval y, my;
4341 interval dx = interval(x), einfachgenau,
4342 error = abs(dx);
4343 real absSupdx = Sup(error), r;
4344
4345 try
4346 {
4347 einfachgenau = asinh(dx);
4348
4349 if (stagprec == 1)
4350 y = asinh(dx);
4351 else if (Inf(dx) == Sup(dx) && Sup(dx) == CXSC_Zero)
4352 y = x;
4353 else // if (absSupdx < MaxReal)
4354 {
4355 if (stagprec < stagmax) stagprec++;
4356 else stagprec = stagmax; // Erhoete Praezision:
4357
4358 if (absSupdx < 2E-108) { // Taylorreihe von ln(x+sqrt(1+x2))
4359 // ln(x+sqrt(1+x2)) = x - x3/6 +- .....
4360 y = x;
4361 dx = interval(absSupdx);
4362 error = dx*dx*dx/6;
4363 r = Sup(error); // r: Oberschranke des absoluten Fehlers
4364 // Jetzt folgt die Addition einer garantierten Einschliessung
4365 // des abs. Fehlers:
4366 y += l_interval(-r,r);
4367 } else
4368 if (Sup(x) < 0.0)
4369 if (absSupdx < 1E+10) y = lnp1( -x+sqrtp1m1(x*x) );
4370 else y = -ln(-x+sqrt1px2(x));
4371 else
4372 if (absSupdx < 1E+10) y = lnp1( x+sqrtp1m1(x*x) );
4373 else y = ln(x+sqrt1px2(x));
4374
4375 stagprec = stagsave;
4376 y = adjust(y);
4377 y = y & einfachgenau;
4378 }
4379 } // try
4380 catch(const ERROR_LINTERVAL_STD_FKT_OUT_OF_DEF &)
4381 {
4382 cxscthrow(ERROR_LINTERVAL_STD_FKT_OUT_OF_DEF("l_interval asinh(const l_interval & x)"));
4383 }
4384 catch(const ERROR_LINTERVAL_FAK_OVERFLOW &)
4385 {
4386 cxscthrow(ERROR_LINTERVAL_FAK_OVERFLOW("l_interval asinh(const l_interval & x)"));
4387 }
4388 return y;
4389} // asinh(x)
4390
4392
4393// ATanh(x), Blomquist 29.12.03;
4394{
4395 int stagsave = stagprec,
4396 stagmax = 19;
4397 l_interval y, my;
4398 interval dx = interval(x),
4399 error = abs(dx),
4400 einfachgenau;
4401 real absSupdx = Sup(error);
4402
4403 einfachgenau = atanh(dx);
4404
4405 if (Inf(x) <= CXSC_MinusOne || Sup(x) >= CXSC_One)
4406 cxscthrow(ERROR_LINTERVAL_STD_FKT_OUT_OF_DEF("l_interval atanh(const l_interval & x)"));
4407 else if (stagprec == 1)
4408 y = atanh(dx);
4409 else if (Inf(dx) == Sup(dx) && Sup(dx) == CXSC_Zero)
4410 y = x;
4411 else
4412 {
4413 if (stagprec < stagmax) stagprec++;
4414 else stagprec = stagmax; // Ab jetzt hoehere Praezision
4415 if (absSupdx < 0.125) {
4416 y = x / (1-x);
4417 times2pown(y,1); // Schnelle Multiplikation mit 2
4418 y = lnp1(y);
4419 times2pown(y,-1); // Schnelle Division durch 2
4420 } else {
4421 y = ln((1.0+x)/(1.0-x));
4422 times2pown(y,-1);
4423 }
4424 stagprec = stagsave;
4425 y = adjust(y);
4426 y = y & einfachgenau;
4427 }
4428
4429 return y;
4430} // atanh
4431
4433
4434// Acoth(x) hier ohne das langsame Newtonverfahren; Fue grosse x gilt:
4435// Acoth = 0.5*ln(1+2/(x-1)); Blomquist, 28.12.03;
4436{
4437 int stagsave = stagprec,
4438 stagmax = 19;
4439 l_interval y, my;
4440 interval dx = interval(x),
4441 einfachgenau,
4442 error = abs(dx);
4443 real absSupdx = Sup(error);
4444
4445 einfachgenau = acoth(dx);
4446
4447 if ((l_interval(Inf(x))) < l_interval(CXSC_MinusOne,CXSC_One)
4448 || (l_interval(Sup(x))) < l_interval(CXSC_MinusOne,CXSC_One))
4449 cxscthrow(ERROR_LINTERVAL_STD_FKT_OUT_OF_DEF
4450 ("l_interval acoth(const l_interval & x)"));
4451 else if (stagprec == 1)
4452 y = acoth(dx);
4453 else
4454 {
4455 if (stagprec < stagmax) stagprec++;
4456 else stagprec = stagmax; // Rechg. in erhoeter Praezision:
4457
4458 if (absSupdx > 1E10) {
4459 y = l_interval(2)/(x-1);
4460 y = lnp1(y);
4461 times2pown(y,-1);
4462 } else {
4463 y = ln((x+1.0)/(x-1.0));
4464 times2pown(y,-1);
4465 }
4466
4467 stagprec = stagsave;
4468 y = adjust(y);
4469 y = y & einfachgenau;
4470 }
4471
4472 return y;
4473} // acoth
4474
4475l_interval lnp1(const l_interval& x) noexcept
4476// ln(1+x) = zeta * P(zeta); zeta := x/(2+x);
4477// P(zeta) := 2 + (2/3)zeta^2 + (2/5)zeta^4 + ...
4478// Approximation of P(zeta) with P_N(zeta), defined by
4479// P_N(zeta) := 2 + (2/3)zeta^(2*1) +
4480// + (2/5)zeta^(2*2) + ... + (2/2N+1)zeta^(2*N);
4481// Absolute approximation error Delta := P(zeta) - P_N(zeta);
4482// With z:=zeta^2 it holds: Delta <= (2*z^(N+1))/{(2N+3)*(1-z)};
4483{
4484 const real c1(0.693147181);
4485 int stagsave(stagprec),
4486 stagmax(19); // For ln-function
4487 if (stagprec>stagmax) stagprec = stagmax; // stagprec now <= 19;
4488 interval dx(x), // dx of type interval inclusing x.
4489 einfachgenau(lnp1(dx)), // einfachgenau inclusing lnp1(x).
4490 ax( abs(dx) ); // ax inclusing the absolute values of dx.
4491 real gr(Sup(ax)), // gr: maximum of the absolute values from x.
4492 lngr,u;
4493 l_interval t(x);
4494 if (gr < 1E-8) // Series expansion, gr = 1E-8 is Ok!
4495 if (sign(gr)==0) t=0;
4496 else { // Calculation of N for the approximation of P(zeta)
4497 // using the polynomial P_N(zeta) of degree N:
4498 int N(0),TwoNp3(3),k;
4499 k = expo(gr);
4500 if (k>-1019) {
4501 k = k - 1 - 53*stagprec;
4502 lngr = ln(gr/2);
4503 if (k < -1074) k = -1074;
4504 k--;
4505 u = k*c1; // u = (k-1)*ln(2);
4506 do {
4507 N++;
4508 TwoNp3 += 2;
4509 } while(TwoNp3*lngr - ln(TwoNp3) > u);
4510 }
4511 // Polynomial degree 0<= N <=17 is now calculated!
4512 // Calculation of the N+1 polynomial coefficients:
4513 l_interval a_lnp1[18], // Declaration of the coeffits.
4514 zeta,z2;
4515 a_lnp1[0] = 2;
4516 for (int i=1; i<=N; ++i) {
4517 a_lnp1[i] = a_lnp1[0]/(2*i + 1);
4518 }
4519 // Polyn.-coeff. a_lnp1[i] are now calculated; 0<=i<=N;
4520 // The calculation of P_N(zeta) follows now:
4521 zeta = t/(2+t);
4522 z2 = sqr(zeta);
4523 dx = z2; // is used for the approximation error!
4524 t = a_lnp1[N]; // Horner-scheme:
4525 for (int i=N-1; i>=0; --i) {
4526 t = t*z2 + a_lnp1[i];
4527 } // t = P_N(zeta) is now calculated
4528 // Calculating now the approximation error
4529 // with dx of type interval:
4530 dx = interval(Sup(dx));
4531 ax = Power(dx,N+1);
4532 times2pown(ax,1); // Multiplication with 2;
4533 dx = ax / ((2*N+3)*(1-dx));
4534 t = t + l_interval(0.0,Sup(dx));
4535 // Approximation error implemented now
4536 t = zeta * t;
4537 }
4538 else
4539 if (gr<1) { // Calculation in higher accuracy:
4540 stagprec++;
4541 if (stagprec>stagmax) stagprec = stagmax;
4542 t = ln(1+t);
4543 } else t = ln(1+t); // Normal accuracy:
4544 stagprec = stagsave;
4545 t = adjust(t);
4546 t = t & einfachgenau; // intersection of t and einfachgenau;
4547 // If x is too wide and therefore t because of the inevitable
4548 // interval overestimations, t is the best possible inclusion.
4549 return t;
4550}
4551
4553// sqrt(x^2-1); Blomquist, 13.04.04;
4554{
4555 int stagsave = stagprec,
4556 stagmax = 19;
4557 l_interval y,z = abs(x);
4558 l_real r1,r2;
4559 interval dx = interval(z),
4560 einfachgenau;
4561 einfachgenau = sqrtx2m1(dx);
4562 if (stagprec>stagmax) stagprec = stagmax;
4563
4564 if (stagprec == 1) y = sqrtx2m1(dx); // simple interval
4565 else if (Inf(z)==1) {
4566 r1=0; z = Sup(z); z = sqrt(z*z-1);
4567 r2 = Sup(z);
4568 y = l_interval(r1,r2);
4569 } else if (expo(Sup(dx))<500) y = sqrt(z*z-1.0); // no overflow!
4570 else { // avoiding overflow using: x-1/(2x) < sqrt(x^2-1) < x
4571 r2 = Sup(z);
4572 z = Inf(z); // z: Point interval
4573 y = 1.0/z;
4574 times2pown(y,-1);
4575 y = z - y;
4576 r1 = Inf(y);
4577 y = l_interval(r1,r2);
4578 }
4579 stagprec = stagsave;
4580 y = adjust(y);
4581 y = y & einfachgenau;
4582 return y;
4583} // sqrtx2m1
4584
4586// sqrt(1-x^2); Blomquist, 13.04.04;
4587{
4588 int stagsave = stagprec,
4589 stagmax = 19;
4590 l_interval y,z = abs(x);
4591 l_real r1,r2;
4592 interval dx = interval(z),
4593 einfachgenau;
4594 einfachgenau = sqrt1mx2(dx);
4595 if (stagprec>stagmax) stagprec = stagmax;
4596
4597 if (stagprec == 1) y = sqrt1mx2(dx); // simple interval
4598 else {
4599 y = comp(0.5,1023); // y = 2^(+1022)
4600 times2pown(z,511);
4601 y = sqrt(y-z*z);
4602 times2pown(y,-511);
4603 }
4604 stagprec = stagsave;
4605 y = adjust(y);
4606 y = y & einfachgenau;
4607 return y;
4608} // sqrt1mx2
4609
4610
4611l_interval ln_sqrtx2y2(const l_interval& x, const l_interval& y) noexcept
4612{ // Inclusion of ln(sqrt{x^2+y^2}) in staggered arithmetic
4613 int stagsave = stagprec;
4614// stagmax = 19;
4615 interval dx = interval(x),
4616 dy = interval(y),
4617 einfachgenau = ln_sqrtx2y2(dx,dy);
4618 dx = abs(dx); dy = abs(dy);
4619 l_interval ax(abs(x)),ay(abs(y)),ar;
4620 int ex_x(expo(Sup(dx))), ex_y(expo(Sup(dy))),
4621 N,N1;
4622 if (ex_y>ex_x) ex_x = ex_y;
4623 if (ex_x>508) { // To avoid overflow:
4624 N = ex_x-500;
4625 times2pown(ax,-N); times2pown(ay,-N);
4626 ar = ax*ax + ay*ay;
4627 ar = ln(ar);
4628 times2pown(ar,-1);
4629 ar += N*li_ln2();
4630 } else
4631 if (ex_x<-20) { // To avoid underflow:
4632 N = 500 - ex_x;
4633 if (N>1023) { // Avoiding an range error with times2pown(...)
4634 N1 = N-1023;
4635 times2pown(ax,1023); times2pown(ax,N1);
4636 times2pown(ay,1023); times2pown(ay,N1);
4637 } else {
4638 times2pown(ax,N);
4639 times2pown(ay,N);
4640 }
4641 ar = ax*ax + ay*ay;
4642 ar = ln(ar);
4643 times2pown(ar,-1);
4644 ar -= N*li_ln2(); // ar = ar - N*ln(2)
4645 } else { // Calculation with function lnp1:
4646 ar = sqr(ax) + sqr(ay) - 1;
4647 ar = lnp1(ar);
4648 times2pown(ar,-1);
4649 }
4650 stagprec = stagsave;
4651 ar = adjust(ar);
4652 ar = ar & einfachgenau;
4653
4654 return ar;
4655}
4656
4658// acoshp1(x) = acosh(1+x); Blomquist, 20.04.05;
4659{
4660 int stagsave = stagprec,ex,
4661 stagmax = 19;
4662 l_interval y,t;
4663 l_real lr;
4664 interval dx = interval(x),
4665 einfachgenau;
4666 einfachgenau = acoshp1(dx);
4667 if (stagprec>stagmax) stagprec = stagmax;
4668
4669 if (stagprec == 1) y = acoshp1(dx); // simple interval
4670 else {
4671 ex = expo(Sup(dx));
4672 if (ex<=-1016) {
4673 t = x;
4674 times2pown(t,1); // fast multiplication with 2 = 2^1;
4675 t = sqrt(t); // t = sqrt(2x);
4676 lr = Inf(t); // Calculating the infimum (Leibniz-series)
4677 y = l_interval(lr)*(1.0 - l_interval(lr)/12.0);
4678 y = l_interval(Inf(y),Sup(t)); // using alternating Leibniz-series
4679 }
4680 else
4681 if (ex<-400) y = lnp1(x*(1.0+sqrt(1+2.0/x)));
4682 else y = acosh(1+x); // x = MaxReal not allowed
4683 }
4684 stagprec = stagsave;
4685 y = adjust(y);
4686 y = y & einfachgenau;
4687 return y;
4688} // acoshp1
4689
4690
4691} // namespace cxsc
4692
The Scalar Type interval.
Definition interval.hpp:55
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
l_interval Sqrt3r_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:2455
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
l_interval SqrtPi_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:1863
const real MinReal
Smallest normalized representable floating-point number.
Definition real.cpp:62
cinterval asinh(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:2718
l_interval Ep2r_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:2899
cinterval coth(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:578
l_interval Pi2_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:1569
l_interval Sqrt2Pir_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:2085
l_interval Sqrt2Pi_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:1937
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
l_interval Sqrt2_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:1135
cinterval power(const cinterval &z, int n) noexcept
Calculates .
Definition cimath.cpp:1941
cinterval log10(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:903
l_interval Ln10r_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:995
l_interval Ep2_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:2825
l_interval EpPi_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:2973
cinterval ln(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:851
cvector diam(const cimatrix_subv &mv) noexcept
Returns the diameter of the matrix.
Definition cimatrix.inl:738
l_interval LnPi_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:2529
cinterval pow(const cinterval &z, const interval &p) noexcept
Calculates .
Definition cimath.cpp:2074
l_interval Ep2Pi_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:3047
l_interval Pir_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:1716
l_interval Pip2_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:2159
interval expmx2(const interval &x)
Calculates .
Definition imath.cpp:192
cinterval sinh(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:231
l_interval li_pi4()
Enclosure-Interval for .
Definition l_imath.hpp:163
l_interval Catalan_l_interval() noexcept
Enclosure-Interval for Catalan Numbers.
Definition l_imath.cpp:3343
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
interval acoshp1(const interval &x)
Calculates .
Definition imath.cpp:617
cinterval tan(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:393
cinterval exp10(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:172
l_interval Sqrt7_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:1276
interval Power(const interval &x, int n)
Calculates .
Definition imath.cpp:536
l_interval E_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:2677
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
l_interval EulerGa_l_interval() noexcept
Enclosure-Interval for Euler Gamma.
Definition l_imath.cpp:3269
cinterval acosh(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:2732
l_interval Pid3_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:1642
l_interval Sqrt3_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:2307
l_interval SqrtPir_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:2011
cinterval cosh(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:223
l_interval Ln2r_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:1349
l_interval EpPid2_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:3121
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
l_interval Sqrt3d2_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:2381
cinterval exp(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:159
l_interval EpPid4_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:3195
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
cinterval expm1(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:177
cinterval cot(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:538
l_interval Er_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:2751
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 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
l_interval Ln2Pi_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:2603
cinterval sqr(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:3342
l_interval Sqrt5_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:1204
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
l_interval li_ln2()
Enclosure-Interval for .
Definition l_imath.hpp:157
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
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
l_interval Pi2r_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:1789