OPAL (Object Oriented Parallel Accelerator Library) 2024.2
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#include "Physics/Physics.h"
27
28#include <algorithm>
29#include <cmath>
30#include <complex>
31#include <iostream>
32#include <type_traits>
33
34// 25. March 2017,
35// http://stackoverflow.com/questions/30736951/templated-class-check-if-complex-at-compile-time
36template<class T> struct is_complex : std::false_type {};
37template<class T> struct is_complex<std::complex<T>> : std::true_type {};
38
39
40// Class FTps; global functions acting on FTps objects.
41// Elementary functions acting on FTps<T,N> objects.
42// ------------------------------------------------------------------------
43
45template <class T, int N>
46FTps<T, N> pow(const FTps<T, N> &x, int y, int trunc = (FTps<T, N>::EXACT));
47
49template <class T, int N>
50FTps<T, N> sqrt(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
51
53template <class T, int N>
54FTps<T, N> sin(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
55
57template <class T, int N>
58FTps<T, N> cos(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
59
61template <class T, int N>
62FTps<T, N> tan(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
63
65template <class T, int N>
66FTps<T, N> cot(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
67
69template <class T, int N>
70FTps<T, N> sec(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
71
73template <class T, int N>
74FTps<T, N> csc(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
75
77template <class T, int N>
78FTps<T, N> exp(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
79
81template <class T, int N>
82FTps<T, N> log(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
83
85template <class T, int N>
86FTps<T, N> sinh(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
87
89template <class T, int N>
90FTps<T, N> cosh(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
91
93template <class T, int N>
94FTps<T, N> tanh(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
95
97template <class T, int N>
98FTps<T, N> coth(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
99
101template <class T, int N>
102FTps<T, N> sech(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
103
105template <class T, int N>
106FTps<T, N> csch(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
107
109template <class T, int N>
110FTps<T, N> erf(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
111
113template <class T, int N>
114FTps<T, N> erfc(const FTps<T, N> &x, int trunc = (FTps<T, N>::EXACT));
115
116
117// Implementation
118// ------------------------------------------------------------------------
119
120template <class T, int N>
121FTps<T, N> pow(const FTps<T, N> &x, int y, int trunc) {
122 // Default: trunc = EXACT
123
124 FTps<T, N> z(T(1));
125
126 if(y > 0) {
127 while(y-- > 0) z = z.multiply(x, trunc);
128 } else if(y < 0) {
129 if(x[0] == T(0) || x.getMinOrder() != 0)
130 throw DomainError("pow(const FTps &,int,int)");
131 FTps<T, N> t = x.inverse(trunc);
132 while(y++ < 0) z = z.multiply(t, trunc);
133 }
134
135 return z;
136}
137
138
139template <class T, int N>
140FTps<T, N> sqrt(const FTps<T, N> &x, int trunc) {
141 // Default: trunc = EXACT
142 int trcOrder = std::min(x.getTruncOrder(), trunc);
143
144 if(trcOrder == FTps<T, N>::EXACT)
145 throw LogicalError("::sqrt(FTps<T,N> &x, int trunc)",
146 "Square-root of EXACT polynomial must be truncated.");
147
148 T aZero = x[0];
149 if( ( std::real(aZero) <= std::real(T(0)) &&
150 std::imag(aZero) <= std::real(T(0))
151 ) || x.getMinOrder() != 0)
152 {
153 std::cerr << "FTps::sqrt(x) called with\nconstant term = " << aZero
154 << "; minOrder = " << x.getMinOrder() << std::endl
155 << "x = " << x << std::endl;
156 throw DomainError("sqrt(const FTps &,int)");
157 }
158
159 T two_aZero = T(2) * aZero;
160 Array1D<T> series(trcOrder + 1);
161 series[0] = sqrt(aZero);
162 for(int i = 1; i <= trcOrder; i++) {
163 series[i] = series[i-1] * double(3 - 2 * i) / (two_aZero * double(i));
164 }
165
166 return x.taylor(series, trcOrder);
167}
168
169
170template <class T, int N>
171FTps<T, N> sin(const FTps<T, N> &x, int trunc) {
172 // Default: trunc = EXACT
173 int trcOrder = std::min(x.getTruncOrder(), trunc);
174
175 if(trcOrder == FTps<T, N>::EXACT)
176 throw LogicalError("::sin(FTps<T,N> &x, int trunc)",
177 "Sine of EXACT polynomial must be truncated.");
178
179 T aZero = x[0];
180 if(x.getMinOrder() != 0) aZero = T(0);
181
182 Array1D<T> series(trcOrder + 1);
183 series[0] = sin(aZero);
184 series[1] = cos(aZero);
185 for(int i = 2; i <= trcOrder; i++) {
186 series[i] = - series[i-2] / double(i * (i - 1));
187 }
188
189 return x.taylor(series, trcOrder);
190}
191
192
193template <class T, int N>
194FTps<T, N> cos(const FTps<T, N> &x, int trunc) {
195 // Default: trunc = EXACT
196 int trcOrder = std::min(x.getTruncOrder(), trunc);
197
198 if(trcOrder == FTps<T, N>::EXACT)
199 throw LogicalError("::cos(FTps<T,N> &x, int trunc)",
200 "Cosine of EXACT polynomial must be truncated.");
201
202 T aZero = x[0];
203 if(x.getMinOrder() != 0) aZero = T(0);
204
205 Array1D<T> series(trcOrder + 1);
206 series[0] = cos(aZero);
207 series[1] = - sin(aZero);
208 for(int i = 2; i <= trcOrder; i++) {
209 series[i] = - series[i-2] / double(i * (i - 1));
210 }
211
212 return x.taylor(series, trcOrder);
213}
214
215
216template <class T, int N>
217FTps<T, N> tan(const FTps<T, N> &x, int trunc) {
218 // Default: trunc = EXACT
219
220 // The direct series expansion requires the Bernoulli numbers to arbitrary order.
221 return sin(x, trunc) / cos(x, trunc);
222}
223
224template <class T, int N>
225FTps<T, N> cot(const FTps<T, N> &x, int trunc) {
226 // Default: trunc = EXACT
227
228 if(x[0] == T(0) || x.getMinOrder() != 0)
229 throw DomainError("cot(const FTps &,int)");
230
231 return cos(x, trunc) / sin(x, trunc);
232}
233
234
235template <class T, int N>
236FTps<T, N> sec(const FTps<T, N> &x, int trunc) {
237 // Default: trunc = EXACT
238
239 return cos(x, trunc).inverse();
240}
241
242
243template <class T, int N>
244FTps<T, N> csc(const FTps<T, N> &x, int trunc) {
245 // Default: trunc = EXACT
246
247 if(x[0] == T(0) || x.getMinOrder() != 0)
248 throw DomainError("csc(const FTps &,int)");
249
250 return sin(x, trunc).inverse();
251}
252
253
254template <class T, int N>
255FTps<T, N> exp(const FTps<T, N> &x, int trunc) {
256 // Default: trunc = EXACT
257 int trcOrder = std::min(x.getTruncOrder(), trunc);
258
259 if(trcOrder == FTps<T, N>::EXACT)
260 throw LogicalError("::exp(FTps<T,N> &x, int trunc)",
261 "Exponential of EXACT polynomial must be truncated.");
262
263 T aZero = x[0];
264 if(x.getMinOrder() != 0) aZero = T(0);
265
266 Array1D<T> series(trcOrder + 1);
267 series[0] = exp(aZero);
268 for(int i = 1; i <= trcOrder; i++) {
269 series[i] = series[i-1] / double(i);
270 }
271
272 return x.taylor(series, trcOrder);
273}
274
275
276template <class T, int N>
277FTps<T, N> log(const FTps<T, N> &x, int trunc) {
278 // Default: trunc = EXACT
279 int trcOrder = std::min(x.getTruncOrder(), trunc);
280
281 if(trcOrder == FTps<T, N>::EXACT)
282 throw LogicalError("::log(FTps<T,N> &x, int trunc)",
283 "Logarithm of EXACT polynomial must be truncated.");
284
285 T aZero = x[0];
286 if(aZero <= T(0) || x.getMinOrder() != 0)
287 throw DomainError("log(const FTps &,int)");
288
289 T a0inv = T(1) / aZero;
290 T ain = a0inv;
291 Array1D<T> series(trcOrder + 1);
292 series[0] = std::log(aZero);
293 series[1] = a0inv;
294 for(int i = 2; i <= trcOrder; i++) {
295 ain *= -a0inv;
296 series[i] = ain / double(i);
297 }
298
299 return x.taylor(series, trcOrder);
300}
301
302
303template <class T, int N>
304FTps<T, N> sinh(const FTps<T, N> &x, int trunc) {
305 // Default: trunc = EXACT
306 int trcOrder = std::min(x.getTruncOrder(), trunc);
307
308 if(trcOrder == FTps<T, N>::EXACT)
309 throw LogicalError("::log(FTps<T,N> &x, int trunc)",
310 "Hyperbolic sine of EXACT polynomial must be truncated.");
311
312 T aZero = x[0];
313 if(x.getMinOrder() != 0) aZero = T(0);
314
315 Array1D<T> series(trcOrder + 1);
316 series[0] = sinh(aZero);
317 series[1] = cosh(aZero);
318 for(int i = 2; i <= trcOrder; i++) {
319 series[i] = series[i-2] / double(i * (i - 1));
320 }
321
322 return x.taylor(series, trcOrder);
323}
324
325
326template <class T, int N>
327FTps<T, N> cosh(const FTps<T, N> &x, int trunc) {
328 // Default: trunc = EXACT
329 int trcOrder = std::min(x.getTruncOrder(), trunc);
330
331 if(trcOrder == FTps<T, N>::EXACT)
332 throw LogicalError("::cosh(FTps<T,N> &x, int trunc)",
333 "Hyperbolic cosine of EXACT polynomial must be truncated.");
334
335 T aZero = x[0];
336 if(x.getMinOrder() != 0) aZero = T(0);
337
338 Array1D<T> series(trcOrder + 1);
339 series[0] = cosh(aZero);
340 series[1] = sinh(aZero);
341 for(int i = 2; i <= trcOrder; i++) {
342 series[i] = series[i-2] / double(i * (i - 1));
343 }
344
345 return x.taylor(series, trcOrder);
346}
347
348
349template <class T, int N>
350FTps<T, N> tanh(const FTps<T, N> &x, int trunc) {
351 // Default: trunc = EXACT
352
353 return sinh(x, trunc) / cosh(x, trunc);
354}
355
356
357template <class T, int N>
358FTps<T, N> coth(const FTps<T, N> &x, int trunc) {
359 // Default: trunc = EXACT
360
361 if(x[0] == T(0) || x.getMinOrder() != 0)
362 throw DomainError("coth(const FTps &,int)");
363
364 return cosh(x, trunc) / sinh(x, trunc);
365}
366
367
368template <class T, int N>
369FTps<T, N> sech(const FTps<T, N> &x, int trunc) {
370 // Default: trunc = EXACT
371
372 return cosh(x, trunc).inverse();
373}
374
375
376template <class T, int N>
377FTps<T, N> csch(const FTps<T, N> &x, int trunc) {
378 // Default: trunc = EXACT
379
380 if(x[0] == T(0) || x.getMinOrder() != 0)
381 throw DomainError("csch(const FTps &,int)");
382
383 return sinh(x, trunc).inverse();
384}
385
387template <class T, int N>
388FTps<T, N> erf(const FTps<T, N> &x, int trunc) {
389
390 if ( is_complex<T>::value )
391 throw LogicalError("::erf(FTps<T,N> &x, int trunc)",
392 "Error function does not support complex numbers.");
393
394 // Default: trunc = EXACT
395 int trcOrder = std::min(x.getTruncOrder(), trunc);
396
397 if(trcOrder == FTps<T, N>::EXACT)
398 throw LogicalError("::erf(FTps<T,N> &x, int trunc)",
399 "Error function of EXACT polynomial must be truncated.");
400
401 T aZero = x[0];
402 if(x.getMinOrder() != 0) aZero = T(0);
403
404 Array1D<T> series(trcOrder + 1);
405 series[0] = std::erf(std::real(aZero));
406 series[1] = 2.0 / std::sqrt(Physics::pi) * std::exp(-aZero*aZero);
407
408 for(int i = 2; i <= trcOrder; ++i) {
409 series[i] = - 2.0 / double(i-1) * double((i-2)) * series[i-2] / double(i);
410 }
411
412 return x.taylor(series, trcOrder);
413}
414
416template <class T, int N>
417FTps<T, N> erfc(const FTps<T, N> &x, int trunc) {
418 // Default: trunc = EXACT
419
420 return T(1) - erf(x, trunc);
421}
422
423#endif // CLASSIC_FTpsMath_HH
csc
FTps< T, N > csc(const FTps< T, N > &x, int trunc=(FTps< T, N >::EXACT))
Cosecant.
Definition FTpsMath.h:244
exp
FTps< T, N > exp(const FTps< T, N > &x, int trunc=(FTps< T, N >::EXACT))
Exponential.
Definition FTpsMath.h:255
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:121
cot
FTps< T, N > cot(const FTps< T, N > &x, int trunc=(FTps< T, N >::EXACT))
Cotangent.
Definition FTpsMath.h:225
sqrt
FTps< T, N > sqrt(const FTps< T, N > &x, int trunc=(FTps< T, N >::EXACT))
Square root.
Definition FTpsMath.h:140
sech
FTps< T, N > sech(const FTps< T, N > &x, int trunc=(FTps< T, N >::EXACT))
Hyperbolic secant.
Definition FTpsMath.h:369
erf
FTps< T, N > erf(const FTps< T, N > &x, int trunc=(FTps< T, N >::EXACT))
Error function.
Definition FTpsMath.h:388
sin
FTps< T, N > sin(const FTps< T, N > &x, int trunc=(FTps< T, N >::EXACT))
Sine.
Definition FTpsMath.h:171
erfc
FTps< T, N > erfc(const FTps< T, N > &x, int trunc=(FTps< T, N >::EXACT))
Complementary error function.
Definition FTpsMath.h:417
tan
FTps< T, N > tan(const FTps< T, N > &x, int trunc=(FTps< T, N >::EXACT))
Tangent.
Definition FTpsMath.h:217
log
FTps< T, N > log(const FTps< T, N > &x, int trunc=(FTps< T, N >::EXACT))
Natural logarithm.
Definition FTpsMath.h:277
csch
FTps< T, N > csch(const FTps< T, N > &x, int trunc=(FTps< T, N >::EXACT))
Hyperbolic cosecant.
Definition FTpsMath.h:377
tanh
FTps< T, N > tanh(const FTps< T, N > &x, int trunc=(FTps< T, N >::EXACT))
Hyperbolic tangent.
Definition FTpsMath.h:350
sinh
FTps< T, N > sinh(const FTps< T, N > &x, int trunc=(FTps< T, N >::EXACT))
Hyperbolic sine.
Definition FTpsMath.h:304
sec
FTps< T, N > sec(const FTps< T, N > &x, int trunc=(FTps< T, N >::EXACT))
Secant.
Definition FTpsMath.h:236
cos
FTps< T, N > cos(const FTps< T, N > &x, int trunc=(FTps< T, N >::EXACT))
Cosine.
Definition FTpsMath.h:194
cosh
FTps< T, N > cosh(const FTps< T, N > &x, int trunc=(FTps< T, N >::EXACT))
Hyperbolic cosine.
Definition FTpsMath.h:327
coth
FTps< T, N > coth(const FTps< T, N > &x, int trunc=(FTps< T, N >::EXACT))
Hyperbolic cotangent.
Definition FTpsMath.h:358
Physics::pi
constexpr double pi
The value of.
Definition Physics.h:30
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 Wed Sep 30 2026 07:14:51 for OPAL (Object Oriented Parallel Accelerator Library) by doxygen 1.9.8