OPAL (Object Oriented Parallel Accelerator Library)  2024.1
OPAL
FTpsMath.h
Go to the documentation of this file.
1 #ifndef CLASSIC_FTpsMath_HH
2 #define CLASSIC_FTpsMath_HH
3 
4 // ------------------------------------------------------------------------
5 // $RCSfile: FTpsMath.h,v $
6 // ------------------------------------------------------------------------
7 // $Revision: 1.1.1.1.2.4 $
8 // ------------------------------------------------------------------------
9 // Copyright: see Copyright.readme
10 // ------------------------------------------------------------------------
11 // Description:
12 //
13 // Declared template functions:
14 //
15 // ------------------------------------------------------------------------
16 // Class category: FixedAlgebra
17 // ------------------------------------------------------------------------
18 //
19 // $Date: 2003/11/07 18:04:21 $
20 // $Author: dbruhwil $
21 //
22 // ------------------------------------------------------------------------
23 
24 #include "FixedAlgebra/FTps.h"
25 #include "Utilities/DomainError.h"
26 
27 #include <complex> // std::real
28 
29 #include <type_traits>
30 
31 // 25. March 2017,
32 // http://stackoverflow.com/questions/30736951/templated-class-check-if-complex-at-compile-time
33 template<class T> struct is_complex : std::false_type {};
34 template<class T> struct is_complex<std::complex<T>> : std::true_type {};
35 
36 
37 // Class FTps; global functions acting on FTps objects.
38 // Elementary functions acting on FTps<T,N> objects.
39 // ------------------------------------------------------------------------
40 
42 template <class T, int N>
43 FTps<T, N> pow(const FTps<T, N> &x, int y, int trunc = (FTps<T, N>::EXACT));
44 
46 template <class T, int N>
47 FTps<T, N> sqrt(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
48 
50 template <class T, int N>
51 FTps<T, N> sin(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
52 
54 template <class T, int N>
55 FTps<T, N> cos(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
56 
58 template <class T, int N>
59 FTps<T, N> tan(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
60 
62 template <class T, int N>
63 FTps<T, N> cot(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
64 
66 template <class T, int N>
67 FTps<T, N> sec(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
68 
70 template <class T, int N>
71 FTps<T, N> csc(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
72 
74 template <class T, int N>
75 FTps<T, N> exp(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
76 
78 template <class T, int N>
79 FTps<T, N> log(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
80 
82 template <class T, int N>
83 FTps<T, N> sinh(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
84 
86 template <class T, int N>
87 FTps<T, N> cosh(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
88 
90 template <class T, int N>
91 FTps<T, N> tanh(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
92 
94 template <class T, int N>
95 FTps<T, N> coth(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
96 
98 template <class T, int N>
99 FTps<T, N> sech(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
100 
102 template <class T, int N>
103 FTps<T, N> csch(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
104 
106 template <class T, int N>
107 FTps<T, N> erf(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
108 
110 template <class T, int N>
111 FTps<T, N> erfc(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
112 
113 
114 // Implementation
115 // ------------------------------------------------------------------------
116 
117 template <class T, int N>
118 FTps<T, N> pow(const FTps<T, N> &x, int y, int trunc) {
119  // Default: trunc = EXACT
120 
121  FTps<T, N> z(T(1));
122 
123  if(y > 0) {
124  while(y-- > 0) z = z.multiply(x, trunc);
125  } else if(y < 0) {
126  if(x[0] == T(0) || x.getMinOrder() != 0)
127  throw DomainError("pow(const FTps &,int,int)");
128  FTps<T, N> t = x.inverse(trunc);
129  while(y++ < 0) z = z.multiply(t, trunc);
130  }
131 
132  return z;
133 }
134 
135 
136 template <class T, int N>
137 FTps<T, N> sqrt(const FTps<T, N> &x, int trunc) {
138  // Default: trunc = EXACT
139  int trcOrder = std::min(x.getTruncOrder(), trunc);
140 
141  if(trcOrder == FTps<T, N>::EXACT)
142  throw LogicalError("::sqrt(FTps<T,N> &x, int trunc)",
143  "Square-root of EXACT polynomial must be truncated.");
144 
145  T aZero = x[0];
146  if( ( std::real(aZero) <= std::real(T(0)) &&
147  std::imag(aZero) <= std::real(T(0))
148  ) || x.getMinOrder() != 0)
149  {
150  std::cerr << "FTps::sqrt(x) called with\nconstant term = " << aZero
151  << "; minOrder = " << x.getMinOrder() << std::endl
152  << "x = " << x << std::endl;
153  throw DomainError("sqrt(const FTps &,int)");
154  }
155 
156  T two_aZero = T(2) * aZero;
157  Array1D<T> series(trcOrder + 1);
158  series[0] = sqrt(aZero);
159  for(int i = 1; i <= trcOrder; i++) {
160  series[i] = series[i-1] * double(3 - 2 * i) / (two_aZero * double(i));
161  }
162 
163  return x.taylor(series, trcOrder);
164 }
165 
166 
167 template <class T, int N>
168 FTps<T, N> sin(const FTps<T, N> &x, int trunc) {
169  // Default: trunc = EXACT
170  int trcOrder = std::min(x.getTruncOrder(), trunc);
171 
172  if(trcOrder == FTps<T, N>::EXACT)
173  throw LogicalError("::sin(FTps<T,N> &x, int trunc)",
174  "Sine of EXACT polynomial must be truncated.");
175 
176  T aZero = x[0];
177  if(x.getMinOrder() != 0) aZero = T(0);
178 
179  Array1D<T> series(trcOrder + 1);
180  series[0] = sin(aZero);
181  series[1] = cos(aZero);
182  for(int i = 2; i <= trcOrder; i++) {
183  series[i] = - series[i-2] / double(i * (i - 1));
184  }
185 
186  return x.taylor(series, trcOrder);
187 }
188 
189 
190 template <class T, int N>
191 FTps<T, N> cos(const FTps<T, N> &x, int trunc) {
192  // Default: trunc = EXACT
193  int trcOrder = std::min(x.getTruncOrder(), trunc);
194 
195  if(trcOrder == FTps<T, N>::EXACT)
196  throw LogicalError("::cos(FTps<T,N> &x, int trunc)",
197  "Cosine of EXACT polynomial must be truncated.");
198 
199  T aZero = x[0];
200  if(x.getMinOrder() != 0) aZero = T(0);
201 
202  Array1D<T> series(trcOrder + 1);
203  series[0] = cos(aZero);
204  series[1] = - sin(aZero);
205  for(int i = 2; i <= trcOrder; i++) {
206  series[i] = - series[i-2] / double(i * (i - 1));
207  }
208 
209  return x.taylor(series, trcOrder);
210 }
211 
212 
213 template <class T, int N>
214 FTps<T, N> tan(const FTps<T, N> &x, int trunc) {
215  // Default: trunc = EXACT
216 
217  // The direct series expansion requires the Bernoulli numbers to arbitrary order.
218  return sin(x, trunc) / cos(x, trunc);
219 }
220 
221 template <class T, int N>
222 FTps<T, N> cot(const FTps<T, N> &x, int trunc) {
223  // Default: trunc = EXACT
224 
225  if(x[0] == T(0) || x.getMinOrder() != 0)
226  throw DomainError("cot(const FTps &,int)");
227 
228  return cos(x, trunc) / sin(x, trunc);
229 }
230 
231 
232 template <class T, int N>
233 FTps<T, N> sec(const FTps<T, N> &x, int trunc) {
234  // Default: trunc = EXACT
235 
236  return cos(x, trunc).inverse();
237 }
238 
239 
240 template <class T, int N>
241 FTps<T, N> csc(const FTps<T, N> &x, int trunc) {
242  // Default: trunc = EXACT
243 
244  if(x[0] == T(0) || x.getMinOrder() != 0)
245  throw DomainError("csc(const FTps &,int)");
246 
247  return sin(x, trunc).inverse();
248 }
249 
250 
251 template <class T, int N>
252 FTps<T, N> exp(const FTps<T, N> &x, int trunc) {
253  // Default: trunc = EXACT
254  int trcOrder = std::min(x.getTruncOrder(), trunc);
255 
256  if(trcOrder == FTps<T, N>::EXACT)
257  throw LogicalError("::exp(FTps<T,N> &x, int trunc)",
258  "Exponential of EXACT polynomial must be truncated.");
259 
260  T aZero = x[0];
261  if(x.getMinOrder() != 0) aZero = T(0);
262 
263  Array1D<T> series(trcOrder + 1);
264  series[0] = exp(aZero);
265  for(int i = 1; i <= trcOrder; i++) {
266  series[i] = series[i-1] / double(i);
267  }
268 
269  return x.taylor(series, trcOrder);
270 }
271 
272 
273 template <class T, int N>
274 FTps<T, N> log(const FTps<T, N> &x, int trunc) {
275  // Default: trunc = EXACT
276  int trcOrder = std::min(x.getTruncOrder(), trunc);
277 
278  if(trcOrder == FTps<T, N>::EXACT)
279  throw LogicalError("::log(FTps<T,N> &x, int trunc)",
280  "Logarithm of EXACT polynomial must be truncated.");
281 
282  T aZero = x[0];
283  if(aZero <= T(0) || x.getMinOrder() != 0)
284  throw DomainError("log(const FTps &,int)");
285 
286  T a0inv = T(1) / aZero;
287  T ain = a0inv;
288  Array1D<T> series(trcOrder + 1);
289  series[0] = log(aZero);
290  series[1] = a0inv;
291  for(int i = 2; i <= trcOrder; i++) {
292  ain *= -a0inv;
293  series[i] = ain / double(i);
294  }
295 
296  return x.taylor(series, trcOrder);
297 }
298 
299 
300 template <class T, int N>
301 FTps<T, N> sinh(const FTps<T, N> &x, int trunc) {
302  // Default: trunc = EXACT
303  int trcOrder = std::min(x.getTruncOrder(), trunc);
304 
305  if(trcOrder == FTps<T, N>::EXACT)
306  throw LogicalError("::log(FTps<T,N> &x, int trunc)",
307  "Hyperbolic sine of EXACT polynomial must be truncated.");
308 
309  T aZero = x[0];
310  if(x.getMinOrder() != 0) aZero = T(0);
311 
312  Array1D<T> series(trcOrder + 1);
313  series[0] = sinh(aZero);
314  series[1] = cosh(aZero);
315  for(int i = 2; i <= trcOrder; i++) {
316  series[i] = series[i-2] / double(i * (i - 1));
317  }
318 
319  return x.taylor(series, trcOrder);
320 }
321 
322 
323 template <class T, int N>
324 FTps<T, N> cosh(const FTps<T, N> &x, int trunc) {
325  // Default: trunc = EXACT
326  int trcOrder = std::min(x.getTruncOrder(), trunc);
327 
328  if(trcOrder == FTps<T, N>::EXACT)
329  throw LogicalError("::cosh(FTps<T,N> &x, int trunc)",
330  "Hyperbolic cosine of EXACT polynomial must be truncated.");
331 
332  T aZero = x[0];
333  if(x.getMinOrder() != 0) aZero = T(0);
334 
335  Array1D<T> series(trcOrder + 1);
336  series[0] = cosh(aZero);
337  series[1] = sinh(aZero);
338  for(int i = 2; i <= trcOrder; i++) {
339  series[i] = series[i-2] / double(i * (i - 1));
340  }
341 
342  return x.taylor(series, trcOrder);
343 }
344 
345 
346 template <class T, int N>
347 FTps<T, N> tanh(const FTps<T, N> &x, int trunc) {
348  // Default: trunc = EXACT
349 
350  return sinh(x, trunc) / cosh(x, trunc);
351 }
352 
353 
354 template <class T, int N>
355 FTps<T, N> coth(const FTps<T, N> &x, int trunc) {
356  // Default: trunc = EXACT
357 
358  if(x[0] == T(0) || x.getMinOrder() != 0)
359  throw DomainError("coth(const FTps &,int)");
360 
361  return cosh(x, trunc) / sinh(x, trunc);
362 }
363 
364 
365 template <class T, int N>
366 FTps<T, N> sech(const FTps<T, N> &x, int trunc) {
367  // Default: trunc = EXACT
368 
369  return cosh(x, trunc).inverse();
370 }
371 
372 
373 template <class T, int N>
374 FTps<T, N> csch(const FTps<T, N> &x, int trunc) {
375  // Default: trunc = EXACT
376 
377  if(x[0] == T(0) || x.getMinOrder() != 0)
378  throw DomainError("csch(const FTps &,int)");
379 
380  return sinh(x, trunc).inverse();
381 }
382 
384 template <class T, int N>
385 FTps<T, N> erf(const FTps<T, N> &x, int trunc) {
386 
387  if ( is_complex<T>::value )
388  throw LogicalError("::erf(FTps<T,N> &x, int trunc)",
389  "Error function does not support complex numbers.");
390 
391  // Default: trunc = EXACT
392  int trcOrder = std::min(x.getTruncOrder(), trunc);
393 
394  if(trcOrder == FTps<T, N>::EXACT)
395  throw LogicalError("::erf(FTps<T,N> &x, int trunc)",
396  "Error function of EXACT polynomial must be truncated.");
397 
398  T aZero = x[0];
399  if(x.getMinOrder() != 0) aZero = T(0);
400 
401  Array1D<T> series(trcOrder + 1);
402  series[0] = std::erf(std::real(aZero));
403  series[1] = 2.0 / std::sqrt(M_PI) * std::exp(-aZero*aZero);
404 
405  for(int i = 2; i <= trcOrder; ++i) {
406  series[i] = - 2.0 / double(i-1) * double((i-2)) * series[i-2] / double(i);
407  }
408 
409  return x.taylor(series, trcOrder);
410 }
411 
413 template <class T, int N>
414 FTps<T, N> erfc(const FTps<T, N> &x, int trunc) {
415  // Default: trunc = EXACT
416 
417  return T(1) - erf(x, trunc);
418 }
419 
420 #endif // CLASSIC_FTpsMath_HH
sinh
Tps< T > sinh(const Tps< T > &x)
Hyperbolic sine.
Definition: TpsMath.h:204
sqrt
Tps< T > sqrt(const Tps< T > &x)
Square root.
Definition: TpsMath.h:91
tan
Tps< T > tan(const Tps< T > &x)
Tangent.
Definition: TpsMath.h:147
FTps
Truncated power series in N variables of type T.
Definition: FLieGenerator.h:29
erf
FTps< T, N > erf(const FTps< T, N > &x, int trunc=(FTps< T, N >::EXACT))
Error function.
Definition: FTpsMath.h:385
imag
FLieGenerator< T, N > imag(const FLieGenerator< std::complex< T >, N > &)
Take imaginary part of a complex generator.
Definition: FLieGenerator.hpp:437
FTps::multiply
FTps multiply(const FTps &y, int trunc=EXACT) const
Multiplication.
Definition: FTps.hpp:650
exp
Tps< T > exp(const Tps< T > &x)
Exponential.
Definition: TpsMath.h:165
FTps::getTruncOrder
int getTruncOrder() const
Get truncation order.
Definition: FTps.h:183
endl
Inform & endl(Inform &inf)
Definition: Inform.cpp:42
FTps::inverse
FTps inverse(int trunc=EXACT) const
Reciprocal, 1/(*this).
Definition: FTps.hpp:707
min
T::PETE_Expr_t::PETE_Return_t min(const PETE_Expr< T > &expr, NDIndex< D > &loc)
Definition: ReductionLoc.h:76
DomainError
Domain error exception.
Definition: DomainError.h:32
tanh
Tps< T > tanh(const Tps< T > &x)
Hyperbolic tangent.
Definition: TpsMath.h:240
FTps::getMinOrder
int getMinOrder() const
Get minimum order.
Definition: FTps.h:165
sech
FTps< T, N > sech(const FTps< T, N > &x, int trunc=(FTps< T, N >::EXACT))
Hyperbolic secant.
Definition: FTpsMath.h:366
coth
FTps< T, N > coth(const FTps< T, N > &x, int trunc=(FTps< T, N >::EXACT))
Hyperbolic cotangent.
Definition: FTpsMath.h:355
cos
Tps< T > cos(const Tps< T > &x)
Cosine.
Definition: TpsMath.h:129
erfc
FTps< T, N > erfc(const FTps< T, N > &x, int trunc=(FTps< T, N >::EXACT))
Complementary error function.
Definition: FTpsMath.h:414
Array1D
One-dimensional array.
Definition: Array1D.h:36
log
Tps< T > log(const Tps< T > &x)
Natural logarithm.
Definition: TpsMath.h:182
FTps::taylor
FTps taylor(const Array1D< T > &series, int order) const
Taylor series.
Definition: FTps.hpp:1481
cot
Tps< T > cot(const Tps< T > &x)
Cotangent.
sec
Tps< T > sec(const Tps< T > &x)
Secant.
Definition: TpsMath.h:153
csc
Tps< T > csc(const Tps< T > &x)
Cosecant.
Definition: TpsMath.h:159
pow
Tps< T > pow(const Tps< T > &x, int y)
Integer power.
Definition: TpsMath.h:76
sin
Tps< T > sin(const Tps< T > &x)
Sine.
Definition: TpsMath.h:111
cosh
Tps< T > cosh(const Tps< T > &x)
Hyperbolic cosine.
Definition: TpsMath.h:222
LogicalError
Logical error exception.
Definition: LogicalError.h:33
interpolation::complex
MMatrix< m_complex > complex(MMatrix< double > real)
Definition: MMatrix.cpp:396
real
FLieGenerator< T, N > real(const FLieGenerator< std::complex< T >, N > &)
Take real part of a complex generator.
Definition: FLieGenerator.hpp:422
csch
FTps< T, N > csch(const FTps< T, N > &x, int trunc=(FTps< T, N >::EXACT))
Hyperbolic cosecant.
Definition: FTpsMath.h:374
Generated on Tue Jul 2 2024 15:05:22 for OPAL (Object Oriented Parallel Accelerator Library) by   doxygen 1.8.5