PLaSK library
Loading...
Searching...
No Matches
math.hpp
Go to the documentation of this file.
1/*
2 * This file is part of PLaSK (https://plask.app) by Photonics Group at TUL
3 * Copyright (c) 2022 Lodz University of Technology
4 *
5 * This program is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, version 3.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 */
14#ifndef PLASK__MATH_H
15#define PLASK__MATH_H
16
17#include <plask/config.hpp>
18
19#ifdef _MSC_VER
20# define _USE_MATH_DEFINES // gives M_PI, etc. in <cmath>
21#endif
22
23#include <cmath>
24#include <limits>
25#include <algorithm>
26
27#include "exceptions.hpp"
28#include <sstream>
29
30#ifdef PLASK_MATH_STD
31# include <complex>
32#endif // PLASK_MATH_STD
33
34namespace plask {
35
36// std. libraries of some compilers, but potentially not all, defines this in cmath (as non-standard extension),
37// for other compilers, numbers are copy/pasted from gcc's math.h:
38#ifndef M_PI /* pi */
39 #define M_PI 3.14159265358979323846
40#endif
41constexpr double PI = M_PI;
42
43
44#ifndef M_E /* e */
45 #define M_E 2.7182818284590452354
46#endif
47constexpr double E = M_E;
48
49
50#ifndef M_LOG2E /* log_2 e */
51 #define M_LOG2E 1.4426950408889634074
52#endif
53constexpr double LOG2E = M_LOG2E;
54
55
56#ifndef M_LOG10E /* log_10 e */
57 #define M_LOG10E 0.43429448190325182765
58#endif
59constexpr double LOG10E = M_LOG10E;
60
61
62#ifndef M_LN2 /* log_e 2 */
63 #define M_LN2 0.69314718055994530942
64#endif
65constexpr double LN2 = M_LN2;
66
67
68#ifndef M_LN10 /* log_e 10 */
69 #define M_LN10 2.30258509299404568402
70#endif
71constexpr double LN10 = M_LN10;
72
73
74#ifndef M_PI_2 /* pi/2 */
75 #define M_PI_2 1.57079632679489661923
76#endif
77constexpr double PI_2 = M_PI_2;
78
79
80#ifndef M_PI_4 /* pi/4 */
81 #define M_PI_4 0.78539816339744830962
82#endif
83constexpr double PI_4 = M_PI_4;
84
85
86#ifndef M_1_PI /* 1/pi */
87 #define M_1_PI 0.31830988618379067154
88#endif
89constexpr double _1_PI = M_1_PI;
90
91
92#ifndef M_2_PI /* 2/pi */
93 #define M_2_PI 0.63661977236758134308
94#endif
95constexpr double _2_PI = M_2_PI;
96
97
98#ifndef M_2_SQRTPI /* 2/sqrt(pi) */
99 #define M_2_SQRTPI 1.12837916709551257390
100#endif
101constexpr double _2_SQRTPI = M_2_SQRTPI;
102
103
104#ifndef M_SQRT2 /* sqrt(2) */
105 #define M_SQRT2 1.41421356237309504880
106#endif
107constexpr double SQRT2 = M_SQRT2;
108
109
110#ifndef M_SQRT1_2 /* 1/sqrt(2) */
111 #define M_SQRT1_2 0.70710678118654752440
112#endif
113constexpr double SQRT1_2 = M_SQRT1_2;
114
115#define _MATH_DEFINES_DEFINED // to shut up MSVC warnings
116
117constexpr double PI_DOUBLED = 6.28318530717958647692;
118
120template <typename T>
121struct NaNImpl {
122 static constexpr T get() { return std::numeric_limits<T>::quiet_NaN(); }
123};
124
129template <typename T> inline constexpr
130typename std::remove_cv<typename std::remove_reference<T>::type>::type NaN() {
132}
133
135template <typename T>
136struct NaNImpl<std::complex<T>> {
137 static constexpr std::complex<T> get() { return std::complex<T>(NaN<T>(), NaN<T>()); }
138};
139
140
142template <typename T>
143struct ZeroImpl {
144 static constexpr T get() { return 0.; }
145};
146
151template <typename T> inline constexpr
152typename std::remove_cv<typename std::remove_reference<T>::type>::type Zero() {
154}
155
156
157// size_t is preferred for array indexing
158using std::size_t;
159using std::ptrdiff_t;
160
161// Complex numbers library
162#ifdef PLASK_MATH_STD
163 using std::complex;
164 using std::conj;
165 using std::sqrt;
166 using std::abs; using std::real; using std::imag;
167 using std::log; using std::exp;
168 using std::sin; using std::cos; using std::tan;
169 using std::sinh; using std::cosh; using std::tanh;
170 using std::asin; using std::acos; using std::atan; using std::atan2;
171 using std::asinh; using std::acosh; using std::atanh;
172 using std::isnan; using std::isinf;
173 typedef complex<double> dcomplex;
174 const dcomplex I(0.,1.);
175#endif // PLASK_MATH_STD
176
177inline long double conj(long double x) { return x; }
178inline double conj(double x) { return x; }
179inline float conj(float x) { return x; }
180
181// Limits for comparing approximate numbers with zero
182constexpr double SMALL = std::numeric_limits<double>::epsilon();
183constexpr double SMALL2 = SMALL*SMALL;
184
189inline bool is_zero(double v, double abs_supremum = SMALL) {
190 return abs(v) < abs_supremum;
191}
192
195inline bool is_zero(dcomplex v) {
196 return real(v)*real(v) + imag(v)*imag(v) < SMALL2;
197}
198
200inline bool isnan(dcomplex v) {
201 return isnan(v.real()) || isnan(v.imag());
202}
203
204
211template <typename T>
212inline T remove_nan(T val, const T nan=Zero<T>()) {
213 return isnan(val)? nan : val;
214}
215
216
217// C++ is lacking some operators
218inline plask::dcomplex operator*(int a, const plask::dcomplex& b) { return double(a) * b; }
219inline plask::dcomplex operator*(const plask::dcomplex& a, int b) { return a * double(b); }
220inline plask::dcomplex operator*(unsigned a, const plask::dcomplex& b) { return double(a) * b; }
221inline plask::dcomplex operator*(const plask::dcomplex& a, unsigned b) { return a * double(b); }
222inline plask::dcomplex operator/(int a, const plask::dcomplex& b) { return double(a) / b; }
223inline plask::dcomplex operator/(const plask::dcomplex& a, int b) { return a / double(b); }
224inline plask::dcomplex operator/(unsigned a, const plask::dcomplex& b) { return double(a) / b; }
225inline plask::dcomplex operator/(const plask::dcomplex& a, unsigned b) { return a / double(b); }
226
227
228// Useful functions
229using std::max; using std::min;
230
231inline double abs2(const dcomplex& x) {
232 return real(x)*real(x) + imag(x)*imag(x);
233}
234
241template <typename T>
242const T& clamp(const T& v, const T& min, const T& max) {
243 if (v < min) return min;
244 if (v > max) return max;
245 return v;
246}
247
254template <typename T>
255bool in_range(const T& v, const T& beg, const T& end) {
256 return beg <= v && v < end;
257}
258
265template <typename T1, typename T2, typename T3>
266auto inline fma(T1 to_mult_1, T2 to_mult_2, T3 to_sum) -> decltype(to_mult_1*to_mult_2+to_sum) {
267 return to_mult_1 * to_mult_2 + to_sum;
268}
269
270inline float fma(float to_mult_1, float to_mult_2, float to_sum) {
271 return std::fma(to_mult_1, to_mult_2, to_sum);
272}
273
274inline double fma(double to_mult_1, double to_mult_2, double to_sum) {
275 return std::fma(to_mult_1, to_mult_2, to_sum);
276}
277
278inline long double fma(long double to_mult_1, long double to_mult_2, long double to_sum) {
279 return std::fma(to_mult_1, to_mult_2, to_sum);
280}
281
282#ifdef __GNUC__
287typedef double v2double __attribute__((vector_size(16)));
289typedef union { v2double simd; double v[2]; } v2double_view;
290#endif
291
292
293// Total order double comparison with NaN greater than all other numbers:
294
302inline bool dbl_compare_eq(double x, double y) {
303 if (std::isnan(x)) return std::isnan(y);
304 return x == y;
305}
306
315inline bool dbl_compare_lt(double x, double y) {
316 if (std::isnan(y)) return !std::isnan(x);
317 return x < y; // if x is NaN it is grater than non-NaN y and std. < operator returns false
318}
319
328inline bool dbl_compare_gt(double x, double y) { return dbl_compare_lt(y, x); }
329
337inline bool dbl_compare_lteq(double x, double y) { return !dbl_compare_gt(x, y); }
338
346inline bool dbl_compare_gteq(double x, double y) { return !dbl_compare_lt(x, y); }
347
352
357 IllFormatedComplex(const std::string& str_to_parse): Exception("Ill-formatted complex number \"{0}\". Allowed formats: 'R+Ij', 'R,Ij', '(R, I)', where R and I are floating point numbers.", str_to_parse)
358 {}
359
360};
361
368template <typename T>
369std::complex<T> parse_complex(std::string str_to_parse) {
370 boost::trim(str_to_parse);
371 if (!str_to_parse.empty() && str_to_parse[0] == '(' && str_to_parse[str_to_parse.length()-1] == ')' &&
372 str_to_parse.find(',') == std::string::npos)
373 str_to_parse = str_to_parse.substr(1, str_to_parse.length()-2);
374 std::istringstream to_parse(str_to_parse);
375 auto check_eof = [&] () {
376 if (!to_parse.eof()) { // we require end-of stream here
377 char c;
378 to_parse >> c; // we check if there is non-white character, this operation should fail
380 }
381 };
382 T real, imag;
383 to_parse >> real;
384 if (to_parse.fail()) { // we will try standard >> operator
385 to_parse.clear(); to_parse.str(str_to_parse);
386 std::complex<T> res;
387 to_parse >> res;
388 if (to_parse.fail()) throw IllFormatedComplex(str_to_parse);
389 check_eof();
390 return res;
391 }
392 if (to_parse.eof()) return std::complex<T>(real);
393 char c;
394 to_parse >> c;
395 if (to_parse.fail()) throw IllFormatedComplex(str_to_parse);
396 if (to_parse.eof()) return std::complex<T>(real);
397 if (c == 'i' || c == 'j') { // only imaginary part is given
398 imag = real;
399 real = 0.0;
400 } else if (c == '+' || c == '-') {
401 char c_ij;
402 to_parse >> imag >> c_ij;
403 if (to_parse.fail() || (c_ij != 'i' && c_ij != 'j')) throw IllFormatedComplex(str_to_parse);
404 if (c == '-') imag = -imag;
405 } else
407 check_eof();
408 return std::complex<T>(real, imag);
409}
410
411extern template PLASK_API std::complex<double> parse_complex<double>(std::string str_to_parse);
412
413
423
424 double s, c; // c is a running compensation for lost low-order bits.
425
426public:
427
428 AccurateSum(double initial = 0.0): s(initial), c(0.0) {}
429
430 AccurateSum& operator=(double v);
431
432 AccurateSum(const AccurateSum& initial) = default;
433 AccurateSum& operator=(const AccurateSum& v) = default;
434
435 AccurateSum& operator+=(double v);
436 AccurateSum& operator-=(double v) { return *this += -v; }
437
438 operator double() const { return s; }
439
440 AccurateSum& operator+=(const AccurateSum& other);
441
442};
443
444
445} // namespace plask
446
447#endif // PLASK__MATH_H