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 
31 namespace cxsc {
32 
33 #define CXSC_One 1.
34 #define CXSC_Zero 0.
35 #define CXSC_MinusOne -1.
36 
37 l_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 
91 l_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 
153 l_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 
232 l_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 
310 l_interval sqrt1px2(const l_interval& x) noexcept
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 
343 l_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 
397 l_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 
441 l_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 
613 l_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 
657 l_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 
705 l_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 
767 l_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 
850 static real CXSC_ln2[21]; // CXSC_ln2[0], ... CXSC_ln2[20]
851 static 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 
921 static real CXSC_ln10[21]; // CXSC_ln10[0], ... CXSC_ln10[20]
922 static 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 
991 static real CXSC_Rln10[21]; // CXSC_Rln10[0], ... CXSC_Rln10[20]
992 static 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 
1061 static real CXSC_pi4[21];
1062 static 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 
1131 static real CXSC_sqrt2[21];
1132 static 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 
1201 static real CXSC_sqrt5[21]; // CXSC_sqrt5[0], ... CXSC_sqrt5[20]
1202 static 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 
1273 static real CXSC_sqrt7[21]; // CXSC_sqrt7[0], ... CXSC_sqrt7[20]
1274 static 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 
1346 static real CXSC_ln2r[21]; // CXSC_ln2r[0], ... CXSC_ln2r[20]
1347 static 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 
1420 static real CXSC_Pi[21]; // CXSC_Pi[0], ... CXSC_Pi[20]
1421 static 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 
1493 static real CXSC_Pid2[21]; // CXSC_Pid2[0], ... CXSC_Pid2[20]
1494 static 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 
1566 static real CXSC_Pi2[21]; // CXSC_Pi2[0], ... CXSC_Pi2[20]
1567 static 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 
1639 static real CXSC_Pid3[21]; // CXSC_Pid3[0], ... CXSC_Pid3[20]
1640 static 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 
1713 static real CXSC_Pir[21]; // CXSC_Pir[0], ... CXSC_Pir[20]
1714 static 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 
1786 static real CXSC_Pi2r[21]; // CXSC_Pi2r[0], ... CXSC_Pi2r[20]
1787 static 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 
1860 static real CXSC_SqrtPi[21]; // CXSC_SqrtPi[0], ... CXSC_SqrtPi[20]
1861 static 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 
1934 static real CXSC_Sqrt2Pi[21]; // CXSC_Sqrt2Pi[0], ... CXSC_Sqrt2Pi[20]
1935 static 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 
2008 static real CXSC_SqrtPir[21]; // CXSC_SqrtPir[0], ... CXSC_SqrtPir[20]
2009 static 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 
2082 static real CXSC_Sqrt2Pir[21]; // CXSC_Sqrt2Pir[0], ... CXSC_Sqrt2Pir[20]
2083 static 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 
2156 static real CXSC_Pip2[21]; // CXSC_Pip2[0], ... CXSC_Pip2[20]
2157 static 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 
2230 static real CXSC_Sqrt2r[21]; // CXSC_Sqrt2r[0], ... CXSC_Sqrt2r[20]
2231 static 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 
2304 static real CXSC_Sqrt3[21]; // CXSC_Sqrt3[0], ... CXSC_Sqrt3[20]
2305 static 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 
2378 static real CXSC_Sqrt3d2[21]; // CXSC_Sqrt3d2[0], ... CXSC_Sqrt3d2[20]
2379 static 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 
2452 static real CXSC_Sqrt3r[21]; // CXSC_Sqrt3r[0], ... CXSC_Sqrt3r[20]
2453 static 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 
2526 static real CXSC_LnPi[21]; // CXSC_LnPi[0], ... CXSC_LnPi[20]
2527 static 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 
2600 static real CXSC_Ln2Pi[21]; // CXSC_Ln2Pi[0], ... CXSC_Ln2Pi[20]
2601 static 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 
2674 static real CXSC_E[21]; // CXSC_E[0], ... CXSC_E[20]
2675 static 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 
2748 static real CXSC_Er[21]; // CXSC_Er[0], ... CXSC_Er[20]
2749 static 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 
2822 static real CXSC_Ep2[21]; // CXSC_Ep2[0], ... CXSC_Ep2[20]
2823 static 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 
2896 static real CXSC_Ep2r[21]; // CXSC_Ep2r[0], ... CXSC_Ep2r[20]
2897 static 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 
2970 static real CXSC_EpPi[21]; // CXSC_EpPi[0], ... CXSC_EpPi[20]
2971 static 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 
3044 static real CXSC_Ep2Pi[21]; // CXSC_Ep2Pi[0], ... CXSC_Ep2Pi[20]
3045 static 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 
3118 static real CXSC_EpPid2[21]; // CXSC_EpPid2[0], ... CXSC_EpPid2[20]
3119 static 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 
3192 static real CXSC_EpPid4[21]; // CXSC_EpPid4[0], ... CXSC_EpPid4[20]
3193 static 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 
3266 static real CXSC_EulerGa[21]; // CXSC_EulerGa[0], ... CXSC_EulerGa[20]
3267 static 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 
3340 static real CXSC_Catalan[21]; // CXSC_Catalan[0], ... CXSC_Catalan[20]
3341 static 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 
3412 l_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 
3510 l_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 
3538 l_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 
3640 l_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 
3656 l_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 
3672 l_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 
3758 int cxsc_zweihoch(int) noexcept;
3759 
3760 l_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 
3910 int 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 
3944 l_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 
3960 l_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 
3976 l_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 
4111 l_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
4157 int ex_tanh[8] = { 1,-1,-2,-4,-5,-6,-8,-9 };
4158 int c_tanhN[8] = { 1,-1,2,-17,62,-1382,21844,-929569 };
4159 int c_tanhD[8] = { 1,3,15,315,2835,155925,6081075,638512875 };
4160 
4161 l_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
4228 int ex_coth[8] = { -1,-5,-8,-12,-15,-18,-22,-25 };
4229 int c_cothN[8] = { 1,-1,2,-1,2,-1382,4,-3617 };
4230 real c_cothD[8] = { 1,15,315,1575,31185,212837625,
4231  6081075,54273594375.0 }; // components exactly stored!
4232 // ************************************************************************
4233 
4234 l_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 
4295 l_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 
4475 l_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 
4611 l_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.
Definition: l_interval.hpp:72
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 Sqrt2_l_interval() noexcept
Enclosure-Interval for .
Definition: l_imath.cpp:1135
l_interval Pid3_l_interval() noexcept
Enclosure-Interval for .
Definition: l_imath.cpp:1642
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 Sqrt3_l_interval() noexcept
Enclosure-Interval for .
Definition: l_imath.cpp:2307
l_interval Pi2_l_interval() noexcept
Enclosure-Interval for .
Definition: l_imath.cpp:1569
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 Pid4_l_interval() noexcept
Enclosure-Interval for .
Definition: l_imath.cpp:1065
l_interval Sqrt5_l_interval() noexcept
Enclosure-Interval for .
Definition: l_imath.cpp:1204
l_interval Sqrt3d2_l_interval() noexcept
Enclosure-Interval for .
Definition: l_imath.cpp:2381
cinterval coth(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:578
l_interval Ln2r_l_interval() noexcept
Enclosure-Interval for .
Definition: l_imath.cpp:1349
l_interval Pid2_l_interval() noexcept
Enclosure-Interval for .
Definition: l_imath.cpp:1496
cinterval log2(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:898
l_interval SqrtPir_l_interval() noexcept
Enclosure-Interval for .
Definition: l_imath.cpp:2011
cinterval power(const cinterval &z, int n) noexcept
Calculates .
Definition: cimath.cpp:1941
l_interval Pip2_l_interval() noexcept
Enclosure-Interval for .
Definition: l_imath.cpp:2159
cinterval log10(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:903
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 SqrtPi_l_interval() noexcept
Enclosure-Interval for .
Definition: l_imath.cpp:1863
cinterval pow(const cinterval &z, const interval &p) noexcept
Calculates .
Definition: cimath.cpp:2074
l_interval EpPid2_l_interval() noexcept
Enclosure-Interval for .
Definition: l_imath.cpp:3121
l_interval Catalan_l_interval() noexcept
Enclosure-Interval for Catalan Numbers.
Definition: l_imath.cpp:3343
l_interval EpPi_l_interval() noexcept
Enclosure-Interval for .
Definition: l_imath.cpp:2973
interval expmx2(const interval &x)
Calculates .
Definition: imath.cpp:192
cinterval sinh(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:231
l_interval EulerGa_l_interval() noexcept
Enclosure-Interval for Euler Gamma.
Definition: l_imath.cpp:3269
l_interval li_pi4()
Enclosure-Interval for .
Definition: l_imath.hpp:163
l_interval Ln2_l_interval() noexcept
Enclosure-Interval for .
Definition: l_imath.cpp:854
cinterval asin(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:2311
interval acoshp1(const interval &x)
Calculates .
Definition: imath.cpp:617
cinterval tan(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:393
cinterval exp10(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:172
interval Power(const interval &x, int n)
Calculates .
Definition: imath.cpp:536
l_interval Sqrt2Pi_l_interval() noexcept
Enclosure-Interval for .
Definition: l_imath.cpp:1937
l_interval EpPid4_l_interval() noexcept
Enclosure-Interval for .
Definition: l_imath.cpp:3195
cinterval acos(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:2553
cinterval sqrtx2m1(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:1109
cinterval acosh(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:2732
l_interval Pi2r_l_interval() noexcept
Enclosure-Interval for .
Definition: l_imath.cpp:1789
l_interval E_l_interval() noexcept
Enclosure-Interval for .
Definition: l_imath.cpp:2677
cinterval cosh(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:223
l_interval Ln10_l_interval() noexcept
Enclosure-Interval for .
Definition: l_imath.cpp:925
cinterval cos(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:207
l_interval Ln2Pi_l_interval() noexcept
Enclosure-Interval for .
Definition: l_imath.cpp:2603
cinterval sqrt1px2(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:1071
l_interval Ep2Pi_l_interval() noexcept
Enclosure-Interval for .
Definition: l_imath.cpp:3047
l_interval LnPi_l_interval() noexcept
Enclosure-Interval for .
Definition: l_imath.cpp:2529
cinterval exp(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:159
l_interval Ep2r_l_interval() noexcept
Enclosure-Interval for .
Definition: l_imath.cpp:2899
l_interval Sqrt2r_l_interval() noexcept
Enclosure-Interval for .
Definition: l_imath.cpp:2233
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 Ep2_l_interval() noexcept
Enclosure-Interval for .
Definition: l_imath.cpp:2825
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
cinterval sqr(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:3342
cinterval lnp1(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:867
l_interval Er_l_interval() noexcept
Enclosure-Interval for .
Definition: l_imath.cpp:2751
cvector mid(const cimatrix_subv &mv) noexcept
Returns the middle of the matrix.
Definition: cimatrix.inl:739
l_interval Pi_l_interval() noexcept
Enclosure-Interval for .
Definition: l_imath.cpp:1423
l_interval li_ln2()
Enclosure-Interval for .
Definition: l_imath.hpp:157
l_interval Pir_l_interval() noexcept
Enclosure-Interval for .
Definition: l_imath.cpp:1716
cinterval atan(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:2938
cinterval atanh(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:3317
l_interval Sqrt2Pir_l_interval() noexcept
Enclosure-Interval for .
Definition: l_imath.cpp:2085
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 Sqrt3r_l_interval() noexcept
Enclosure-Interval for .
Definition: l_imath.cpp:2455
l_interval Ln10r_l_interval() noexcept
Enclosure-Interval for .
Definition: l_imath.cpp:995
l_interval Sqrt7_l_interval() noexcept
Enclosure-Interval for .
Definition: l_imath.cpp:1276