<|start_filename|>include/xorshift.hpp<|end_filename|>
#ifndef XORSHIFT_HPP
#define XORSHIFT_HPP

#include <array>
#include <cstdint>
#include <limits>

// Xorshift implementation from [triSYCL](https://github.com/triSYCL/triSYCL)
template <int Nbits = 32> struct xorshift {
  /// Bit size of the random generator data type
  static auto constexpr bit_size = Nbits;

  /// The default initial internal state of the pseudo random generator
  static auto constexpr initial_state = [] {
    if constexpr (bit_size == 32) {
      /* Default value for xor32 from Marsaglia,
         p. 4 Section "3 Application to Xorshift RNGs" */
      return std::uint32_t { 2463534242 };
    } else if constexpr (bit_size == 64) {
      /* Default value for xor32 from Marsaglia,
         p. 4 Section "3 Application to Xorshift RNGs" */
      return std::uint64_t { 88172645463325252 };
    } else if constexpr (bit_size == 128) {
      /// Default value for xor128 from Marsaglia, p. 5 Section "4 Summary"
      return std::array<std::uint32_t, 4> { 123456789, 362436069, 521288629,
                                            88675123 };
    } else
      return nullptr;
  }();

  /// The type of the internal state of the generator
  using value_type = std::remove_const_t<decltype(initial_state)>;

  /// The type of the result
  using result_type = value_type;

  static_assert(!std::is_same_v<value_type, std::nullptr_t>,
                "Bit size not implemented");

  /// The internal state of the pseudo random generator
  value_type state = initial_state;

  /// The minimum returned value
  static auto constexpr min() {
    // It cannot return 0
    return std::numeric_limits<result_type>::min() + 1;
  };

  /// The maximum returned value
  static auto constexpr max() {
    return std::numeric_limits<result_type>::max();
  };

  /** Initialize the internal state from a user-given value

      Do not use the value 0 or the output is always 0
  */
  xorshift(const value_type& s)
      : state { s } {}

  xorshift() = default;

  /// Compute a new pseudo random integer
  const result_type& operator()() {
    if constexpr (bit_size == 32) {
      /* Pick the one of type "I" with best bit equidistribution from
         Panneton & L'Écuyer, Section "5.1 Equidistribution
         properties"

         X4 = (I + Ra )(I + Lb )(I + Rc )
      */
      state ^= state >> 7;
      state ^= state << 1;
      state ^= state >> 9;
    } else if constexpr (bit_size == 64) {
      /* The xor64 from Marsaglia, p. 4 Section "3 Application to
         Xorshift RNGs" */
      state ^= state << 13;
      state ^= state >> 7;
      state ^= state << 17;
    } else if constexpr (bit_size == 128) {
      // The xor128 from Marsaglia, p. 5 Section "4 Summary"
      auto t = state[0] ^ (state[0] << 11);
      state[0] = state[1];
      state[1] = state[2];
      state[2] = state[3];
      state[3] = state[3] ^ (state[3] >> 19) ^ (t ^ (t >> 8));
    } else
      // Just to error in that case
      return nullptr;

    return state;
  }
};
#endif // XORSHIFT_HPP


<|start_filename|>include/build_parameters.hpp<|end_filename|>
#ifndef BUILD_PARAMETERS_HPP
#define BUILD_PARAMETERS_HPP
namespace buildparams {

#ifdef USE_SINGLE_TASK
constexpr bool use_single_task = true;
#else
constexpr bool use_single_task = false;
#endif

#ifdef USE_SYCL_COMPILER
constexpr bool use_sycl_compiler = USE_SYCL_COMPILER;
#else
constexpr bool use_sycl_compiler = false;
#endif

constexpr int output_width = OUTPUT_WIDTH;
constexpr int output_height = OUTPUT_HEIGHT;
} // namespace buildparams

#endif // BUILD_PARAMETERS_HPP


<|start_filename|>include/box.hpp<|end_filename|>
#ifndef BOX_HPP
#define BOX_HPP

#include "rectangle.hpp"
#include "rtweekend.hpp"
#include "visit.hpp"

