math.hpp

Go to the documentation of this file.

/*
 * This file is part of PLaSK (https://plask.app) by Photonics Group at TUL
 * Copyright (c) 2022 Lodz University of Technology
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, version 3.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 * GNU General Public License for more details.
 */
#ifndef PLASK__MATH_H
#define PLASK__MATH_H
 
#include <plask/config.hpp>
 
#ifdef _MSC_VER
#   define _USE_MATH_DEFINES   // gives M_PI, etc. in <cmath>
#endif
 
#include <cmath>
#include <limits>
#include <algorithm>
 
#include "exceptions.hpp"
#include <sstream>
 
#ifdef PLASK_MATH_STD
#   include <complex>
#endif // PLASK_MATH_STD
 
namespace plask {
 
// std. libraries of some compilers, but potentially not all, defines this in cmath (as non-standard extension),
// for other compilers, numbers are copy/pasted from gcc's math.h:
#ifndef M_PI     /* pi */
    #define M_PI 3.14159265358979323846
#endif
constexpr double PI = M_PI;
 
 
#ifndef M_E      /* e */
    #define M_E 2.7182818284590452354
#endif
constexpr double E = M_E;
 
 
#ifndef M_LOG2E  /* log_2 e */
    #define M_LOG2E 1.4426950408889634074
#endif
constexpr double LOG2E = M_LOG2E;
 
 
#ifndef M_LOG10E /* log_10 e */
    #define M_LOG10E 0.43429448190325182765
#endif
constexpr double LOG10E = M_LOG10E;
 
 
#ifndef M_LN2 /* log_e 2 */
    #define M_LN2 0.69314718055994530942
#endif
constexpr double LN2 = M_LN2;
 
 
#ifndef M_LN10 /* log_e 10 */
    #define M_LN10 2.30258509299404568402
#endif
constexpr double LN10 = M_LN10;
 
 
#ifndef M_PI_2 /* pi/2 */
    #define M_PI_2 1.57079632679489661923
#endif
constexpr double PI_2 = M_PI_2;
 
 
#ifndef M_PI_4 /* pi/4 */
    #define M_PI_4 0.78539816339744830962
#endif
constexpr double PI_4 = M_PI_4;
 
 
#ifndef M_1_PI /* 1/pi */
    #define M_1_PI 0.31830988618379067154
#endif
constexpr double _1_PI = M_1_PI;
 
 
#ifndef M_2_PI /* 2/pi */
    #define M_2_PI 0.63661977236758134308
#endif
constexpr double _2_PI = M_2_PI;
 
 
#ifndef M_2_SQRTPI /* 2/sqrt(pi) */
    #define M_2_SQRTPI 1.12837916709551257390
#endif
constexpr double _2_SQRTPI = M_2_SQRTPI;
 
 
#ifndef M_SQRT2 /* sqrt(2) */
    #define M_SQRT2 1.41421356237309504880
#endif
constexpr double SQRT2 = M_SQRT2;
 
 
#ifndef M_SQRT1_2 /* 1/sqrt(2) */
    #define M_SQRT1_2 0.70710678118654752440
#endif
constexpr double SQRT1_2 = M_SQRT1_2;
 
#define _MATH_DEFINES_DEFINED   // to shut up MSVC warnings
 
constexpr double PI_DOUBLED = 6.28318530717958647692;
 
template <typename T>
struct NaNImpl {
    static constexpr T get() { return std::numeric_limits<T>::quiet_NaN(); }
};
 
template <typename T> inline constexpr
typename std::remove_cv<typename std::remove_reference<T>::type>::type NaN() {
    return NaNImpl<typename std::remove_cv<typename std::remove_reference<T>::type>::type>::get();
}
 
template <typename T>
struct NaNImpl<std::complex<T>> {
    static constexpr std::complex<T> get() { return std::complex<T>(NaN<T>(), NaN<T>()); }
};
 
 
template <typename T>
struct ZeroImpl {
    static constexpr T get() { return 0.; }
};
 
template <typename T> inline constexpr
typename std::remove_cv<typename std::remove_reference<T>::type>::type Zero() {
    return ZeroImpl<typename std::remove_cv<typename std::remove_reference<T>::type>::type>::get();
}
 
 
// size_t is preferred for array indexing
using std::size_t;
using std::ptrdiff_t;
 
// Complex numbers library
#ifdef PLASK_MATH_STD
    using std::complex;
    using std::conj;
    using std::sqrt;
    using std::abs; using std::real; using std::imag;
    using std::log; using std::exp;
    using std::sin; using std::cos; using std::tan;
    using std::sinh; using std::cosh; using std::tanh;
    using std::asin; using std::acos; using std::atan; using std::atan2;
    using std::asinh; using std::acosh; using std::atanh;
    using std::isnan; using std::isinf;
    typedef complex<double> dcomplex;
    const dcomplex I(0.,1.);
#endif // PLASK_MATH_STD
 
inline long double conj(long double x) { return x; }
inline double conj(double x) { return x; }
inline float conj(float x) { return x; }
 
// Limits for comparing approximate numbers with zero
constexpr double SMALL = std::numeric_limits<double>::epsilon(); 
constexpr double SMALL2 = SMALL*SMALL; 
 
inline bool is_zero(double v, double abs_supremum = SMALL) {
    return abs(v) < abs_supremum;
}
 
inline bool is_zero(dcomplex v) {
    return real(v)*real(v) + imag(v)*imag(v) < SMALL2;
}
 
inline bool isnan(dcomplex v) {
    return isnan(v.real()) || isnan(v.imag());
}
 
 
template <typename T>
inline T remove_nan(T val, const T nan=Zero<T>()) {
    return isnan(val)? nan : val;
}
 
 
// C++ is lacking some operators
inline plask::dcomplex operator*(int a, const plask::dcomplex& b) { return double(a) * b; }
inline plask::dcomplex operator*(const plask::dcomplex& a, int b) { return a * double(b); }
inline plask::dcomplex operator*(unsigned a, const plask::dcomplex& b) { return double(a) * b; }
inline plask::dcomplex operator*(const plask::dcomplex& a, unsigned b) { return a * double(b); }
inline plask::dcomplex operator/(int a, const plask::dcomplex& b) { return double(a) / b; }
inline plask::dcomplex operator/(const plask::dcomplex& a, int b) { return a / double(b); }
inline plask::dcomplex operator/(unsigned a, const plask::dcomplex& b) { return double(a) / b; }
inline plask::dcomplex operator/(const plask::dcomplex& a, unsigned b) { return a / double(b); }
 
 
// Useful functions
using std::max; using std::min;
 
inline double abs2(const dcomplex& x) {
    return real(x)*real(x) + imag(x)*imag(x);
}
 
template <typename T>
const T& clamp(const T& v, const T& min, const T& max) {
    if (v < min) return min;
    if (v > max) return max;
    return v;
}
 
template <typename T>
bool in_range(const T& v, const T& beg, const T& end) {
    return beg <= v && v < end;
}
 
template <typename T1, typename T2, typename T3>
auto inline fma(T1 to_mult_1, T2 to_mult_2, T3 to_sum) -> decltype(to_mult_1*to_mult_2+to_sum) {
    return to_mult_1 * to_mult_2 + to_sum;
}
 
inline float fma(float to_mult_1, float to_mult_2, float to_sum) {
    return std::fma(to_mult_1, to_mult_2, to_sum);
}
 
inline double fma(double to_mult_1, double to_mult_2, double to_sum) {
    return std::fma(to_mult_1, to_mult_2, to_sum);
}
 
inline long double fma(long double to_mult_1, long double to_mult_2, long double to_sum) {
    return std::fma(to_mult_1, to_mult_2, to_sum);
}
 
#ifdef __GNUC__
typedef double v2double __attribute__((vector_size(16)));
typedef union { v2double simd; double v[2]; } v2double_view;
#endif
 
 
// Total order double comparison with NaN greater than all other numbers:
 
inline bool dbl_compare_eq(double x, double y) {
    if (std::isnan(x)) return std::isnan(y);
    return x == y;
}
 
inline bool dbl_compare_lt(double x, double y) {
      if (std::isnan(y)) return !std::isnan(x);
      return x < y; // if x is NaN it is grater than non-NaN y and std. < operator returns false
}
 
inline bool dbl_compare_gt(double x, double y) { return dbl_compare_lt(y, x); }
 
inline bool dbl_compare_lteq(double x, double y) { return !dbl_compare_gt(x, y); }
 
inline bool dbl_compare_gteq(double x, double y) { return !dbl_compare_lt(x, y); }
 
struct PLASK_API IllFormatedComplex: public Exception {
 
    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)
    {}
 
};
 
