template<typename T> struct quat_base {
  static constexpr int Size = 4;
  static constexpr int type_length = Size;

  using base_type = T;
  using uint_type = quat_base<as_uint_type<T>>;

  T w, x, y, z;

  quat_base() = default;

  explicit quat_base(uint value)
  {
    for (int i = 0; i < Size; i++) {
      (*this)[i] = static_cast<T>(value);
    }
  }

  explicit quat_base(int value)
  {
    for (int i = 0; i < Size; i++) {
      (*this)[i] = static_cast<T>(value);
    }
  }

  explicit quat_base(float value)
  {
    for (int i = 0; i < Size; i++) {
      (*this)[i] = static_cast<T>(value);
    }
  }

  explicit quat_base(double value)
  {
    for (int i = 0; i < Size; i++) {
      (*this)[i] = static_cast<T>(value);
    }
  }

  quat_base(T _w, T _x, T _y, T _z)
  {
    (*this)[0] = _w;
    (*this)[1] = _x;
    (*this)[2] = _y;
    (*this)[3] = _z;
  }

  /** Mixed scalar-vector constructors. */

  template<typename U>
  quat_base(T re, vec_base<U, 3> im)
      : quat_base(
            static_cast<T>(re), static_cast<T>(im.x), static_cast<T>(im.y), static_cast<T>(im.z))
  {
  }

  /** Conversion from pointers (from C-style vectors). */

  quat_base(const T *ptr)
  {
    for (int i = 0; i < Size; i++) {
      (*this)[i] = ptr[i];
    }
  }

  template<typename U, BLI_ENABLE_IF((std::is_convertible_v<U, T>))>
  explicit quat_base(const U *ptr)
  {
    for (int i = 0; i < Size; i++) {
      (*this)[i] = ptr[i];
    }
  }

  quat_base(const T (*ptr)[Size]) : quat_base(static_cast<const T *>(ptr[0]))
  {
  }

  /** Conversion from other vector types. */

  template<typename U> explicit quat_base(const vec_base<U, Size> &vec)
  {
    for (int i = 0; i < Size; i++) {
      (*this)[i] = static_cast<T>(vec[i]);
    }
  }

  /** C-style pointer dereference. */

  operator const T *() const
  {
    return reinterpret_cast<const T *>(this);
  }

  operator T *()
  {
    return reinterpret_cast<T *>(this);
  }

  /** Array access. */

  const T &operator[](int index) const
  {
    BLI_assert(index >= 0);
    BLI_assert(index < Size);
    return reinterpret_cast<const T *>(this)[index];
  }

  T &operator[](int index)
  {
    BLI_assert(index >= 0);
    BLI_assert(index < Size);
    return reinterpret_cast<T *>(this)[index];
  }

  /** Internal Operators Macro. */

#define BLI_VEC_OP_IMPL(_result, _i, _op) \
  quat_base _result; \
  for (int _i = 0; _i < Size; _i++) { \
    _op; \
  } \
  return _result;

#define BLI_VEC_OP_IMPL_SELF(_i, _op) \
  for (int _i = 0; _i < Size; _i++) { \
    _op; \
  } \
  return *this;

  /** Arithmetic operators. */

  friend quat_base operator+(const quat_base &a, const quat_base &b)
  {
    BLI_VEC_OP_IMPL(ret, i, ret[i] = a[i] + b[i]);
  }

  friend quat_base operator+(const quat_base &a, const T &b)
  {
    BLI_VEC_OP_IMPL(ret, i, ret[i] = a[i] + b);
  }

  friend quat_base operator+(const T &a, const quat_base &b)
  {
    return b + a;
  }

  quat_base &operator+=(const quat_base &b)
  {
    BLI_VEC_OP_IMPL_SELF(i, (*this)[i] += b[i]);
  }

  quat_base &operator+=(const T &b)
  {
    BLI_VEC_OP_IMPL_SELF(i, (*this)[i] += b);
  }

  friend quat_base operator-(const quat_base &a)
  {
    BLI_VEC_OP_IMPL(ret, i, ret[i] = -a[i]);
  }

  friend quat_base operator-(const quat_base &a, const quat_base &b)
  {
    BLI_VEC_OP_IMPL(ret, i, ret[i] = a[i] - b[i]);
  }

  friend quat_base operator-(const quat_base &a, const T &b)
  {
    BLI_VEC_OP_IMPL(ret, i, ret[i] = a[i] - b);
  }

  friend quat_base operator-(const T &a, const quat_base &b)
  {
    BLI_VEC_OP_IMPL(ret, i, ret[i] = a - b[i]);
  }

  quat_base &operator-=(const quat_base &b)
  {
    BLI_VEC_OP_IMPL_SELF(i, (*this)[i] -= b[i]);
  }

  quat_base &operator-=(const T &b)
  {
    BLI_VEC_OP_IMPL_SELF(i, (*this)[i] -= b);
  }

  friend quat_base operator*(const quat_base &a, const quat_base &b)
  {
    /* Quaternion multiplication is not component-wise */
    return quat_base{a[0] * b[0] - a[1] * b[1] - a[2] * b[2] - a[3] * b[3],
                     a[0] * b[1] + a[1] * b[0] + a[2] * b[3] - a[3] * b[2],
                     a[0] * b[2] + a[2] * b[0] + a[3] * b[1] - a[1] * b[3],
                     a[0] * b[3] + a[3] * b[0] + a[1] * b[2] - a[2] * b[1]};
  }

  template<typename FactorT> friend quat_base operator*(const quat_base &a, FactorT b)
  {
    /* Multiplying by scalar is equivalent to multiplying by (b, 0, 0, 0) */
    BLI_VEC_OP_IMPL(ret, i, ret[i] = a[i] * b);
  }

  friend quat_base operator*(T a, const quat_base &b)
  {
    return b * a;
  }

  friend vec_base<T, 3> operator*(const quat_base &q, const vec_base<T, 3> &r)
  {
    quat_base t{-q[1] * r[0] - q[2] * r[1] - q[3] * r[2],
                q[0] * r[0] + q[2] * r[2] - q[3] * r[1],
                q[0] * r[1] + q[3] * r[0] - q[1] * r[2],
                q[0] * r[2] + q[1] * r[1] - q[2] * r[0]};
    return float3{t[0] * -q[1] + t[1] * q[0] - t[2] * q[3] + t[3] * q[2],
                  t[0] * -q[2] + t[2] * q[0] - t[3] * q[1] + t[1] * q[3],
                  t[0] * -q[3] + t[3] * q[0] - t[1] * q[2] + t[2] * q[1]};
  }

  quat_base &operator*=(const quat_base &b)
  {
    BLI_VEC_OP_IMPL_SELF(i, (*this)[i] *= b[i]);
  }

  quat_base &operator*=(T b)
  {
    BLI_VEC_OP_IMPL_SELF(i, (*this)[i] *= b);
  }

  friend quat_base operator/(const quat_base &a, T b)
  {
    BLI_assert(b != T(0));
    BLI_VEC_OP_IMPL(ret, i, ret[i] = a[i] / b);
  }

  friend quat_base operator/(T a, const quat_base &b)
  {
    for (int i = 0; i < Size; i++) {
      BLI_assert(b[i] != T(0));
    }
    BLI_VEC_OP_IMPL(ret, i, ret[i] = a / b[i]);
  }

  quat_base &operator/=(T b)
  {
    BLI_assert(b != T(0));
    BLI_VEC_OP_IMPL_SELF(i, (*this)[i] /= b);
  }

#undef BLI_VEC_OP_IMPL
#undef BLI_VEC_OP_IMPL_SELF

  /** Compare. */

  friend bool operator==(const quat_base &a, const quat_base &b)
  {
    for (int i = 0; i < Size; i++) {
      if (a[i] != b[i]) {
        return false;
      }
    }
    return true;
  }

  friend bool operator!=(const quat_base &a, const quat_base &b)
  {
    return !(a == b);
  }

  /** Misc. */

  uint64_t hash() const
  {
    return math::vector_hash(*this);
  }

  friend std::ostream &operator<<(std::ostream &stream, const quat_base &v)
  {
    stream << "(";
    for (int i = 0; i < Size; i++) {
      stream << v[i];
      if (i != Size - 1) {
        stream << ", ";
      }
    }
    stream << ")";
    return stream;
  }

  static quat_base unit()
  {
    return quat_base(1.0f, 0.0f, 0.0f, 0.0f);
  }

  inline T real() const
  {
    return w;
  }

  inline vec_base<T, 3> imag() const
  {
    return vec_base<T, 3>(x, y, z);
  }
};

