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
30namespace cxsc {
31
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
133l_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
153l_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):
208l_real Ln2_l_real() noexcept // ln(2)
209{ return mid( Ln2_l_interval() ); }
210l_real Ln10_l_real() noexcept // ln(10)
211{ return mid( Ln10_l_interval() ); }
212l_real Ln10r_l_real() noexcept // 1/ln(10)
213{ return mid( Ln10r_l_interval()); }
214l_real Pid4_l_real() noexcept // Pi/4
215{ return mid( Pid4_l_interval() ); }
216l_real Sqrt2_l_real() noexcept // sqrt(2)
217{ return mid( Sqrt2_l_interval() ); }
218l_real Sqrt5_l_real() noexcept // sqrt(5)
219{ return mid( Sqrt5_l_interval() ); }
220l_real Sqrt7_l_real() noexcept // sqrt(7)
221{ return mid( Sqrt7_l_interval() ); }
222l_real Ln2r_l_real() noexcept // 1/ln(2)
223{ return mid( Ln2r_l_interval() ); }
224l_real Pi_l_real() noexcept // Pi
225{ return mid( Pi_l_interval() ); }
226l_real Pid2_l_real() noexcept // Pi/2
227{ return mid( Pid2_l_interval() ); }
228l_real Pi2_l_real() noexcept // 2*Pi
229{ return mid( Pi2_l_interval() ); }
230l_real Pid3_l_real() noexcept // Pi/3
231{ return mid( Pid3_l_interval() ); }
232l_real Pir_l_real() noexcept // 1/Pi
233{ return mid( Pir_l_interval() ); }
234l_real Pi2r_l_real() noexcept // 1/(2*Pi)
235{ return mid( Pi2r_l_interval() ); }
236l_real SqrtPi_l_real() noexcept // sqrt(Pi)
237{ return mid( SqrtPi_l_interval() ); }
238l_real Sqrt2Pi_l_real() noexcept // sqrt(2*Pi)
239{ return mid( Sqrt2Pi_l_interval() ); }
240l_real SqrtPir_l_real() noexcept // 1/sqrt(Pi)
241{ return mid( SqrtPir_l_interval() ); }
242l_real Sqrt2Pir_l_real() noexcept // 1/sqrt(2*Pi)
243{ return mid( Sqrt2Pir_l_interval() ); }
244l_real Pip2_l_real() noexcept // Pi^2
245{ return mid( Pip2_l_interval() ); }
246l_real Sqrt2r_l_real() noexcept // 1/sqrt(2)
247{ return mid( Sqrt2r_l_interval() ); }
248l_real Sqrt3_l_real() noexcept // sqrt(3)
249{ return mid( Sqrt3_l_interval() ); }
250l_real Sqrt3d2_l_real() noexcept // sqrt(3)/2
251{ return mid( Sqrt3d2_l_interval() ); }
252l_real Sqrt3r_l_real() noexcept // 1/sqrt(3)
253{ return mid( Sqrt3r_l_interval() ); }
254l_real LnPi_l_real() noexcept // ln(Pi)
255{ return mid( LnPi_l_interval() ); }
256l_real Ln2Pi_l_real() noexcept // ln(2*Pi)
257{ return mid( Ln2Pi_l_interval() ); }
258l_real E_l_real() noexcept // e = exp(1)
259{ return mid( E_l_interval() ); }
260l_real Er_l_real() noexcept // 1/e
261{ return mid( Er_l_interval() ); }
262l_real Ep2_l_real() noexcept // e^2
263{ return mid( Ep2_l_interval() ); }
264l_real Ep2r_l_real() noexcept // 1/e^2
265{ return mid( Ep2r_l_interval() ); }
266l_real EpPi_l_real() noexcept // e^Pi
267{ return mid( EpPi_l_interval() ); }
268l_real Ep2Pi_l_real() noexcept // e^(2*Pi)
269{ return mid( Ep2Pi_l_interval() ); }
270l_real EpPid2_l_real() noexcept // e^(Pi/2)
271{ return mid( EpPid2_l_interval() ); }
272l_real EpPid4_l_real() noexcept // e^(Pi/4)
273{ return mid( EpPid4_l_interval() ); }
274l_real EulerGa_l_real() noexcept // EulerGamma
275{ return mid(EulerGa_l_interval() ); }
276l_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 Sqrt3r_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:2455
l_interval SqrtPi_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:1863
l_real Sqrt3r_l_real() noexcept
Approximation of .
Definition l_rmath.cpp:252
l_interval Ep2r_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:2899
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_real E_l_real() noexcept
Approximation of .
Definition l_rmath.cpp:258
l_interval Pi2_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:1569
l_interval Sqrt2Pir_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:2085
l_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 Sqrt2Pi_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:1937
l_interval Pi_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:1423
l_interval Sqrt2_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:1135
cinterval power(const cinterval &z, int n) noexcept
Calculates .
Definition cimath.cpp:1941
l_interval Ln10r_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:995
l_interval Ep2_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:2825
l_interval EpPi_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:2973
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 LnPi_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:2529
l_real Catalan_l_real() noexcept
Approximation of Catalan Numbers.
Definition l_rmath.cpp:276
l_interval Ep2Pi_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:3047
l_real EpPid2_l_real() noexcept
Approximation of .
Definition l_rmath.cpp:270
l_interval Pir_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:1716
l_real SqrtPi_l_real() noexcept
Approximation of .
Definition l_rmath.cpp:236
l_interval Pip2_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:2159
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 Catalan_l_interval() noexcept
Enclosure-Interval for Catalan Numbers.
Definition l_imath.cpp:3343
l_interval Pid2_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:1496
l_interval Sqrt7_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:1276
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_real EpPi_l_real() noexcept
Approximation of .
Definition l_rmath.cpp:266
l_interval E_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:2677
l_interval Sqrt2r_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:2233
l_real Sqrt2Pi_l_real() noexcept
Approximation of .
Definition l_rmath.cpp:238
l_interval EulerGa_l_interval() noexcept
Enclosure-Interval for Euler Gamma.
Definition l_imath.cpp:3269
l_real Sqrt3_l_real() noexcept
Approximation of .
Definition l_rmath.cpp:248
l_interval Pid3_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:1642
l_real Ep2_l_real() noexcept
Approximation of .
Definition l_rmath.cpp:262
l_interval Sqrt3_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:2307
l_interval SqrtPir_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:2011
l_real Pir_l_real() noexcept
Approximation of .
Definition l_rmath.cpp:232
l_interval Ln2r_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:1349
l_interval EpPid2_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:3121
l_interval Pid4_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:1065
cinterval sqrt1px2(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:1071
l_interval Sqrt3d2_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:2381
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 EpPid4_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:3195
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 Er_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:2751
ivector abs(const cimatrix_subv &mv) noexcept
Returns the absolute value of the matrix.
Definition cimatrix.inl:737
cinterval sqrt(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:1007
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_interval Ln2Pi_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:2603
l_interval Sqrt5_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:1204
l_interval Ln2_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:854
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
cvector mid(const cimatrix_subv &mv) noexcept
Returns the middle of the matrix.
Definition cimatrix.inl:739
l_real Sqrt3d2_l_real() noexcept
Approximation of .
Definition l_rmath.cpp:250
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 Ln10_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:925
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 Pi2r_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:1789