/// This class implements a axis aligned cuboid using 6 rectangles
class box {
 public:
  box() = default;

  /// p0 = { x0, y0, z0 } and p1 = { x1, y1. z1 }
  /// where x0 <= x1, y0 <= y1 and z0 <= z1
  box(const point& p0, const point& p1, const material_t& mat_type)
      : box_min { p0 }
      , box_max { p1 }
      , material_type { mat_type } {
    /// Add six sides of the box based on box_min and box_max to sides
    sides[0] = xy_rect(p0.x(), p1.x(), p0.y(), p1.y(), p1.z(), mat_type);
    sides[1] = xy_rect(p0.x(), p1.x(), p0.y(), p1.y(), p0.z(), mat_type);
    sides[2] = xz_rect(p0.x(), p1.x(), p0.z(), p1.z(), p1.y(), mat_type);
    sides[3] = xz_rect(p0.x(), p1.x(), p0.z(), p1.z(), p0.y(), mat_type);
    sides[4] = yz_rect(p0.y(), p1.y(), p0.z(), p1.z(), p1.x(), mat_type);
    sides[5] = yz_rect(p0.y(), p1.y(), p0.z(), p1.z(), p0.x(), mat_type);
  }

  /// Compute ray interaction with the box
  bool hit(auto& ctx, const ray& r, real_t min, real_t max, hit_record& rec,
           material_t& hit_material_type) const {
    hit_record temp_rec;
    material_t temp_material_type;
    auto hit_anything = false;
    auto closest_so_far = max;
    // Checking if the ray hits any of the sides
    for (const auto& side : sides) {
      if (dev_visit(
              [&](auto&& arg) {
                return arg.hit(ctx, r, min, closest_so_far, temp_rec,
                               temp_material_type);
              },
              side)) {
        hit_anything = true;
        closest_so_far = temp_rec.t;
        rec = temp_rec;
        hit_material_type = temp_material_type;
      }
    }
    return hit_anything;
  }

  point box_min;
  point box_max;
  material_t material_type;
  std::array<rectangle_t, 6> sides;
};

#endif


<|start_filename|>include/rectangle.hpp<|end_filename|>
#ifndef RECT_HPP
#define RECT_HPP

#include "material.hpp"
#include "ray.hpp"
#include "rtweekend.hpp"
#include "texture.hpp"
#include "vec.hpp"

/** The Following classes implement:

        -
   https://raytracing.github.io/books/RayTracingTheNextWeek.html#rectanglesandlights/creatingrectangleobjectsa
*/

class xy_rect {
 public:
  xy_rect() = default;

  /// x0 <= x1 and y0 <= y1
  xy_rect(real_t _x0, real_t _x1, real_t _y0, real_t _y1, real_t _k,
          const material_t& mat_type)
      : x0 { _x0 }
      , x1 { _x1 }
      , y0 { _y0 }
      , y1 { _y1 }
      , k { _k }
      , material_type { mat_type } {}

  /// Compute ray interaction with rectangle
  bool hit(auto&, const ray& r, real_t min, real_t max, hit_record& rec,
           material_t& hit_material_type) const {
    hit_material_type = material_type;

    auto t = (k - r.origin().z()) / r.direction().z();
    if (t < min || t > max)
      return false;
    auto x = r.origin().x() + t * r.direction().x();
    auto y = r.origin().y() + t * r.direction().y();
    if (x < x0 || x > x1 || y < y0 || y > y1)
      return false;
    rec.u = (x - x0) / (x1 - x0);
    rec.v = (y - y0) / (y1 - y0);
    rec.t = t;
    rec.p = r.at(rec.t);
    vec outward_normal = vec(0, 0, 1);
    rec.set_face_normal(r, outward_normal);
    return true;
  }
  real_t x0, x1, y0, y1, k;
  material_t material_type;
};

class xz_rect {
 public:
  xz_rect() = default;

