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
34
namespace
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
41
constexpr
double
PI
=
M_PI
;
42
43
44
#ifndef M_E
/* e */
45
#define M_E 2.7182818284590452354
46
#endif
47
constexpr
double
E
=
M_E
;
48
49
50
#ifndef M_LOG2E
/* log_2 e */
51
#define M_LOG2E 1.4426950408889634074
52
#endif
53
constexpr
double
LOG2E
=
M_LOG2E
;
54
55
56
#ifndef M_LOG10E
/* log_10 e */
57
#define M_LOG10E 0.43429448190325182765
58
#endif
59
constexpr
double
LOG10E
=
M_LOG10E
;
60
61
62
#ifndef M_LN2
/* log_e 2 */
63
#define M_LN2 0.69314718055994530942
64
#endif
65
constexpr
double
LN2
=
M_LN2
;
66
67
68
#ifndef M_LN10
/* log_e 10 */
69
#define M_LN10 2.30258509299404568402
70
#endif
71
constexpr
double
LN10
=
M_LN10
;
72
73
74
#ifndef M_PI_2
/* pi/2 */
75
#define M_PI_2 1.57079632679489661923
76
#endif
77
constexpr
double
PI_2
=
M_PI_2
;
78
79
80
#ifndef M_PI_4
/* pi/4 */
81
#define M_PI_4 0.78539816339744830962
82
#endif
83
constexpr
double
PI_4
=
M_PI_4
;
84
85
86
#ifndef M_1_PI
/* 1/pi */
87
#define M_1_PI 0.31830988618379067154
88
#endif
89
constexpr
double
_1_PI
=
M_1_PI
;
90
91
92
#ifndef M_2_PI
/* 2/pi */
93
#define M_2_PI 0.63661977236758134308
94
#endif
95
constexpr
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
101
constexpr
double
_2_SQRTPI
=
M_2_SQRTPI
;
102
103
104
#ifndef M_SQRT2
/* sqrt(2) */
105
#define M_SQRT2 1.41421356237309504880
106
#endif
107
constexpr
double
SQRT2
=
M_SQRT2
;
108
109
110
#ifndef M_SQRT1_2
/* 1/sqrt(2) */
111
#define M_SQRT1_2 0.70710678118654752440
112
#endif
113
constexpr
double
SQRT1_2
=
M_SQRT1_2
;
114
115
#define _MATH_DEFINES_DEFINED
// to shut up MSVC warnings
116
117
constexpr
double
PI_DOUBLED
= 6.28318530717958647692;
118
120
template
<
typename
T>
121
struct
NaNImpl
{
122
static
constexpr
T
get
() {
return
std::numeric_limits<T>::quiet_NaN(); }
123
};
124
129
template
<
typename
T>
inline
constexpr
130
typename
std::remove_cv<typename std::remove_reference<T>::type>::type
NaN
() {
131
return
NaNImpl<typename std::remove_cv<typename std::remove_reference<T>::type
>::type>::get();
132
}
133
135
template
<
typename
T>
136
struct
NaNImpl
<
std
::
complex
<T>> {
137
static
constexpr
std::complex<T>
get
() {
return
std::complex<T>(
NaN<T>
(),
NaN<T>
()); }
138
};
139
140
142
template
<
typename
T>
143
struct
ZeroImpl
{
144
static
constexpr
T
get
() {
return
0.; }
145
};
146
151
template
<
typename
T>
inline
constexpr
152
typename
std::remove_cv<typename std::remove_reference<T>::type>::type
Zero
() {
153
return
ZeroImpl<typename std::remove_cv<typename std::remove_reference<T>::type
>::type>::get();
154
}
155
156
157
// size_t is preferred for array indexing
158
using
std::size_t;
159
using
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
177
inline
long
double
conj
(
long
double
x) {
return
x; }
178
inline
double
conj
(
double
x) {
return
x; }
179
inline
float
conj
(
float
x) {
return
x; }
180
181
// Limits for comparing approximate numbers with zero
182
constexpr
double
SMALL
= std::numeric_limits<double>::epsilon();
183
constexpr
double
SMALL2
=
SMALL
*
SMALL
;
184
189
inline
bool
is_zero
(
double
v
,
double
abs_supremum
=
SMALL
) {
190
return
abs(
v
) <
abs_supremum
;
191
}
192
195
inline
bool
is_zero
(dcomplex
v
) {
196
return
real
(
v
)*
real
(
v
) +
imag
(
v
)*
imag
(
v
) <
SMALL2
;
197
}
198
200
inline
bool
isnan
(dcomplex
v
) {
201
return
isnan
(
v
.real()) ||
isnan
(
v
.imag());
202
}
203
204
211
template
<
typename
T>
212
inline
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
218
inline
plask::dcomplex
operator*
(
int
a
,
const
plask::dcomplex&
b
) {
return
double
(
a
) *
b
; }
219
inline
plask::dcomplex
operator*
(
const
plask::dcomplex&
a
,
int
b
) {
return
a
*
double
(
b
); }
220
inline
plask::dcomplex
operator*
(
unsigned
a
,
const
plask::dcomplex&
b
) {
return
double
(
a
) *
b
; }
221
inline
plask::dcomplex
operator*
(
const
plask::dcomplex&
a
,
unsigned
b
) {
return
a
*
double
(
b
); }
222
inline
plask::dcomplex
operator/
(
int
a
,
const
plask::dcomplex&
b
) {
return
double
(
a
) /
b
; }
223
inline
plask::dcomplex
operator/
(
const
plask::dcomplex&
a
,
int
b
) {
return
a
/
double
(
b
); }
224
inline
plask::dcomplex
operator/
(
unsigned
a
,
const
plask::dcomplex&
b
) {
return
double
(
a
) /
b
; }
225
inline
plask::dcomplex
operator/
(
const
plask::dcomplex&
a
,
unsigned
b
) {
return
a
/
double
(
b
); }
226
227
228
// Useful functions
229
using
std::max;
using
std::min;
230
231
inline
double
abs2
(
const
dcomplex& x) {
232
return
real
(x)*
real
(x) +
imag
(x)*
imag
(x);
233
}
234
241
template
<
typename
T>
242
const
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
254
template
<
typename
T>
255
bool
in_range
(
const
T&
v
,
const
T&
beg
,
const
T& end) {
256
return
beg
<=
v
&&
v
< end;
257
}
258
265
template
<
typename
T1,
typename
T2,
typename
T3>
266
auto
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
270
inline
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
274
inline
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
278
inline
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__
287
typedef
double
v2double
__attribute__
((
vector_size
(16)));
289
typedef
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
302
inline
bool
dbl_compare_eq
(
double
x,
double
y) {
303
if
(
std::isnan
(x))
return
std::isnan
(y);
304
return
x == y;
305
}
306
315
inline
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
328
inline
bool
dbl_compare_gt
(
double
x,
double
y) {
return
dbl_compare_lt
(y, x); }
329
337
inline
bool
dbl_compare_lteq
(
double
x,
double
y) {
return
!
dbl_compare_gt
(x, y); }
338
346
inline
bool
dbl_compare_gteq
(
double
x,
double
y) {
return
!
dbl_compare_lt
(x, y); }
347
351
struct
PLASK_API
IllFormatedComplex
:
public
Exception
{
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
368
template
<
typename
T>
369
std::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
379
if
(
to_parse
)
throw
IllFormatedComplex
(
str_to_parse
);
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
406
throw
IllFormatedComplex
(
str_to_parse
);
407
check_eof
();
408
return
std::complex<T>(
real
,
imag
);
409
}
410
411
extern
template
PLASK_API
std::complex<double>
parse_complex<double>
(std::string
str_to_parse
);
412
413
422
class
AccurateSum
{
423
424
double
s, c;
// c is a running compensation for lost low-order bits.
425
426
public
:
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
plask
math.hpp
Generated by
1.9.8