namespace math
{

/** Quaternion functions */
template<typename T, BLI_ENABLE_IF((is_math_float_type<T>))>
inline T length_squared(const quat_base<T> &q)
{
  return q[0] * q[0] + q[1] * q[1] + q[2] * q[2] + q[3] * q[3];
}

template<typename T, BLI_ENABLE_IF((is_math_float_type<T>))>
inline T length(const quat_base<T> &q)
{
  return std::sqrt(length_squared(q));
}

template<typename T, BLI_ENABLE_IF((is_math_float_type<T>))>
quat_base<T> conjugate(const quat_base<T> &q)
{
  return quat_base<T>(q[0], -q[1], -q[2], -q[3]);
};

template<typename T, BLI_ENABLE_IF((is_math_float_type<T>))>
quat_base<T> inverse(const quat_base<T> &q)
{
  const T len_sq = length_squared(q);
  if (len_sq == 0.0f) {
    return q;
  }
  return conjugate(q) / len_sq;
};

template<typename T, BLI_ENABLE_IF((is_math_float_type<T>))>
quat_base<T> axis_angle_normalized_to_quat(const vec_base<T, 3> &axis, const float angle)
{
  const float phi = 0.5f * angle;
  const float si = sinf(phi);
  const float co = cosf(phi);
  BLI_ASSERT_UNIT_V3(axis);
  return Quat(co, axis * si);
}

template<typename T, BLI_ENABLE_IF((is_math_float_type<T>))>
quat_base<T> expmap_to_quat(const vec_base<T, 3> &expmap)
{
  float angle;
  const float3 axis = math::normalize_and_get_length(expmap, angle);

  /* Obtain axis/angle representation. */
  if (LIKELY(angle != 0.0f)) {
    return axis_angle_normalized_to_quat(axis, angle_wrap_rad(angle));
  }
  else {
    return Quat::unit();
  }
}


template<typename T, BLI_ENABLE_IF((is_math_float_type<T>))>
inline quat_base<T> normalize_and_get_length(const quat_base<T> &v, T &out_length)
{
  out_length = length_squared(v);
  /* A larger value causes normalize errors in a scaled down models with camera extreme close. */
  constexpr T threshold = std::is_same_v<T, double> ? 1.0e-70 : 1.0e-35f;
  if (out_length > threshold) {
    out_length = sqrt(out_length);
    return v / out_length;
  }
  /* Either the vector is small or one of it's values contained `nan`. */
  out_length = 0.0;
  return quat_base<T>(0.0);
}

template<typename T, BLI_ENABLE_IF((is_math_float_type<T>))>
inline quat_base<T> normalize(const quat_base<T> &v)
{
  T len;
  return normalize_and_get_length(v, len);
}

}  // namespace math

using Quat = quat_base<float>;