C-XSC - A C++ Class Library for Extended Scientific Computing  2.5.4
l_rmath.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_rmath.cpp,v 1.28 2014/01/30 17:23:46 cxsc Exp $ */
25 
26 #include <l_rmath.hpp>
27 #include <rmath.hpp>
28 #include <l_interval.hpp>
29 
30 namespace cxsc {
31 
32 l_real sqrt(const l_real &x)
33 // Blomquist, additional scaling, 10.12.02
34 {
35  int stagsave = stagprec, stagmax = 19, stagcalc;
36  l_real y;
37  if (sign(x[1])<0) // Blomquist: is faster than x < 0.0;
38  cxscthrow(ERROR_LREAL_STD_FKT_OUT_OF_DEF("l_real sqrt(const l_real &x)"));
39  else if (zero_(x) || x == 1.0 ) // Blomquist, faster than x == 0;
40  y = x;
41  else
42  {
43  l_real x1 = x;
44  int ex = expo(x1[1]);
45  ex = 1021-ex; // ex: optimal scaling factor
46  if (ex%2) ex--; // ex is now even;
47  times2pown(x1,ex); // scaling x1 with 2^ex
48 
49  // Einsatz des Newton-Verfahrens y = y-f(y)/f'(y)
50  if (stagprec < stagmax) stagcalc = stagprec+1;
51  else stagcalc = stagmax+1; // Blomquist with +1; 30.11.02;
52  y = sqrt(_real(x1));
53  stagprec = 1;
54  while (stagprec < stagcalc)
55  {
56  stagprec += stagprec; // now calculation in double precision:
57  if (stagprec > stagcalc) stagprec = stagcalc;
58  // y = 0.5*((x1/y)+y);
59  y += x1/y; // Blomquist, 30.11.02;
60  times2pown(y,-1); // Blomquist, this multiplication is faster!!
61  }
62  times2pown(y,-ex/2); // backscaling
63  stagprec = stagsave; // restore old stagprec
64  adjust(y);
65  }
66  return y;
67 }
68 
70  const l_real& y) noexcept // Blomquist 6.12.02
71 // Calculation of an approximation of sqrt(x^2+y^2).
72 // In general the maximum precision is stagprec=19, predifined by the used
73 // sqrt-function declared in l_rmath.hpp.
74 // If the difference |exa-exb| of the exponents (base 2) is sufficient high,
75 // precision and accuracy can be choosen greater 19.
76 {
77  int stagsave = stagprec, stagcalc, stagmax = 19;
78  l_real a,b,r,r1;
79  int exa,exb,ex;
80  a = x; b = y;
81  exa = expo(a[1]);
82  exb = expo(b[1]);
83  if (exb > exa)
84  { // Permutation of a,b:
85  r = a; a = b; b = r;
86  ex = exa; exa = exb; exb = ex;
87  } // |a| >= |b|
88  if (sign(a[1]) < 0) a = -a; // a == abs(a);
89  if (sign(b[1]==0)) return a;
90 
91  // |a| = a >= |b| > 0:
92  ex = 0; // initialization
93  if (6*exb < 5*exa-1071)
94  { // approximation with a + 0.5*b^2/a - b^4/(8*a^3), if the next
95  // Taylor term b^6/(16*a^5) is less than 2^-1074==minreal;
96  // minreal is the smallest positive number in IEEE-format.
97  r1 = (b/a);
98  r = r1*r1;
99  times2pown(r,-2); // r = r/4;
100  r = 1 - r;
101  r1 *= b;
102  times2pown(r1,-1); // r = r/2;
103  r *= r1;
104  r += a; // r = a + 0.5*b^2/a - b^4/(8*a^3);
105  } else
106  { // Scaling ubwards or downwards to achiev optimal accuracy and
107  // to avoid overflow by accu-reading:
108  stagcalc = stagprec;
109  if (stagcalc > stagmax) stagcalc = stagmax; // using sqrt(...) a
110  stagprec = stagcalc; // higher precision than stagmax makes no sense!
111 
112  ex = 511 - exa; // scaling factor to avoide overflow by accu-reading
113  if (ex < 0)
114  { // sqrt(a^2+b^2) = a*sqrt( 1+(b/a)^2 )
115  r = b/a;
116  r = a*sqrt(1+r*r); // sqrt(...) declared in l_rmath.hpp
117  } else
118  {
119  times2pown(a,ex); // exact scaling with ex >= 0
120  times2pown(b,ex); // exact scaling eith ex >= 0
121  dotprecision dot(0.0);
122  accumulate(dot,a,a);
123  accumulate(dot,b,b);
124  r = dot; // r with no overflow!
125  r = sqrt(r); // sqrt(...) declared in l_rmath.hpp
126  times2pown(r,-ex); // back-scaling
127  }
128  stagprec = stagsave; // restore old stagprec
129  }
130  return r;
131 } // sqrtx2y2(...)
132 
133 l_real sqrt1px2(const l_real& x) noexcept
134 // Inclusion of sqrt(1+x^2); Blomquist, 13.12.02;
135 // With stagmax=19 we get about 16*19=304 exact decimal digits.
136 {
137  l_real y,t;
138  int stagsave, stagmax=19;
139  stagsave = stagprec;
140  if (stagprec > stagmax) stagprec = stagmax;
141  if (expo(x[1]) > 260)
142  { // sqrt(1+x^2) = |x| + 1/(2*|x|)
143  y = abs(x);
144  t = 1/y;
145  times2pown(t,-1);
146  y += t;
147  } else y = sqrt(1+x*x);
148  stagprec = stagsave;
149  y = adjust(y);
150  return y;
151 }
152 
153 l_real power(const l_real& x, int n)
154 {
155  int stagsave = stagprec, stagmax = 19;
156  long int zhi = 2;
157  l_real y, neu;
158  bool neg=false;
159 
160  if (x == 1.0)
161  y = x;
162  else if (n == 0)
163  y = adjust(_l_real(1.0));
164  else
165  {
166  if (stagprec < stagmax)
167  stagprec++;
168  else
169  stagprec = stagmax;
170 
171  if (n == 1)
172  y = x;
173  else if (n == 2)
174  y = x*x;
175  else
176  {
177  if (n < 0)
178  {
179  neg = true;
180  n = -n;
181  }
182 
183  // Initialisierung
184  if (n%2)
185  y = x;
186  else
187  y = 1.0; // Praezision wird bei 1 Mult. auf
188  // aktuellen Wert gesetzt;
189  // Berechnung durch binaere Darstellung der n
190  neu = x*x;
191  do {
192  if ((n/zhi)%2)
193  y *= neu;
194  zhi += zhi;
195  if (zhi <= n) // letzte Mult. entfaellt --> schneller
196  neu *= neu;
197  } while (zhi <= n);
198  if (neg)
199  y = 1/(y);
200  }
201  stagprec = stagsave;
202  y = adjust(y);
203  }
204  return y;
205 }
206 
207 // real staggered constants (the same as in l_interval.hpp):
208 l_real Ln2_l_real() noexcept // ln(2)
209 { return mid( Ln2_l_interval() ); }
210 l_real Ln10_l_real() noexcept // ln(10)
211 { return mid( Ln10_l_interval() ); }
212 l_real Ln10r_l_real() noexcept // 1/ln(10)
213 { return mid( Ln10r_l_interval()); }
214 l_real Pid4_l_real() noexcept // Pi/4
215 { return mid( Pid4_l_interval() ); }
216 l_real Sqrt2_l_real() noexcept // sqrt(2)
217 { return mid( Sqrt2_l_interval() ); }
218 l_real Sqrt5_l_real() noexcept // sqrt(5)
219 { return mid( Sqrt5_l_interval() ); }
220 l_real Sqrt7_l_real() noexcept // sqrt(7)
221 { return mid( Sqrt7_l_interval() ); }
222 l_real Ln2r_l_real() noexcept // 1/ln(2)
223 { return mid( Ln2r_l_interval() ); }
224 l_real Pi_l_real() noexcept // Pi
225 { return mid( Pi_l_interval() ); }
226 l_real Pid2_l_real() noexcept // Pi/2
227 { return mid( Pid2_l_interval() ); }
228 l_real Pi2_l_real() noexcept // 2*Pi
229 { return mid( Pi2_l_interval() ); }
230 l_real Pid3_l_real() noexcept // Pi/3
231 { return mid( Pid3_l_interval() ); }
232 l_real Pir_l_real() noexcept // 1/Pi
233 { return mid( Pir_l_interval() ); }
234 l_real Pi2r_l_real() noexcept // 1/(2*Pi)
235 { return mid( Pi2r_l_interval() ); }
236 l_real SqrtPi_l_real() noexcept // sqrt(Pi)
237 { return mid( SqrtPi_l_interval() ); }
238 l_real Sqrt2Pi_l_real() noexcept // sqrt(2*Pi)
239 { return mid( Sqrt2Pi_l_interval() ); }
240 l_real SqrtPir_l_real() noexcept // 1/sqrt(Pi)
241 { return mid( SqrtPir_l_interval() ); }
242 l_real Sqrt2Pir_l_real() noexcept // 1/sqrt(2*Pi)
243 { return mid( Sqrt2Pir_l_interval() ); }
244 l_real Pip2_l_real() noexcept // Pi^2
245 { return mid( Pip2_l_interval() ); }
246 l_real Sqrt2r_l_real() noexcept // 1/sqrt(2)
247 { return mid( Sqrt2r_l_interval() ); }
248 l_real Sqrt3_l_real() noexcept // sqrt(3)
249 { return mid( Sqrt3_l_interval() ); }
250 l_real Sqrt3d2_l_real() noexcept // sqrt(3)/2
251 { return mid( Sqrt3d2_l_interval() ); }
252 l_real Sqrt3r_l_real() noexcept // 1/sqrt(3)
253 { return mid( Sqrt3r_l_interval() ); }
254 l_real LnPi_l_real() noexcept // ln(Pi)
255 { return mid( LnPi_l_interval() ); }
256 l_real Ln2Pi_l_real() noexcept // ln(2*Pi)
257 { return mid( Ln2Pi_l_interval() ); }
258 l_real E_l_real() noexcept // e = exp(1)
259 { return mid( E_l_interval() ); }
260 l_real Er_l_real() noexcept // 1/e
261 { return mid( Er_l_interval() ); }
262 l_real Ep2_l_real() noexcept // e^2
263 { return mid( Ep2_l_interval() ); }
264 l_real Ep2r_l_real() noexcept // 1/e^2
265 { return mid( Ep2r_l_interval() ); }
266 l_real EpPi_l_real() noexcept // e^Pi
267 { return mid( EpPi_l_interval() ); }
268 l_real Ep2Pi_l_real() noexcept // e^(2*Pi)
269 { return mid( Ep2Pi_l_interval() ); }
270 l_real EpPid2_l_real() noexcept // e^(Pi/2)
271 { return mid( EpPid2_l_interval() ); }
272 l_real EpPid4_l_real() noexcept // e^(Pi/4)
273 { return mid( EpPid4_l_interval() ); }
274 l_real EulerGa_l_real() noexcept // EulerGamma
275 { return mid(EulerGa_l_interval() ); }
276 l_real Catalan_l_real() noexcept // Catalan
277 { return mid( Catalan_l_interval() ); }
278 
279 } // namespace cxsc
280 
The Data Type dotprecision.
Definition: dot.hpp:112
The Multiple-Precision Data Type l_real.
Definition: l_real.hpp:78
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
l_interval Sqrt3_l_interval() noexcept
Enclosure-Interval for .
Definition: l_imath.cpp:2307
l_real Sqrt3r_l_real() noexcept
Approximation of .
Definition: l_rmath.cpp:252
l_interval Pi2_l_interval() noexcept
Enclosure-Interval for .
Definition: l_imath.cpp:1569
l_real Ln2_l_real() noexcept
Approximation of .
Definition: l_rmath.cpp:208
l_real Pid3_l_real() noexcept
Approximation of .
Definition: l_rmath.cpp:230
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
l_real E_l_real() noexcept
Approximation of .
Definition: l_rmath.cpp:258
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
l_real Er_l_real() noexcept
Approximation of .
Definition: l_rmath.cpp:260
l_real SqrtPir_l_real() noexcept
Approximation of .
Definition: l_rmath.cpp:240
l_real LnPi_l_real() noexcept
Approximation of .
Definition: l_rmath.cpp:254
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
l_real Pip2_l_real() noexcept
Approximation of .
Definition: l_rmath.cpp:244
l_real Ln10r_l_real() noexcept
Approximation of .
Definition: l_rmath.cpp:212
l_real Ln10_l_real() noexcept
Approximation of .
Definition: l_rmath.cpp:210
l_interval SqrtPi_l_interval() noexcept
Enclosure-Interval for .
Definition: l_imath.cpp:1863
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
l_real Catalan_l_real() noexcept
Approximation of Catalan Numbers.
Definition: l_rmath.cpp:276
l_real EpPid2_l_real() noexcept
Approximation of .
Definition: l_rmath.cpp:270
l_real SqrtPi_l_real() noexcept
Approximation of .
Definition: l_rmath.cpp:236
l_real EulerGa_l_real() noexcept
Approximation of Euler Gamma.
Definition: l_rmath.cpp:274
l_real Pi2r_l_real() noexcept
Approximation of .
Definition: l_rmath.cpp:234
l_interval EulerGa_l_interval() noexcept
Enclosure-Interval for Euler Gamma.
Definition: l_imath.cpp:3269
l_interval Ln2_l_interval() noexcept
Enclosure-Interval for .
Definition: l_imath.cpp:854
l_real Ln2Pi_l_real() noexcept
Approximation of .
Definition: l_rmath.cpp:256
l_real Sqrt2Pir_l_real() noexcept
Approximation of .
Definition: l_rmath.cpp:242
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
l_real EpPi_l_real() noexcept
Approximation of .
Definition: l_rmath.cpp:266
l_real Sqrt2Pi_l_real() noexcept
Approximation of .
Definition: l_rmath.cpp:238
l_real Sqrt3_l_real() noexcept
Approximation of .
Definition: l_rmath.cpp:248
l_real Ep2_l_real() noexcept
Approximation of .
Definition: l_rmath.cpp:262
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
l_real Pir_l_real() noexcept
Approximation of .
Definition: l_rmath.cpp:232
l_interval Ln10_l_interval() noexcept
Enclosure-Interval for .
Definition: l_imath.cpp:925
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
l_real Sqrt5_l_real() noexcept
Approximation of .
Definition: l_rmath.cpp:218
l_real Ep2Pi_l_real() noexcept
Approximation of .
Definition: l_rmath.cpp:268
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
l_real Sqrt2_l_real() noexcept
Approximation of .
Definition: l_rmath.cpp:216
l_real Ep2r_l_real() noexcept
Approximation of .
Definition: l_rmath.cpp:264
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
void times2pown(cinterval &x, int n) noexcept
Fast multiplication of reference parameter [z] with .
Definition: cimath.cpp:2059
l_real Sqrt2r_l_real() noexcept
Approximation of .
Definition: l_rmath.cpp:246
l_real Pi2_l_real() noexcept
Approximation of .
Definition: l_rmath.cpp:228
l_real Sqrt7_l_real() noexcept
Approximation of .
Definition: l_rmath.cpp:220
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_real Sqrt3d2_l_real() noexcept
Approximation of .
Definition: l_rmath.cpp:250
l_interval Pir_l_interval() noexcept
Enclosure-Interval for .
Definition: l_imath.cpp:1716
l_real Pi_l_real() noexcept
Approximation of .
Definition: l_rmath.cpp:224
l_real Pid4_l_real() noexcept
Approximation of .
Definition: l_rmath.cpp:214
l_interval Sqrt2Pir_l_interval() noexcept
Enclosure-Interval for .
Definition: l_imath.cpp:2085
l_real Pid2_l_real() noexcept
Approximation of .
Definition: l_rmath.cpp:226
interval sqrtx2y2(const interval &x, const interval &y) noexcept
Calculates .
Definition: imath.cpp:80
l_real EpPid4_l_real() noexcept
Approximation of .
Definition: l_rmath.cpp:272
l_real Ln2r_l_real() noexcept
Approximation of .
Definition: l_rmath.cpp:222
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