C-XSC - A C++ Class Library for Extended Scientific Computing 2.5.4
l_cmath.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_cmath.cpp,v 1.11 2014/01/30 17:23:46 cxsc Exp $ */
25
26#include "l_cmath.hpp"
27
28namespace cxsc {
29
30 l_complex sqrt(const l_complex& z) noexcept
31// Computation of sqrt(z); stagprec <= stagmax=19 determines the maximum
32// accuracy of about 16*19 = 304 decimal digits.
33// The branch cut of this sqrt-function is the negative imaginary axis.
34// stagmax=19 is predetermined by the internal used function
35// l_real sqrt(const l_real&), which has its own maximum precision of
36// stagprec=19;
37 {
38 l_real x=Re(z), y=Im(z), w;
39 if ( zero_(x) && zero_(y) ) return l_complex(0,0);
40 // Now z != 0
41 int stagsave = stagprec, stagmax = 19;
42 if (stagprec > stagmax) stagprec = stagmax;
43 // stagprec>19 makes no sense, because stagmax(sqrt(...))=19;
44 int ex = expo(x[1]), exy = expo(y[1]);
45 if (exy > ex) ex = exy; // ex: maximum of the exponents;
46 ex = 400-ex; // ex: optimal scaling factor
47 if (ex%2) ex--; // ex is now even;
48 times2pown(x,ex); times2pown(y,ex); // scaling with 2^ex
49 bool neg = sign(x[1]) < 0;
50 if (neg) x = -x;
51 w = abs( l_complex(x,y) ) + x;
52 times2pown(w,1);
53 w = sqrt(w);
54 if (neg)
55 {
56 x = abs(y)/w;
57 y = sign(y[1]) < 0 ? -w : w;
58 times2pown(y,-1);
59 } else
60 {
61 x = w;
62 times2pown(x,-1);
63 y /= w;
64 }
65 ex /= 2; // Backscaling of the current result with 2^(-ex):
66 times2pown(x,-ex); times2pown(y,-ex);
67 stagprec = stagsave; // restore old value of the stagprec variable
68 return l_complex(x,y);
69 } // sqrt
70
71 l_complex sqrtp1m1(const l_complex& z) noexcept
72 { return mid(sqrtp1m1(l_cinterval(z))); }
73
74 l_complex sqrt1px2(const l_complex& z) noexcept
75 { return mid(sqrt1px2(l_cinterval(z))); }
76
77 l_complex sqrtx2m1(const l_complex& z) noexcept
78 { return mid(sqrtx2m1(l_cinterval(z))); }
79
80 l_complex sqrt1mx2(const l_complex& z) noexcept
81 { return mid(sqrt1mx2(l_cinterval(z))); }
82
83 l_complex exp(const l_complex& z) noexcept
84 { return mid(exp(l_cinterval(z))); }
85
86 l_complex expm1(const l_complex& z) noexcept
87 { return mid(expm1(l_cinterval(z))); }
88
89 l_complex exp2(const l_complex& z) noexcept
90 { return mid(exp2(l_cinterval(z))); }
91
92 l_complex exp10(const l_complex& z) noexcept
93 { return mid(exp10(l_cinterval(z))); }
94
95 l_complex sin(const l_complex& z) noexcept
96 { return mid(sin(l_cinterval(z))); }
97
98 l_complex cos(const l_complex& z) noexcept
99 { return mid(cos(l_cinterval(z))); }
100
101 l_complex tan(const l_complex& z) noexcept
102 { return mid(tan(l_cinterval(z))); }
103
104 l_complex cot(const l_complex& z) noexcept
105 { return mid(cot(l_cinterval(z))); }
106
107 l_complex asin(const l_complex& z) noexcept
108 { return mid(asin(l_cinterval(z))); }
109
110 l_complex acos(const l_complex& z) noexcept
111 { return mid(acos(l_cinterval(z))); }
112
113 l_complex atan(const l_complex& z) noexcept
114 { return mid(atan(l_cinterval(z))); }
115
116 l_complex acot(const l_complex& z) noexcept
117 { return mid(acot(l_cinterval(z))); }
118
119 l_complex sinh(const l_complex& z) noexcept
120 { return mid(sinh(l_cinterval(z))); }
121
122 l_complex cosh(const l_complex& z) noexcept
123 { return mid(cosh(l_cinterval(z))); }
124
125 l_complex tanh(const l_complex& z) noexcept
126 { return mid(tanh(l_cinterval(z))); }
127
128 l_complex coth(const l_complex& z) noexcept
129 { return mid(coth(l_cinterval(z))); }
130
131 l_complex asinh(const l_complex& z) noexcept
132 { return mid(asinh(l_cinterval(z))); }
133
134 l_complex acosh(const l_complex& z) noexcept
135 { return mid(acosh(l_cinterval(z))); }
136
137 l_complex atanh(const l_complex& z) noexcept
138 { return mid(atanh(l_cinterval(z))); }
139
140 l_complex acoth(const l_complex& z) noexcept
141 { return mid(acoth(l_cinterval(z))); }
142
143 // sqrt_all(c) computes a list of 2 values for all square roots of c
144 std::list<l_complex> sqrt_all( const l_complex& c )
145 {
146 l_complex lc;
147 lc = sqrt(c);
148
149 std::list<l_complex> res;
150 res.push_back( lc );
151 res.push_back( -lc );
152
153 return res;
154 } // end sqrt_all
155
156 l_complex sqrt(const l_complex& z, int n) noexcept
157 { return mid(sqrt(l_cinterval(z),n)); }
158
159 l_real arg(const l_complex& z) noexcept
160 { return mid(arg(l_cinterval(z))); }
161
162 l_real Arg(const l_complex& z) noexcept
163 { return mid(Arg(l_cinterval(z))); }
164
165 std::list<l_complex> sqrt_all( const l_complex& z, int n )
166//
167// sqrt_all(z,n) computes a list of n values approximating all n-th roots of z
168//
169// For n >=3, computing the optimal approximations is very expensive
170// and thus not considered cost-effective.
171//
172// Hence, the polar form is used to calculate sqrt_all(z,n)
173//
174 {
175 std::list<l_complex> res;
176
177 if( n == 0 )
178 {
179 res.push_back( l_complex(1,0) );
180 return res;
181 }
182 else if( n == 1 )
183 {
184 res.push_back(z);
185 return res;
186 }
187 else if( n == 2 ) return sqrt_all( z );
188 else
189 {
190 l_real
191 arg_z = arg( z ), root_abs_z = sqrt( abs( z ), n );
192
193 for(int k = 0; k < n; k++)
194 {
195 l_real arg_k = ( arg_z + 2 * k * mid(Pi_l_interval()) ) / n;
196
197 res.push_back( l_complex( root_abs_z * cos( arg_k ),
198 root_abs_z * sin( arg_k ) ) );
199 }
200 return res;
201 }
202 }
203 //
204//-- end sqrt_all -------------------------------------------------------------
205
206 l_complex ln(const l_complex& z) noexcept
207 { return mid(ln(l_cinterval(z))); }
208
209 l_complex lnp1(const l_complex& z) noexcept
210 { return mid(lnp1(l_cinterval(z))); }
211
212 l_complex power_fast(const l_complex& z, int n) noexcept
213 {
214 if( n == 0 )
215 return l_complex(1,0);
216 else
217 if( n == 1 ) return z;
218 else
219 if( n == -1 ) return 1 / z;
220 else
221 if( n == 2 ) return sqr(z);
222 else
223 {
224 l_real abs_z = abs(z);
225
226 if( n < 0 && abs_z == 0.0 )
227 // z contains 0
228 cxscthrow (STD_FKT_OUT_OF_DEF
229 ("l_complex power_fast(const l_complex& z, int n ); z == 0."));
230 if( abs_z == 0.0 )
231 return l_complex(0,0);
232 else
233 {
234 l_real arg_z = arg(z);
235 l_real abs_z_n = exp( n * ln( abs_z ) );
236
237 return l_complex( abs_z_n * cos( n * arg_z ),
238 abs_z_n * sin( n * arg_z ) );
239 }
240 }
241}
242
243l_complex power(const l_complex& z, int n) noexcept
244{ return mid( power(l_cinterval(z),n) ); }
245
246l_complex log2(const l_complex& z) noexcept
247{ return mid(log2(l_cinterval(z))); }
248
249l_complex log10(const l_complex& z) noexcept
250{ return mid(log10(l_cinterval(z))); }
251
252l_complex pow(const l_complex& z, const l_real& p) noexcept
253{ return mid( pow( l_cinterval(z) , l_interval(p) ) ); }
254
255l_complex pow(const l_complex& z, const l_complex& p) noexcept
256{ return mid( pow( l_cinterval(z) , l_cinterval(p) ) ); }
257
258
259
260} // namespace cxsc
The Multiple-Precision Data Type l_cinterval.
The Multiple-Precision Data Type l_complex.
Definition l_complex.hpp:46
The Multiple-Precision Data Type l_interval.
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
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
cinterval asinh(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:2718
cinterval coth(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:578
cinterval log2(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:898
l_interval Pi_l_interval() noexcept
Enclosure-Interval for .
Definition l_imath.cpp:1423
cinterval power(const cinterval &z, int n) noexcept
Calculates .
Definition cimath.cpp:1941
cinterval log10(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:903
cinterval ln(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:851
cinterval pow(const cinterval &z, const interval &p) noexcept
Calculates .
Definition cimath.cpp:2074
cinterval sinh(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:231
cinterval asin(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:2311
cinterval tan(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:393
cinterval exp10(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:172
interval arg(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:741
std::list< cinterval > sqrt_all(const cinterval &z)
Calculates and returns all possible solutions.
Definition cimath.cpp:1176
cinterval acos(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:2553
cinterval sqrtx2m1(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:1109
cinterval acosh(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:2732
cinterval cosh(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:223
cinterval cos(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:207
cinterval sqrt1px2(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:1071
cinterval exp(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:159
cinterval tanh(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:565
cinterval expm1(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:177
cinterval cot(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:538
ivector abs(const cimatrix_subv &mv) noexcept
Returns the absolute value of the matrix.
Definition cimatrix.inl:737
cinterval sqrt(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:1007
cinterval power_fast(const cinterval &z, int n) noexcept
Calculates .
Definition cimath.cpp:1520
cinterval acot(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:3130
void times2pown(cinterval &x, int n) noexcept
Fast multiplication of reference parameter [z] with .
Definition cimath.cpp:2059
cinterval sqr(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:3342
cinterval lnp1(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:867
cvector mid(const cimatrix_subv &mv) noexcept
Returns the middle of the matrix.
Definition cimatrix.inl:739
cinterval atan(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:2938
cinterval atanh(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:3317
interval Arg(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:654
cinterval acoth(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:3330
cinterval sin(const cinterval &z) noexcept
Calculates .
Definition cimath.cpp:215