misc.h

/*
 * ########################################################################
 * The contents of this file is free and unencumbered software released into the
 * public domain. For more information, please refer to <http://unlicense.org/>
 * ########################################################################
 */
 
#ifndef BSPLINE_MISC_H
#define BSPLINE_MISC_H
 
#include <array>
 
#ifndef BSPLINE_DOXYGEN_IGNORE
namespace bspline::internal {
 
template <typename T, size_t size>
std::array<T, size> make_array(T val) {
  std::array<T, size> ret;
  ret.fill(val);
  return ret;
}
 
template <typename T, size_t sizea, size_t sizeb>
std::array<T, std::max(sizea, sizeb)> add(const std::array<T, sizea> &a,
                                          const std::array<T, sizeb> &b) {
  if constexpr (sizeb > sizea) {
    return add(b, a);
  } else {
    std::array<T, sizea> ret = a;
    for (size_t i = 0; i < sizeb; i++) {
      ret[i] += b[i];
    }
    return ret;
  }
}
 
template <typename T, size_t sizein, size_t sizeout>
std::array<T, sizeout> changearraysize(const std::array<T, sizein> &in) {
  static_assert(sizeout >= sizein,
                "sizeout must be bigger or equal to sizein.");
  if constexpr (sizeout == sizein)
    return in;
  else {
    std::array<T, sizeout> ret;
    for (size_t i = 0; i < sizein; i++) ret[i] = in[i];
    for (size_t i = sizein; i < sizeout; i++) ret[i] = static_cast<T>(0);
    return ret;
  }
}
 
template <typename T, size_t size>
T evaluateInterval(const T &x, const std::array<T, size> &coeffs, const T &xm) {
  // Use Horner's scheme to evaluate.
  const T dx = x - xm;
  T result = coeffs.back();
  for (auto it = coeffs.rbegin() + 1; it != coeffs.rend(); it++) {
    result = dx * result + (*it);
  }
  return result;
}
 
template <typename T>
T faculty(size_t n) {
  T retVal = static_cast<T>(1);
  for (size_t i = 2; i <= n; i++) {
    retVal *= static_cast<T>(i);
  }
  return retVal;
}
 
template <typename T>
T facultyRatio(size_t counter, size_t denominator) {
  if (denominator > counter) {
    return static_cast<T>(1) / facultyRatio<T>(denominator, counter);
  } else {
    T retVal = static_cast<T>(1);
    for (size_t i = denominator + 1; i <= counter; i++) {
      retVal *= static_cast<T>(i);
    }
    return retVal;
  }
}
 
template <typename T>
T binomialCoefficient(size_t n, size_t k) {
  if (k > n) {
    return static_cast<T>(0);
  }
 
  const size_t smaller = std::min(k, n - k);
  const size_t larger = std::max(k, n - k);
 
  return facultyRatio<T>(n, larger) / faculty<T>(smaller);
}
 
}  // end namespace bspline::internal
 
#endif  // BSPLINE_DOXYGEN_IGNORE
#endif  // BSPLINE_MISC_H