Quaternion (tools/quaternion.hpp)

View this file on GitHub
View document part on GitHub
Last update: 2024-08-31 13:46:12+09:00
Include: #include "tools/quaternion.hpp"

It is a quaternion.

License

Author

anqooqie

Constructor

(1) quaternion<T> q();
(2) quaternion<T> q(T x, T y, T z, T w);

(1)
- It creates a quaternion $0$.
(2)
- It creates a quaternion $w + xi + yj + zk$.

Constraints

std::is_floating_point_v<T> is true.

Time Complexity

$O(1)$

x

T q.x;

It is the X component of $q$.

y

T q.y;

It is the Y component of $q$.

z

T q.z;

It is the Z component of $q$.

w

T q.w;

It is the W component of $q$.

real

(1) T q.real();
(2) void q.real(T val);

(1)
- It returns the W component of $q$.
(2)
- It updates the W component of $q$ to val.

Constraints

None

Time Complexity

$O(1)$

imag

(1) tools::vector3<T> q.imag();
(2) void q.imag(tools::vector3<T> val);

(1)
- It returns the X component, the Y component and the Z component of $q$.
(2)
- It updates the X component of $q$ to val.x, the Y component of $q$ to val.y and the Z component of $q$ to val.z.

Constraints

None

Time Complexity

$O(1)$

abs

T q.abs();

It returns $|q| = \sqrt{w^2 + x^2 + y^2 + z^2}$.

Constraints

None

Time Complexity

$O(1)$

norm

T q.norm();

It returns $q q^\ast = w^2 + x^2 + y^2 + z^2$.

Constraints

None

Time Complexity

$O(1)$

conj

quaternion<T> q.conj();

It returns $q^\ast = w - xi - yj - zk$.

Constraints

None

Time Complexity

$O(1)$

inv

quaternion<T> q.inv();

It returns $q^{-1} = \frac{q^\ast}{|q|^2}$.

Constraints

$q \neq 0$

Time Complexity

$O(1)$

angle

T q.angle();

It returns the angle of the rotation which $q$ represents in radians. In other words, it returns $2 \cos^{-1}(w)$.

Constraints

$|q| \approx 1$

Time Complexity

$O(1)$

axis

tools::vector3<T> q.axis();

It returns the axis of the rotation which $q$ represents. In other words, it returns the following unit vector where $\mathbb{v} = (x, y, z)$.

