numutils.hpp

Go to the documentation of this file.

#ifndef MEDUSA_BITS_UTILS_NUMUTILS_HPP_
#define MEDUSA_BITS_UTILS_NUMUTILS_HPP_
 
#include <medusa/Config.hpp>
#include <medusa/bits/utils/assert.hpp>
#include <medusa/bits/types/Vec.hpp>
#include <medusa/bits/types/Range.hpp>
#include <cmath>
 
namespace mm {
 
template<typename T>
int iceil(T x) { return static_cast<int>(std::ceil(x)); }
 
template<typename T>
int ifloor(T x) { return static_cast<int>(std::floor(x)); }
 
template <unsigned int exponent>
inline double ipow(double base) {
    return ipow<exponent - 1>(base) * base;
}
 
template <>
inline double ipow<0>(double) { return 1; }
 
template <typename T>
inline T ipow(T b, int e) {
    T r = 1.0;
    while (e > 0) {
        r *= b;
        --e;
    }
    return r;
}
 
template <typename T>
inline T ipowneg(T b, int e) {
    if (e < 0) return 1.0 / ipow(b, -e);
    else return ipow(b, e);
}
 
template <typename T>
inline constexpr int signum(T x, std::false_type) {
    return T(0) < x;
}
template <typename T>
inline constexpr int signum(T x, std::true_type) {
    return (T(0) < x) - (x < T(0));
}
template <typename T>
inline constexpr int signum(T x) {
    return signum(x, std::is_signed<T>());
}
 
template <int dim>
bool incrementCounter(Vec<int, dim>& counter, const Vec<int, dim>& limit) {
    for (int i = dim - 1; i >= 0; --i) {
        if (counter[i] >= limit[i] - 1) {
            counter[i] = 0;
        } else {
            counter[i]++;
            return true;
        }
    }
    return false;
}
 
template <int dim>
bool incrementCounter(Vec<int, dim>& counter, const Vec<int, dim>& low, const Vec<int, dim>& high) {
    for (int i = dim - 1; i >= 0; --i) {
        if (counter[i] >= high[i] - 1) {
            counter[i] = low[i];
        } else {
            counter[i]++;
            return true;
        }
    }
    return false;
}
 
template <int dim>
bool incrementCyclicCounter(Vec<int, dim>& counter, const Vec<int, dim>& low,
                            const Vec<int, dim>& high, const Vec<int, dim>& size) {
    for (int i = dim - 1; i >= 0; --i) {
        if (counter[i] == high[i] - 1 || (high[i] == 0 && counter[i] == size[i] - 1)) {
            counter[i] = low[i];
        } else {
            counter[i]++;
            if (counter[i] >= size[i]) counter[i] = 0;
            return true;
        }
    }
    return false;
}
 
template <class scalar_t, int dim>
Range<Vec<scalar_t, dim>> linspace(const Vec<scalar_t, dim>& beg, const Vec<scalar_t, dim>& end,
                                   const Vec<int, dim>& counts,
                                   const Vec<bool, dim> include_boundary = true) {
    Range<Vec<scalar_t, dim>> ret;
    for (int i = 0; i < dim; ++i) {
        if (include_boundary[i]) assert(counts[i] >= 2);
        else assert(counts[i] >= 0);
        if (counts[i] == 0) return ret;
    }
    Vec<scalar_t, dim> step;
    for (int i = 0; i < dim; ++i) {
        if (include_boundary[i]) step[i] = (end[i] - beg[i]) / (counts[i] - 1);
        else step[i] = (end[i] - beg[i]) / (counts[i] + 1);
    }
    Vec<int, dim> counter = 0;
    do {
        Vec<scalar_t, dim> tmp;
        for (int i = 0; i < dim; ++i) {
            if (include_boundary[i]) tmp[i] = beg[i] + counter[i] * step[i];
            else tmp[i] = beg[i] + (counter[i] + 1) * step[i];
        }
        ret.push_back(tmp);
    } while (incrementCounter(counter, counts));
    return ret;
}
 
template <class scalar_t, int dim>
Range<Vec<scalar_t, dim>> linspace(const Vec<scalar_t, dim>& beg, const Vec<scalar_t, dim>& end,
                                   const Vec<int, dim>& counts, bool include_boundary) {
    return linspace(beg, end, counts, Vec<bool, dim>(include_boundary));
}
 
template <typename scalar_t>
Range<scalar_t> linspace(scalar_t beg, scalar_t end, int count, bool include_boundary = true) {
    Range<scalar_t> ret;
    if (include_boundary) {
        assert_msg(count >= 2, "Count must be >= 2, got %d.", count);
        scalar_t step = (end - beg) / (count - 1);
        for (int i = 0; i < count; ++i) {
            ret.push_back(beg + step * i);
        }
    } else {
        assert_msg(count >= 0, "Count must be >= 0, got %d.", count);
        scalar_t step = (end - beg) / (count + 1);
        for (int i = 1; i <= count; ++i) {
            ret.push_back(beg + step * i);
        }
    }
    return ret;
}
 
template<typename function_t, typename input_t, typename output_t, bool verbose = false>
input_t bisection(const function_t& f, input_t lo, input_t hi, output_t target = 0.0,
                  input_t tolerance = 1e-4, int max_iter = 40) {
    input_t a = lo, d = hi - lo;
    output_t fa = f(lo) - target;
    output_t fb = f(hi) - target;
 
    assert_msg(signum(fa) != signum(fb),
               "Function values have the same sign, f(a) = %s, f(b) = %s.", fa, fb);
 
    while (d > tolerance) {
        d /= static_cast<input_t>(2.);
        output_t fc = f(a+d) - target;
 
        if (verbose) std::cout << "f(" << a+d << ")=" << fc + target << "\n";
        if (signum(fa) == signum(fc)) {
            a += d;
            fa = fc;
        }
        --max_iter;
        if (max_iter < 0) {
            std::cerr << "\nBisection ended without finding the solution after max_num"
                    " iterations." << std::endl;
            break;
        }
    }
    return a+d;
}
 
}  // namespace mm
 
#endif  // MEDUSA_BITS_UTILS_NUMUTILS_HPP_