template <typename T>
std::complex<T> parse_complex(std::string str_to_parse) {
    boost::trim(str_to_parse);
    if (!str_to_parse.empty() && str_to_parse[0] == '(' && str_to_parse[str_to_parse.length()-1] == ')' &&
        str_to_parse.find(',') == std::string::npos)
        str_to_parse = str_to_parse.substr(1, str_to_parse.length()-2);
    std::istringstream to_parse(str_to_parse);
    auto check_eof = [&] () {
        if (!to_parse.eof()) {  // we require end-of stream here
            char c;
            to_parse >> c;  // we check if there is non-white character, this operation should fail
            if (to_parse) throw IllFormatedComplex(str_to_parse);
        }
    };
    T real, imag;
    to_parse >> real;
    if (to_parse.fail()) {  // we will try standard >> operator
        to_parse.clear(); to_parse.str(str_to_parse);
        std::complex<T> res;
        to_parse >> res;
        if (to_parse.fail()) throw IllFormatedComplex(str_to_parse);
        check_eof();
        return res;
    }
    if (to_parse.eof()) return std::complex<T>(real);
    char c;
    to_parse >> c;
    if (to_parse.fail()) throw IllFormatedComplex(str_to_parse);
    if (to_parse.eof()) return std::complex<T>(real);
    if (c == 'i' || c == 'j') { // only imaginary part is given
        imag = real;
        real = 0.0;
    } else if (c == '+' || c == '-') {
        char c_ij;
        to_parse >> imag >> c_ij;
        if (to_parse.fail() || (c_ij != 'i' && c_ij != 'j')) throw IllFormatedComplex(str_to_parse);
        if (c == '-') imag = -imag;
    } else
         throw IllFormatedComplex(str_to_parse);
    check_eof();
    return std::complex<T>(real, imag);
}
 
extern template PLASK_API std::complex<double> parse_complex<double>(std::string str_to_parse);
 
 
class AccurateSum {
 
    double s, c;    // c is a running compensation for lost low-order bits.
 
public:
 
    AccurateSum(double initial = 0.0): s(initial), c(0.0) {}
 
    AccurateSum& operator=(double v);
 
    AccurateSum(const AccurateSum& initial) = default;
    AccurateSum& operator=(const AccurateSum& v) = default;
 
    AccurateSum& operator+=(double v);
    AccurateSum& operator-=(double v) { return *this += -v; }
 
    operator double() const { return s; }
 
    AccurateSum& operator+=(const AccurateSum& other);
 
};
 
 
} // namespace plask
 
#endif // PLASK__MATH_H