  /// x0 <= x1 and z0 <= z1
  xz_rect(real_t _x0, real_t _x1, real_t _z0, real_t _z1, real_t _k,
          const material_t& mat_type)
      : x0 { _x0 }
      , x1 { _x1 }
      , z0 { _z0 }
      , z1 { _z1 }
      , k { _k }
      , material_type { mat_type } {}

  /// Compute ray interaction with rectangle
  bool hit(auto&, const ray& r, real_t min, real_t max, hit_record& rec,
           material_t& hit_material_type) const {
    hit_material_type = material_type;

    auto t = (k - r.origin().y()) / r.direction().y();
    if (t < min || t > max)
      return false;
    auto x = r.origin().x() + t * r.direction().x();
    auto z = r.origin().z() + t * r.direction().z();
    if (x < x0 || x > x1 || z < z0 || z > z1)
      return false;
    rec.u = (x - x0) / (x1 - x0);
    rec.v = (z - z0) / (z1 - z0);
    rec.t = t;
    rec.p = r.at(rec.t);
    vec outward_normal = vec(0, 1, 0);
    rec.set_face_normal(r, outward_normal);
    return true;
  }
  real_t x0, x1, z0, z1, k;
  material_t material_type;
};

class yz_rect {
 public:
  yz_rect() = default;

  /// y0 <= y1 and z0 <= z1
  yz_rect(real_t _y0, real_t _y1, real_t _z0, real_t _z1, real_t _k,
          const material_t& mat_type)
      : y0 { _y0 }
      , y1 { _y1 }
      , z0 { _z0 }
      , z1 { _z1 }
      , k { _k }
      , material_type { mat_type } {}

  /// Compute ray interaction with rectangle
  bool hit(auto&, const ray& r, real_t min, real_t max, hit_record& rec,
           material_t& hit_material_type) const {
    hit_material_type = material_type;

    auto t = (k - r.origin().x()) / r.direction().x();
    if (t < min || t > max)
      return false;
    auto y = r.origin().y() + t * r.direction().y();
    auto z = r.origin().z() + t * r.direction().z();
    if (y < y0 || y > y1 || z < z0 || z > z1)
      return false;
    rec.u = (y - y0) / (y1 - y0);
    rec.v = (z - z0) / (z1 - z0);
    rec.t = t;
    rec.p = r.at(rec.t);
    vec outward_normal = vec(1, 0, 0);
    rec.set_face_normal(r, outward_normal);
    return true;
  }
  real_t y0, y1, z0, z1, k;
  material_t material_type;
};

using rectangle_t = std::variant<xy_rect, xz_rect, yz_rect>;

#endif


<|start_filename|>include/rtweekend.hpp<|end_filename|>
#ifndef RT_SYCL_RTWEEKEND_HPP
#define RT_SYCL_RTWEEKEND_HPP

#include <algorithm>
#include <cmath>
#include <cstdlib>
#include <limits>
#include <memory>
#include <random>
#include <type_traits>
#include <variant>
#include <vector>

#include <sycl.hpp>

#include <build_parameters.hpp>
#include <xorshift.hpp>

// Constants

constexpr float infinity = std::numeric_limits<float>::infinity();
constexpr float pi = 3.1415926535897932385f;

// type aliases for float3 - vec, point and color
using point = sycl::float3;
using color = sycl::float3;
using vec = sycl::float3;

// Utility Functions

inline float degrees_to_radians(float degrees) { return degrees * pi / 180.0f; }

class LocalPseudoRNG {
 public:
  inline LocalPseudoRNG(std::uint32_t init_state = xorshift<>::initial_state)
      : generator { init_state } {}

  // Returns a random float in 0.f 1.
  inline float float_t() {
    constexpr float scale = 1.f / (uint64_t { 1 } << 32);
    return generator() * scale;
  }

  // Returns a random float in min, max
  inline float float_t(float min, float max) {
    // TODO use FMA ?
    return min + (max - min) * float_t();
  }

  // Returns a random vector with coordinates in 0.f 1.
  inline vec vec_t() { return { float_t(), float_t(), float_t() }; }

  // Returns a random vec with coordinates in min, max
  inline vec vec_t(float min, float max) {
    auto scale = max - min;
    return vec_t() * scale + min;
  }

