commit a4b72e2f8595a99ec974bfd4935631b894514afc
Author: Zack Buhman <zack@buhman.org>
Date:   Mon Jan 23 21:08:13 2023 -0800

    initial: draw raytraced spheres on the sega saturn

diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..24318d3
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,7 @@
+*.gch
+*.o
+*.elf
+*.bin
+*.iso
+*.ppm
+*.png
diff --git a/Makefile b/Makefile
new file mode 100644
index 0000000..20e45b3
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,10 @@
+CFLAGS = -Isaturn -Imath
+OPT = -O3
+
+all: raytracing.iso
+
+LIB = ./saturn
+include $(LIB)/common.mk
+
+LIBGCC = $(shell $(CC) -print-file-name=libgcc.a)
+raytracing.elf: main-saturn.o raytracing.o $(LIBGCC)
diff --git a/fp.hpp b/fp.hpp
new file mode 100644
index 0000000..07914d8
--- /dev/null
+++ b/fp.hpp
@@ -0,0 +1,144 @@
+#pragma once
+
+#include <stdint.h>
+
+struct fp_raw_tag {};
+
+template <typename T, typename I, int B>
+struct fp
+{
+  T value;
+
+  constexpr inline fp(T n) noexcept
+    : value(n * (1 << B))
+  {}
+
+  constexpr inline explicit fp(T n, struct fp_raw_tag) noexcept
+    : value(n)
+  {}
+
+  constexpr inline fp<T, I, B> operator-() const noexcept
+  {
+    return fp(-value, fp_raw_tag{});
+  }
+};
+
+template <typename T, typename I, int B>
+constexpr inline fp<T, I, B> operator+(const fp<T, I, B>& a, const fp<T, I, B>& b) noexcept
+{
+  return fp<T, I, B>(a.value + b.value, fp_raw_tag{});
+}
+
+template <typename T, typename I, int B>
+constexpr inline fp<T, I, B> operator-(const fp<T, I, B>& a, const fp<T, I, B>& b) noexcept
+{
+  return fp<T, I, B>(a.value - b.value, fp_raw_tag{});
+}
+
+template <typename T, typename I, int B>
+constexpr inline fp<T, I, B> operator*(const fp<T, I, B>& a, const fp<T, I, B>& b) noexcept
+{
+  I p = (static_cast<I>(a.value) * static_cast<I>(b.value));
+  return fp<T, I, B>(static_cast<T>(p >> B), fp_raw_tag{});
+}
+
+template <typename T, typename I, int B>
+constexpr inline fp<T, I, B> operator*(const fp<T, I, B>& a, T b) noexcept
+{
+  I p = (static_cast<I>(a.value) * static_cast<I>(b));
+  return fp<T, I, B>(static_cast<T>(p), fp_raw_tag{});
+}
+
+template <typename T, typename I, int B>
+constexpr inline fp<T, I, B> operator*(T b, const fp<T, I, B>& a) noexcept
+{
+  I p = (static_cast<I>(a.value) * static_cast<I>(b));
+  return fp<T, I, B>(static_cast<T>(p), fp_raw_tag{});
+}
+
+template <typename T, typename I, int B>
+constexpr inline fp<T, I, B> operator/(const fp<T, I, B>& a, const fp<T, I, B>& b) noexcept
+{
+  I p = (static_cast<I>(a.value) * (static_cast<I>(1) << B)) / static_cast<I>(b.value);
+  return fp<T, I, B>(static_cast<T>(p), fp_raw_tag{});
+}
+
+// comparison
+
+template <typename T, typename I, int B>
+constexpr inline bool operator==(const fp<T, I, B>& a, const fp<T, I, B>& b) noexcept
+{
+    return a.value == b.value;
+}
+
+template <typename T, typename I, int B>
+constexpr inline bool operator!=(const fp<T, I, B>& a, const fp<T, I, B>& b) noexcept
+{
+    return a.value != b.value;
+}
+
+template <typename T, typename I, int B>
+constexpr inline bool operator<(const fp<T, I, B>& a, const fp<T, I, B>& b) noexcept
+{
+    return a.value < b.value;
+}
+
+template <typename T, typename I, int B>
+constexpr inline bool operator>(const fp<T, I, B>& a, const fp<T, I, B>& b) noexcept
+{
+    return a.value > b.value;
+}
+
+template <typename T, typename I, int B>
+constexpr inline bool operator<=(const fp<T, I, B>& a, const fp<T, I, B>& b) noexcept
+{
+    return a.value <= b.value;
+}
+
+template <typename T, typename I, int B>
+constexpr inline bool operator>=(const fp<T, I, B>& a, const fp<T, I, B>& b) noexcept
+{
+    return a.value >= b.value;
+}
+
+// limits
+
+template <typename F>
+struct fp_limits;
+
+template <typename T, typename I, int B>
+struct fp_limits<fp<T, I, B>>
+{
+  static constexpr fp<T, I, B> min() noexcept
+  {
+    return fp<T, I, B>(-(1 << (2 * B - 1)), fp_raw_tag{});
+  }
+
+  static constexpr fp<T, I, B> max() noexcept
+  {
+    return fp<T, I, B>((static_cast<I>(1) << (2 * B - 1)) - 1, fp_raw_tag{});
+  }
+};
+
+// specializations
+
+using fp16_16 = fp<int32_t, int64_t, 16>;
+
+constexpr fp16_16 sqrt(fp16_16 n) noexcept
+{
+  int32_t x = n.value;
+  int32_t c = 0;
+  int32_t d = 1 << 30;
+  while (d > (1 << 6))
+  {
+    int32_t t = c + d;
+    if (x >= t)
+    {
+      x -= t;
+      c = t + d;
+    }
+    x <<= 1;
+    d >>= 1;
+  }
+  return fp16_16(c >> 8, fp_raw_tag{});
+}
diff --git a/main-hosted.cpp b/main-hosted.cpp
new file mode 100644
index 0000000..4d8ffde
--- /dev/null
+++ b/main-hosted.cpp
@@ -0,0 +1,54 @@
+#include <iostream>
+#include <cstdint>
+
+using namespace std;
+
+#include "vec.hpp"
+#include "fp.hpp"
+#include "raytracing.hpp"
+
+typedef vec<3, uint8_t> pixel;
+static pixel frame[canvas::height][canvas::width] = { 0 };
+
+fp16_16 clamp(fp16_16 const& n)
+{
+  return (n > fp16_16(1) ? fp16_16(1) : (n < fp16_16(0) ? fp16_16(0) : n));
+};
+
+uint8_t to_uint8_t(fp16_16 const& v)
+{
+  return static_cast<uint8_t>(v.value >> 16);
+}
+
+void put_pixel(int32_t x, int32_t y, const vec3& color)
+{
+  using namespace canvas;
+
+  int sx = width / 2 + x;
+  int sy = height / 2 - y;
+
+  if (!(sx >= 0 && sx < width && sy >= 0 && sy < height)) {
+    return;
+  }
+
+  vec3 px255 = functor1(clamp, color) * fp16_16(255);
+  frame[sy][sx] = functor1(to_uint8_t, px255);
+}
+
+void render_ppm(ostream& out)
+{
+  using namespace canvas;
+
+  out << "P3 " << width << ' ' << height << " 255\n";
+  for (int sy = 0; sy < height; sy++) {
+    for (int sx = 0; sx < width; sx++) {
+      const pixel& px = frame[sy][sx];
+      out << +px.r << ' ' << +px.g << ' ' << +px.b << '\n';
+    }
+  }
+}
+
+int main()
+{
+  render(put_pixel);
+}
diff --git a/main-saturn.cpp b/main-saturn.cpp
new file mode 100644
index 0000000..5e4295f
--- /dev/null
+++ b/main-saturn.cpp
@@ -0,0 +1,58 @@
+#include "vdp2.h"
+#include "vdp1.h"
+#include "scu.h"
+#include "smpc.h"
+#include "sh2.h"
+
+#include "vec.hpp"
+#include "fp.hpp"
+#include "raytracing.hpp"
+
+fp16_16 clamp(fp16_16 const& n)
+{
+  return (n > fp16_16(1) ? fp16_16(1) : (n < fp16_16(0) ? fp16_16(0) : n));
+};
+
+uint16_t rgb15(const vec3& color)
+{
+  vec3 c = functor1(clamp, color) * fp16_16(255);
+
+  uint8_t red = (c.r.value >> 16) & 0xff;
+  uint8_t green = (c.g.value >> 16) & 0xff;
+  uint8_t blue = (c.b.value >> 16) & 0xff;
+
+  return (blue << 10) | (green << 5) | (red << 0);
+}
+
+void put_pixel(int32_t x, int32_t y, const vec3& color)
+{
+  int sx = 320 / 2 + x;
+  int sy = 240 / 2 - y;
+
+  vdp2.vram.u16[512 * sy + sx] = (1 << 15) | rgb15(color);
+}
+
+void main_asdf()
+{
+  // DISP: Please make sure to change this bit from 0 to 1 during V blank.
+  vdp2.reg.TVMD = ( TVMD__DISP | TVMD__LSMD__NON_INTERLACE
+                  | TVMD__VRESO__240 | TVMD__HRESO__NORMAL_320);
+
+  vdp2.reg.BGON = BGON__N0ON;
+
+  vdp2.reg.CHCTLA = ( CHCTLA__N0CHCN__32K_COLOR     // 15 bits per pixel, RGB
+                    | CHCTLA__N0BMSZ__512x256_DOT
+                    | CHCTLA__N0BMEN__BITMAP_FORMAT
+                    );
+
+  vdp2.reg.MPOFN = MPOFN__N0MP(0);
+
+  render(put_pixel);
+}
+
+extern "C"
+void start(void)
+{
+  main_asdf();
+  while (1) {}
+}
diff --git a/math.hpp b/math.hpp
new file mode 100644
index 0000000..889d296
--- /dev/null
+++ b/math.hpp
@@ -0,0 +1,4 @@
+#pragma once
+
+template <typename T>
+constexpr T sqrt(T n) noexcept;
diff --git a/raytracing.cpp b/raytracing.cpp
new file mode 100644
index 0000000..fe1478f
--- /dev/null
+++ b/raytracing.cpp
@@ -0,0 +1,124 @@
+#include <stdint.h>
+
+#include "vec.hpp"
+#include "fp.hpp"
+#include "raytracing.hpp"
+
+namespace viewport {
+  constexpr int width = 1;
+  constexpr int height = 1;
+}
+
+vec3 canvas_to_viewport(int cx, int cy)
+{
+  return vec3(
+    fp16_16(((cx * viewport::width) * (1 << 16)) >> canvas::bit_width, fp_raw_tag{}),
+    fp16_16(((cy * viewport::height) * (1 << 16)) >> canvas::bit_height, fp_raw_tag{}),
+    fp16_16(1)
+  );
+}
+
+struct sphere {
+  vec3 center;
+  fp16_16 radius;
+  vec3 color;
+};
+
+struct scene {
+  sphere spheres[3];
+};
+
+constexpr scene scene {
+  { // spheres
+    {
+      {0, -1, 3}, // center
+      fp16_16(1), // radius
+      {1, 0, 0},  // color
+    },
+    {
+      {2, 0, 4},
+      fp16_16(1),
+      {0, 0, 1},
+    },
+    {
+      {-2, 0, 4},
+      fp16_16(1),
+      {0, 1, 0},
+    }
+  }
+};
+
+static_assert(scene.spheres[0].center.z.value == (3 << 16));
+
+struct t1_t2 {
+  fp16_16 t1;
+  fp16_16 t2;
+};
+
+t1_t2 intersect_ray_sphere(const vec3& origin, const vec3& direction, const sphere& sphere)
+{
+  fp16_16 r = sphere.radius;
+  vec3 CO = origin - sphere.center;
+
+  auto a = dot(direction, direction);
+  auto b = dot(CO, direction) * static_cast<int32_t>(2);
+  auto c = dot(CO, CO) - r*r;
+
+  auto discriminant = b*b - static_cast<int32_t>(4)*a*c;
+
+  if (discriminant < fp16_16(0)) {
+    return {fp_limits<fp16_16>::max(), fp_limits<fp16_16>::max()};
+  } else {
+    auto sqrt_d = sqrt(discriminant);
+    auto a2 = fp16_16(a*static_cast<int32_t>(2));
+    auto t1 = (-b + sqrt_d) / a2;
+    auto t2 = (-b - sqrt_d) / a2;
+    return {t1, t2};
+  }
+}
+
+static vec3 trace_ray
+(
+  const vec3& origin,
+  const vec3& direction,
+  const fp16_16 t_min,
+  const fp16_16 t_max
+)
+{
+  fp16_16 closest_t = fp_limits<fp16_16>::max();
+  const sphere * closest_sphere = nullptr;
+  for (int i = 0; i < 3; i++) {
+    auto& sphere = scene.spheres[i];
+    auto [t1, t2] = intersect_ray_sphere(origin, direction, sphere);
+    if (t1 >= t_min && t1 < t_max && t1 < closest_t) {
+      closest_t = t1;
+      closest_sphere = &sphere;
+    }
+    if (t2 >= t_min && t2 < t_max && t2 < closest_t) {
+      closest_t = t2;
+      closest_sphere = &sphere;
+    }
+  }
+  if (closest_sphere == nullptr) {
+    return vec3(0, 0, 0);
+  } else {
+    return closest_sphere->color;
+  }
+}
+
+void render(void (&put_pixel) (int32_t x, int32_t y, const vec3& c))
+{
+  using namespace canvas;
+
+  vec3 origin = vec3(0, 0, 0);
+
+  for (int x = -(width/2); x < (width/2); x++) {
+    for (int y = -(height/2 + 1); y < (height/2 + 1); y++) {
+      vec3 direction = canvas_to_viewport(x, y);
+      vec3 color = trace_ray(origin, direction,
+                             fp16_16(1),
+                             fp_limits<fp16_16>::max());
+      put_pixel(x, y, color);
+    }
+  }
+}
diff --git a/raytracing.hpp b/raytracing.hpp
new file mode 100644
index 0000000..b8eeb4e
--- /dev/null
+++ b/raytracing.hpp
@@ -0,0 +1,15 @@
+#pragma once
+
+#include "fp.hpp"
+#include "raytracing.hpp"
+
+using vec3 = vec<3, fp16_16>;
+
+namespace canvas {
+  constexpr int bit_width = 8;
+  constexpr int bit_height = 8;
+  constexpr int width = (1 << bit_width);
+  constexpr int height = (1 << bit_height);
+}
+
+void render(void (&put_pixel) (int32_t x, int32_t y, const vec3& c));
diff --git a/saturn b/saturn
new file mode 120000
index 0000000..73286b8
--- /dev/null
+++ b/saturn
@@ -0,0 +1 @@
+../saturn-c
\ No newline at end of file
diff --git a/vec.hpp b/vec.hpp
new file mode 100644
index 0000000..b43afc8
--- /dev/null
+++ b/vec.hpp
@@ -0,0 +1,147 @@
+#pragma once
+
+#include "math.hpp"
+
+template<int L, typename T>
+struct vec;
+
+//
+// vec3
+//
+
+template<typename T>
+struct vec<3, T>
+{
+  union
+  {
+    struct { T x, y, z; };
+    struct { T r, g, b; };
+    struct { T s, t, p; };
+  };
+
+  inline constexpr vec();
+  inline constexpr vec(T scalar);
+  inline constexpr vec(T _x, T _y, T _z);
+
+  inline constexpr T const& operator[](int i) const;
+
+  template<typename U>
+  inline constexpr vec<3, T>& operator=(vec<3, U> const& v);
+
+  template<typename U>
+  inline constexpr vec<3, T>& operator+=(vec<3, U> const& v);
+
+  template<typename U>
+  inline constexpr vec<3, T>& operator-=(vec<3, U> const& v);
+};
+
+template <typename T>
+inline constexpr vec<3, T>::vec()
+  : x(0), y(0), z(0)
+{}
+
+template <typename T>
+inline constexpr vec<3, T>::vec(T scalar)
+  : x(scalar), y(scalar), z(scalar)
+{}
+
+template <typename T>
+inline constexpr vec<3, T>::vec(T _x, T _y, T _z)
+  : x(_x), y(_y), z(_z)
+{}
+
+template <typename T>
+inline constexpr T const& vec<3, T>::operator[](int i) const
+{
+  switch(i)
+  {
+  default:
+  case 0:
+    return x;
+  case 1:
+    return y;
+  case 2:
+    return z;
+  }
+}
+
+template<typename T>
+template<typename U>
+inline constexpr vec<3, T>& vec<3, T>::operator=(vec<3, U> const& v)
+{
+  this->x = static_cast<T>(v.x);
+  this->y = static_cast<T>(v.y);
+  this->z = static_cast<T>(v.z);
+  return *this;
+}
+
+template<typename T>
+template<typename U>
+inline constexpr vec<3, T>& vec<3, T>::operator+=(vec<3, U> const& v)
+{
+  *this = *this + vec<3, T>(v);
+  return *this;
+}
+
+template<typename T>
+template<typename U>
+inline constexpr vec<3, T>& vec<3, T>::operator-=(vec<3, U> const& v)
+{
+  *this = *this + vec<3, T>(v);
+  return *this;
+}
+
+template <typename T>
+inline constexpr vec<3, T> operator+(vec<3, T> const& v1, vec<3, T> const& v2)
+{
+  return vec<3, T>(v1.x + v2.x,
+                   v1.y + v2.y,
+                   v1.z + v2.z);
+}
+
+template <typename T>
+inline constexpr vec<3, T> operator-(vec<3, T> const& v1, vec<3, T> const& v2)
+{
+  return vec<3, T>(v1.x - v2.x,
+                   v1.y - v2.y,
+                   v1.z - v2.z);
+}
+
+template <typename T>
+inline constexpr vec<3, T> operator*(vec<3, T> const& v1, vec<3, T> const& v2)
+{
+  return vec<3, T>(v1.x * v2.x,
+                   v1.y * v2.y,
+                   v1.z * v2.z);
+}
+
+template <typename T>
+inline constexpr vec<3, T> operator*(vec<3, T> const& v1, T const& scalar)
+{
+  return v1 * vec<3, T>(scalar);
+}
+
+template <typename T>
+inline constexpr T dot(vec<3, T> const& v1, vec<3, T> const& v2)
+{
+  vec<3, T> tmp(v1 * v2);
+  return tmp.x + tmp.y + tmp.z;
+}
+
+template <typename T>
+inline constexpr vec<3, T> functor1(T (&func) (T const& x), vec<3, T> const& v)
+{
+  return vec<3, T>(func(v.x), func(v.y), func(v.z));
+}
+
+template <typename T, typename U>
+inline constexpr vec<3, U> functor1(U (&func) (T const& x), vec<3, T> const& v)
+{
+  return vec<3, U>(func(v.x), func(v.y), func(v.z));
+}
+
+template <typename T>
+inline constexpr T length(vec<3, T> const& v)
+{
+  return sqrt(dot(v, v));
+}