\[\begin{align*} \left\{\begin{array}{ll} \frac{\mathbb{v}}{|\mathbb{v}|} & \text{(if $|\mathbb{v}| > 0$)}\\ (1, 0, 0) & \text{(if $|\mathbb{v}| = 0$)} \end{array}\right.& \end{align*}\]

Constraints

$|q| \approx 1$

Time Complexity

$O(1)$

Unary plus operator

quaternion<T> q.operator+();

It returns the copy of $q$.

Constraints

None

Time Complexity

$O(1)$

Unary negation operator

quaternion<T> q.operator-();

It returns $-w - xi - yj - zk$.

Constraints

None

Time Complexity

$O(1)$

Addition operators

(1) quaternion<T>& q1.operator+=(quaternion<T> q2);
(2) quaternion<T> operator+(quaternion<T> q1, quaternion<T> q2);

(1)
- It updates $q_1$ to $(w_1 + w_2) + (x_1 + x_2)i + (y_1 + y_2)j + (z_1 + z_2)k$, and returns $q_1$.
(2)
- It returns $(w_1 + w_2) + (x_1 + x_2)i + (y_1 + y_2)j + (z_1 + z_2)k$.

Constraints

None

Time Complexity

$O(1)$

Subtraction operators

(1) quaternion<T>& q1.operator-=(quaternion<T> q2);
(2) quaternion<T> operator-(quaternion<T> q1, quaternion<T> q2);

(1)
- It updates $q_1$ to $(w_1 - w_2) + (x_1 - x_2)i + (y_1 - y_2)j + (z_1 - z_2)k$, and returns $q_1$.
(2)
- It returns $(w_1 - w_2) + (x_1 - x_2)i + (y_1 - y_2)j + (z_1 - z_2)k$.

Constraints

None

Time Complexity

$O(1)$

Multiplication operators

(1) quaternion<T>& q.operator*=(T c);
(2) quaternion<T> operator*(quaternion<T> q, T c);
(3) quaternion<T> operator*(T c, quaternion<T> q);
(4) quaternion<T>& q1.operator*=(quaternion<T> q2);
(5) quaternion<T> operator*(quaternion<T> q1, quaternion<T> q2);
(6) tools::vector3<T> operator*(quaternion<T> q, tools::vector3<T> p);

(1)
- It updates $q$ to $cw + cxi + cyj + czk$, and returns $q$.
(2), (3)
- It returns $cw + cxi + cyj + czk$.
(4)
- It updates $q_1$ to $(w_1 w_2 - \mathbb{v}_1 \cdot \mathbb{v}_2) + (w_1 \mathbb{v}_2 + w_2 \mathbb{v}_1 + \mathbb{v}_1 \times \mathbb{v}_2)$ where $\mathbb{v}_1 = (x_1, y_1, z_1)$ and $\mathbb{v}_2 = (x_2, y_2, z_2)$, and returns $q_1$.
(5)
- It returns $(w_1 w_2 - \mathbb{v}_1 \cdot \mathbb{v}_2) + (w_1 \mathbb{v}_2 + w_2 \mathbb{v}_1 + \mathbb{v}_1 \times \mathbb{v}_2)$ where $\mathbb{v}_1 = (x_1, y_1, z_1)$ and $\mathbb{v}_2 = (x_2, y_2, z_2)$.
(6)
- It applies the rotation which $q$ represents to $\mathbb{p}$. In other words, it returns the imaginary part of $q (\mathbb{p}_x i + \mathbb{p}_y j + \mathbb{p}_z k) q^\ast$, that is $(\cos \theta) \mathbb{p} + (1 - \cos \theta) (\hat{\mathbb{v}} \cdot \mathbb{p}) \hat{\mathbb{v}} + \sin \theta (\hat{\mathbb{v}} \times \mathbb{p})$ where $\theta$ = q.angle() and $\hat{\mathbb{v}}$ = q.axis().

Constraints

None

Time Complexity

$O(1)$

Division operators

(1) quaternion<T>& q.operator/=(T c);
(2) quaternion<T> operator/(quaternion<T> q, T c);
(3) quaternion<T>& q1.operator/=(quaternion<T> q2);
(4) quaternion<T> operator/(quaternion<T> q1, quaternion<T> q2);

(1)
- It updates $q$ to $\frac{w}{c} + \frac{x}{c}i + \frac{y}{c}j + \frac{z}{c}k$, and returns $q$.
(2)
- It returns $\frac{w}{c} + \frac{x}{c}i + \frac{y}{c}j + \frac{z}{c}k$.
(3)
- It updates $q_1$ to $q_1 q_2^{-1}$, and returns $q_1$.
(4)
- It returns $q_1 q_2^{-1}$.

Constraints

(1), (2)
- $c \neq 0$
(3), (4)
- $q_2 \neq 0$

Time Complexity

$O(1)$

Comparison operators

(1) bool operator==(quaternion<T> q1, quaternion<T> q2);
(2) bool operator!=(quaternion<T> q1, quaternion<T> q2);

(1)
- It returns whether $w_1 = w_2 \land x_1 = x_2 \land y_1 = y_2 \land z_1 = z_2$.
(2)
- It returns whether $w_1 \neq w_2 \lor x_1 \neq x_2 \lor y_1 \neq y_2 \lor z_1 \neq z_2$.

Constraints

None

Time Complexity

$O(1)$

operator<<

std::ostream& operator<<(std::ostream& os, const quaternion<T>& self);

It is equivalent to the following code.

std::ostringstream s;
s.flags(os.flags());
s.imbue(os.getloc());
s.precision(os.precision());
s << '(' << self.x << ',' << self.y << ',' << self.z << ',' << self.w << ')';
return os << s.str();

Constraints

None

Time Complexity

$O(1)$

angle_axis

quaternion<T> quaternion<T>::angle_axis(T theta, tools::vector3<T> v);

It returns a rotation which rotates $\theta$ radians around $\mathbb{v}$. In other words, it returns $\cos\left(\frac{\theta}{2}\right) + \sin\left(\frac{\theta}{2}\right) \frac{\mathbb{v}}{|\mathbb{v}|}$.

Constraints

$\mathbb{v} \neq (0, 0, 0)$

Time Complexity

$O(1)$

from_to_rotation

quaternion<T> quaternion<T>::from_to_rotation(tools::vector3<T> u, tools::vector3<T> v);

It returns a rotation which rotates from $\mathbb{u}$ to $\mathbb{v}$.

Constraints

$\mathbb{u} \neq (0, 0, 0)$
$\mathbb{v} \neq (0, 0, 0)$

Time Complexity

$O(1)$

look_rotation

quaternion<T> quaternion<T>::look_rotation(tools::vector3<T> f, tools::vector3<T> u);

It returns a rotation with the specified forward direction $\mathbb{f}$ and upwards direction $\mathbb{u}$. Z axis will be aligned with $\mathbb{f}$, X axis aligned with $\mathbb{f} \times \mathbb{u}$, and Y axis aligned with cross product between Z and X.

Constraints

$\mathbb{u} \neq (0, 0, 0)$
$\mathbb{v} \neq (0, 0, 0)$

Time Complexity

$O(1)$

slerp

quaternion<T> quaternion<T>::slerp(quaternion<T> q0, quaternion<T> q1, T t);

It spherically interpolates between quaternions $q_0$ and $q_1$ by ratio $t$. In other words, it returns $q_0 \left(q_0^{-1} q_1\right)^t$.

Constraints

$\mathbb{q_0} \neq 0$

Time Complexity

$O(1)$

identity

quaternion<T> quaternion<T>::identity();

It returns $1$.

Constraints

None

Time Complexity

$O(1)$

random

template <typename R>
quaternion<T> quaternion<T>::random(R& engine);

It returns a random rotation.

Constraints

engine is a pseudorandom number generator.

Time Complexity

$O(1)$

dice_rotations

std::array<quaternion<T>, 24> quaternion<T>::dice_rotations();

It returns all the states which a dice can have..

Constraints

None

Time Complexity

$O(1)$

tools::exp

quaternion<T> tools::exp(quaternion<T> q);

It returns $e^q = \sum_{n = 0}^\infty \frac{q^n}{n!}$.

Note

$e^q$ is equal to the following quaternion where $\mathbb{v} = xi + yj + zk$.

\[\begin{align*} \left\{\begin{array}{ll} e^w \left(\cos|\mathbb{v}| + \sin|\mathbb{v}| \frac{\mathbb{v}}{|\mathbb{v}|}\right) & \text{(if $|\mathbb{v}| > 0$)}\\ e^w & \text{(if $|\mathbb{v}| = 0$)} \end{array}\right.& \end{align*}\]

Constraints

None

Time Complexity

$O(1)$

tools::log

quaternion<T> tools::log(quaternion<T> q);

It returns $\mathrm{Log}(q)$. $\mathrm{Log}(q)$ is defined as follows where $\mathbb{v} = xi + yj + zk$.

\[\begin{align*} \mathrm{Log}(q) &= \left\{\begin{array}{ll} \log |q| + \mathrm{atan2}(|\mathbb{v}|, w) \frac{\mathbb{v}}{|\mathbb{v}|} & \text{(if $|\mathbb{v}| > 0$)}\\ \log |q| + \mathrm{atan2}(|\mathbb{v}|, w) i & \text{(if $|\mathbb{v}| = 0$ and $w \neq 0$)} \end{array}\right.& \end{align*}\]

Note

$\mathrm{Log}(q)$ is one of the quaternions $q’$ such that $e^{q’} = q$ holds. All the quaternions $q’$ such that $e^{q’} = q$ holds can be written in the following form for any unit three-dimensional vector $\hat{\mathbb{u}}$ and any integer $k$.

\[\begin{align*} \left\{\begin{array}{ll} \log |q| + (\mathrm{atan2}(|\mathbb{v}|, w) + 2k \pi) \frac{\mathbb{v}}{|\mathbb{v}|} & \text{(if $|\mathbb{v}| > 0$)}\\ \log |q| + (\mathrm{atan2}(|\mathbb{v}|, w) + 2k \pi) \hat{\mathbb{u}} & \text{(if $|\mathbb{v}| = 0$ and $w \neq 0$)} \end{array}\right.& \end{align*}\]

Constraints

$q \neq 0$

Time Complexity

$O(1)$

tools::pow

quaternion<T> tools::pow(quaternion<T> q, T t);

It returns $q^t$. $q^t$ is defined as $e^{t \mathrm{Log}(q)}$.

Constraints

$q \neq 0$ or $t \neq 0$

Time Complexity

$O(1)$

Depends on

Verified with

Code

#ifndef TOOLS_QUATERNION_HPP
#define TOOLS_QUATERNION_HPP

#include <type_traits>
#include <cmath>
#include <algorithm>
#include <cassert>
#include <sstream>
#include <random>
#include "tools/exp.hpp"
#include "tools/log.hpp"
#include "tools/pow.hpp"
#include "tools/vector4.hpp"
#include "tools/vector3.hpp"

namespace tools {
  template <typename T>
  class quaternion;

  template <typename T>
  ::tools::quaternion<T> exp(const ::tools::quaternion<T>& q);
  template <typename T>
  ::tools::quaternion<T> log(const ::tools::quaternion<T>& q);
  template <typename T>
  ::tools::quaternion<T> pow(const ::tools::quaternion<T>& base, T exponent);

  template <typename T>
  class quaternion {
    static_assert(::std::is_floating_point_v<T>);

  private:
    static constexpr T eps = 1e-5;
    static constexpr T pi = 3.14159265358979323846264338327950288419716939937510L;

  public:
    T x;
    T y;
    T z;
    T w;

    quaternion() = default;
    quaternion(const ::tools::quaternion<T>&) = default;
    quaternion(::tools::quaternion<T>&&) = default;
    ~quaternion() = default;
    ::tools::quaternion<T>& operator=(const ::tools::quaternion<T>&) = default;
    ::tools::quaternion<T>& operator=(::tools::quaternion<T>&&) = default;

    quaternion(const T x, const T y, const T z, const T w) : x(x), y(y), z(z), w(w) {}

    T real() const {
      return this->w;
    }
    void real(const T val) {
      this->w = val;
    }
    ::tools::vector3<T> imag() const {
      return ::tools::vector3<T>(this->x, this->y, this->z);
    }
    void imag(const ::tools::vector3<T>& val) {
      this->x = val.x;
      this->y = val.y;
      this->z = val.z;
    }

    T abs() const {
      return ::tools::vector4<T>(this->x, this->y, this->z, this->w).l2_norm();
    }
    T norm() const {
      return ::tools::vector4<T>(this->x, this->y, this->z, this->w).squared_l2_norm();
    }
    ::tools::quaternion<T> conj() const {
      return ::tools::quaternion<T>(-this->x, -this->y, -this->z, this->w);
    }

    T angle() const {
      assert(::std::abs(this->norm() - 1) <= eps);
      return ::std::acos(::std::clamp<T>(this->real(), -1, 1)) * 2;
    }
    ::tools::vector3<T> axis() const {
      assert(::std::abs(this->norm() - 1) <= eps);
      return ::std::abs(this->real()) >= 1.0 ? ::tools::vector3<T>(1, 0, 0) : this->imag().normalized();
    }

    ::tools::quaternion<T> operator+() const {
      return *this;
    }
    ::tools::quaternion<T> operator-() const {
      return ::tools::quaternion<T>(-this->x, -this->y, -this->z, -this->w);
    }

    ::tools::quaternion<T>& operator+=(const ::tools::quaternion<T>& other) {
      this->x += other.x;
      this->y += other.y;
      this->z += other.z;
      this->w += other.w;
      return *this;
    }
    friend ::tools::quaternion<T> operator+(const ::tools::quaternion<T>& lhs, const ::tools::quaternion<T>& rhs) {
      return ::tools::quaternion<T>(lhs) += rhs;
    }

    ::tools::quaternion<T>& operator-=(const ::tools::quaternion<T>& other) {
      this->x -= other.x;
      this->y -= other.y;
      this->z -= other.z;
      this->w -= other.w;
      return *this;
    }
    friend ::tools::quaternion<T> operator-(const ::tools::quaternion<T>& lhs, const ::tools::quaternion<T>& rhs) {
      return ::tools::quaternion<T>(lhs) -= rhs;
    }

    ::tools::quaternion<T>& operator*=(const T c) {
      this->x *= c;
      this->y *= c;
      this->z *= c;
      this->w *= c;
      return *this;
    }
    friend ::tools::quaternion<T> operator*(const ::tools::quaternion<T>& self, const T c) {
      return ::tools::quaternion<T>(self) *= c;
    }
    friend ::tools::quaternion<T> operator*(const T c, const ::tools::quaternion<T>& self) {
      return ::tools::quaternion<T>(self) *= c;
    }

    friend ::tools::quaternion<T> operator*(const ::tools::quaternion<T>& lhs, const ::tools::quaternion<T>& rhs) {
      const auto real = lhs.real() * rhs.real() - lhs.imag().inner_product(rhs.imag());
      const auto imag = lhs.real() * rhs.imag() + rhs.real() * lhs.imag() + lhs.imag().outer_product(rhs.imag());
      return ::tools::quaternion<T>(imag.x, imag.y, imag.z, real);
    }
    ::tools::quaternion<T>& operator*=(const ::tools::quaternion<T>& other) {
      return *this = *this * other;
    }
    friend ::tools::vector3<T> operator*(const ::tools::quaternion<T>& q, const ::tools::vector3<T>& v) {
      assert(::std::abs(q.norm() - 1) <= eps);
      return (q * ::tools::quaternion<T>(v.x, v.y, v.z, 0) * q.conj()).imag();
    }

    ::tools::quaternion<T>& operator/=(const T c) {
      this->x /= c;
      this->y /= c;
      this->z /= c;
      this->w /= c;
      return *this;
    }
    friend ::tools::quaternion<T> operator/(const ::tools::quaternion<T>& self, const T c) {
      return ::tools::quaternion<T>(self) /= c;
    }

    ::tools::quaternion<T> inv() const {
      const auto norm = this->norm();
      assert(norm != 0);
      return this->conj() / norm;
    }
    friend ::tools::quaternion<T> operator/(const ::tools::quaternion<T>& lhs, const ::tools::quaternion<T>& rhs) {
      return lhs * rhs.inv();
    }
    ::tools::quaternion<T>& operator/=(const ::tools::quaternion<T>& other) {
      return *this = *this / other;
    }

    friend bool operator==(const ::tools::quaternion<T>& lhs, const ::tools::quaternion<T>& rhs) {
      return lhs.x == rhs.x && lhs.y == rhs.y && lhs.z == rhs.z && lhs.w == rhs.w;
    }
    friend bool operator!=(const ::tools::quaternion<T>& lhs, const ::tools::quaternion<T>& rhs) {
      return !(lhs == rhs);
    }

    friend ::std::ostream& operator<<(::std::ostream& os, const ::tools::quaternion<T>& self) {
      ::std::ostringstream s;
      s.flags(os.flags());
      s.imbue(os.getloc());
      s.precision(os.precision());
      s << '(' << self.x << ',' << self.y << ',' << self.z << ',' << self.w << ')';
      return os << s.str();
    }

    static ::tools::quaternion<T> angle_axis(const T angle, const ::tools::vector3<T>& axis) {
      assert(axis != ::tools::vector3<T>(0, 0, 0));
      const auto real = ::std::cos(angle / 2);
      const auto imag = ::std::sin(angle / 2) * axis.normalized();
      return ::tools::quaternion<T>(imag.x, imag.y, imag.z, real);
    }

    static ::tools::quaternion<T> from_to_rotation(::tools::vector3<T> from_direction, ::tools::vector3<T> to_direction) {
      assert(from_direction != ::tools::vector3<T>(0, 0, 0));
      assert(to_direction != ::tools::vector3<T>(0, 0, 0));
      from_direction = from_direction.normalized();
      to_direction = to_direction.normalized();
      if (from_direction.inner_product(to_direction) <= -1 + eps) {
        return ::tools::quaternion<T>::angle_axis(pi, from_direction.orthonormal_basis()[1]);
      }
      const auto bisect = (from_direction + to_direction).normalized();
      const auto real = from_direction.inner_product(bisect);
      const auto imag = from_direction.outer_product(bisect);
      return ::tools::quaternion<T>(imag.x, imag.y, imag.z, real);
    }

    static ::tools::quaternion<T> look_rotation(::tools::vector3<T> forward, ::tools::vector3<T> upwards = ::tools::vector3<T>(0, 1, 0)) {
      assert(forward != ::tools::vector3<T>(0, 0, 0));
      assert(upwards != ::tools::vector3<T>(0, 0, 0));

      auto side = upwards.outer_product(forward);
      if (side == ::tools::vector3<T>(0, 0, 0)) {
        upwards = ::tools::vector3<T>(0, 1, 0);
        side = upwards.outer_product(forward);
      }
      if (side == ::tools::vector3<T>(0, 0, 0)) {
        side = forward.orthonormal_basis()[1];
      }
      upwards = forward.outer_product(side);

      forward = forward.normalized();
      upwards = upwards.normalized();

      const auto q1 = ::tools::quaternion<T>::from_to_rotation(::tools::vector3<T>(0, 0, 1), forward);
      const auto theta = ::std::atan2((q1 * ::tools::vector3<T>(0, 1, 0)).outer_product(upwards).inner_product(forward), (q1 * ::tools::vector3<T>(0, 1, 0)).inner_product(upwards));
      const auto q2 = ::tools::quaternion<T>::angle_axis(theta, forward);
      return q2 * q1;
    }

    static ::tools::quaternion<T> slerp(const ::tools::quaternion<T>& q0, const ::tools::quaternion<T>& q1, const T t) {
      return q0 * ::tools::pow(q0.inv() * q1, t);
    }

    static ::tools::quaternion<T> identity() {
      return ::tools::quaternion<T>(0, 0, 0, 1);
    }

    template <typename R>
    static ::tools::quaternion<T> random(R& engine) {
      static ::std::uniform_real_distribution<T> dist(0, 1);
      const auto u1 = dist(engine);
      const auto u2 = dist(engine);
      const auto u3 = dist(engine);
      return ::tools::quaternion<T>(
        ::std::sqrt(1 - u1) * ::std::sin(2 * pi * u2),
        ::std::sqrt(1 - u1) * ::std::cos(2 * pi * u2),
        ::std::sqrt(u1) * ::std::sin(2 * pi * u3),
        ::std::sqrt(u1) * ::std::cos(2 * pi * u3)
      );
    }

    static ::std::array<::tools::quaternion<T>, 24> dice_rotations() {
      return ::std::array<::tools::quaternion<T>, 24>({
        ::tools::quaternion<double>::identity(),
        ::tools::quaternion<double>::angle_axis(-0.5 * pi, ::tools::vector3<double>(1, 0, 0)),
        ::tools::quaternion<double>::angle_axis(0.5 * pi, ::tools::vector3<double>(1, 0, 0)),
        ::tools::quaternion<double>::angle_axis(pi, ::tools::vector3<double>(1, 0, 0)),
        ::tools::quaternion<double>::angle_axis(-0.5 * pi, ::tools::vector3<double>(0, 1, 0)),
        ::tools::quaternion<double>::angle_axis(0.5 * pi, ::tools::vector3<double>(0, 1, 0)),
        ::tools::quaternion<double>::angle_axis(pi, ::tools::vector3<double>(0, 1, 0)),
        ::tools::quaternion<double>::angle_axis(-0.5 * pi, ::tools::vector3<double>(0, 0, 1)),
        ::tools::quaternion<double>::angle_axis(0.5 * pi, ::tools::vector3<double>(0, 0, 1)),
        ::tools::quaternion<double>::angle_axis(pi, ::tools::vector3<double>(0, 0, 1)),
        ::tools::quaternion<double>::angle_axis(-2.0 / 3.0 * pi, ::tools::vector3<double>(-1, -1, 1)),
        ::tools::quaternion<double>::angle_axis(2.0 / 3.0 * pi, ::tools::vector3<double>(-1, -1, 1)),
        ::tools::quaternion<double>::angle_axis(-2.0 / 3.0 * pi, ::tools::vector3<double>(-1, 1, 1)),
        ::tools::quaternion<double>::angle_axis(2.0 / 3.0 * pi, ::tools::vector3<double>(-1, 1, 1)),
        ::tools::quaternion<double>::angle_axis(-2.0 / 3.0 * pi, ::tools::vector3<double>(1, -1, 1)),
        ::tools::quaternion<double>::angle_axis(2.0 / 3.0 * pi, ::tools::vector3<double>(1, -1, 1)),
        ::tools::quaternion<double>::angle_axis(-2.0 / 3.0 * pi, ::tools::vector3<double>(1, 1, 1)),
        ::tools::quaternion<double>::angle_axis(2.0 / 3.0 * pi, ::tools::vector3<double>(1, 1, 1)),
        ::tools::quaternion<double>::angle_axis(pi, ::tools::vector3<double>(0, -1, 1)),
        ::tools::quaternion<double>::angle_axis(pi, ::tools::vector3<double>(0, 1, 1)),
        ::tools::quaternion<double>::angle_axis(pi, ::tools::vector3<double>(-1, 0, 1)),
        ::tools::quaternion<double>::angle_axis(pi, ::tools::vector3<double>(1, 0, 1)),
        ::tools::quaternion<double>::angle_axis(pi, ::tools::vector3<double>(-1, 1, 0)),
        ::tools::quaternion<double>::angle_axis(pi, ::tools::vector3<double>(1, 1, 0))
      });
    }
  };

  template <typename T>
  ::tools::quaternion<T> exp(const ::tools::quaternion<T>& q) {
    const auto inorm = q.imag().l2_norm();
    if (inorm == 0) {
      return ::tools::quaternion<T>(0, 0, 0, ::std::exp(q.real()));
    }

    const auto rexp = ::std::exp(q.real());
    const auto real = rexp * ::std::cos(inorm);
    const auto imag = rexp * ::std::sin(inorm) / inorm * q.imag();
    return ::tools::quaternion<T>(imag.x, imag.y, imag.z, real);
  }

  template <typename T>
  ::tools::quaternion<T> log(const ::tools::quaternion<T>& q) {
    assert(q != ::tools::quaternion<T>(0, 0, 0, 0));
    const auto inorm = q.imag().l2_norm();
    const auto uimag = inorm == 0 ? ::tools::vector3<T>(1, 0, 0) : q.imag() / inorm;
    const auto real = ::std::log(q.abs());
    const auto imag = ::std::atan2(inorm, q.real()) * uimag;
    return ::tools::quaternion<T>(imag.x, imag.y, imag.z, real);
  }

  template <typename T>
  ::tools::quaternion<T> pow(const ::tools::quaternion<T>& base, const T exponent) {
    const auto norm = base.norm();
    if (norm == 0) {
      assert(exponent != 0);
      return ::tools::quaternion<T>(0, 0, 0, 0);
    }

    return ::tools::exp(exponent * ::tools::log(base));
  }
}

#endif

#line 1 "tools/quaternion.hpp"



#include <type_traits>
#include <cmath>
#include <algorithm>
#include <cassert>
#include <sstream>
#include <random>
#line 1 "tools/exp.hpp"



#line 5 "tools/exp.hpp"

namespace tools {
  template <typename T>
  auto exp(const T x) {
    return ::std::exp(x);
  }
}


#line 1 "tools/log.hpp"



#line 5 "tools/log.hpp"

namespace tools {
  template <typename T>
  auto log(const T x) {
    return ::std::log(x);
  }
}


#line 1 "tools/pow.hpp"



#line 1 "tools/monoid.hpp"



#line 6 "tools/monoid.hpp"
#include <limits>
#line 1 "tools/gcd.hpp"



#line 5 "tools/gcd.hpp"
#include <numeric>

namespace tools {
  template <typename M, typename N>
  constexpr ::std::common_type_t<M, N> gcd(const M m, const N n) {
    return ::std::gcd(m, n);
  }
}


#line 9 "tools/monoid.hpp"

namespace tools {
  namespace monoid {
    template <typename M, M ...dummy>
    struct max;

    template <typename M>
    struct max<M> {
      static_assert(::std::is_arithmetic_v<M>, "M must be a built-in arithmetic type.");

      using T = M;
      static T op(const T lhs, const T rhs) {
        return ::std::max(lhs, rhs);
      }
      static T e() {
        if constexpr (::std::is_integral_v<M>) {
          return ::std::numeric_limits<M>::min();
        } else {
          return -::std::numeric_limits<M>::infinity();
        }
      }
    };

    template <typename M, M E>
    struct max<M, E> {
      static_assert(::std::is_integral_v<M>, "M must be a built-in integral type.");

      using T = M;
      static T op(const T lhs, const T rhs) {
        assert(E <= lhs);
        assert(E <= rhs);
        return ::std::max(lhs, rhs);
      }
      static T e() {
        return E;
      }
    };

    template <typename M, M ...dummy>
    struct min;

    template <typename M>
    struct min<M> {
      static_assert(::std::is_arithmetic_v<M>, "M must be a built-in arithmetic type.");

      using T = M;
      static T op(const T lhs, const T rhs) {
        return ::std::min(lhs, rhs);
      }
      static T e() {
        if constexpr (::std::is_integral_v<M>) {
          return ::std::numeric_limits<M>::max();
        } else {
          return ::std::numeric_limits<M>::infinity();
        }
      }
    };

    template <typename M, M E>
    struct min<M, E> {
      static_assert(::std::is_integral_v<M>, "M must be a built-in integral type.");

      using T = M;
      static T op(const T lhs, const T rhs) {
        assert(lhs <= E);
        assert(rhs <= E);
        return ::std::min(lhs, rhs);
      }
      static T e() {
        return E;
      }
    };

    template <typename M>
    struct multiplies {
    private:
      using VR = ::std::conditional_t<::std::is_arithmetic_v<M>, const M, const M&>;

    public:
      using T = M;
      static T op(VR lhs, VR rhs) {
        return lhs * rhs;
      }
      static T e() {
        return T(1);
      }
    };

    template <>
    struct multiplies<bool> {
      using T = bool;
      static T op(const bool lhs, const bool rhs) {
        return lhs && rhs;
      }
      static T e() {
        return true;
      }
    };

    template <typename M>
    struct gcd {
    private:
      static_assert(!::std::is_arithmetic_v<M> || (::std::is_integral_v<M> && !::std::is_same_v<M, bool>), "If M is a built-in arithmetic type, it must be integral except for bool.");
      using VR = ::std::conditional_t<::std::is_arithmetic_v<M>, const M, const M&>;

    public:
      using T = M;
      static T op(VR lhs, VR rhs) {
        return ::tools::gcd(lhs, rhs);
      }
      static T e() {
        return T(0);
      }
    };

    template <typename M, M E>
    struct update {
      static_assert(::std::is_integral_v<M>, "M must be a built-in integral type.");

      using T = M;
      static T op(const T lhs, const T rhs) {
        return lhs == E ? rhs : lhs;
      }
      static T e() {
        return E;
      }
    };
  }
}


#line 1 "tools/square.hpp"



#line 1 "tools/is_monoid.hpp"



#line 5 "tools/is_monoid.hpp"
#include <utility>

namespace tools {

  template <typename M, typename = void>
  struct is_monoid : ::std::false_type {};

  template <typename M>
  struct is_monoid<M, ::std::enable_if_t<
    ::std::is_same_v<typename M::T, decltype(M::op(::std::declval<typename M::T>(), ::std::declval<typename M::T>()))> &&
    ::std::is_same_v<typename M::T, decltype(M::e())>
  , void>> : ::std::true_type {};

  template <typename M>
  inline constexpr bool is_monoid_v = ::tools::is_monoid<M>::value;
}


#line 6 "tools/square.hpp"

namespace tools {

  template <typename M>
  ::std::enable_if_t<::tools::is_monoid_v<M>, typename M::T> square(const typename M::T& x) {
    return M::op(x, x);
  }

  template <typename T>
  ::std::enable_if_t<!::tools::is_monoid_v<T>, T> square(const T& x) {
    return x * x;
  }
}


#line 9 "tools/pow.hpp"

namespace tools {

  template <typename M, typename E>
  ::std::enable_if_t<::std::is_integral_v<E>, typename M::T> pow(const typename M::T& base, const E exponent) {
    assert(exponent >= 0);
    return exponent == 0
      ? M::e()
      : exponent % 2 == 0
        ? ::tools::square<M>(::tools::pow<M>(base, exponent / 2))
        : M::op(::tools::pow<M>(base, exponent - 1), base);
  }

  template <typename T, typename E>
  ::std::enable_if_t<::std::is_integral_v<E>, T> pow(const T& base, const E exponent) {
    assert(exponent >= 0);
    return ::tools::pow<::tools::monoid::multiplies<T>>(base, exponent);
  }

  template <typename T, typename E>
  auto pow(const T base, const E exponent) -> ::std::enable_if_t<!::std::is_integral_v<E>, decltype(::std::pow(base, exponent))> {
    return ::std::pow(base, exponent);
  }
}


#line 1 "tools/vector4.hpp"



#line 1 "tools/vector.hpp"



#include <cstddef>
#include <array>
#include <initializer_list>
#line 9 "tools/vector.hpp"
#include <vector>
#line 11 "tools/vector.hpp"
#include <iterator>
#line 15 "tools/vector.hpp"
#include <iostream>
#include <string>
#include <functional>
#include <tuple>
#line 1 "tools/abs.hpp"



namespace tools {
  constexpr float abs(const float x) {
    return x < 0 ? -x : x;
  }
  constexpr double abs(const double x) {
    return x < 0 ? -x : x;
  }
  constexpr long double abs(const long double x) {
    return x < 0 ? -x : x;
  }
  constexpr int abs(const int x) {
    return x < 0 ? -x : x;
  }
  constexpr long abs(const long x) {
    return x < 0 ? -x : x;
  }
  constexpr long long abs(const long long x) {
    return x < 0 ? -x : x;
  }
  constexpr unsigned int abs(const unsigned int x) {
    return x;
  }
  constexpr unsigned long abs(const unsigned long x) {
    return x;
  }
  constexpr unsigned long long abs(const unsigned long long x) {
    return x;
  }
}


#line 1 "tools/tuple_hash.hpp"



#line 1 "tools/now.hpp"



#include <chrono>

namespace tools {
  inline long long now() {
    return ::std::chrono::duration_cast<::std::chrono::nanoseconds>(::std::chrono::high_resolution_clock::now().time_since_epoch()).count();
  }
}


#line 1 "tools/hash_combine.hpp"



#line 6 "tools/hash_combine.hpp"

// Source: https://github.com/google/cityhash/blob/f5dc54147fcce12cefd16548c8e760d68ac04226/src/city.h
// License: MIT
// Author: Google Inc.

// Copyright (c) 2011 Google, Inc.
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in
// all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
// THE SOFTWARE.

namespace tools {
  template <typename T>
  void hash_combine(::std::size_t& seed, const T& v) {
    static const ::std::hash<T> hasher;
    static constexpr ::std::size_t k_mul = 0x9ddfea08eb382d69ULL;
    ::std::size_t a = (hasher(v) ^ seed) * k_mul;
    a ^= (a >> 47);
    ::std::size_t b = (seed ^ a) * k_mul;
    b ^= (b >> 47);
    seed = b * k_mul;
  }
}


#line 11 "tools/tuple_hash.hpp"

namespace tools {
  template <typename... Ts>
  struct tuple_hash {
    template <::std::size_t I = sizeof...(Ts) - 1>
    ::std::size_t operator()(const ::std::tuple<Ts...>& key) const {
      if constexpr (I == ::std::numeric_limits<::std::size_t>::max()) {
        static const ::std::size_t seed = ::tools::now();
        return seed;
      } else {
        ::std::size_t seed = this->operator()<I - 1>(key);
        ::tools::hash_combine(seed, ::std::get<I>(key));
        return seed;
      }
    }
  };
}


#line 21 "tools/vector.hpp"

namespace tools {
  namespace detail {
    namespace vector {
      template <typename T, ::std::size_t N>
      class members {
      protected:
        constexpr static bool variable_sized = false;
        constexpr static bool has_aliases = false;
        ::std::array<T, N> m_values;
        members() : m_values() {}
        members(const ::std::initializer_list<T> il) : m_values(il) {
          assert(il.size() == N);
        }
      };

      template <typename T>
      class members<T, 2> {
      protected:
        constexpr static bool variable_sized = false;
        constexpr static bool has_aliases = true;
        ::std::array<T*, 2> m_values;
        members() : m_values{&this->x, &this->y} {}
        members(const T& x, const T& y) : m_values{&this->x, &this->y}, x(x), y(y) {}
        members(const ::std::initializer_list<T> il) : m_values{&this->x, &this->y}, x(il.begin()[0]), y(il.begin()[1]) {
          assert(il.size() == 2);
        }

      public:
        T x;
        T y;
      };

      template <typename T>
      class members<T, 3> {
      protected:
        constexpr static bool variable_sized = false;
        constexpr static bool has_aliases = true;
        ::std::array<T*, 3> m_values;
        members() : m_values{&this->x, &this->y, &this->z} {}
        members(const T& x, const T& y, const T& z) : m_values{&this->x, &this->y, &this->z}, x(x), y(y), z(z) {}
        members(const ::std::initializer_list<T> il) : m_values{&this->x, &this->y, &this->z}, x(il.begin()[0]), y(il.begin()[1]), z(il.begin()[2]) {
          assert(il.size() == 3);
        }

      public:
        T x;
        T y;
        T z;
      };

      template <typename T>
      class members<T, 4> {
      protected:
        constexpr static bool variable_sized = false;
        constexpr static bool has_aliases = true;
        ::std::array<T*, 4> m_values;
        members() : m_values{&this->x, &this->y, &this->z, &this->w} {}
        members(const T& x, const T& y, const T& z, const T& w) : m_values{&this->x, &this->y, &this->z, &this->w}, x(x), y(y), z(z), w(w) {}
        members(const ::std::initializer_list<T> il) : m_values{&this->x, &this->y, &this->z, &this->w}, x(il.begin()[0]), y(il.begin()[1]), z(il.begin()[2]), w(il.begin()[3]) {
          assert(il.size() == 4);
        }

      public:
        T x;
        T y;
        T z;
        T w;
      };

      template <typename T>
      class members<T, ::std::numeric_limits<::std::size_t>::max()> {
      protected:
        constexpr static bool variable_sized = true;
        constexpr static bool has_aliases = false;
        ::std::vector<T> m_values;
        members() = default;
        members(const ::std::size_t n) : m_values(n) {}
        members(const ::std::size_t n, const T& value) : m_values(n, value) {}
        template <typename InputIter>
        members(const InputIter first, const InputIter last) : m_values(first, last) {}
        members(const ::std::initializer_list<T> il) : m_values(il) {}
      };
    }
  }

  template <typename T, ::std::size_t N = ::std::numeric_limits<::std::size_t>::max()>
  class vector : public ::tools::detail::vector::members<T, N> {
  private:
    using Base = ::tools::detail::vector::members<T, N>;
    using F = ::std::conditional_t<::std::is_floating_point_v<T>, T, double>;
    using V = ::tools::vector<T, N>;
    constexpr static bool variable_sized = Base::variable_sized;
    constexpr static bool has_aliases = Base::has_aliases;

  public:
    using reference = T&;
    using const_reference = const T&;
    using size_type = ::std::size_t;
    using difference_type = ::std::ptrdiff_t;
    using pointer = T*;
    using const_pointer = const T*;
    using value_type = T;
    class iterator {
    private:
      V* m_parent;
      size_type m_i;

    public:
      using difference_type = ::std::ptrdiff_t;
      using value_type = T;
      using reference = T&;
      using pointer = T*;
      using iterator_category = ::std::random_access_iterator_tag;

      iterator() = default;
      iterator(const iterator&) = default;
      iterator(iterator&&) = default;
      ~iterator() = default;
      iterator& operator=(const iterator&) = default;
      iterator& operator=(iterator&&) = default;

      iterator(V * const parent, const size_type i) : m_parent(parent), m_i(i) {}

      reference operator*() const {
        return (*this->m_parent)[this->m_i];
      }
      pointer operator->() const {
        return &(*(*this));
      }

      iterator& operator++() {
        ++this->m_i;
        return *this;
      }
      iterator operator++(int) {
        const auto self = *this;
        ++*this;
        return self;
      }
      iterator& operator--() {
        --this->m_i;
        return *this;
      }
      iterator operator--(int) {
        const auto self = *this;
        --*this;
        return self;
      }
      iterator& operator+=(const difference_type n) {
        this->m_i += n;
        return *this;
      }
      iterator& operator-=(const difference_type n) {
        this->m_i -= n;
        return *this;
      }
      friend iterator operator+(const iterator& self, const difference_type n) {
        return iterator(self) += n;
      }
      friend iterator operator+(const difference_type n, const iterator& self) {
        return iterator(self) += n;
      }
      friend iterator operator-(const iterator& self, const difference_type n) {
        return iterator(self) -= n;
      }
      friend difference_type operator-(const iterator& lhs, const iterator& rhs) {
        assert(lhs.m_parent == rhs.m_parent);
        return static_cast<difference_type>(lhs.m_i) - static_cast<difference_type>(rhs.m_i);
      }
      reference operator[](const difference_type n) const {
        return *(*this + n);
      }

      friend bool operator==(const iterator& lhs, const iterator& rhs) {
        assert(lhs.m_parent == rhs.m_parent);
        return lhs.m_i == rhs.m_i;
      }
      friend bool operator!=(const iterator& lhs, const iterator& rhs) {
        assert(lhs.m_parent == rhs.m_parent);
        return lhs.m_i != rhs.m_i;
      }
      friend bool operator<(const iterator& lhs, const iterator& rhs) {
        assert(lhs.m_parent == rhs.m_parent);
        return lhs.m_i < rhs.m_i;
      }
      friend bool operator<=(const iterator& lhs, const iterator& rhs) {
        assert(lhs.m_parent == rhs.m_parent);
        return lhs.m_i <= rhs.m_i;
      }
      friend bool operator>(const iterator& lhs, const iterator& rhs) {
        assert(lhs.m_parent == rhs.m_parent);
        return lhs.m_i > rhs.m_i;
      }
      friend bool operator>=(const iterator& lhs, const iterator& rhs) {
        assert(lhs.m_parent == rhs.m_parent);
        return lhs.m_i >= rhs.m_i;
      }
    };
    class const_iterator {
    private:
      const V *m_parent;
      size_type m_i;

    public:
      using difference_type = ::std::ptrdiff_t;
      using value_type = T;
      using reference = const T&;
      using pointer = const T*;
      using iterator_category = ::std::random_access_iterator_tag;

      const_iterator() = default;
      const_iterator(const const_iterator&) = default;
      const_iterator(const_iterator&&) = default;
      ~const_iterator() = default;
      const_iterator& operator=(const const_iterator&) = default;
      const_iterator& operator=(const_iterator&&) = default;

      const_iterator(const V * const parent, const size_type i) : m_parent(parent), m_i(i) {}

      reference operator*() const {
        return (*this->m_parent)[this->m_i];
      }
      pointer operator->() const {
        return &(*(*this));
      }

      const_iterator& operator++() {
        ++this->m_i;
        return *this;
      }
      const_iterator operator++(int) {
        const auto self = *this;
        ++*this;
        return self;
      }
      const_iterator& operator--() {
        --this->m_i;
        return *this;
      }
      const_iterator operator--(int) {
        const auto self = *this;
        --*this;
        return self;
      }
      const_iterator& operator+=(const difference_type n) {
        this->m_i += n;
        return *this;
      }
      const_iterator& operator-=(const difference_type n) {
        this->m_i -= n;
        return *this;
      }
      friend const_iterator operator+(const const_iterator& self, const difference_type n) {
        return const_iterator(self) += n;
      }
      friend const_iterator operator+(const difference_type n, const const_iterator& self) {
        return const_iterator(self) += n;
      }
      friend const_iterator operator-(const const_iterator& self, const difference_type n) {
        return const_iterator(self) -= n;
      }
      friend difference_type operator-(const const_iterator& lhs, const const_iterator& rhs) {
        assert(lhs.m_parent == rhs.m_parent);
        return static_cast<difference_type>(lhs.m_i) - static_cast<difference_type>(rhs.m_i);
      }
      reference operator[](const difference_type n) const {
        return *(*this + n);
      }

      friend bool operator==(const const_iterator& lhs, const const_iterator& rhs) {
        assert(lhs.m_parent == rhs.m_parent);
        return lhs.m_i == rhs.m_i;
      }
      friend bool operator!=(const const_iterator& lhs, const const_iterator& rhs) {
        assert(lhs.m_parent == rhs.m_parent);
        return lhs.m_i != rhs.m_i;
      }
      friend bool operator<(const const_iterator& lhs, const const_iterator& rhs) {
        assert(lhs.m_parent == rhs.m_parent);
        return lhs.m_i < rhs.m_i;
      }
      friend bool operator<=(const const_iterator& lhs, const const_iterator& rhs) {
        assert(lhs.m_parent == rhs.m_parent);
        return lhs.m_i <= rhs.m_i;
      }
      friend bool operator>(const const_iterator& lhs, const const_iterator& rhs) {
        assert(lhs.m_parent == rhs.m_parent);
        return lhs.m_i > rhs.m_i;
      }
      friend bool operator>=(const const_iterator& lhs, const const_iterator& rhs) {
        assert(lhs.m_parent == rhs.m_parent);
        return lhs.m_i >= rhs.m_i;
      }
    };
    using reverse_iterator = ::std::reverse_iterator<iterator>;
    using const_reverse_iterator = ::std::reverse_iterator<const_iterator>;

    vector() = default;
    vector(const V& other) : Base() {
      if constexpr (has_aliases) {
        ::std::copy(other.begin(), other.end(), this->begin());
      } else {
        this->m_values = other.m_values;
      }
    }
    vector(V&& other) noexcept {
      if constexpr (has_aliases) {
        ::std::copy(other.begin(), other.end(), this->begin());
      } else {
        this->m_values = ::std::move(other.m_values);
      }
    }
    ~vector() = default;
    V& operator=(const V& other) {
      if constexpr (has_aliases) {
        ::std::copy(other.begin(), other.end(), this->begin());
      } else {
        this->m_values = other.m_values;
      }
      return *this;
    }
    V& operator=(V&& other) noexcept {
      if constexpr (has_aliases) {
        ::std::copy(other.begin(), other.end(), this->begin());
      } else {
        this->m_values = ::std::move(other.m_values);
      }
      return *this;
    }

    template <bool SFINAE = variable_sized, ::std::enable_if_t<SFINAE, ::std::nullptr_t> = nullptr>
    explicit vector(size_type n) : Base(n) {}
    template <bool SFINAE = variable_sized, ::std::enable_if_t<SFINAE, ::std::nullptr_t> = nullptr>
    vector(size_type n, const_reference value) : Base(n, value) {}
    template <typename InputIter, bool SFINAE = variable_sized, ::std::enable_if_t<SFINAE, ::std::nullptr_t> = nullptr>
    vector(const InputIter first, const InputIter last) : Base(first, last) {}
    template <bool SFINAE = N == 2, ::std::enable_if_t<SFINAE, ::std::nullptr_t> = nullptr>
    vector(const T& x, const T& y) : Base(x, y) {}
    template <bool SFINAE = N == 3, ::std::enable_if_t<SFINAE, ::std::nullptr_t> = nullptr>
    vector(const T& x, const T& y, const T& z) : Base(x, y, z) {}
    template <bool SFINAE = N == 4, ::std::enable_if_t<SFINAE, ::std::nullptr_t> = nullptr>
    vector(const T& x, const T& y, const T& z, const T& w) : Base(x, y, z, w) {}
    vector(const ::std::initializer_list<T> il) : Base(il) {}

    iterator begin() noexcept { return iterator(this, 0); }
    const_iterator begin() const noexcept { return const_iterator(this, 0); }
    const_iterator cbegin() const noexcept { return const_iterator(this, 0); }
    iterator end() noexcept { return iterator(this, this->size()); }
    const_iterator end() const noexcept { return const_iterator(this, this->size()); }
    const_iterator cend() const noexcept { return const_iterator(this, this->size()); }
    reverse_iterator rbegin() noexcept { return ::std::make_reverse_iterator(this->end()); }
    const_reverse_iterator rbegin() const noexcept { return ::std::make_reverse_iterator(this->end()); }
    const_reverse_iterator crbegin() const noexcept { return ::std::make_reverse_iterator(this->cend()); }
    reverse_iterator rend() noexcept { return ::std::make_reverse_iterator(this->begin()); }
    const_reverse_iterator rend() const noexcept { return ::std::make_reverse_iterator(this->begin()); }
    const_reverse_iterator crend() const noexcept { return ::std::make_reverse_iterator(this->cbegin()); }

    size_type size() const noexcept { return this->m_values.size(); }
    bool empty() const noexcept { return this->m_values.empty(); }

    reference operator[](const size_type n) {
      if constexpr (has_aliases) {
        return *this->m_values[n];
      } else {
        return this->m_values[n];
      }
    }
    const_reference operator[](const size_type n) const {
      if constexpr (has_aliases) {
        return *this->m_values[n];
      } else {
        return this->m_values[n];
      }
    }
    reference front() { return *this->begin(); }
    const_reference front() const { return *this->begin(); }
    reference back() { return *this->rbegin(); }
    const_reference back() const { return *this->rbegin(); }

    V operator+() const {
      return *this;
    }
    V operator-() const {
      V res = *this;
      for (auto& v : res) v = -v;
      return res;
    }
    V& operator+=(const V& other) {
      assert(this->size() == other.size());
      for (::std::size_t i = 0; i < this->size(); ++i) {
        (*this)[i] += other[i];
      }
      return *this;
    }
    friend V operator+(const V& lhs, const V& rhs) {
      return V(lhs) += rhs;
    }
    V& operator-=(const V& other) {
      assert(this->size() == other.size());
      for (::std::size_t i = 0; i < this->size(); ++i) {
        (*this)[i] -= other[i];
      }
      return *this;
    }
    friend V operator-(const V& lhs, const V& rhs) {
      return V(lhs) -= rhs;
    }
    V& operator*=(const T& c) {
      for (auto& v : *this) v *= c;
      return *this;
    }
    friend V operator*(const T& lhs, const V& rhs) {
      return V(rhs) *= lhs;
    }
    friend V operator*(const V& lhs, const T& rhs) {
      return V(lhs) *= rhs;
    }
    V& operator/=(const T& c) {
      for (auto& v : *this) v /= c;
      return *this;
    }
    friend V operator/(const V& lhs, const T& rhs) {
      return V(lhs) /= rhs;
    }

    friend bool operator==(const V& lhs, const V& rhs) {
      return ::std::equal(lhs.begin(), lhs.end(), rhs.begin(), rhs.end());
    }
    friend bool operator!=(const V& lhs, const V& rhs) {
      return !(lhs == rhs);
    }

    T inner_product(const V& other) const {
      assert(this->size() == other.size());
      T res{};
      for (::std::size_t i = 0; i < this->size(); ++i) {
        res += (*this)[i] * other[i];
      }
      return res;
    }
    T l1_norm() const {
      T res{};
      for (const auto& v : *this) {
        res += ::tools::abs(v);
      }
      return res;
    }
    T squared_l2_norm() const {
      return this->inner_product(*this);
    }
    F l2_norm() const {
      return ::std::sqrt(static_cast<F>(this->squared_l2_norm()));
    }
    template <bool SFINAE = ::std::is_floating_point_v<T>, ::std::enable_if_t<SFINAE, ::std::nullptr_t> = nullptr>
    V normalized() const {
      return *this / this->l2_norm();
    }

    friend ::std::ostream& operator<<(::std::ostream& os, const V& self) {
      os << '(';
      ::std::string delimiter = "";
      for (const auto& v : self) {
        os << delimiter << v;
        delimiter = ", ";
      }
      return os << ')';
    }
    friend ::std::istream& operator>>(::std::istream& is, V& self) {
      for (auto& v : self) {
        is >> v;
      }
      return is;
    }

    template <bool SFINAE = N == 2, ::std::enable_if_t<SFINAE, ::std::nullptr_t> = nullptr>
    T outer_product(const V& other) const {
      return this->x * other.y - this->y * other.x;
    }
    template <bool SFINAE = N == 2, ::std::enable_if_t<SFINAE, ::std::nullptr_t> = nullptr>
    V turned90() const {
      return V{-this->y, this->x};
    }
    template <bool SFINAE = N == 2, ::std::enable_if_t<SFINAE, ::std::nullptr_t> = nullptr>
    V turned270() const {
      return V{this->y, -this->x};
    }

    template <bool SFINAE = N == 2, ::std::enable_if_t<SFINAE, ::std::nullptr_t> = nullptr>
    static const ::std::array<V, 4>& four_directions() {
      static const ::std::array<V, 4> res = {
        V{T(1), T(0)},
        V{T(0), T(1)},
        V{T(-1), T(0)},
        V{T(0), T(-1)}
      };
      return res;
    }
    template <bool SFINAE = N == 2, ::std::enable_if_t<SFINAE, ::std::nullptr_t> = nullptr>
    static const ::std::array<V, 8>& eight_directions() {
      static const ::std::array<V, 8> res = {
        V{T(1), T(0)},
        V{T(1), T(1)},
        V{T(0), T(1)},
        V{T(-1), T(1)},
        V{T(-1), T(0)},
        V{T(-1), T(-1)},
        V{T(0), T(-1)},
        V{T(1), T(-1)}
      };
      return res;
    }

    template <bool SFINAE = N == 3, ::std::enable_if_t<SFINAE, ::std::nullptr_t> = nullptr>
    V outer_product(const V& other) const {
      return V{this->y * other.z - this->z * other.y, this->z * other.x - this->x * other.z, this->x * other.y - this->y * other.x};
    }
    template <bool SFINAE = N == 3 && ::std::is_floating_point_v<T>, ::std::enable_if_t<SFINAE, ::std::nullptr_t> = nullptr>
    ::std::array<V, 3> orthonormal_basis() const {
      assert((*this != V{0, 0, 0}));

      ::std::array<V, 3> v;
      v[0] = *this;
      v[1] = V{0, this->z, -this->y};
      if (v[1] == V{0, 0, 0}) {
        v[1] = V{-this->z, 0, this->x};
      }
      v[1] -= v[0].inner_product(v[1]) / v[0].inner_product(v[0]) * v[0];

      v[0] = v[0].normalized();
      v[1] = v[1].normalized();
      v[2] = v[0].outer_product(v[1]);

      return v;
    }
  };
}

namespace std {
  template <typename T>
  struct hash<::tools::vector<T, 2>> {
    using result_type = ::std::size_t;
    using argument_type = ::tools::vector<T, 2>;
    result_type operator()(const argument_type& key) const {
      static const ::tools::tuple_hash<T, T> hasher;
      return hasher(::std::make_tuple(key.x, key.y));
    }
  };
  template <typename T>
  struct hash<::tools::vector<T, 3>> {
    using result_type = ::std::size_t;
    using argument_type = ::tools::vector<T, 3>;
    result_type operator()(const argument_type& key) const {
      static const ::tools::tuple_hash<T, T, T> hasher;
      return hasher(::std::make_tuple(key.x, key.y, key.z));
    }
  };
  template <typename T>
  struct hash<::tools::vector<T, 4>> {
    using result_type = ::std::size_t;
    using argument_type = ::tools::vector<T, 4>;
    result_type operator()(const argument_type& key) const {
      static const ::tools::tuple_hash<T, T, T, T> hasher;
      return hasher(::std::make_tuple(key.x, key.y, key.z, key.w));
    }
  };
}


#line 5 "tools/vector4.hpp"

namespace tools {
  template <typename T>
  using vector4 = ::tools::vector<T, 4>;
}


#line 1 "tools/vector3.hpp"



#line 5 "tools/vector3.hpp"

namespace tools {
  template <typename T>
  using vector3 = ::tools::vector<T, 3>;
}


#line 15 "tools/quaternion.hpp"

namespace tools {
  template <typename T>
  class quaternion;

  template <typename T>
  ::tools::quaternion<T> exp(const ::tools::quaternion<T>& q);
  template <typename T>
  ::tools::quaternion<T> log(const ::tools::quaternion<T>& q);
  template <typename T>
  ::tools::quaternion<T> pow(const ::tools::quaternion<T>& base, T exponent);

  template <typename T>
  class quaternion {
    static_assert(::std::is_floating_point_v<T>);

  private:
    static constexpr T eps = 1e-5;
    static constexpr T pi = 3.14159265358979323846264338327950288419716939937510L;

  public:
    T x;
    T y;
    T z;
    T w;

    quaternion() = default;
    quaternion(const ::tools::quaternion<T>&) = default;
    quaternion(::tools::quaternion<T>&&) = default;
    ~quaternion() = default;
    ::tools::quaternion<T>& operator=(const ::tools::quaternion<T>&) = default;
    ::tools::quaternion<T>& operator=(::tools::quaternion<T>&&) = default;

    quaternion(const T x, const T y, const T z, const T w) : x(x), y(y), z(z), w(w) {}

    T real() const {
      return this->w;
    }
    void real(const T val) {
      this->w = val;
    }
    ::tools::vector3<T> imag() const {
      return ::tools::vector3<T>(this->x, this->y, this->z);
    }
    void imag(const ::tools::vector3<T>& val) {
      this->x = val.x;
      this->y = val.y;
      this->z = val.z;
    }

    T abs() const {
      return ::tools::vector4<T>(this->x, this->y, this->z, this->w).l2_norm();
    }
    T norm() const {
      return ::tools::vector4<T>(this->x, this->y, this->z, this->w).squared_l2_norm();
    }
    ::tools::quaternion<T> conj() const {
      return ::tools::quaternion<T>(-this->x, -this->y, -this->z, this->w);
    }

    T angle() const {
      assert(::std::abs(this->norm() - 1) <= eps);
      return ::std::acos(::std::clamp<T>(this->real(), -1, 1)) * 2;
    }
    ::tools::vector3<T> axis() const {
      assert(::std::abs(this->norm() - 1) <= eps);
      return ::std::abs(this->real()) >= 1.0 ? ::tools::vector3<T>(1, 0, 0) : this->imag().normalized();
    }

    ::tools::quaternion<T> operator+() const {
      return *this;
    }
    ::tools::quaternion<T> operator-() const {
      return ::tools::quaternion<T>(-this->x, -this->y, -this->z, -this->w);
    }

    ::tools::quaternion<T>& operator+=(const ::tools::quaternion<T>& other) {
      this->x += other.x;
      this->y += other.y;
      this->z += other.z;
      this->w += other.w;
      return *this;
    }
    friend ::tools::quaternion<T> operator+(const ::tools::quaternion<T>& lhs, const ::tools::quaternion<T>& rhs) {
      return ::tools::quaternion<T>(lhs) += rhs;
    }

    ::tools::quaternion<T>& operator-=(const ::tools::quaternion<T>& other) {
      this->x -= other.x;
      this->y -= other.y;
      this->z -= other.z;
      this->w -= other.w;
      return *this;
    }
    friend ::tools::quaternion<T> operator-(const ::tools::quaternion<T>& lhs, const ::tools::quaternion<T>& rhs) {
      return ::tools::quaternion<T>(lhs) -= rhs;
    }

    ::tools::quaternion<T>& operator*=(const T c) {
      this->x *= c;
      this->y *= c;
      this->z *= c;
      this->w *= c;
      return *this;
    }
    friend ::tools::quaternion<T> operator*(const ::tools::quaternion<T>& self, const T c) {
      return ::tools::quaternion<T>(self) *= c;
    }
    friend ::tools::quaternion<T> operator*(const T c, const ::tools::quaternion<T>& self) {
      return ::tools::quaternion<T>(self) *= c;
    }

    friend ::tools::quaternion<T> operator*(const ::tools::quaternion<T>& lhs, const ::tools::quaternion<T>& rhs) {
      const auto real = lhs.real() * rhs.real() - lhs.imag().inner_product(rhs.imag());
      const auto imag = lhs.real() * rhs.imag() + rhs.real() * lhs.imag() + lhs.imag().outer_product(rhs.imag());
      return ::tools::quaternion<T>(imag.x, imag.y, imag.z, real);
    }
    ::tools::quaternion<T>& operator*=(const ::tools::quaternion<T>& other) {
      return *this = *this * other;
    }
    friend ::tools::vector3<T> operator*(const ::tools::quaternion<T>& q, const ::tools::vector3<T>& v) {
      assert(::std::abs(q.norm() - 1) <= eps);
      return (q * ::tools::quaternion<T>(v.x, v.y, v.z, 0) * q.conj()).imag();
    }

    ::tools::quaternion<T>& operator/=(const T c) {
      this->x /= c;
      this->y /= c;
      this->z /= c;
      this->w /= c;
      return *this;
    }
    friend ::tools::quaternion<T> operator/(const ::tools::quaternion<T>& self, const T c) {
      return ::tools::quaternion<T>(self) /= c;
    }

    ::tools::quaternion<T> inv() const {
      const auto norm = this->norm();
      assert(norm != 0);
      return this->conj() / norm;
    }
    friend ::tools::quaternion<T> operator/(const ::tools::quaternion<T>& lhs, const ::tools::quaternion<T>& rhs) {
      return lhs * rhs.inv();
    }
    ::tools::quaternion<T>& operator/=(const ::tools::quaternion<T>& other) {
      return *this = *this / other;
    }

    friend bool operator==(const ::tools::quaternion<T>& lhs, const ::tools::quaternion<T>& rhs) {
      return lhs.x == rhs.x && lhs.y == rhs.y && lhs.z == rhs.z && lhs.w == rhs.w;
    }
    friend bool operator!=(const ::tools::quaternion<T>& lhs, const ::tools::quaternion<T>& rhs) {
      return !(lhs == rhs);
    }

    friend ::std::ostream& operator<<(::std::ostream& os, const ::tools::quaternion<T>& self) {
      ::std::ostringstream s;
      s.flags(os.flags());
      s.imbue(os.getloc());
      s.precision(os.precision());
      s << '(' << self.x << ',' << self.y << ',' << self.z << ',' << self.w << ')';
      return os << s.str();
    }

    static ::tools::quaternion<T> angle_axis(const T angle, const ::tools::vector3<T>& axis) {
      assert(axis != ::tools::vector3<T>(0, 0, 0));
      const auto real = ::std::cos(angle / 2);
      const auto imag = ::std::sin(angle / 2) * axis.normalized();
      return ::tools::quaternion<T>(imag.x, imag.y, imag.z, real);
    }

    static ::tools::quaternion<T> from_to_rotation(::tools::vector3<T> from_direction, ::tools::vector3<T> to_direction) {
      assert(from_direction != ::tools::vector3<T>(0, 0, 0));
      assert(to_direction != ::tools::vector3<T>(0, 0, 0));
      from_direction = from_direction.normalized();
      to_direction = to_direction.normalized();
      if (from_direction.inner_product(to_direction) <= -1 + eps) {
        return ::tools::quaternion<T>::angle_axis(pi, from_direction.orthonormal_basis()[1]);
      }
      const auto bisect = (from_direction + to_direction).normalized();
      const auto real = from_direction.inner_product(bisect);
      const auto imag = from_direction.outer_product(bisect);
      return ::tools::quaternion<T>(imag.x, imag.y, imag.z, real);
    }

    static ::tools::quaternion<T> look_rotation(::tools::vector3<T> forward, ::tools::vector3<T> upwards = ::tools::vector3<T>(0, 1, 0)) {
      assert(forward != ::tools::vector3<T>(0, 0, 0));
      assert(upwards != ::tools::vector3<T>(0, 0, 0));

      auto side = upwards.outer_product(forward);
      if (side == ::tools::vector3<T>(0, 0, 0)) {
        upwards = ::tools::vector3<T>(0, 1, 0);
        side = upwards.outer_product(forward);
      }
      if (side == ::tools::vector3<T>(0, 0, 0)) {
        side = forward.orthonormal_basis()[1];
      }
      upwards = forward.outer_product(side);

      forward = forward.normalized();
      upwards = upwards.normalized();

      const auto q1 = ::tools::quaternion<T>::from_to_rotation(::tools::vector3<T>(0, 0, 1), forward);
      const auto theta = ::std::atan2((q1 * ::tools::vector3<T>(0, 1, 0)).outer_product(upwards).inner_product(forward), (q1 * ::tools::vector3<T>(0, 1, 0)).inner_product(upwards));
      const auto q2 = ::tools::quaternion<T>::angle_axis(theta, forward);
      return q2 * q1;
    }

    static ::tools::quaternion<T> slerp(const ::tools::quaternion<T>& q0, const ::tools::quaternion<T>& q1, const T t) {
      return q0 * ::tools::pow(q0.inv() * q1, t);
    }

    static ::tools::quaternion<T> identity() {
      return ::tools::quaternion<T>(0, 0, 0, 1);
    }

    template <typename R>
    static ::tools::quaternion<T> random(R& engine) {
      static ::std::uniform_real_distribution<T> dist(0, 1);
      const auto u1 = dist(engine);
      const auto u2 = dist(engine);
      const auto u3 = dist(engine);
      return ::tools::quaternion<T>(
        ::std::sqrt(1 - u1) * ::std::sin(2 * pi * u2),
        ::std::sqrt(1 - u1) * ::std::cos(2 * pi * u2),
        ::std::sqrt(u1) * ::std::sin(2 * pi * u3),
        ::std::sqrt(u1) * ::std::cos(2 * pi * u3)
      );
    }

    static ::std::array<::tools::quaternion<T>, 24> dice_rotations() {
      return ::std::array<::tools::quaternion<T>, 24>({
        ::tools::quaternion<double>::identity(),
        ::tools::quaternion<double>::angle_axis(-0.5 * pi, ::tools::vector3<double>(1, 0, 0)),
        ::tools::quaternion<double>::angle_axis(0.5 * pi, ::tools::vector3<double>(1, 0, 0)),
        ::tools::quaternion<double>::angle_axis(pi, ::tools::vector3<double>(1, 0, 0)),
        ::tools::quaternion<double>::angle_axis(-0.5 * pi, ::tools::vector3<double>(0, 1, 0)),
        ::tools::quaternion<double>::angle_axis(0.5 * pi, ::tools::vector3<double>(0, 1, 0)),
        ::tools::quaternion<double>::angle_axis(pi, ::tools::vector3<double>(0, 1, 0)),
        ::tools::quaternion<double>::angle_axis(-0.5 * pi, ::tools::vector3<double>(0, 0, 1)),
        ::tools::quaternion<double>::angle_axis(0.5 * pi, ::tools::vector3<double>(0, 0, 1)),
        ::tools::quaternion<double>::angle_axis(pi, ::tools::vector3<double>(0, 0, 1)),
        ::tools::quaternion<double>::angle_axis(-2.0 / 3.0 * pi, ::tools::vector3<double>(-1, -1, 1)),
        ::tools::quaternion<double>::angle_axis(2.0 / 3.0 * pi, ::tools::vector3<double>(-1, -1, 1)),
        ::tools::quaternion<double>::angle_axis(-2.0 / 3.0 * pi, ::tools::vector3<double>(-1, 1, 1)),
        ::tools::quaternion<double>::angle_axis(2.0 / 3.0 * pi, ::tools::vector3<double>(-1, 1, 1)),
        ::tools::quaternion<double>::angle_axis(-2.0 / 3.0 * pi, ::tools::vector3<double>(1, -1, 1)),
        ::tools::quaternion<double>::angle_axis(2.0 / 3.0 * pi, ::tools::vector3<double>(1, -1, 1)),
        ::tools::quaternion<double>::angle_axis(-2.0 / 3.0 * pi, ::tools::vector3<double>(1, 1, 1)),
        ::tools::quaternion<double>::angle_axis(2.0 / 3.0 * pi, ::tools::vector3<double>(1, 1, 1)),
        ::tools::quaternion<double>::angle_axis(pi, ::tools::vector3<double>(0, -1, 1)),
        ::tools::quaternion<double>::angle_axis(pi, ::tools::vector3<double>(0, 1, 1)),
        ::tools::quaternion<double>::angle_axis(pi, ::tools::vector3<double>(-1, 0, 1)),
        ::tools::quaternion<double>::angle_axis(pi, ::tools::vector3<double>(1, 0, 1)),
        ::tools::quaternion<double>::angle_axis(pi, ::tools::vector3<double>(-1, 1, 0)),
        ::tools::quaternion<double>::angle_axis(pi, ::tools::vector3<double>(1, 1, 0))
      });
    }
  };

  template <typename T>
  ::tools::quaternion<T> exp(const ::tools::quaternion<T>& q) {
    const auto inorm = q.imag().l2_norm();
    if (inorm == 0) {
      return ::tools::quaternion<T>(0, 0, 0, ::std::exp(q.real()));
    }

    const auto rexp = ::std::exp(q.real());
    const auto real = rexp * ::std::cos(inorm);
    const auto imag = rexp * ::std::sin(inorm) / inorm * q.imag();
    return ::tools::quaternion<T>(imag.x, imag.y, imag.z, real);
  }

  template <typename T>
  ::tools::quaternion<T> log(const ::tools::quaternion<T>& q) {
    assert(q != ::tools::quaternion<T>(0, 0, 0, 0));
    const auto inorm = q.imag().l2_norm();
    const auto uimag = inorm == 0 ? ::tools::vector3<T>(1, 0, 0) : q.imag() / inorm;
    const auto real = ::std::log(q.abs());
    const auto imag = ::std::atan2(inorm, q.real()) * uimag;
    return ::tools::quaternion<T>(imag.x, imag.y, imag.z, real);
  }

  template <typename T>
  ::tools::quaternion<T> pow(const ::tools::quaternion<T>& base, const T exponent) {
    const auto norm = base.norm();
    if (norm == 0) {
      assert(exponent != 0);
      return ::tools::quaternion<T>(0, 0, 0, 0);
    }

    return ::tools::exp(exponent * ::tools::log(base));
  }
}