#include <cmath>
#include <concepts>

#define Float double

namespace fmath {

constexpr Float pi = 3.141592653589793;

template <std::floating_point T>
constexpr T clamp(T v, T mi, T ma) {
  return v <= mi ? mi : v >= ma ? ma : v;
}

template <std::floating_point T>
constexpr T saturate(T v) {
  return clamp(v, 0.0, 1.0);
}

template <std::floating_point T>
constexpr T smoothstep(T a, T b, T x) {
  T t = saturate((x - a) / (b - a));
  return t * t * (3.0 - (2.0 * t));
}

template <std::floating_point T>
constexpr T mix(T x, T y, T a) {
  return x * (1 - a) + y * a;
}

struct vec2 {
  vec2() = default;
  vec2(const vec2& v) = default;
  vec2& operator=(const vec2& v) = default;
  vec2(vec2&& v) = default;
  vec2& operator=(vec2&& v) = default;

  explicit vec2(Float v) : x{v}, y{v} {}
  vec2(Float x_, Float y_) : x{x_}, y{y_} {}

  vec2 operator+(const vec2& v) const { return vec2{x + v.x, y + v.y}; }
  vec2& operator+=(const vec2& v) {
    x += v.x;
    y += v.y;
    return *this;
  }
  vec2 operator+(Float c) const { return vec2{x + c, y + c}; }
  vec2& operator+=(Float c) {
    x += c;
    y += c;
    return *this;
  }
  vec2 operator-(const vec2& v) const { return vec2{x - v.x, y - v.y}; }
  vec2& operator-=(const vec2& v) {
    x -= v.x;
    y -= v.y;
    return *this;
  }
  vec2 operator-(Float c) const { return vec2{x - c, y - c}; }
  vec2& operator-=(Float c) {
    x -= c;
    y -= c;
    return *this;
  }
  vec2 operator*(const vec2& v) const { return vec2{x * v.x, y * v.y}; }
  vec2& operator*=(const vec2& v) {
    x *= v.x;
    y *= v.y;
    return *this;
  }
  vec2 operator*(Float c) const { return vec2{x * c, y * c}; }
  vec2& operator*=(Float c) {
    x *= c;
    y *= c;
    return *this;
  }
  vec2 operator/(const vec2& v) const { return vec2{x / v.x, y / v.y}; }
  vec2& operator/=(const vec2& v) {
    x /= v.x;
    y /= v.y;
    return *this;
  }
  vec2 operator/(Float c) const { return vec2{x / c, y / c}; }
  vec2& operator/=(Float c) {
    x /= c;
    y /= c;
    return *this;
  }

  vec2 operator-() { return vec2{-x, -y}; }

  friend vec2 operator+(Float c, vec2 v);
  friend vec2 operator-(Float c, vec2 v);
  friend vec2 operator*(Float c, vec2 v);
  friend vec2 operator/(Float c, vec2 v);

  Float x{0}, y{0};
};

struct vec3 {
  vec3() = default;
  vec3(const vec3& v) = default;
  vec3& operator=(const vec3& v) = default;
  vec3(vec3&& v) = default;
  vec3& operator=(vec3&& v) = default;

  explicit vec3(Float v) : x{v}, y{v}, z{v} {}
  vec3(Float x_, Float y_, Float z_) : x{x_}, y{y_}, z{z_} {}