  // Returns a random unit vector
  inline vec unit_vec() {
    auto x = float_t(-1.f, 1.f);
    auto maxy = sycl::sqrt(1 - x * x);
    auto y = float_t(-maxy, maxy);
    auto absz = sycl::sqrt(maxy * maxy - y * y);
    auto z = (float_t() > 0.5) ? absz : -absz;
    return vec(x, y, z);
  }

  // Returns a random vector in the unit ball of usual norm
  inline vec in_unit_ball() {
    // Polar coordinates r, theta, phi
    auto r = float_t();
    auto theta = float_t(0, 2 * pi);
    auto phi = float_t(0, pi);

    auto plan_seed = r * sycl::sin(phi);
    auto z = r * sycl::cos(phi);

    return { plan_seed * sycl::cos(theta), plan_seed * sycl::sin(theta), z };
  }

  // Return a random vector in the unit disk of usual norm in plane x, y
  inline vec in_unit_disk() {
    auto x = float_t(-1.f, 1.f);
    auto maxy = sycl::sqrt(1 - x * x);
    auto y = float_t(-maxy, maxy);
    return { x, y, 0.f };
  }

 private:
  xorshift<> generator;
};

/**
 @brief Used as a poorman's cooperative ersatz of device global variable
        The task context is (manually) passed through the call stack to all
        kernel callees
 */
struct task_context {
  LocalPseudoRNG rng;
  // See image_texture in texture.hpp for more details
  sycl::global_ptr<uint8_t> texture_data;
};

// Common Headers
#include "ray.hpp"
#include "vec.hpp"

#endif


<|start_filename|>include/triangle.hpp<|end_filename|>
#ifndef TRIANGLE_HPP
#define TRIANGLE_HPP

#include "material.hpp"
#include "ray.hpp"
#include "rtweekend.hpp"
#include "texture.hpp"
#include "vec.hpp"

struct _triangle_coord {
  point v0, v1, v2;
};

inline bool badouel_ray_triangle_intersec(const ray& r,
                                          _triangle_coord const& tri,
                                          real_t min, real_t max,
                                          hit_record& rec) {
  // Get triangle edge vectors and plane normal
  auto u = tri.v1 - tri.v0;
  auto v = tri.v2 - tri.v0;
  vec outward_normal = sycl::cross(u, v);

  auto w0 = r.origin() - tri.v0;
  auto a = -sycl::dot(outward_normal, w0);
  auto b = sycl::dot(outward_normal, r.direction());

  // ray is parallel to triangle plane
  if (sycl::fabs(b) < 0.000001f)
    return false;

  // intersection point of ray with triangle
  real_t length = a / b;
  if (length < 0)
    return false;
  else if (length < min || length > max)
    return false;

  vec hit_pt = r.at(length);
  auto uu = sycl::dot(u, u);
  auto uv = sycl::dot(u, v);
  auto vv = sycl::dot(v, v);
  auto w = hit_pt - tri.v0;
  auto wu = sycl::dot(w, u);
  auto wv = sycl::dot(w, v);
  auto D = uv * uv - uu * vv;

  auto s = (uv * wv - vv * wu) / D;
  auto t = (uv * wu - uu * wv) / D;
  if (s < 0.0f || s > 1.0f || t < 0.0f || (s + t) > 1.0f)
    return false;

  rec.set_face_normal(r, outward_normal);
  rec.t = length;
  rec.p = hit_pt;
  return true;
};

