#ifndef VEC3_H
#define VEC3_H

#include <cmath>
#include <limits>

#ifdef __CUDACC__
#define CUDA_CALLABLE __host__ __device__
#else
#define CUDA_CALLABLE
#endif

template <typename T> struct Vec3 {
  T x;
  T y;
  T z;

  CUDA_CALLABLE inline Vec3<T> operator+(Vec3<T> other) const {
    return {x + other.x, y + other.y, z + other.z};
  };

  CUDA_CALLABLE inline Vec3<T> operator-(Vec3<T> other) const {
    return {x - other.x, y - other.y, z - other.z};
  };

  CUDA_CALLABLE inline Vec3 scale(T scalar) {
    return {x * scalar, y * scalar, z * scalar};
  };

  CUDA_CALLABLE inline T dot(Vec3<T> other) const {
    return x * other.x + y * other.y + z * other.z;
  }

  CUDA_CALLABLE inline Vec3<T> cross(Vec3<T> other) const {
    return {y * other.z - z * other.y, z * other.x - x * other.z,
            x * other.y - y * other.x};
  }

  CUDA_CALLABLE inline T squared_norm2() const { return x * x + y * y + z * z; }

  CUDA_CALLABLE inline T norm2() const { return std::sqrt(squared_norm2()); }

  CUDA_CALLABLE inline Vec3<T> normalized() {
    // Add epsilon to the norm for stability when the norm is 0
    T norm = std::max(norm2(), std::numeric_limits<T>::epsilon());
    return {x / norm, y / norm, z / norm};
  }
};

#endif