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 
28 namespace 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 
243 l_complex power(const l_complex& z, int n) noexcept
244 { return mid( power(l_cinterval(z),n) ); }
245 
246 l_complex log2(const l_complex& z) noexcept
247 { return mid(log2(l_cinterval(z))); }
248 
249 l_complex log10(const l_complex& z) noexcept
250 { return mid(log10(l_cinterval(z))); }
251 
252 l_complex pow(const l_complex& z, const l_real& p) noexcept
253 { return mid( pow( l_cinterval(z) , l_interval(p) ) ); }
254 
255 l_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
cxsc::power
cinterval power(const cinterval &z, int n) noexcept
Calculates .
Definition: cimath.cpp:1941
cxsc::mid
cvector mid(const cimatrix_subv &mv) noexcept
Returns the middle of the matrix.
Definition: cimatrix.inl:739
cxsc::lnp1
cinterval lnp1(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:867
cxsc::asinh
cinterval asinh(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:2718
cxsc::sqrt1px2
cinterval sqrt1px2(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:1071
cxsc::coth
cinterval coth(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:578
cxsc::sin
cinterval sin(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:215
cxsc::sqrt_all
std::list< cinterval > sqrt_all(const cinterval &z)
Calculates and returns all possible solutions.
Definition: cimath.cpp:1176
cxsc::cot
cinterval cot(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:538
cxsc::exp10
cinterval exp10(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:172
cxsc::log10
cinterval log10(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:903
cxsc::arg
interval arg(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:741
cxsc::tan
cinterval tan(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:393
cxsc::tanh
cinterval tanh(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:565
cxsc::sqrt
cinterval sqrt(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:1007
cxsc::ln
cinterval ln(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:851
cxsc::log2
cinterval log2(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:898
cxsc::abs
ivector abs(const cimatrix_subv &mv) noexcept
Returns the absolute value of the matrix.
Definition: cimatrix.inl:737
cxsc::sqrtx2m1
cinterval sqrtx2m1(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:1109
cxsc::Arg
interval Arg(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:654
cxsc::exp
cinterval exp(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:159
cxsc::acos
cinterval acos(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:2553
cxsc::l_complex
The Multiple-Precision Data Type l_complex.
Definition: l_complex.hpp:46
cxsc::l_interval
The Multiple-Precision Data Type l_interval.
Definition: l_interval.hpp:72
cxsc::cosh
cinterval cosh(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:223
cxsc::Pi_l_interval
l_interval Pi_l_interval() noexcept
Enclosure-Interval for .
Definition: l_imath.cpp:1423
cxsc::sinh
cinterval sinh(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:231
cxsc::pow
cinterval pow(const cinterval &z, const interval &p) noexcept
Calculates .
Definition: cimath.cpp:2074
cxsc::exp2
cinterval exp2(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:167
cxsc::acoth
cinterval acoth(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:3330
cxsc
The namespace cxsc, providing all functionality of the class library C-XSC.
Definition: cdot.cpp:29
cxsc::l_real
The Multiple-Precision Data Type l_real.
Definition: l_real.hpp:78
cxsc::sqrtp1m1
cinterval sqrtp1m1(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:1054
cxsc::expm1
cinterval expm1(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:177
cxsc::atanh
cinterval atanh(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:3317
cxsc::l_cinterval
The Multiple-Precision Data Type l_cinterval.
Definition: l_cinterval.hpp:54
cxsc::times2pown
void times2pown(cinterval &x, int n) noexcept
Fast multiplication of reference parameter [z] with .
Definition: cimath.cpp:2059
cxsc::cos
cinterval cos(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:207
cxsc::atan
cinterval atan(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:2938
cxsc::acot
cinterval acot(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:3130
cxsc::power_fast
cinterval power_fast(const cinterval &z, int n) noexcept
Calculates .
Definition: cimath.cpp:1520
cxsc::sqrt1mx2
cinterval sqrt1mx2(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:1140
cxsc::acosh
cinterval acosh(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:2732
cxsc::sqr
cinterval sqr(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:3342
cxsc::asin
cinterval asin(const cinterval &z) noexcept
Calculates .
Definition: cimath.cpp:2311