math.h

Go to the documentation of this file.

#ifndef COTILA_SCALAR_MATH_H_
#define COTILA_SCALAR_MATH_H_

#include <cotila/detail/type_traits.h>
#include <cotila/detail/assert.h>
#include <type_traits>

namespace cotila {

constexpr double sqrt(double x) {
  if (x < 0)
    throw "sqrt argument must be positive";
  double prev = 0;
  double est = (1 + x) / 2;
  while (prev != est) {
    prev = est;
    est = (est + x / est) / 2;
  }
  return est;
}

constexpr float sqrt(float x) { return float(sqrt(double(x))); }

template <typename T> constexpr detail::remove_complex_t<T> abs(T x) {
  COTILA_DETAIL_ASSERT_ARITHMETIC(T);
  if constexpr (detail::is_complex_v<T>)
    return sqrt(x.real() * x.real() + x.imag() * x.imag());
  else
    return x > 0 ? x : -x;
}

constexpr double exponentiate(double x, int n) {
  if (n == 0)
    return 1;
  if (n < 0) {
    x = 1. / x;
    n = -n;
  }
  double y = 1.;
  while (n > 1) {
    if (n % 2 == 0) {
      n /= 2;
    } else {
      y *= x;
      n = (n - 1) / 2;
    }
    x *= x;
  }
  return x * y;
}

constexpr double nthroot(double x, int n) {
  if (x < 0)
    throw "nth root argument must be positive";
  double prev = -1;
  double est = 1;
  while (prev != est) {
    prev = est;
    double dxk = 1. / n * (x / exponentiate(prev, n - 1) - prev);
    est = prev + dxk;
  }
  return est;
}

template <typename T> constexpr T conj(T x) {
  COTILA_DETAIL_ASSERT_ARITHMETIC(T);
  if constexpr (detail::is_complex_v<T>)
    return {x.real(), -x.imag()};
  else
    return x;
}

} // namespace cotila

#endif // COTILA_SCALAR_MATH_H_