OPAL (Object Oriented Parallel Accelerator Library) 2022.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
33template<class T> struct is_complex : std::false_type {};
34template<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
42template <class T, int N>
43FTps<T, N> pow(const FTps<T, N> &x, int y, int trunc = (FTps<T, N>::EXACT));
44
46template <class T, int N>
47FTps<T, N> sqrt(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
48
50template <class T, int N>
51FTps<T, N> sin(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
52
54template <class T, int N>
55FTps<T, N> cos(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
56
58template <class T, int N>
59FTps<T, N> tan(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
60
62template <class T, int N>
63FTps<T, N> cot(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
64
66template <class T, int N>
67FTps<T, N> sec(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
68
70template <class T, int N>
71FTps<T, N> csc(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
72
74template <class T, int N>
75FTps<T, N> exp(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
76
78template <class T, int N>
79FTps<T, N> log(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
80
82template <class T, int N>
83FTps<T, N> sinh(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
84
86template <class T, int N>
87FTps<T, N> cosh(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
88
90template <class T, int N>
91FTps<T, N> tanh(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
92
94template <class T, int N>
95FTps<T, N> coth(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
96
98template <class T, int N>
99FTps<T, N> sech(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
100
102template <class T, int N>
103FTps<T, N> csch(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
104
106template <class T, int N>
107FTps<T, N> erf(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
108
110template <class T, int N>
111FTps<T, N> erfc(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
112
113
114// Implementation
115// ------------------------------------------------------------------------
116
117template <class T, int N>
118FTps<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
136template <class T, int N>
137FTps<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
167template <class T, int N>
168FTps<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
190template <class T, int N>
191FTps<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
213template <class T, int N>
214FTps<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
221template <class T, int N>
222FTps<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
232template <class T, int N>
233FTps<T, N> sec(const FTps<T, N> &x, int trunc) {
234 // Default: trunc = EXACT
235
236 return cos(x, trunc).inverse();
237}
238
239
240template <class T, int N>
241FTps<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
251template <class T, int N>
252FTps<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
273template <class T, int N>
274FTps<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
300template <class T, int N>
301FTps<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
323template <class T, int N>
324FTps<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
346template <class T, int N>
347FTps<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
354template <class T, int N>
355FTps<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
365template <class T, int N>
366FTps<T, N> sech(const FTps<T, N> &x, int trunc) {
367 // Default: trunc = EXACT
368
369 return cosh(x, trunc).inverse();
370}
371
372
373template <class T, int N>
374FTps<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
384template <class T, int N>
385FTps<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
413template <class T, int N>
414FTps<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
real
FLieGenerator< T, N > real(const FLieGenerator< std::complex< T >, N > &)
Take real part of a complex generator.
Definition: FLieGenerator.hpp:422
imag
FLieGenerator< T, N > imag(const FLieGenerator< std::complex< T >, N > &)
Take imaginary part of a complex generator.
Definition: FLieGenerator.hpp:437
csc
FTps< T, N > csc(const FTps< T, N > &x, int trunc=(FTps< T, N >::EXACT))
Cosecant.
Definition: FTpsMath.h:241
exp
FTps< T, N > exp(const FTps< T, N > &x, int trunc=(FTps< T, N >::EXACT))
Exponential.
Definition: FTpsMath.h:252
pow
FTps< T, N > pow(const FTps< T, N > &x, int y, int trunc=(FTps< T, N >::EXACT))
Tps x to the power (int y).
Definition: FTpsMath.h:118
cot
FTps< T, N > cot(const FTps< T, N > &x, int trunc=(FTps< T, N >::EXACT))
Cotangent.
Definition: FTpsMath.h:222
sqrt
FTps< T, N > sqrt(const FTps< T, N > &x, int trunc=(FTps< T, N >::EXACT))
Square root.
Definition: FTpsMath.h:137
sech
FTps< T, N > sech(const FTps< T, N > &x, int trunc=(FTps< T, N >::EXACT))
Hyperbolic secant.
Definition: FTpsMath.h:366
erf
FTps< T, N > erf(const FTps< T, N > &x, int trunc=(FTps< T, N >::EXACT))
Error function.
Definition: FTpsMath.h:385
sin
FTps< T, N > sin(const FTps< T, N > &x, int trunc=(FTps< T, N >::EXACT))
Sine.
Definition: FTpsMath.h:168
erfc
FTps< T, N > erfc(const FTps< T, N > &x, int trunc=(FTps< T, N >::EXACT))
Complementary error function.
Definition: FTpsMath.h:414
tan
FTps< T, N > tan(const FTps< T, N > &x, int trunc=(FTps< T, N >::EXACT))
Tangent.
Definition: FTpsMath.h:214
log
FTps< T, N > log(const FTps< T, N > &x, int trunc=(FTps< T, N >::EXACT))
Natural logarithm.
Definition: FTpsMath.h:274
csch
FTps< T, N > csch(const FTps< T, N > &x, int trunc=(FTps< T, N >::EXACT))
Hyperbolic cosecant.
Definition: FTpsMath.h:374
tanh
FTps< T, N > tanh(const FTps< T, N > &x, int trunc=(FTps< T, N >::EXACT))
Hyperbolic tangent.
Definition: FTpsMath.h:347
sinh
FTps< T, N > sinh(const FTps< T, N > &x, int trunc=(FTps< T, N >::EXACT))
Hyperbolic sine.
Definition: FTpsMath.h:301
sec
FTps< T, N > sec(const FTps< T, N > &x, int trunc=(FTps< T, N >::EXACT))
Secant.
Definition: FTpsMath.h:233
cos
FTps< T, N > cos(const FTps< T, N > &x, int trunc=(FTps< T, N >::EXACT))
Cosine.
Definition: FTpsMath.h:191
cosh
FTps< T, N > cosh(const FTps< T, N > &x, int trunc=(FTps< T, N >::EXACT))
Hyperbolic cosine.
Definition: FTpsMath.h:324
coth
FTps< T, N > coth(const FTps< T, N > &x, int trunc=(FTps< T, N >::EXACT))
Hyperbolic cotangent.
Definition: FTpsMath.h:355
min
T::PETE_Expr_t::PETE_Return_t min(const PETE_Expr< T > &expr, NDIndex< D > &loc)
Definition: ReductionLoc.h:76
endl
Inform & endl(Inform &inf)
Definition: Inform.cpp:42
interpolation::complex
MMatrix< m_complex > complex(MMatrix< double > real)
Definition: MMatrix.cpp:396
Array1D
One-dimensional array.
Definition: Array1D.h:36
FTps
Truncated power series in N variables of type T.
Definition: FTps.h:45
FTps::inverse
FTps inverse(int trunc=EXACT) const
Reciprocal, 1/(*this).
Definition: FTps.hpp:707
FTps::taylor
FTps taylor(const Array1D< T > &series, int order) const
Taylor series.
Definition: FTps.hpp:1481
FTps::getMinOrder
int getMinOrder() const
Get minimum order.
Definition: FTps.h:165
FTps::getTruncOrder
int getTruncOrder() const
Get truncation order.
Definition: FTps.h:183
FTps::multiply
FTps multiply(const FTps &y, int trunc=EXACT) const
Multiplication.
Definition: FTps.hpp:650
DomainError
Domain error exception.
Definition: DomainError.h:32
LogicalError
Logical error exception.
Definition: LogicalError.h:33
Generated on Thu Oct 20 2022 17:40:10 for OPAL (Object Oriented Parallel Accelerator Library) by doxygen 1.9.3