Math.h

Go to the documentation of this file.

#ifndef INLCUDE_MOLASSEMBLER_TEMPLE_CONSTEXPR_MATH_H
#define INLCUDE_MOLASSEMBLER_TEMPLE_CONSTEXPR_MATH_H

#include "Molassembler/Temple/Preprocessor.h"

#include <cmath>
#include <limits>
#include <type_traits>

namespace Scine {
namespace Molassembler {
namespace Temple {
namespace Math {
namespace Traits {

template<typename T, typename U>
using enableIfFloatingWithReturn = std::enable_if_t<
  std::is_floating_point<T>::value,
  U
>;

template<typename T, typename U>
using enableIfIntegralWithReturn = std::enable_if_t<
  std::is_integral<T>::value,
  U
>;

template<typename T, typename U>
using enableIfArithmeticWithReturn = std::enable_if_t<
  std::is_arithmetic<T>::value,
  U
>;

template<typename FloatingPoint>
PURITY_STRONG constexpr inline std::enable_if_t<
  std::is_floating_point<FloatingPoint>::value,
  bool
> isnan(const FloatingPoint x) {
  // see notes at http://en.cppreference.com/w/cpp/numeric/math/isnan
  return x != x;
}

} // namespace Traits

/* Logic */
template<typename ... Bools>
constexpr bool XOR(Bools ... bools);

/* Some very basic math functions for arithmetic types */
template<typename T>
inline constexpr Traits::enableIfArithmeticWithReturn<T, T> abs(T x) noexcept;

template<typename T>
constexpr Traits::enableIfArithmeticWithReturn<T, T> max(T a, T b) noexcept;

template<typename T>
constexpr Traits::enableIfArithmeticWithReturn<T, T> min(T a, T b) noexcept;

/* Floating-point math functions */

template<typename T>
constexpr Traits::enableIfFloatingWithReturn<T, T> toRadians(T inDegrees) noexcept;

template<typename T>
constexpr Traits::enableIfFloatingWithReturn<T, T> toDegrees(T inRadians) noexcept;

template<typename T, typename U>
constexpr Traits::enableIfFloatingWithReturn<T, T> fmod(T value, U divider) noexcept;

template<typename T>
constexpr Traits::enableIfFloatingWithReturn<T, long> ceil(T value) noexcept;

template<typename T>
constexpr Traits::enableIfFloatingWithReturn<T, int> floor(T value) noexcept;

template<typename T>
constexpr Traits::enableIfArithmeticWithReturn<T, T> pow(T base, unsigned exponent) noexcept;

template<typename T>
constexpr Traits::enableIfArithmeticWithReturn<T, double> pow(T base, int exponent) noexcept;

template<typename T>
constexpr Traits::enableIfFloatingWithReturn<T, T> sqrt(T x);

template<typename T>
constexpr Traits::enableIfIntegralWithReturn<T, T> factorial(T x);

template<typename T>
constexpr Traits::enableIfFloatingWithReturn<T, T> ln(T x);

template<typename T>
constexpr Traits::enableIfFloatingWithReturn<T, T> log10(T x);

template<typename T>
constexpr Traits::enableIfFloatingWithReturn<T, T> log(T x, T base);

// Inverse trigonometry
template<typename T>
constexpr Traits::enableIfFloatingWithReturn<T, T> asin(T x);

template<typename T>
constexpr Traits::enableIfFloatingWithReturn<T, T> acos(T x);

template<typename T>
constexpr Traits::enableIfFloatingWithReturn<T, T> atan(T x);


/* Implementations begin here ------------------------------------------------*/

namespace Detail { // Implementation helpers

// Specialization of TPPSum for empty parameter pack
constexpr unsigned TPPSum() {
  return 0;
}

// Template parameter pack sum (needed for XOR function)
template<typename T1, typename... T>
constexpr unsigned TPPSum(T1 a, T ... pack) {
  return a + TPPSum(pack ...);
}

/* Based on series expansion of ln x:
 *
 *   y(x) = (x - 1) / (x + 1)
 *   ln x = 2 [ y + y^3/3 + y^5/5 + ...]
 *
 *   for Re x ≥ 0 and x ≠ 0
 */
template<typename T>
PURITY_STRONG constexpr T lnSeries(const T x) {
  if(x <= 0) {
    throw "Ln domain error: x <= 0";
  }

  const T epsilon = std::numeric_limits<T>::epsilon();

  // Since every term, y is multiplied by y², we keep a counting square
  T countingPower { // at n = 1
    (x - 1) / (x + 1)
  };

  // This is just y², which is what we multiply which each term
  const T iterationMultiplier = countingPower * countingPower;

  // Initial value for n = 1 (terms up until y)
  T value {
    2 * countingPower
  };

  // Term n = 0
  T previous {0};

  for(unsigned n = 3; Temple::Math::abs(previous - value) > epsilon; n += 2) {
    // previous iteration
    previous = value;

    // Update the power: y^(n-2) * y^2 becomes y^n
    countingPower *= iterationMultiplier;

    // Add the full term to the running value
    value += 2 * countingPower / n;
  }

  return value;
}

/* Implements an approximation to asin over 0 < x < 1 with accuracy > 2e-8
 * from Abramowitz, M., Stegun, I.: Handbook of Mathematical Functions, p. 64,
 * 1964, from http://people.math.sfu.ca/~cbm/aands/abramowitz_and_stegun.pdf
 */
template<typename T>
PURITY_STRONG constexpr T asinApprox(const T x) {
  if(!(0.0 <= x && x <= 1.0)) {
    throw "Asin approximation domain error: only applicable for 0 < x < 1!";
  }

  const T x2 = x * x;
  const T x4 = x2 * x2;

  return (
    M_PI / 2
    - Math::sqrt(1 - x) * (
      1.5707963050
      - 0.2145988016 * x
      + 0.0889789874 * x2
      - 0.0501743046 * x2 * x
      + 0.0308918810 * x4
      - 0.0170881256 * x4 * x
      + 0.0066700901 * x4 * x2
      - 0.0012624911 * x4 * x2 * x
    )
  );
}

} // namespace Detail

template<typename ... Bools>
constexpr bool XOR(Bools ... bools) {
  return Detail::TPPSum(bools ...) == 1;
}

template<typename T>
PURITY_STRONG inline constexpr Traits::enableIfArithmeticWithReturn<T, T> abs(const T x) noexcept {
  return (x >= 0) ? x : -x;
}

template<typename T>
PURITY_STRONG constexpr Traits::enableIfArithmeticWithReturn<T, T> max(const T a, const T b) noexcept {
  return (a > b) ? a : b;
}

template<typename T>
PURITY_STRONG constexpr Traits::enableIfArithmeticWithReturn<T, T> min(const T a, const T b) noexcept {
  return (a < b) ? a : b;
}

template<typename T>
PURITY_STRONG constexpr Traits::enableIfFloatingWithReturn<T, T> toRadians(const T inDegrees) noexcept {
  return M_PI * inDegrees / 180;
}

template<typename T>
PURITY_STRONG constexpr Traits::enableIfFloatingWithReturn<T, T> toDegrees(const T inRadians) noexcept {
  return 180 * inRadians / M_PI;
}

template<typename T, typename U>
PURITY_STRONG constexpr Traits::enableIfFloatingWithReturn<T, T>  fmod(const T value, const U divider) noexcept {
 T remainder = value;
 while (divider < remainder) {
   remainder -= divider;
 }
 return remainder;
}

template<typename T>
PURITY_STRONG constexpr Traits::enableIfFloatingWithReturn<T, long> ceil(const T value) noexcept {
  // Truncate to an integer
  const auto truncated = static_cast<long>(value);

  // we first check if the given floating point number is actually an integer 
  // this avoids rounding mistakes, occuring with some compilers
  const double eps = 1e-12;
  if (fmod(value, 1) < eps && abs(truncated - value) < eps) { 
    return truncated;
  }

  if(truncated < value) {
    return truncated + 1;
  }

  return truncated;
}

template<typename T>
PURITY_STRONG constexpr Traits::enableIfFloatingWithReturn<T, int> floor(const T value) noexcept {
  // Truncate to an int
  const auto truncated = static_cast<int>(value);

  if(truncated > value) {
    return truncated - 1;
  }

  return truncated;
}

// Really weak first implementation
template<typename T>
PURITY_STRONG constexpr Traits::enableIfArithmeticWithReturn<T, T> pow(const T base, const unsigned exponent) noexcept {
  if(exponent == 0) {
    return 1;
  }

  double value = base;

  for(unsigned n = 1; n < exponent; n++) {
    value *= base;
  }

  return value;
}

template<typename T>
PURITY_STRONG constexpr T recPow(const T base, const unsigned exponent) noexcept {
  if(exponent == 0) {
    return 1;
  }

  if(exponent == 1) {
    return base;
  }

  if(exponent % 2 == 0) {
    auto halfProblem = recPow(base, exponent / 2);
    return halfProblem * halfProblem;
  }

  return base * recPow(base, exponent - 1);
}

template<typename T>
PURITY_STRONG constexpr Traits::enableIfArithmeticWithReturn<T, double> pow(const T base, const int exponent) noexcept {
  if(exponent < 0) {
    return 1.0 / pow(base, static_cast<unsigned>(Temple::Math::abs(exponent)));
  }

  if(exponent == 0) {
    return 1;
  }

  return pow(base, static_cast<unsigned>(exponent));
}

/* Implements Newton's iteration to compute the square root of a positive number
 */
template<typename T>
PURITY_STRONG constexpr Traits::enableIfFloatingWithReturn<T, T> sqrt(const T x) {
  if(x < 0) {
    throw "Square-root domain error: Only real if x >= 0!";
  }

  const T epsilon = std::numeric_limits<T>::epsilon();
  T value = 1;
  T previous = 2;

  while(Temple::Math::abs(previous - value) > epsilon) {
    // store the previous value
    previous = value;

    // compute next iteration
    value = 0.5 * (value + x / value);
  }

  return value;
}

template<typename T>
PURITY_STRONG constexpr Traits::enableIfIntegralWithReturn<T, T> factorial(const T x) {
  if(x < 0) {
    throw "Factorial domain error!";
  }

  if(x == 0) {
    return 1;
  }

  return x * factorial(x - 1);
}

template<typename T>
PURITY_STRONG constexpr Traits::enableIfFloatingWithReturn<T, T> ln(const T x) {
  unsigned decimalReduction = 0;
  T calcX = x;

  while(abs(calcX) > 10) {
    calcX /= 10;
    decimalReduction += 1;
  }

  // Ensure last division leads to value closer to 1
  if(Math::abs(calcX / 10 - 1) < Math::abs(calcX - 1)) {
    calcX /= 10;
    decimalReduction += 1;
  }

  return Detail::lnSeries(calcX) + decimalReduction * M_LN10;
}

template<typename T>
PURITY_STRONG constexpr Traits::enableIfFloatingWithReturn<T, T> log10(const T x) {
  if(x <= 0) {
    throw "Log10 domain error!";
  }

  /* ln(z) = ln(10) * log10(z)
   * -> log10(z) = ln(z) / ln(10)
   */
  return ln(x) / M_LN10;
}

template<typename T>
PURITY_STRONG constexpr Traits::enableIfFloatingWithReturn<T, T> log(const T x, const T base) {
  if(x <= 0) {
    throw "Log domain error!";
  }

  /* ln(z) = ln(b) * log_b(z)
   * -> log_b(z) = ln(z) / ln(b)
   */
  return ln(x) / ln(base);
}

/* Implements the infinite series where the derivative is expanded as a binomial
 * series and every term is integrated. Deviates most strongly from std::asin at
 * values very close to the borders. Perhaps it is best to use the approximate
 * form there?
 */
template<typename T>
PURITY_STRONG constexpr Traits::enableIfFloatingWithReturn<T, T> asin(const T x) {
  if(!(-1 <= x && x <= 1)) {
    throw "Inverse sine domain error: only real if -1 < x < 1!";
  }

  if(Temple::Math::abs(x) > 0.90) {
    return (x >= 0) ? Detail::asinApprox(x) : -Detail::asinApprox(-x);
  }

  const T epsilon = std::numeric_limits<T>::epsilon();

  T value = x;
  T upper_factorial = 1;
  T lower_factorial = 1;
  T term = 1;

  for(unsigned n = 1; Temple::Math::abs(term) > epsilon; ++n) {
    upper_factorial *= 2 * (n - 1) + 1;
    lower_factorial *= 2 * n;

    term = (
      upper_factorial / lower_factorial
    ) * pow(x, 2 * n + 1) / (2 * n + 1);

    if(Traits::isnan(term)) {
      break;
    }

    value += term;
  }

  return value;
}

template<typename T>
PURITY_STRONG constexpr Traits::enableIfFloatingWithReturn<T, T> acos(const T x) {
  if(!(-1 <= x && x <= 1)) {
    throw "Inverse cosine domain error: only real if -1 <= x <= 1!";
  }

  return M_PI / 2 - asin(x);
}

template<typename T>
PURITY_STRONG constexpr Traits::enableIfFloatingWithReturn<T, T> atan(const T x) {
  if(!(-M_PI / 2 < x && x < M_PI / 2)) {
    throw "Inverse cosine domain error: only real if -1 < x < 1!";
  }

  return asin(
    x / sqrt(x * x + 1)
  );
}

} // namespace Math
} // namespace Temple
} // namespace Molassembler
} // namespace Scine

#endif