diff --git a/demos/ComputerGraphics/build.sh b/demos/ComputerGraphics/build.sh
new file mode 100755
index 000000000..1003e1013
--- /dev/null
+++ b/demos/ComputerGraphics/build.sh
@@ -0,0 +1,2 @@
+#!/bin/bash
+clang -Ofast -Xclang -add-plugin -Xclang clad -Xclang -load -Xclang /build/clad/build-clad-tracer/lib/clad.so -I../../include/ -I/usr/lib/gcc/x86_64-linux-gnu/10/include -x c++ -std=c++14 -lstdc++ -lm SmallPT-1.cpp -fopenmp=libiomp5 -o SmallPT-1 "$@"
diff --git a/demos/ComputerGraphics/cpp-smallpt-d/README.md b/demos/ComputerGraphics/cpp-smallpt-d/README.md
new file mode 100644
index 000000000..681b31fee
--- /dev/null
+++ b/demos/ComputerGraphics/cpp-smallpt-d/README.md
@@ -0,0 +1,6 @@
+Demo project 'cpp-smallpt-d' is base on:
+  C++ modification of Kevin Baeson's 99 line C++ path tracer
+  https://github.com/matt77hias/cpp-smallpt
+
+  A Differentiable Renderer supports reverse Parameter and Scene
+  Reconstruction.
diff --git a/demos/ComputerGraphics/cpp-smallpt-d/build.sh b/demos/ComputerGraphics/cpp-smallpt-d/build.sh
new file mode 100755
index 000000000..56eb18c05
--- /dev/null
+++ b/demos/ComputerGraphics/cpp-smallpt-d/build.sh
@@ -0,0 +1,4 @@
+#!/bin/bash
+clang -Ofast -Xclang -add-plugin -Xclang clad -Xclang -load -Xclang /build/clad/build-clad-tracer/lib/clad.so -I../../../include/ -x c++ -std=c++14 -lstdc++ -lm cpp-smallpt.cpp -fopenmp=libiomp5 -o cpp-smallpt "$@"
+
+clang -Ofast -Xclang -add-plugin -Xclang clad -Xclang -load -Xclang /build/clad/build-clad-tracer/lib/clad.so -I../../../include/ -x c++ -std=c++14 -lstdc++ -lm diff-tracer-1.cpp -fopenmp=libiomp5 -o diff-tracer-1 "$@"
diff --git a/demos/ComputerGraphics/cpp-smallpt-d/cpp-smallpt b/demos/ComputerGraphics/cpp-smallpt-d/cpp-smallpt
new file mode 100755
index 000000000..32dc9fdcd
Binary files /dev/null and b/demos/ComputerGraphics/cpp-smallpt-d/cpp-smallpt differ
diff --git a/demos/ComputerGraphics/cpp-smallpt-d/cpp-smallpt.cpp b/demos/ComputerGraphics/cpp-smallpt-d/cpp-smallpt.cpp
new file mode 100644
index 000000000..8887bbb85
--- /dev/null
+++ b/demos/ComputerGraphics/cpp-smallpt-d/cpp-smallpt.cpp
@@ -0,0 +1,222 @@
+//--------------------------------------------------------------------*- C++ -*-
+// clad - The C++ Clang-based Automatic Differentiator
+//
+// A demo, describing how to use clad in simple Differentiable Path tracer.
+//
+// Author:  Alexander Penev <alexander_penev-at-yahoo.com>, 2022
+//          Based on smallpt, a Path Tracer by Kevin Beason, 2008
+// Based on: C++ modification of Kevin Baeson's 99 line C++ path tracer
+//           https://github.com/matt77hias/cpp-smallpt
+//----------------------------------------------------------------------------//
+
+// To compile the demo please type:
+// path/to/clang -O3 -Xclang -add-plugin -Xclang clad -Xclang -load -Xclang \
+// path/to/libclad.so -I../../include/ -x c++ -std=c++11 -lstdc++ -lm \
+// cpp-smallpt-d.cpp -fopenmp=libiomp5 -o SmallPT
+//
+// To run the demo please type:
+// ./cpp-smallpt-d 500 && xv image.ppm
+
+// A typical invocation would be:
+// ../../../../../obj/Debug+Asserts/bin/clang -O3 -Xclang -add-plugin -Xclang clad \
+// -Xclang -load -Xclang ../../../../../obj/Debug+Asserts/lib/libclad.dylib \
+// -I../../include/ -x c++ -std=c++11 -lstdc++ -lm cpp-smalpt-d.cpp -fopenmp=libiomp5 \
+// -o cpp-smallpt-d
+// ./cpp-smallpt-d 500 && xv image.ppm
+
+
+//-----------------------------------------------------------------------------
+// Includes
+//-----------------------------------------------------------------------------
+#pragma region
+
+#include "imageio.hpp"
+#include "sampling.hpp"
+#include "specular.hpp"
+#include "sphere.hpp"
+
+#pragma endregion
+
+//-----------------------------------------------------------------------------
+// System Includes
+//-----------------------------------------------------------------------------
+#pragma region
+
+#include <iterator>
+#include <memory>
+
+#pragma endregion
+
+//-----------------------------------------------------------------------------
+// Defines
+//-----------------------------------------------------------------------------
+#pragma region
+
+#define REFRACTIVE_INDEX_OUT 1.0
+#define REFRACTIVE_INDEX_IN 1.5
+
+#pragma endregion
+
+//-----------------------------------------------------------------------------
+// Declarations and Definitions
+//-----------------------------------------------------------------------------
+namespace smallpt {
+
+  /*constexpr*/ Sphere* scene[] = {
+      new Sphere(1e5,  Vector3(1e5 + 1, 40.8, 81.6),   Vector3(),   Vector3(0.75, 0.25, 0.25), Reflection_t::Diffuse),   // Left
+      new Sphere(1e5,  Vector3(-1e5 + 99, 40.8, 81.6), Vector3(),   Vector3(0.25, 0.25, 0.75), Reflection_t::Diffuse),   // Right
+      new Sphere(1e5,  Vector3(50, 40.8, 1e5),         Vector3(),   Vector3(0.75),             Reflection_t::Diffuse),   // Back
+      new Sphere(1e5,  Vector3(50, 40.8, -1e5 + 170),  Vector3(),   Vector3(),                 Reflection_t::Diffuse),   // Front
+      new Sphere(1e5,  Vector3(50, 1e5, 81.6),         Vector3(),   Vector3(0.75),             Reflection_t::Diffuse),   // Bottom
+      new Sphere(1e5,  Vector3(50, -1e5 + 81.6, 81.6), Vector3(),   Vector3(0.75),             Reflection_t::Diffuse),   // Top
+      new Sphere(16.5, Vector3(27, 16.5, 47),          Vector3(),   Vector3(0.999),            Reflection_t::Refractive),  // Glass
+//      new Sphere(16.5, Vector3(55, 30, 57),          Vector3(),   Vector3(0.999),            Reflection_t::Refractive),  // Glassr
+//      new Sphere(16.5, Vector3(80, 60, 67),          Vector3(),   Vector3(0.999),            Reflection_t::Refractive),  // Glass
+      new Sphere(16.5, Vector3(73, 16.5, 78),          Vector3(),   Vector3(0.999),            Reflection_t::Specular),// Mirror
+      new Sphere(600,  Vector3(50, 681.6 - .27, 81.6), Vector3(12), Vector3(),                 Reflection_t::Diffuse)    // Light
+  };
+
+  [[nodiscard]]
+  /*constexpr*/ size_t Intersect(Sphere* g_scene[], const size_t n_scene, const Ray& ray) noexcept {
+    size_t hit = SIZE_MAX;
+    for (size_t i = 0u; i < n_scene; ++i) {
+      if (g_scene[i]->Intersect(ray)) {
+        hit = i;
+      }
+    }
+
+    return hit;
+  }
+
+  [[nodiscard]]
+  static const Vector3 Radiance(Sphere* g_scene[], const size_t n_scene, const Ray& ray, RNG& rng) noexcept {
+    Ray r = ray;
+    Vector3 L;
+    Vector3 F(1.0);
+
+    while (true) {
+      const auto hit = Intersect(g_scene, n_scene, r);
+      if (hit == SIZE_MAX) {
+        return L;
+      }
+
+      const Sphere* shape = g_scene[hit];
+      const Vector3 p = r(r.m_tmax);
+      const Vector3 n = Normalize(p - shape->m_p);
+
+      L += F * shape->m_e;
+      F *= shape->m_f;
+
+      // Russian roulette
+      if (4u < r.m_depth) {
+        const double continue_probability = shape->m_f.Max();
+        if (rng.Uniform() >= continue_probability) {
+          return L;
+        }
+        F /= continue_probability;
+      }
+
+      // Next path segment
+      switch (shape->m_reflection_t) {
+
+        case Reflection_t::Specular: {
+          const Vector3 d = IdealSpecularReflect(r.m_d, n);
+          r = Ray(p, d, EPSILON_SPHERE, INFINITY, r.m_depth + 1u);
+          break;
+        }
+
+        case Reflection_t::Refractive: {
+          double pr;
+          const Vector3 d = IdealSpecularTransmit(r.m_d, n,
+                                                  REFRACTIVE_INDEX_OUT,
+                                                  REFRACTIVE_INDEX_IN, pr, rng);
+          F *= pr;
+          r = Ray(p, d, EPSILON_SPHERE, INFINITY, r.m_depth + 1u);
+          break;
+        }
+
+        default: {
+          const Vector3 w = (0.0 > n.Dot(r.m_d)) ? n : -n;
+          const Vector3 u = Normalize((std::abs(w.m_x) > 0.1
+                                           ? Vector3(0.0, 1.0, 0.0)
+                                           : Vector3(1.0, 0.0, 0.0))
+                                          .Cross(w));
+          const Vector3 v = w.Cross(u);
+
+          const Vector3
+              sample_d = CosineWeightedSampleOnHemisphere(rng.Uniform(),
+                                                          rng.Uniform());
+          const Vector3 d = Normalize(sample_d.m_x * u + sample_d.m_y * v +
+                                      sample_d.m_z * w);
+          r = Ray(p, d, EPSILON_SPHERE, INFINITY, r.m_depth + 1u);
+          break;
+        }
+      }
+    }
+  }
+
+  static void Render(
+      Sphere* g_scene[], const size_t n_scene, // Geometry, Lights
+      double Vx, double Vy, double Vz, // Params - Center of one sphere // must be Vector3()
+      const std::uint32_t w, const std::uint32_t h, std::uint32_t nb_samples, size_t fr, // Camera
+      Vector3 Ls[], const char* fileName // Result
+    ) noexcept {
+
+    RNG rng;
+
+    // Camera params
+    const Vector3 eye = {50.0, 52.0, 295.6};
+    const Vector3 gaze = Normalize(Vector3(0.0, -0.042612, -1.0));
+    const double fov = 0.5135;
+    const Vector3 cx = {w * fov / h, 0.0, 0.0};
+    const Vector3 cy = Normalize(cx.Cross(gaze)) * fov;
+
+    // Change unfixed geometry center
+    g_scene[6]->m_p = Vector3(Vx, Vy, Vz);
+
+    #pragma omp parallel for schedule(static)
+    for (int y = 0; y < static_cast<int>(h); ++y) { // pixel row
+      for (std::size_t x = 0u; x < w; ++x) { // pixel column
+        for (std::size_t sy = 0u, i = (h - 1u - y) * w + x; sy < 2u; ++sy) { // 2 subpixel row
+          for (std::size_t sx = 0u; sx < 2u; ++sx) { // 2 subpixel column
+            Vector3 L;
+
+            for (std::size_t s = 0u; s < nb_samples; ++s) { // samples per subpixel
+              const double u1 = 2.0 * rng.Uniform();
+              const double u2 = 2.0 * rng.Uniform();
+              const double dx = u1 < 1.0 ? sqrt(u1) - 1.0 : 1.0 - sqrt(2.0 - u1);
+              const double dy = u2 < 1.0 ? sqrt(u2) - 1.0 : 1.0 - sqrt(2.0 - u2);
+              const Vector3 d = cx * (((sx + 0.5 + dx) * 0.5 + x) / w - 0.5) +
+                                cy * (((sy + 0.5 + dy) * 0.5 + y) / h - 0.5) +
+                                gaze;
+
+              L += Radiance(g_scene, n_scene, Ray(eye + d * 130.0, Normalize(d), EPSILON_SPHERE), rng) * (1.0 / nb_samples);
+            }
+
+            Ls[i] += 0.25 * Clamp(L);
+          }
+        }
+      }
+    }
+
+    WritePPM(w, h, Ls, fileName);
+  }
+} // namespace smallpt
+
+#ifndef CPP_EMBED_SMALLPT
+using namespace smallpt;
+int main(int argc, char* argv[]) {
+  const std::uint32_t nb_samples = (2 == argc) ? atoi(argv[1]) / 4 : 1;
+  const std::uint32_t w = 1024;
+  const std::uint32_t h = 768;
+
+  smallpt::Render(
+    scene, *(&scene + 1) - scene, // Geometry, Lights
+    27, 16.5, 47, // Params - Center of one sphere // must be Vector3()
+    w, h, nb_samples, 0, // Camera
+    new Vector3[w*h], "image.ppm" // Result
+  );
+
+  return 0;
+}
+#endif
diff --git a/demos/ComputerGraphics/cpp-smallpt-d/diff-tracer-1.cpp b/demos/ComputerGraphics/cpp-smallpt-d/diff-tracer-1.cpp
new file mode 100644
index 000000000..d682dfe0d
--- /dev/null
+++ b/demos/ComputerGraphics/cpp-smallpt-d/diff-tracer-1.cpp
@@ -0,0 +1,294 @@
+//--------------------------------------------------------------------*- C++ -*-
+// clad - The C++ Clang-based Automatic Differentiator
+//
+// Demo application of a simple gradient descent algorithm to find
+// some params (one sphere position) in ray traced 3D scene.
+// Clad is used to calculate the gradient of the cost function.
+// Cost function is calculated as norm of differece between two images -
+// first is target image and second is interemediate image generated
+// during gradient decent algorithm.
+//
+//----------------------------------------------------------------------------//
+
+// To run the demo please type:
+// path/to/clang  -Xclang -add-plugin -Xclang clad -Xclang -load -Xclang \
+// path/to/libclad.so  -I../include/ -x c++ -lstdc++ -lm diff-tracer-1.cpp
+//
+// A typical invocation would be:
+// ../../../../obj/Debug+Asserts/bin/clang  -Xclang -add-plugin -Xclang clad \
+// -Xclang -load -Xclang ../../../../obj/Debug+Asserts/lib/libclad.dylib     \
+// -I../include/ -x c++ -lstdc++ -lm diff-tracer-1.cpp
+//
+// To plot the results:
+// Results are stored in images: image-target.ppm, and image-i-(0..N).ppm
+
+//#include <fstream>  // For plotting data.
+//#include <iostream> // For std::*
+//#include <vector>   // For std::vector.
+
+// Necessary for clad to work include
+#include "clad/Differentiator/Differentiator.h"
+
+// Include traces as embed "library" of functions
+#define CPP_EMBED_SMALLPT 1
+#include "cpp-smallpt.cpp"
+
+using namespace smallpt;
+
+// Structure to represent input dataset
+// w and h - represetn width and height of images
+// target, current - represents the input/target and
+// output data/images respectively
+struct Dataset {
+  const std::uint32_t w = 1024u;
+  const std::uint32_t h = 768u;
+  const std::uint32_t nb_samples;
+  Vector3* target;
+  Vector3* current;
+  double learning_rate = 1e-2;
+  // result
+  double Vx, Vy, Vz;
+  int step = 0;
+
+  // Populate the data with target image (render example + random noice)
+  Dataset(const std::uint32_t nb_samples = 1, const std::uint32_t w = 1024u, const std::uint32_t h = 768u) :
+    nb_samples(nb_samples), w(w), h(h) {
+    // Create Target data/image
+    char fileName[4096];
+    snprintf(fileName, sizeof(fileName), "image-target.ppm");
+    target = new Vector3[w*h];
+    current = new Vector3[w*h];
+    Render(
+      scene, *(&scene + 1) - scene, // Geometry, Lights
+      27, 16.5, 47, // Params - Center of one sphere // must be Vector3()
+      w, h, nb_samples, 0, // Camera
+      target, fileName // Result
+    );
+
+    // For initial position of sphere we use randomized values
+    Vx = 80.0 * (((double)std::rand())/RAND_MAX) + 10.0;
+    Vy = 80.0 * (((double)std::rand())/RAND_MAX) + 10.0;
+    Vz = 80.0 * (((double)std::rand())/RAND_MAX) + 10.0;
+  }
+};
+
+// Function to perform a minimization step
+// theta_x are the hypothesis parameters, dt is the generated dataset and
+// clad_grad is the gradient function generated by Clad
+template <typename T>
+void performStep(double& theta_0, double& theta_1, double& theta_2, Dataset dt, T clad_grad) {
+  double J_theta[3+2];
+  double result[3] = {0, 0, 0};
+  // current?
+  for (size_t i = 0; i < dt.w*dt.h; i++) {
+    J_theta[0] = J_theta[1] = J_theta[2] = J_theta[3] = 0;
+//    clad_grad.execute(theta_0, theta_1, theta_2, dt.target[i], dt.current[i],
+//        &J_theta[0], &J_theta[1], &J_theta[2], &J_theta[3], &J_theta[4]
+//    );
+
+    result[0] += J_theta[0];
+    result[1] += J_theta[1];
+    result[2] += J_theta[2];
+  }
+
+  theta_0 -= dt.learning_rate * result[0] / (2 * dt.w*dt.h);
+  theta_1 -= dt.learning_rate * result[1] / (2 * dt.w*dt.h);
+  theta_2 -= dt.learning_rate * result[2] / (2 * dt.w*dt.h);
+}
+
+// The cost function to minimize using gradient descent
+// theta_x are the parameters to learn; x, y are the inputs and outputs of f
+double cost(double theta_0, double theta_1, double theta_2, Dataset dt) {
+  char fileName[4096];
+  snprintf(fileName, sizeof(fileName), "image-%d.ppm", dt.step++);
+  Render(
+      scene, *(&scene + 1) - scene, // Geometry, Lights
+      dt.Vx, dt.Vy, dt.Vz, // Params - Center of one sphere // must be Vector3()
+      dt.w, dt.h, dt.nb_samples, 0, // Camera
+      dt.current, fileName // Result
+    );
+  for (int i=0; i<=dt.h; i++) {
+    for (int j=0; j<=dt.w; j++) {
+//      sum += 
+    }
+  }
+
+//  double f_x = f(theta_0, theta_1, theta_2, x);
+//  return (f_x - y) * (f_x - y);
+  return 0;
+}
+
+double f(double x, double y) {
+  double t = 0;
+  for (int i=0; i<10; i++) {
+    t += x*y;
+  }
+  return t;
+}
+
+double f1(double x[], int cnt) {
+  double sum = 0;
+  for (int i=0; i<cnt; i++) {
+    sum += x[i];
+  }
+  return sum;
+}
+
+// Function to optimize the cost function of interest
+// theta is the hypothesis parameter list and maxSteps is the maximum steps to
+// perform
+void optimize(double theta[3], Dataset dt, unsigned int maxSteps, double eps) {
+  auto diff = theta;
+  bool hasConverged = false;
+  int currentStep = 0;
+
+  // Call for Clad to differentiate the cost function specified before
+  auto clad_grad = clad::gradient(cost);
+
+  do {
+    performStep(theta[0], theta[1], theta[2], dt, clad_grad);
+
+    std::cout << "Steps #" << currentStep << " Theta 0: " << theta[0]
+              << " Theta 1: " << theta[1] << " Theta 2: " << theta[2] << std::endl;
+
+    hasConverged = abs(diff[0] - theta[0]) <= eps &&
+                   abs(diff[1] - theta[1]) <= eps &&
+                   abs(diff[2] - theta[2]) <= eps;
+
+    diff = theta;
+  } while (currentStep++ < maxSteps && !hasConverged);
+}
+
+int main(int argc, char* argv[]) {
+  const std::uint32_t nb_samples = (2 == argc) ? atoi(argv[1]) / 4 : 1;
+
+auto f_dx = clad::differentiate(f, "x");
+auto f_dy = clad::differentiate(f, "y");
+auto f_g = clad::gradient(f);
+
+auto f1_dx = clad::differentiate(f1, "x[0]");
+auto f1_g = clad::gradient(f1);
+
+  Dataset dt(nb_samples, 1024u, 768u);
+  double theta[] = {dt.Vx, dt.Vy, dt.Vz};
+  optimize(theta, dt, 10000, 1e-6);
+
+  std::cout << "Result: "
+            << "(" << theta[0] << ", " << theta[1] << ", " << theta[2] << ")" << std::endl;
+
+  return 0;
+}
+/*
+double f(double x, double y) {
+  double t = 0;
+  for (int i=0; i<10; i++) {
+    t += x*y;
+  }
+  return t;
+}
+
+double f_dx(double x, double y) {
+  double t = 0;
+  for (int i=0; i<cnt; i++) {
+    t += y;
+  }
+  return t;
+}
+
+double f_dy(double x, double y) {
+  double t = 0;
+  for (int i=0; i<cnt; i++) {
+    t += x;
+  }
+  return t;
+}
+
+double f_darg0(double x, double y) {
+    double _d_x = 1;
+    double _d_y = 0;
+    double _d_t = 0;
+    double t = 0;
+    {
+        int _d_i = 0;
+        for (int i = 0; i < 10; i++) {
+            _d_t += _d_x * y + x * _d_y;
+            t += x * y;
+        }
+    }
+    return _d_t;
+}
+
+double f_darg1(double x, double y) {
+    double _d_x = 0;
+    double _d_y = 1;
+    double _d_t = 0;
+    double t = 0;
+    {
+        int _d_i = 0;
+        for (int i = 0; i < 10; i++) {
+            _d_t += _d_x * y + x * _d_y;
+            t += x * y;
+        }
+    }
+    return _d_t;
+}
+
+void f_grad(double x, double y, clad::array_ref<double> _d_x, clad::array_ref<double> _d_y) {
+    double _d_t = 0;
+    unsigned long _t0;
+    int _d_i = 0;
+    clad::tape<double> _t1 = {};
+    clad::tape<double> _t2 = {};
+    double t = 0;
+    _t0 = 0;
+    for (int i = 0; i < 10; i++) {
+        _t0++;
+        t += clad::push(_t2, x) * clad::push(_t1, y);
+    }
+    double f_return = t;
+    goto _label0;
+  _label0:
+    _d_t += 1;
+    for (; _t0; _t0--) {
+        {
+            double _r_d0 = _d_t;
+            _d_t += _r_d0;
+            double _r0 = _r_d0 * clad::pop(_t1);
+            * _d_x += _r0;
+            double _r1 = clad::pop(_t2) * _r_d0;
+            * _d_y += _r1;
+            _d_t -= _r_d0;
+        }
+    }
+}
+
+void f_g(double x, double y, double result[2]) {
+  double t_x = 0;
+  double t_y = 0;
+  for (int i=0; i<cnt; i++) {
+    t_x += x;
+    t_y += y;
+  }
+  result[0] = t_x;
+  result[1] = t_y;
+}
+
+void f1(double x[], int cnt, double dx[]) {
+  double sum = 0;
+  for (int i=0; i<cnt; i++) {
+    sum += x[i]*x[i];
+  }
+  return sum;
+
+  double sum_0 = 0;
+  double sum_1 = 0;
+  for (int i=0; i<cnt; i++) {
+    sum_0 += (i==0) * 2*x[i];
+    sum_1 += (i==1) * 2*x[i];
+  }
+  dx[0] = sum_0;
+  dx[1] = sum_1;
+}
+
+
+*/
\ No newline at end of file
diff --git a/demos/ComputerGraphics/cpp-smallpt-d/geometry.hpp b/demos/ComputerGraphics/cpp-smallpt-d/geometry.hpp
new file mode 100644
index 000000000..e0d00cddb
--- /dev/null
+++ b/demos/ComputerGraphics/cpp-smallpt-d/geometry.hpp
@@ -0,0 +1,63 @@
+#pragma once
+
+//-----------------------------------------------------------------------------
+// Includes
+//-----------------------------------------------------------------------------
+#pragma region
+
+#include "vector.hpp"
+
+#pragma endregion
+
+//-----------------------------------------------------------------------------
+// Declarations and Definitions
+//-----------------------------------------------------------------------------
+namespace smallpt {
+
+  struct Ray {
+
+    //---------------------------------------------------------------------
+    // Constructors and Destructors
+    //---------------------------------------------------------------------
+
+    constexpr explicit Ray(
+        Vector3 o, Vector3 d, double tmin = 0.0,
+        double tmax = std::numeric_limits<double>::infinity(),
+        std::uint32_t depth = 0u) noexcept
+        : m_o(std::move(o)), m_d(std::move(d)), m_tmin(tmin), m_tmax(tmax),
+          m_depth(depth){};
+    constexpr Ray(const Ray& ray) noexcept = default;
+    constexpr Ray(Ray&& ray) noexcept = default;
+    ~Ray() = default;
+
+    //---------------------------------------------------------------------
+    // Assignment Operators
+    //---------------------------------------------------------------------
+
+    Ray& operator=(const Ray& ray) = default;
+    Ray& operator=(Ray&& ray) = default;
+
+    //---------------------------------------------------------------------
+    // Member Methods
+    //---------------------------------------------------------------------
+
+    [[nodiscard]]
+    constexpr const Vector3 operator()(double t) const noexcept {
+      return m_o + m_d * t;
+    }
+
+    //---------------------------------------------------------------------
+    // Member Variables
+    //---------------------------------------------------------------------
+
+    Vector3 m_o, m_d;
+    mutable double m_tmin, m_tmax;
+    std::uint32_t m_depth;
+  };
+
+  inline std::ostream& operator<<(std::ostream& os, const Ray& r) {
+    os << "o: " << r.m_o << std::endl;
+    os << "d: " << r.m_d << std::endl;
+    return os;
+  }
+} // namespace smallpt
diff --git a/demos/ComputerGraphics/cpp-smallpt-d/imageio.hpp b/demos/ComputerGraphics/cpp-smallpt-d/imageio.hpp
new file mode 100644
index 000000000..35cd280d4
--- /dev/null
+++ b/demos/ComputerGraphics/cpp-smallpt-d/imageio.hpp
@@ -0,0 +1,38 @@
+#pragma once
+
+//-----------------------------------------------------------------------------
+// Includes
+//-----------------------------------------------------------------------------
+#pragma region
+
+#include "vector.hpp"
+
+#pragma endregion
+
+//-----------------------------------------------------------------------------
+// System Includes
+//-----------------------------------------------------------------------------
+#pragma region
+
+#include <cstdio>
+
+#pragma endregion
+
+//-----------------------------------------------------------------------------
+// Declarations and Definitions
+//-----------------------------------------------------------------------------
+namespace smallpt {
+
+  inline void WritePPM(std::uint32_t w, std::uint32_t h, const Vector3* Ls,
+                       const char* fname = "cpp-smallpt-d.ppm") noexcept {
+
+    FILE* fp = fopen(fname, "w");
+
+    std::fprintf(fp, "P3\n%u %u\n%u\n", w, h, 255u);
+    for (std::size_t i = 0; i < w * h; ++i) {
+      std::fprintf(fp, "%u %u %u ", ToByte(Ls[i].m_x), ToByte(Ls[i].m_y), ToByte(Ls[i].m_z));
+    }
+
+    std::fclose(fp);
+  }
+} // namespace smallpt
diff --git a/demos/ComputerGraphics/cpp-smallpt-d/math.hpp b/demos/ComputerGraphics/cpp-smallpt-d/math.hpp
new file mode 100644
index 000000000..592bcc504
--- /dev/null
+++ b/demos/ComputerGraphics/cpp-smallpt-d/math.hpp
@@ -0,0 +1,47 @@
+#pragma once
+
+//-----------------------------------------------------------------------------
+// System Includes
+//-----------------------------------------------------------------------------
+#pragma region
+
+#include <algorithm>
+#include <cmath>
+#include <cstdint>
+#include <limits>
+
+#pragma endregion
+
+//-----------------------------------------------------------------------------
+// Declarations and Definitions
+//-----------------------------------------------------------------------------
+namespace smallpt {
+
+  //-------------------------------------------------------------------------
+  // Constants
+  //-------------------------------------------------------------------------
+
+  constexpr double g_pi = 3.14159265358979323846;
+
+  //-------------------------------------------------------------------------
+  // Utilities
+  //-------------------------------------------------------------------------
+
+  template <class T>
+  [[nodiscard]]
+  constexpr const T& clamp(const T& v, const T& lo, const T& hi) noexcept {
+    return clamp(v, lo, hi, [](const T& a, const T& b) noexcept { return (a < b); });
+  }
+
+  template <class T, class Compare>
+  [[nodiscard]]
+  constexpr const T& clamp(const T& v, const T& lo, const T& hi, Compare comp) noexcept {
+    return comp(v, lo) ? lo : comp(hi, v) ? hi : v;
+  }
+
+  [[nodiscard]]
+  inline std::uint8_t ToByte(double color, double gamma = 2.2) noexcept {
+    const double gcolor = std::pow(color, 1.0 / gamma);
+    return static_cast<std::uint8_t>(clamp(255.0 * gcolor, 0.0, 255.0));
+  }
+} // namespace smallpt
diff --git a/demos/ComputerGraphics/cpp-smallpt-d/rng.hpp b/demos/ComputerGraphics/cpp-smallpt-d/rng.hpp
new file mode 100644
index 000000000..48546c1e1
--- /dev/null
+++ b/demos/ComputerGraphics/cpp-smallpt-d/rng.hpp
@@ -0,0 +1,56 @@
+#pragma once
+
+//-----------------------------------------------------------------------------
+// System Includes
+//-----------------------------------------------------------------------------
+#pragma region
+
+#include <cstdint>
+#include <random>
+#include <limits>
+
+#pragma endregion
+
+//-----------------------------------------------------------------------------
+// Declarations and Definitions
+//-----------------------------------------------------------------------------
+namespace smallpt {
+
+  class RNG {
+  public:
+    //---------------------------------------------------------------------
+    // Constructors and Destructors
+    //---------------------------------------------------------------------
+
+    explicit RNG(std::uint32_t seed = 606418532u) noexcept
+        : m_generator(), m_distribution() {
+      Seed(seed);
+    }
+    RNG(const RNG& rng) noexcept = default;
+    RNG(RNG&& rng) noexcept = default;
+    ~RNG() = default;
+
+    //---------------------------------------------------------------------
+    // Assignment Operators
+    //---------------------------------------------------------------------
+
+    RNG& operator=(const RNG& rng) = delete;
+    RNG& operator=(RNG&& rng) = delete;
+
+    //---------------------------------------------------------------------
+    // Member Methods
+    //---------------------------------------------------------------------
+
+    void Seed(uint32_t seed) noexcept { /*m_generator.seed(seed);*/ }
+
+    double Uniform() noexcept { return m_distribution(m_generator); }
+
+  private:
+    //---------------------------------------------------------------------
+    // Member Variables
+    //---------------------------------------------------------------------
+
+    std::default_random_engine m_generator;
+    std::uniform_real_distribution<double> m_distribution;
+  };
+} // namespace smallpt
diff --git a/demos/ComputerGraphics/cpp-smallpt-d/sampling.hpp b/demos/ComputerGraphics/cpp-smallpt-d/sampling.hpp
new file mode 100644
index 000000000..9a85dbe18
--- /dev/null
+++ b/demos/ComputerGraphics/cpp-smallpt-d/sampling.hpp
@@ -0,0 +1,40 @@
+#pragma once
+
+//-----------------------------------------------------------------------------
+// Includes
+//-----------------------------------------------------------------------------
+#pragma region
+
+#include "vector.hpp"
+
+#pragma endregion
+
+//-----------------------------------------------------------------------------
+// Declarations and Definitions
+//-----------------------------------------------------------------------------
+namespace smallpt {
+
+  [[nodiscard]]
+  inline const Vector3 UniformSampleOnSphere(double u1, double u2) noexcept {
+    const double cos_theta = 1.0 - 2.0 * u1;
+    const double sin_theta = std::sqrt(std::max(0.0, 1.0 - cos_theta * cos_theta));
+    const double phi = 2.0 * g_pi * u2;
+    return {std::cos(phi) * sin_theta, std::sin(phi) * sin_theta, cos_theta};
+  }
+
+  [[nodiscard]]
+  inline const Vector3 UniformSampleOnHemisphere(double u1, double u2) noexcept {
+    // u1 := cos_theta
+    const double sin_theta = std::sqrt(std::max(0.0, 1.0 - u1 * u1));
+    const double phi = 2.0 * g_pi * u2;
+    return {std::cos(phi) * sin_theta, std::sin(phi) * sin_theta, u1};
+  }
+
+  [[nodiscard]]
+  inline const Vector3 CosineWeightedSampleOnHemisphere(double u1, double u2) noexcept {
+    const double cos_theta = sqrt(1.0 - u1);
+    const double sin_theta = sqrt(u1);
+    const double phi = 2.0 * g_pi * u2;
+    return {std::cos(phi) * sin_theta, std::sin(phi) * sin_theta, cos_theta};
+  }
+} // namespace smallpt
diff --git a/demos/ComputerGraphics/cpp-smallpt-d/specular.hpp b/demos/ComputerGraphics/cpp-smallpt-d/specular.hpp
new file mode 100644
index 000000000..56d4bcb27
--- /dev/null
+++ b/demos/ComputerGraphics/cpp-smallpt-d/specular.hpp
@@ -0,0 +1,67 @@
+#pragma once
+
+//-----------------------------------------------------------------------------
+// Includes
+//-----------------------------------------------------------------------------
+#pragma region
+
+#include "rng.hpp"
+#include "vector.hpp"
+
+#pragma endregion
+
+//-----------------------------------------------------------------------------
+// Declarations and Definitions
+//-----------------------------------------------------------------------------
+namespace smallpt {
+
+  [[nodiscard]]
+  constexpr double Reflectance0(double n1, double n2) noexcept {
+    const double sqrt_R0 = (n1 - n2) / (n1 + n2);
+    return sqrt_R0 * sqrt_R0;
+  }
+
+  [[nodiscard]]
+  constexpr double SchlickReflectance(double n1, double n2, double c) noexcept {
+    const double R0 = Reflectance0(n1, n2);
+    return R0 + (1.0 - R0) * c * c * c * c * c;
+  }
+
+  [[nodiscard]]
+  constexpr const Vector3 IdealSpecularReflect(const Vector3& d, const Vector3& n) noexcept {
+    return d - 2.0 * n.Dot(d) * n;
+  }
+
+  [[nodiscard]]
+  inline const Vector3 IdealSpecularTransmit(const Vector3& d, const Vector3& n, double n_out,
+                        double n_in, double& pr, RNG& rng) noexcept {
+    const Vector3 d_Re = IdealSpecularReflect(d, n);
+
+    const bool out_to_in = (0.0 > n.Dot(d));
+    const Vector3 nl = out_to_in ? n : -n;
+    const double nn = out_to_in ? n_out / n_in : n_in / n_out;
+    const double cos_theta = d.Dot(nl);
+    const double cos2_phi = 1.0 - nn * nn * (1.0 - cos_theta * cos_theta);
+
+    // Total Internal Reflection
+    if (0.0 > cos2_phi) {
+      pr = 1.0;
+      return d_Re;
+    }
+
+    const Vector3 d_Tr = Normalize(nn * d - nl * (nn * cos_theta + sqrt(cos2_phi)));
+    const double c = 1.0 - (out_to_in ? -cos_theta : d_Tr.Dot(n));
+
+    const double Re = SchlickReflectance(n_out, n_in, c);
+    const double p_Re = 0.25 + 0.5 * Re;
+    if (rng.Uniform() < p_Re) {
+      pr = (Re / p_Re);
+      return d_Re;
+    } else {
+      const double Tr = 1.0 - Re;
+      const double p_Tr = 1.0 - p_Re;
+      pr = (Tr / p_Tr);
+      return d_Tr;
+    }
+  }
+} // namespace smallpt
diff --git a/demos/ComputerGraphics/cpp-smallpt-d/sphere.hpp b/demos/ComputerGraphics/cpp-smallpt-d/sphere.hpp
new file mode 100644
index 000000000..b89fc6009
--- /dev/null
+++ b/demos/ComputerGraphics/cpp-smallpt-d/sphere.hpp
@@ -0,0 +1,116 @@
+#pragma once
+
+//-----------------------------------------------------------------------------
+// Includes
+//-----------------------------------------------------------------------------
+#pragma region
+
+#include "geometry.hpp"
+
+#pragma endregion
+
+//-----------------------------------------------------------------------------
+// Defines
+//-----------------------------------------------------------------------------
+#pragma region
+
+#define EPSILON_SPHERE 1e-4
+
+#pragma endregion
+
+//-----------------------------------------------------------------------------
+// Declarations and Definitions
+//-----------------------------------------------------------------------------
+namespace smallpt {
+
+  //-------------------------------------------------------------------------
+  // Declarations and Definitions: Reflection_t
+  //-------------------------------------------------------------------------
+
+  enum struct Reflection_t : std::uint8_t {
+    Diffuse = 0u,
+    Specular,
+    Refractive
+  };
+
+  //-------------------------------------------------------------------------
+  // Declarations and Definitions: Sphere
+  //-------------------------------------------------------------------------
+
+  struct Sphere {
+
+    //---------------------------------------------------------------------
+    // Constructors and Destructors
+    //---------------------------------------------------------------------
+
+    constexpr explicit Sphere(double r, Vector3 p, Vector3 e, Vector3 f,
+                              Reflection_t reflection_t) noexcept
+        : m_r(r), m_p(std::move(p)), m_e(std::move(e)), m_f(std::move(f)),
+          m_reflection_t(reflection_t) {}
+    constexpr Sphere(const Sphere& sphere) noexcept = default;
+    constexpr Sphere(Sphere&& sphere) noexcept = default;
+    ~Sphere() = default;
+
+    //---------------------------------------------------------------------
+    // Assignment Operators
+    //---------------------------------------------------------------------
+
+    Sphere& operator=(const Sphere& sphere) = default;
+    Sphere& operator=(Sphere&& sphere) = default;
+
+    //---------------------------------------------------------------------
+    // Member Methods
+    //---------------------------------------------------------------------
+
+    [[nodiscard]]
+    constexpr bool Intersect(const Ray& ray) const noexcept {
+      // (o + t*d - p) . (o + t*d - p) - r*r = 0
+      // <=> (d . d) * t^2 + 2 * d . (o - p) * t + (o - p) . (o - p) - r*r = 0
+      //
+      // Discriminant check
+      // (2 * d . (o - p))^2 - 4 * (d . d) * ((o - p) . (o - p) - r*r) <? 0
+      // <=> (d . (o - p))^2 - (d . d) * ((o - p) . (o - p) - r*r) <? 0
+      // <=> (d . op)^2 - 1 * (op . op - r*r) <? 0
+      // <=> b^2 - (op . op) + r*r <? 0
+      // <=> D <? 0
+      //
+      // Solutions
+      // t = (- 2 * d . (o - p) +- 2 * sqrt(D)) / (2 * (d . d))
+      // <=> t = dop +- sqrt(D)
+
+      const Vector3 op = m_p - ray.m_o;
+      const double dop = ray.m_d.Dot(op);
+      const double D = dop * dop - op.Dot(op) + m_r * m_r;
+
+      if (0.0 > D) {
+        return false;
+      }
+
+      const double sqrtD = sqrt(D);
+
+      const double tmin = dop - sqrtD;
+      if (ray.m_tmin < tmin && tmin < ray.m_tmax) {
+        ray.m_tmax = tmin;
+        return true;
+      }
+
+      const double tmax = dop + sqrtD;
+      if (ray.m_tmin < tmax && tmax < ray.m_tmax) {
+        ray.m_tmax = tmax;
+        return true;
+      }
+
+      return false;
+    }
+
+    //---------------------------------------------------------------------
+    // Member Variables
+    //---------------------------------------------------------------------
+
+    double m_r;
+    Vector3 m_p; // position
+    Vector3 m_e; // emission
+    Vector3 m_f; // reflection
+    Reflection_t m_reflection_t;
+  };
+} // namespace smallpt
diff --git a/demos/ComputerGraphics/cpp-smallpt-d/vector.hpp b/demos/ComputerGraphics/cpp-smallpt-d/vector.hpp
new file mode 100644
index 000000000..4b3b46c17
--- /dev/null
+++ b/demos/ComputerGraphics/cpp-smallpt-d/vector.hpp
@@ -0,0 +1,310 @@
+#pragma once
+
+//-----------------------------------------------------------------------------
+// Includes
+//-----------------------------------------------------------------------------
+#pragma region
+
+#include "math.hpp"
+
+#pragma endregion
+
+//-----------------------------------------------------------------------------
+// System Includes
+//-----------------------------------------------------------------------------
+#pragma region
+
+#include <iostream>
+
+#pragma endregion
+
+//-----------------------------------------------------------------------------
+// Declarations and Definitions
+//-----------------------------------------------------------------------------
+namespace smallpt {
+
+  //-------------------------------------------------------------------------
+  // Vector3
+  //-------------------------------------------------------------------------
+
+  struct Vector3 {
+
+  public:
+    //---------------------------------------------------------------------
+    // Constructors and Destructors
+    //---------------------------------------------------------------------
+
+    constexpr explicit Vector3(double xyz = 0.0) noexcept
+        : Vector3(xyz, xyz, xyz) {}
+    constexpr Vector3(double x, double y, double z) noexcept
+        : m_x(x), m_y(y), m_z(z) {}
+    constexpr Vector3(const Vector3& v) noexcept = default;
+    constexpr Vector3(Vector3&& v) noexcept = default;
+    ~Vector3() = default;
+
+    //---------------------------------------------------------------------
+    // Assignment Operators
+    //---------------------------------------------------------------------
+
+    Vector3& operator=(const Vector3& v) = default;
+    Vector3& operator=(Vector3&& v) = default;
+
+    //---------------------------------------------------------------------
+    // Member Methods
+    //---------------------------------------------------------------------
+
+    [[nodiscard]] bool HasNaNs() const noexcept {
+      return std::isnan(m_x) || std::isnan(m_y) || std::isnan(m_z);
+    }
+
+    [[nodiscard]] constexpr const Vector3 operator-() const noexcept {
+      return {-m_x, -m_y, -m_z};
+    }
+
+    [[nodiscard]] constexpr const Vector3
+    operator+(const Vector3& v) const noexcept {
+      return {m_x + v.m_x, m_y + v.m_y, m_z + v.m_z};
+    }
+
+    [[nodiscard]] constexpr const Vector3
+    operator-(const Vector3& v) const noexcept {
+      return {m_x - v.m_x, m_y - v.m_y, m_z - v.m_z};
+    }
+
+    [[nodiscard]] constexpr const Vector3
+    operator*(const Vector3& v) const noexcept {
+      return {m_x * v.m_x, m_y * v.m_y, m_z * v.m_z};
+    }
+
+    [[nodiscard]] constexpr const Vector3
+    operator/(const Vector3& v) const noexcept {
+      return {m_x / v.m_x, m_y / v.m_y, m_z / v.m_z};
+    }
+
+    [[nodiscard]] constexpr const Vector3 operator+(double a) const noexcept {
+      return {m_x + a, m_y + a, m_z + a};
+    }
+
+    [[nodiscard]] constexpr const Vector3 operator-(double a) const noexcept {
+      return {m_x - a, m_y - a, m_z - a};
+    }
+
+    [[nodiscard]] constexpr const Vector3 operator*(double a) const noexcept {
+      return {m_x * a, m_y * a, m_z * a};
+    }
+
+    [[nodiscard]] constexpr const Vector3 operator/(double a) const noexcept {
+      const double inv_a = 1.0 / a;
+      return {m_x * inv_a, m_y * inv_a, m_z * inv_a};
+    }
+
+    Vector3& operator+=(const Vector3& v) noexcept {
+      m_x += v.m_x;
+      m_y += v.m_y;
+      m_z += v.m_z;
+      return *this;
+    }
+
+    Vector3& operator-=(const Vector3& v) noexcept {
+      m_x -= v.m_x;
+      m_y -= v.m_y;
+      m_z -= v.m_z;
+      return *this;
+    }
+
+    Vector3& operator*=(const Vector3& v) noexcept {
+      m_x *= v.m_x;
+      m_y *= v.m_y;
+      m_z *= v.m_z;
+      return *this;
+    }
+
+    Vector3& operator/=(const Vector3& v) noexcept {
+      m_x /= v.m_x;
+      m_y /= v.m_y;
+      m_z /= v.m_z;
+      return *this;
+    }
+
+    Vector3& operator+=(double a) noexcept {
+      m_x += a;
+      m_y += a;
+      m_z += a;
+      return *this;
+    }
+
+    Vector3& operator-=(double a) noexcept {
+      m_x -= a;
+      m_y -= a;
+      m_z -= a;
+      return *this;
+    }
+
+    Vector3& operator*=(double a) noexcept {
+      m_x *= a;
+      m_y *= a;
+      m_z *= a;
+      return *this;
+    }
+
+    Vector3& operator/=(double a) noexcept {
+      const double inv_a = 1.0 / a;
+      m_x *= inv_a;
+      m_y *= inv_a;
+      m_z *= inv_a;
+      return *this;
+    }
+
+    [[nodiscard]] constexpr double Dot(const Vector3& v) const noexcept {
+      return m_x * v.m_x + m_y * v.m_y + m_z * v.m_z;
+    }
+
+    [[nodiscard]] constexpr const Vector3
+    Cross(const Vector3& v) const noexcept {
+      return {m_y * v.m_z - m_z * v.m_y, m_z * v.m_x - m_x * v.m_z,
+              m_x * v.m_y - m_y * v.m_x};
+    }
+
+    [[nodiscard]] constexpr bool operator==(const Vector3& rhs) const {
+      return m_x == rhs.m_x && m_y == rhs.m_y && m_z == rhs.m_z;
+    }
+
+    [[nodiscard]] constexpr bool operator!=(const Vector3& rhs) const {
+      return !(*this == rhs);
+    }
+
+    [[nodiscard]] double& operator[](std::size_t i) noexcept {
+      return (&m_x)[i];
+    }
+
+    [[nodiscard]] constexpr double operator[](std::size_t i) const noexcept {
+      return (&m_x)[i];
+    }
+
+    [[nodiscard]] constexpr std::size_t MinDimension() const noexcept {
+      return (m_x < m_y && m_x < m_z) ? 0u : ((m_y < m_z) ? 1u : 2u);
+    }
+
+    [[nodiscard]] constexpr std::size_t MaxDimension() const noexcept {
+      return (m_x > m_y && m_x > m_z) ? 0u : ((m_y > m_z) ? 1u : 2u);
+    }
+
+    [[nodiscard]] constexpr double Min() const noexcept {
+      return std::min(m_x, std::min(m_y, m_z));
+    }
+    [[nodiscard]] constexpr double Max() const noexcept {
+      return std::max(m_x, std::max(m_y, m_z));
+    }
+
+    [[nodiscard]] constexpr double Norm2_squared() const noexcept {
+      return m_x * m_x + m_y * m_y + m_z * m_z;
+    }
+
+    [[nodiscard]] double Norm2() const noexcept {
+      return std::sqrt(Norm2_squared());
+    }
+
+    void Normalize() noexcept {
+      const double a = 1.0 / Norm2();
+      m_x *= a;
+      m_y *= a;
+      m_z *= a;
+    }
+
+    //---------------------------------------------------------------------
+    // Member Variables
+    //---------------------------------------------------------------------
+
+    double m_x, m_y, m_z;
+  };
+
+  //-------------------------------------------------------------------------
+  // Vector3 Utilities
+  //-------------------------------------------------------------------------
+
+  std::ostream& operator<<(std::ostream& os, const Vector3& v) {
+    os << '[' << v.m_x << ' ' << v.m_y << ' ' << v.m_z << ']';
+    return os;
+  }
+
+  [[nodiscard]] constexpr const Vector3 operator+(double a,
+                                                  const Vector3& v) noexcept {
+    return {a + v.m_x, a + v.m_y, a + v.m_z};
+  }
+
+  [[nodiscard]] constexpr const Vector3 operator-(double a,
+                                                  const Vector3& v) noexcept {
+    return {a - v.m_x, a - v.m_y, a - v.m_z};
+  }
+
+  [[nodiscard]] constexpr const Vector3 operator*(double a,
+                                                  const Vector3& v) noexcept {
+    return {a * v.m_x, a * v.m_y, a * v.m_z};
+  }
+
+  [[nodiscard]] constexpr const Vector3 operator/(double a,
+                                                  const Vector3& v) noexcept {
+    return {a / v.m_x, a / v.m_y, a / v.m_z};
+  }
+
+  [[nodiscard]] inline const Vector3 Sqrt(const Vector3& v) noexcept {
+    return {std::sqrt(v.m_x), std::sqrt(v.m_y), std::sqrt(v.m_z)};
+  }
+
+  [[nodiscard]] inline const Vector3 Pow(const Vector3& v, double a) noexcept {
+    return {std::pow(v.m_x, a), std::pow(v.m_y, a), std::pow(v.m_z, a)};
+  }
+
+  [[nodiscard]] inline const Vector3 Abs(const Vector3& v) noexcept {
+    return {std::abs(v.m_x), std::abs(v.m_y), std::abs(v.m_z)};
+  }
+
+  [[nodiscard]] constexpr const Vector3 Min(const Vector3& v1,
+                                            const Vector3& v2) noexcept {
+    return {std::min(v1.m_x, v2.m_x), std::min(v1.m_y, v2.m_y),
+            std::min(v1.m_z, v2.m_z)};
+  }
+
+  [[nodiscard]] constexpr const Vector3 Max(const Vector3& v1,
+                                            const Vector3& v2) noexcept {
+    return {std::max(v1.m_x, v2.m_x), std::max(v1.m_y, v2.m_y),
+            std::max(v1.m_z, v2.m_z)};
+  }
+
+  [[nodiscard]] inline const Vector3 Round(const Vector3& v) noexcept {
+    return {std::round(v.m_x), std::round(v.m_y), std::round(v.m_z)};
+  }
+
+  [[nodiscard]] inline const Vector3 Floor(const Vector3& v) noexcept {
+    return {std::floor(v.m_x), std::floor(v.m_y), std::floor(v.m_z)};
+  }
+
+  [[nodiscard]] inline const Vector3 Ceil(const Vector3& v) noexcept {
+    return {std::ceil(v.m_x), std::ceil(v.m_y), std::ceil(v.m_z)};
+  }
+
+  [[nodiscard]] inline const Vector3 Trunc(const Vector3& v) noexcept {
+    return {std::trunc(v.m_x), std::trunc(v.m_y), std::trunc(v.m_z)};
+  }
+
+  [[nodiscard]] constexpr const Vector3
+  Clamp(const Vector3& v, double low = 0.0, double high = 1.0) noexcept {
+
+    return {clamp(v.m_x, low, high), clamp(v.m_y, low, high),
+            clamp(v.m_z, low, high)};
+  }
+  [[nodiscard]] constexpr const Vector3 Lerp(double a, const Vector3& v1,
+                                             const Vector3& v2) noexcept {
+    return v1 + a * (v2 - v1);
+  }
+
+  template <std::size_t X, std::size_t Y, std::size_t Z>
+  [[nodiscard]] constexpr const Vector3 Permute(const Vector3& v) noexcept {
+    return {v[X], v[Y], v[Z]};
+  }
+
+  [[nodiscard]] inline const Vector3 Normalize(const Vector3& v) noexcept {
+    const double a = 1.0 / v.Norm2();
+    return a * v;
+  }
+} // namespace smallpt