  vec3 operator+(const vec3& v) const {
    return vec3{x + v.x, y + v.y, z + v.z};
  }
  vec3& operator+=(const vec3& v) {
    x += v.x;
    y += v.y;
    z += v.z;
    return *this;
  }
  vec3 operator+(Float c) const { return vec3{x + c, y + c, z + c}; }
  vec3& operator+=(Float c) {
    x += c;
    y += c;
    z += c;
    return *this;
  }
  vec3 operator-(const vec3& v) const {
    return vec3{x - v.x, y - v.y, z - v.z};
  }
  vec3& operator-=(const vec3& v) {
    x -= v.x;
    y -= v.y;
    z -= v.z;
    return *this;
  }
  vec3 operator-(Float c) const { return vec3{x - c, y - c, z - c}; }
  vec3& operator-=(Float c) {
    x -= c;
    y -= c;
    z -= c;
    return *this;
  }
  vec3 operator*(const vec3& v) const {
    return vec3{x * v.x, y * v.y, z * v.z};
  }
  vec3& operator*=(const vec3& v) {
    x *= v.x;
    y *= v.y;
    z *= v.z;
    return *this;
  }
  vec3 operator*(Float c) const { return vec3{x * c, y * c, z * c}; }
  vec3& operator*=(Float c) {
    x *= c;
    y *= c;
    z *= c;
    return *this;
  }
  vec3 operator/(const vec3& v) const {
    return vec3{x / v.x, y / v.y, z / v.z};
  }
  vec3& operator/=(const vec3& v) {
    x /= v.x;
    y /= v.y;
    z /= v.z;
    return *this;
  }
  vec3 operator/(Float c) const { return vec3{x / c, y / c, z / c}; }
  vec3& operator/=(Float c) {
    x /= c;
    y /= c;
    z /= c;
    return *this;
  }

  vec3 operator-() { return vec3{-x, -y, -z}; }

  friend vec3 operator+(Float c, vec3 v);
  friend vec3 operator-(Float c, vec3 v);
  friend vec3 operator*(Float c, vec3 v);
  friend vec3 operator/(Float c, vec3 v);

  Float x{0.0}, y{0.0}, z{0.0};
};

vec2 operator+(Float c, vec2 v) { return vec2{c + v.x, c + v.y}; }
vec2 operator-(Float c, vec2 v) { return vec2{c - v.x, c - v.y}; }
vec2 operator*(Float c, vec2 v) { return vec2{c * v.x, c * v.y}; }
vec2 operator/(Float c, vec2 v) { return vec2{c / v.x, c / v.y}; }

vec3 operator+(Float c, vec3 v) { return vec3{c + v.x, c + v.y, c + v.z}; }
vec3 operator-(Float c, vec3 v) { return vec3{c - v.x, c - v.y, c - v.z}; }
vec3 operator*(Float c, vec3 v) { return vec3{c * v.x, c * v.y, c * v.z}; }
vec3 operator/(Float c, vec3 v) { return vec3{c / v.x, c / v.y, c / v.z}; }

Float dot(const vec3& v0, const vec3& v1) {
  return v0.x * v1.x + v0.y * v1.y + v0.z * v1.z;
}

Float luminance(const vec3& v) {
  return v.x * 0.299 + v.y * 0.587 + v.z * 0.114;
}

Float pow2(Float x) { return x * x; }

Float pow3(Float x) { return x * x * x; }

Float pow1_4(Float x) { return std::sqrt(std::sqrt(x)); }

Float pow4(Float x) { return x * x * x * x; }

Float pow5(Float x) { return x * x * x * x * x; }

vec3 pow(const vec3& b, const vec3& e) {
  return vec3{std::pow(b.x, e.x), std::pow(b.y, e.y), std::pow(b.z, e.z)};
}

vec3 sin(const vec3& v) {
  return vec3{std::sin(v.x), std::sin(v.y), std::sin(v.z)};
}

vec3 cos(const vec3& v) {
  return vec3{std::cos(v.x), std::cos(v.y), std::cos(v.z)};
}

vec3 mix(const vec3& x, const vec3& y, const vec3& a) {
  return vec3{
      mix(x.x, y.x, a.x),
      mix(x.y, y.y, a.y),
      mix(x.z, y.z, a.z),
  };
}

vec3 cross(const vec3& v0, const vec3& v1) {
  vec3 v0_yzx = vec3(v0.y, v0.z, v0.x);
  vec3 v1_yzx = vec3(v1.y, v1.z, v1.x);
  vec3 v2 = v0 * v1_yzx - v0_yzx * v1;
  return vec3{v2.y, v2.z, v2.x};
}

vec3 sqrt(const vec3& v) {
  return vec3{std::sqrt(v.x), std::sqrt(v.y), std::sqrt(v.z)};
}

vec3 rsqrt(const vec3& v) { return 1.0 / fmath::sqrt(v); }

vec3 normalize(const vec3& v) { return rsqrt(vec3(dot(v, v))) * v; }

Float rcp(Float x) { return (Float)1.0 / x; }

}  // namespace fmath