inline bool moller_trumbore_triangle_intersec(const ray& r,
                                              _triangle_coord const& tri,
                                              real_t min, real_t max,
                                              hit_record& rec) {
  constexpr auto epsilon = 0.0000001f;

  // Get triangle edge vectors and plane normal
  auto edge1 = tri.v1 - tri.v0;
  auto edge2 = tri.v2 - tri.v0;
  vec h = sycl::cross(r.direction(), edge2);
  auto a = sycl::dot(edge1, h);
  auto a_abs = sycl::fabs(a);

  if (a_abs < epsilon)
    return false; // This ray is parallel to this triangle

  auto a_pos = a > 0.f;

  auto s = r.origin() - tri.v0;
  auto u = sycl::dot(s, h);
  auto u_pos = u > 0.f;

  if ((u_pos xor a_pos) || sycl::fabs(u) > a_abs)
    return false;

  auto q = sycl::cross(s, edge1);
  auto v = sycl::dot(r.direction(), q);
  auto v_pos = v > 0.f;
  if ((v_pos xor a_pos) || (sycl::fabs(u + v) > a_abs))
    return false;

  auto length = sycl::dot(edge2, q) / a;

  if (length < min || length > max)
    return false;

  vec hit_pt = r.at(length);

  rec.set_face_normal(r, sycl::cross(edge1, edge2));
  rec.t = length;
  rec.p = hit_pt;
  return true;
};

// A triangle based on 3 points
template <auto IntersectionStrategy = moller_trumbore_triangle_intersec>
class _triangle : public _triangle_coord {
 public:
  _triangle() = default;
  _triangle(const point& _v0, const point& _v1, const point& _v2,
            const material_t& mat_type)
      : _triangle_coord { _v0, _v1, _v2 }
      , material_type { mat_type } {}

  /// Compute ray interaction with triangle
  bool hit(auto&, const ray& r, real_t min, real_t max, hit_record& rec,
           material_t& hit_material_type) const {
    hit_material_type = material_type;
    return IntersectionStrategy(r, *this, min, max, rec);
  }

  material_t material_type;
};

using triangle = _triangle<>;
#endif


<|start_filename|>include/sycl.hpp<|end_filename|>
#ifdef USE_SYCL_COMPILER
#include <CL/sycl.hpp>
using namespace cl;
#else
#include <sycl/sycl.hpp>
#endif


<|start_filename|>include/ray.hpp<|end_filename|>
#ifndef RT_SYCL_RAY_HPP
#define RT_SYCL_RAY_HPP

#include "vec.hpp"

class ray {
 public:
  ray() = default;

  ray(const point& origin, const vec& direction, real_t time = 0.0f)
      : orig { origin }
      , dir { direction }
      , tm { time } {}

  point origin() const { return orig; }
  vec direction() const { return dir; }
  real_t time() const { return tm; }

  // returns point along the ray at distance t from ray's origin
  // the ray P(t) = Origin + t*direction
  point at(float t) const { return orig + t * dir; }

 public:
  // To store the origin and direction of the ray
  point orig;
  vec dir;
  real_t tm;
};

#endif


<|start_filename|>include/vec.hpp<|end_filename|>
#ifndef RT_SYCL_VEC_HPP
#define RT_SYCL_VEC_HPP

#include "rtweekend.hpp"
#include <cmath>
#include <iostream>

using real_t = float;

// vec Utility Functions
inline float length_squared(const vec& v) {
  return sycl::fma(v.x(), v.x(), sycl::fma(v.y(), v.y(), v.z() * v.z()));
}

inline std::ostream& operator<<(std::ostream& out, const vec& v) {
  return out << v.x() << ' ' << v.y() << ' ' << v.z();
}

// Missing operator from the SYCL specification for now
vec operator-(const vec& u) { return vec(-u.x(), -u.y(), -u.z()); }

// Compute a unit vector from a non-null vector
inline vec unit_vector(const vec& v) { return v / sycl::length(v); }

// Compute reflected ray's direction
vec reflect(const vec& v, const vec& n) { return v - 2 * sycl::dot(v, n) * n; }

// Computes refracted ray's direction based on refractive index
vec refract(const vec& uv, const vec& n, float etai_over_etat) {
  auto cos_theta = sycl::fmin(-sycl::dot(uv, n), 1.0f);
  vec r_out_perp = etai_over_etat * (uv + cos_theta * n);
  vec r_out_parallel =
      -sycl::sqrt(sycl::fabs(1.0f - length_squared(r_out_perp))) * n;
  return r_out_perp + r_out_parallel;
}

#endif
