diff --git a/common.mk b/common.mk
index 8bbee2b..caf702e 100644
--- a/common.mk
+++ b/common.mk
@@ -2,7 +2,7 @@ MAKEFILE_PATH := $(abspath $(lastword $(MAKEFILE_LIST)))
 DIR := $(dir $(MAKEFILE_PATH))
 
 LIB ?= .
-OPT ?= -O3
+OPT ?= -O1
 DEBUG ?= -g -gdwarf-4
 GENERATED ?=
 
diff --git a/example/cube.cpp b/example/cube.cpp
new file mode 100644
index 0000000..05b4b3f
--- /dev/null
+++ b/example/cube.cpp
@@ -0,0 +1,208 @@
+#include <cstdint>
+
+#include "align.hpp"
+
+#include "vga.hpp"
+#include "holly.hpp"
+#include "holly/core.hpp"
+#include "holly/core_bits.hpp"
+#include "holly/ta_fifo_polygon_converter.hpp"
+#include "holly/ta_parameter.hpp"
+#include "holly/ta_bits.hpp"
+#include "holly/region_array.hpp"
+#include "holly/background.hpp"
+#include "holly/texture_memory_alloc.hpp"
+#include "memorymap.hpp"
+#include "serial.hpp"
+
+#include "geometry/cube.hpp"
+#include "math/vec4.hpp"
+
+using vec4 = vec<4, float>;
+
+constexpr float half_degree = 0.01745329f / 2;
+
+vec3 rotate(const vec3& vertex, float theta)
+{
+  float x = vertex.x;
+  float y = vertex.y;
+  float z = vertex.z;
+  float t;
+
+  t  = y * cos(theta) - z * sin(theta);
+  z  = y * sin(theta) + z * cos(theta);
+  y  = t;
+
+  float theta2 = 3.14 * sin(theta / 2);
+
+  t  = x * cos(theta2) - z * sin(theta2);
+  z  = x * sin(theta2) + z * cos(theta2);
+  x  = t;
+
+  return vec3(x, y, z);
+}
+
+void transform(ta_parameter_writer& parameter,
+               const uint32_t face_ix,
+               const float theta,
+               const vec3 lights[2])
+{
+  const uint32_t parameter_control_word = para_control::para_type::polygon_or_modifier_volume
+                                        | para_control::list_type::opaque
+                                        | obj_control::col_type::floating_color
+                                        | obj_control::gouraud;
+
+  const uint32_t isp_tsp_instruction_word = isp_tsp_instruction_word::depth_compare_mode::greater
+                                          | isp_tsp_instruction_word::culling_mode::cull_if_positive;
+
+  const uint32_t tsp_instruction_word = tsp_instruction_word::src_alpha_instr::one
+                                      | tsp_instruction_word::dst_alpha_instr::zero
+                                      | tsp_instruction_word::fog_control::no_fog;
+
+  parameter.append<global_polygon_type_0>() = global_polygon_type_0(parameter_control_word,
+                                                                    isp_tsp_instruction_word,
+                                                                    tsp_instruction_word,
+                                                                    0);
+  auto& face = cube::faces[face_ix];
+
+  constexpr uint32_t strip_length = 3;
+  for (uint32_t i = 0; i < strip_length; i++) {
+    bool end_of_strip = i == strip_length - 1;
+
+    // world transform
+    uint32_t vertex_ix = face[i].vertex - 1;
+    auto& vertex = cube::vertices[vertex_ix];
+    auto point = rotate(vertex, theta);
+
+    // lighting transform
+    uint32_t normal_ix = face[i].normal - 1;
+    auto& normal = cube::normals[normal_ix];
+    auto n = rotate(normal, theta);
+
+    vec4 color = {0.3, 0.3, 0.3, 1.0};
+
+    // intensity calculation
+    {
+      auto l = lights[0] - point;
+      auto n_dot_l = dot(n, l);
+      if (n_dot_l > 0) {
+        color.x += 0.35 * n_dot_l / (length(n) * length(l));
+      }
+    }
+
+    {
+      auto l = lights[1] - point;
+      auto n_dot_l = dot(n, l);
+      if (n_dot_l > 0) {
+        color.y += 0.35 * n_dot_l / (length(n) * length(l));
+      }
+    }
+
+    float x = point.x;
+    float y = point.y;
+    float z = point.z;
+
+    // camera transform
+    z += 3;
+
+    // perspective
+    x = x / z;
+    y = y / z;
+
+    // screen space transform
+    x *= 240.f;
+    y *= 240.f;
+    x += 320.f;
+    y += 240.f;
+    z = 1 / z;
+
+    parameter.append<vertex_polygon_type_1>() =
+      vertex_polygon_type_1(x, y, z,
+			    color.w, // alpha
+			    color.x, // r
+			    color.y, // g
+			    color.z, // b
+			    end_of_strip);
+  }
+}
+
+void init_texture_memory(const struct opb_size& opb_size)
+{
+  auto mem = reinterpret_cast<volatile texture_memory_alloc *>(texture_memory32);
+
+  background_parameter(mem->background, 0xff220000);
+
+  region_array2(mem->region_array,
+	        (offsetof (struct texture_memory_alloc, object_list)),
+		640 / 32, // width
+		480 / 32, // height
+		opb_size
+		);
+}
+
+uint32_t _ta_parameter_buf[((32 * (5 * 6 + 1)) + 32) / 4];
+
+void main()
+{
+  vga();
+
+  // The address of `ta_parameter_buf` must be a multiple of 32 bytes.
+  // This is mandatory for ch2-dma to the ta fifo polygon converter.
+  uint32_t * ta_parameter_buf = align_32byte(_ta_parameter_buf);
+
+  constexpr uint32_t ta_alloc = ta_alloc_ctrl::pt_opb::no_list
+			      | ta_alloc_ctrl::tm_opb::no_list
+			      | ta_alloc_ctrl::t_opb::no_list
+			      | ta_alloc_ctrl::om_opb::no_list
+                              | ta_alloc_ctrl::o_opb::_16x4byte;
+
+  constexpr struct opb_size opb_size = { .opaque = 16 * 4
+				       , .opaque_modifier = 0
+				       , .translucent = 0
+				       , .translucent_modifier = 0
+				       , .punch_through = 0
+				       };
+
+  constexpr uint32_t tiles = (640 / 32) * (320 / 32);
+
+  holly.SOFTRESET = softreset::pipeline_soft_reset
+		  | softreset::ta_soft_reset;
+  holly.SOFTRESET = 0;
+
+  core_init();
+  init_texture_memory(opb_size);
+
+  uint32_t frame_ix = 0;
+  constexpr uint32_t num_frames = 1;
+
+  float theta = 0;
+  vec3 lights[2] = {
+    {0.f, 0.f, -4.f},
+    {0.f, 0.f, 0.f},
+  };
+
+  while (1) {
+    ta_polygon_converter_init(opb_size.total() * tiles, ta_alloc,
+                              640, 480);
+
+    //lights[0].x = cos(theta) * 10;
+    //lights[0].z = sin(theta) * 10;
+
+    //lights[1].x = cos(theta + half_degree * 90.f) * 10;
+    //lights[1].z = sin(theta + half_degree * 90.f) * 10;
+
+    auto parameter = ta_parameter_writer(ta_parameter_buf);
+    for (uint32_t i = 0; i < 12; i++) {
+      transform(parameter, i, theta, lights);
+    }
+    parameter.append<global_end_of_list>() = global_end_of_list();
+    ta_polygon_converter_transfer(ta_parameter_buf, parameter.offset);
+    ta_wait_opaque_list();
+    core_start_render(frame_ix, num_frames);
+
+    v_sync_in();
+    core_wait_end_of_render_video(frame_ix, num_frames);
+    theta += half_degree;
+    frame_ix += 1;
+  }
+}
diff --git a/example/example.mk b/example/example.mk
index 948292b..6231139 100644
--- a/example/example.mk
+++ b/example/example.mk
@@ -96,6 +96,17 @@ TRANSLUCENCY_OBJ = \
 example/translucency.elf: LDSCRIPT = $(LIB)/alt.lds
 example/translucency.elf: $(START_OBJ) $(TRANSLUCENCY_OBJ)
 
+CUBE_OBJ = \
+	example/cube.o \
+	vga.o \
+	holly/core.o \
+	holly/region_array.o \
+	holly/background.o \
+	holly/ta_fifo_polygon_converter.o
+
+example/cube.elf: LDSCRIPT = $(LIB)/alt.lds
+example/cube.elf: $(START_OBJ) $(CUBE_OBJ)
+
 MACAW_CUBE_OBJ = \
 	example/macaw_cube.o \
 	vga.o \
diff --git a/example/macaw_cube.cpp b/example/macaw_cube.cpp
index 7781a5d..06954b8 100644
--- a/example/macaw_cube.cpp
+++ b/example/macaw_cube.cpp
@@ -54,7 +54,7 @@ vertex cube_faces[][4] = {
 };
 constexpr uint32_t num_faces = (sizeof (cube_faces)) / (sizeof (cube_faces[0]));
 
-constexpr uint32_t color = 0xffff00ff;
+constexpr uint32_t color = 0xff0000ff;
 
 static float theta = 0;
 constexpr float half_degree = 0.01745329f / 2.f;
@@ -64,7 +64,7 @@ void transform(ta_parameter_writer& parameter,
 		   const uint32_t strip_length)
 {
   uint32_t texture_address = (offsetof (struct texture_memory_alloc, texture));
-  parameter.append<global_polygon_type_0>() = global_polygon_type_0(texture_address);
+  parameter.append<global_polygon_type_0>() = global_polygon_type_0();
 
   for (uint32_t i = 0; i < strip_length; i++) {
     bool end_of_strip = i == strip_length - 1;
@@ -98,10 +98,10 @@ void transform(ta_parameter_writer& parameter,
 
     z = 1 / z;
 
-    parameter.append<vertex_polygon_type_3>() =
-      vertex_polygon_type_3(x, y, z,
-			    strip_vertices[i].u,
-			    strip_vertices[i].v,
+    parameter.append<vertex_polygon_type_0>() =
+      vertex_polygon_type_0(x, y, z,
+			    //strip_vertices[i].u,
+			    //strip_vertices[i].v,
 			    color,
 			    end_of_strip);
   }
@@ -111,7 +111,7 @@ void init_texture_memory(const struct opb_size& opb_size)
 {
   auto mem = reinterpret_cast<volatile texture_memory_alloc *>(texture_memory32);
 
-  background_parameter(mem->background);
+  background_parameter(mem->background, 0xffff00ff);
 
   region_array2(mem->region_array,
 	        (offsetof (struct texture_memory_alloc, object_list)),
@@ -173,7 +173,7 @@ void main()
   constexpr uint32_t num_frames = 1;
 
   while (1) {
-    ta_polygon_converter_init(opb_size.total() * tiles, ta_alloc);
+    ta_polygon_converter_init(opb_size.total() * tiles, ta_alloc, 640, 480);
     auto parameter = ta_parameter_writer(ta_parameter_buf);
     for (uint32_t i = 0; i < num_faces; i++) {
       transform(parameter, cube_faces[i], 4);
diff --git a/geometry/cube.hpp b/geometry/cube.hpp
new file mode 100644
index 0000000..adb8c62
--- /dev/null
+++ b/geometry/cube.hpp
@@ -0,0 +1,54 @@
+#pragma once
+
+#include <cstdint>
+
+#include "math/vec3.hpp"
+
+using vec3 = vec<3, float>;
+
+struct vertex__normal {
+  uint8_t vertex;
+  uint8_t normal;
+};
+
+using face = vertex__normal[3];
+
+namespace cube {
+
+  constexpr vec3 vertices[] = {
+    { 1.f,  1.f, -1.f},
+    { 1.f, -1.f, -1.f},
+    { 1.f,  1.f,  1.f},
+    { 1.f, -1.f,  1.f},
+    {-1.f,  1.f, -1.f},
+    {-1.f, -1.f, -1.f},
+    {-1.f,  1.f,  1.f},
+    {-1.f, -1.f,  1.f},
+  };
+
+  constexpr vec3 normals[] = {
+    {-0.f,  1.f, -0.f},
+    {-0.f, -0.f,  1.f},
+    {-1.f, -0.f, -0.f},
+    {-0.f, -1.f, -0.f},
+    { 1.f, -0.f, -0.f},
+    {-0.f, -0.f, -1.f},
+  };
+
+  constexpr face faces[] = {
+    {{5, 1}, {3, 1}, {1, 1}},
+    {{3, 2}, {8, 2}, {4, 2}},
+    {{7, 3}, {6, 3}, {8, 3}},
+    {{2, 4}, {8, 4}, {6, 4}},
+    {{1, 5}, {4, 5}, {2, 5}},
+    {{5, 6}, {2, 6}, {6, 6}},
+    {{5, 1}, {7, 1}, {3, 1}},
+    {{3, 2}, {7, 2}, {8, 2}},
+    {{7, 3}, {5, 3}, {6, 3}},
+    {{2, 4}, {4, 4}, {8, 4}},
+    {{1, 5}, {3, 5}, {4, 5}},
+    {{5, 6}, {1, 6}, {2, 6}},
+  };
+
+  constexpr int num_faces = (sizeof (faces)) / (sizeof (face));
+}
diff --git a/holly/ta_parameter.hpp b/holly/ta_parameter.hpp
index be739ca..6548e20 100644
--- a/holly/ta_parameter.hpp
+++ b/holly/ta_parameter.hpp
@@ -90,6 +90,67 @@ struct vertex_polygon_type_0 {
     { }
 };
 
+struct vertex_polygon_type_1 {
+  uint32_t parameter_control_word;
+  float x;
+  float y;
+  float z;
+  float base_color_alpha;
+  float base_color_r;
+  float base_color_g;
+  float base_color_b;
+
+  vertex_polygon_type_1(const float x,
+                        const float y,
+                        const float z,
+                        const float base_color_alpha,
+                        const float base_color_r,
+                        const float base_color_g,
+                        const float base_color_b,
+                        const bool end_of_strip
+                        )
+    : parameter_control_word( para_control::para_type::vertex_parameter
+			    | (end_of_strip ? para_control::end_of_strip : 0)
+                            )
+    , x(x)
+    , y(y)
+    , z(z)
+    , base_color_alpha(base_color_alpha)
+    , base_color_r(base_color_r)
+    , base_color_g(base_color_g)
+    , base_color_b(base_color_b)
+    { }
+};
+
+struct vertex_polygon_type_2 {
+  uint32_t parameter_control_word;
+  float x;
+  float y;
+  float z;
+  uint32_t _res0;
+  uint32_t _res1;
+  float base_intensity;
+  uint32_t _res2;
+
+  vertex_polygon_type_2(const float x,
+                        const float y,
+                        const float z,
+                        const float base_intensity,
+                        const bool end_of_strip
+                        )
+    : parameter_control_word( para_control::para_type::vertex_parameter
+			    | (end_of_strip ? para_control::end_of_strip : 0)
+                            )
+    , x(x)
+    , y(y)
+    , z(z)
+    , _res0(0)
+    , _res1(0)
+    , base_intensity(base_intensity)
+    , _res2(0)
+    { }
+};
+
 struct vertex_polygon_type_3 {
   uint32_t parameter_control_word;
   float x;
@@ -171,17 +232,31 @@ struct global_polygon_type_0 {
   uint32_t data_size_for_sort_dma;
   uint32_t next_address_for_sort_dma;
 
+  global_polygon_type_0(const uint32_t parameter_control_word,
+                        const uint32_t isp_tsp_instruction_word,
+                        const uint32_t tsp_instruction_word,
+                        const uint32_t texture_control_word)
+    : parameter_control_word(parameter_control_word)
+    , isp_tsp_instruction_word(isp_tsp_instruction_word)
+    , tsp_instruction_word(tsp_instruction_word)
+    , texture_control_word(texture_control_word)
+    , _res0(0)
+    , _res1(0)
+    , data_size_for_sort_dma(0)
+    , next_address_for_sort_dma(0)
+  { }
+
   // untextured
   global_polygon_type_0()
     : parameter_control_word( para_control::para_type::polygon_or_modifier_volume
 			    | para_control::list_type::opaque
                             | obj_control::col_type::packed_color )
 
-    , isp_tsp_instruction_word( isp_tsp_instruction_word::depth_compare_mode::always
+    , isp_tsp_instruction_word( isp_tsp_instruction_word::depth_compare_mode::greater
 			      | isp_tsp_instruction_word::culling_mode::no_culling )
 
-    , tsp_instruction_word( tsp_instruction_word::src_alpha_instr::src_alpha
-			  | tsp_instruction_word::dst_alpha_instr::inverse_src_alpha
+    , tsp_instruction_word( tsp_instruction_word::src_alpha_instr::one
+			  | tsp_instruction_word::dst_alpha_instr::zero
 			  | tsp_instruction_word::fog_control::no_fog )
 
     , texture_control_word( 0 )
@@ -234,6 +309,35 @@ static_assert((offsetof (struct global_polygon_type_0, _res1))
 static_assert((offsetof (struct global_polygon_type_0, data_size_for_sort_dma))    == 0x18);
 static_assert((offsetof (struct global_polygon_type_0, next_address_for_sort_dma)) == 0x1c);
 
+struct global_polygon_type_1 {
+  uint32_t parameter_control_word;
+  uint32_t isp_tsp_instruction_word;
+  uint32_t tsp_instruction_word;
+  uint32_t texture_control_word;
+  float face_color_alpha;
+  float face_color_r;
+  float face_color_g;
+  float face_color_b;
+
+  global_polygon_type_1(const uint32_t parameter_control_word,
+                        const uint32_t isp_tsp_instruction_word,
+                        const uint32_t tsp_instruction_word,
+                        const uint32_t texture_control_word,
+                        const float face_color_alpha,
+                        const float face_color_r,
+                        const float face_color_g,
+                        const float face_color_b)
+    : parameter_control_word(parameter_control_word)
+    , isp_tsp_instruction_word(isp_tsp_instruction_word)
+    , tsp_instruction_word(tsp_instruction_word)
+    , texture_control_word(texture_control_word)
+    , face_color_alpha(face_color_alpha)
+    , face_color_r(face_color_r)
+    , face_color_g(face_color_g)
+    , face_color_b(face_color_b)
+  { }
+};
+
 struct global_sprite {
   uint32_t parameter_control_word;
   uint32_t isp_tsp_instruction_word;
diff --git a/math/div.hpp b/math/div.hpp
new file mode 100644
index 0000000..f0d4e55
--- /dev/null
+++ b/math/div.hpp
@@ -0,0 +1,120 @@
+#pragma once
+
+#include <stdint.h>
+
+#ifndef USE_SH2_DVSR
+inline constexpr uint32_t
+__udiv32(uint32_t n, uint32_t d)
+{
+  uint32_t q = 0;
+  uint32_t r = 0;
+
+  for (int i = 31; i >= 0; --i) {
+    q = q << 1;
+    r = r << 1;
+
+    r |= (n >> 31) & 1;
+    n = n << 1;
+
+    if (d <= r) {
+      r = r - d;
+      q = q | 1;
+    }
+  }
+
+  return q;
+}
+
+inline constexpr uint32_t
+__udiv64_32(uint64_t n, uint32_t base)
+{
+  uint64_t rem = n;
+  uint64_t b = base;
+  uint64_t res = 0, d = 1;
+  uint32_t high = rem >> 32;
+
+  if (high >= base) {
+    high = __udiv32(high, base);
+    res = (uint64_t)high << 32;
+    rem -= (uint64_t)(high*base) << 32;
+  }
+
+  while ((int64_t)b > 0 && b < rem) {
+    b = b+b;
+    d = d+d;
+  }
+
+  do {
+    if (rem >= b) {
+      rem -= b;
+      res += d;
+    }
+    b >>= 1;
+    d >>= 1;
+  } while (d);
+
+  return res;
+}
+#else
+#include "sh2.h"
+inline uint32_t
+__udiv64_32(uint64_t n, uint32_t d)
+{
+  sh2.reg.DVSR = d;
+  sh2.reg.DVDNTH = (uint32_t)(n >> 32);
+  sh2.reg.DVDNTL = (uint32_t)(n);
+
+  // 39 cycles
+  asm volatile ("nop");
+  asm volatile ("nop");
+  asm volatile ("nop");
+  asm volatile ("nop");
+  asm volatile ("nop");
+  asm volatile ("nop");
+  asm volatile ("nop");
+  asm volatile ("nop");
+  asm volatile ("nop");
+  asm volatile ("nop");
+  asm volatile ("nop");
+  asm volatile ("nop");
+  asm volatile ("nop");
+  asm volatile ("nop");
+  asm volatile ("nop");
+  asm volatile ("nop");
+  asm volatile ("nop");
+  asm volatile ("nop");
+  asm volatile ("nop");
+  asm volatile ("nop");
+  asm volatile ("nop");
+  asm volatile ("nop");
+  asm volatile ("nop");
+  asm volatile ("nop");
+  asm volatile ("nop");
+  asm volatile ("nop");
+  asm volatile ("nop");
+  asm volatile ("nop");
+  asm volatile ("nop");
+  asm volatile ("nop");
+  asm volatile ("nop");
+  asm volatile ("nop");
+  asm volatile ("nop");
+  asm volatile ("nop");
+  asm volatile ("nop");
+  asm volatile ("nop");
+  asm volatile ("nop");
+  asm volatile ("nop");
+  asm volatile ("nop");
+
+  return sh2.reg.DVDNTL;
+}
+#endif
+
+inline int32_t
+__div64_32(int64_t n, int32_t d)
+{
+  uint64_t n_abs = n >= 0 ? (uint64_t)n : -(uint64_t)n;
+  uint32_t d_abs = d >= 0 ? (uint32_t)d : -(uint32_t)d;
+  uint32_t q_abs = __udiv64_32(n_abs, d_abs);
+
+  return (n < 0) == (d < 0) ? (int32_t)q_abs : -(int32_t)q_abs;
+}
diff --git a/math/mat4x4.hpp b/math/mat4x4.hpp
new file mode 100644
index 0000000..1a5c25e
--- /dev/null
+++ b/math/mat4x4.hpp
@@ -0,0 +1,124 @@
+#pragma once
+
+#include <utility>
+
+template <int R, int C, typename T>
+struct mat;
+
+//
+// mat4x4
+//
+
+template <typename T>
+struct mat<4, 4, T>
+{
+  typedef vec<4, T> row_type;
+  typedef vec<4, T> col_type;
+
+private:
+  row_type value[4];
+
+public:
+  inline constexpr mat();
+
+  inline constexpr mat
+  (
+    T const& a00, T const& a01, T const& a02, T const& a03,
+    T const& a10, T const& a11, T const& a12, T const& a13,
+    T const& a20, T const& a21, T const& a22, T const& a23,
+    T const& a30, T const& a31, T const& a32, T const& a33
+  );
+
+  inline static constexpr int length() { return 4; }
+
+  inline constexpr typename mat<4, 4, T>::row_type const &
+  operator[](int i) const;
+
+  void operator=(const mat<4, 4, T>&) = delete;
+
+};
+
+
+template<typename T>
+inline constexpr mat<4, 4, T>::mat()
+  : value{std::move(row_type(1, 0, 0, 0)),
+	  std::move(row_type(0, 1, 0, 0)),
+          std::move(row_type(0, 0, 1, 0)),
+          std::move(row_type(0, 0, 0, 1))}
+{ }
+
+template<typename T>
+inline constexpr mat<4, 4, T>::mat
+(
+  T const& a00, T const& a01, T const& a02, T const& a03,
+  T const& a10, T const& a11, T const& a12, T const& a13,
+  T const& a20, T const& a21, T const& a22, T const& a23,
+  T const& a30, T const& a31, T const& a32, T const& a33
+)
+  : value{std::move(row_type(a00, a01, a02, a03)),
+	  std::move(row_type(a10, a11, a12, a13)),
+          std::move(row_type(a20, a21, a22, a23)),
+          std::move(row_type(a30, a31, a32, a33))}
+{ }
+
+template <typename T>
+inline constexpr typename mat<4, 4, T>::row_type const &
+mat<4, 4, T>::operator[](int i) const
+{
+  switch (i)
+  {
+  default: [[fallthrough]];
+  case 0:
+    return value[0];
+  case 1:
+    return value[1];
+  case 2:
+    return value[2];
+  case 3:
+    return value[3];
+  }
+}
+
+template<typename T>
+inline constexpr mat<4, 4, T> operator*(mat<4, 4, T> const& m1, mat<4, 4, T> const& m2)
+{
+#define c(i, j) (                               \
+    m1[i][0] * m2[0][j]                         \
+  + m1[i][1] * m2[1][j]                         \
+  + m1[i][2] * m2[2][j]                         \
+  + m1[i][3] * m2[3][j] )
+
+  return mat<4, 4, T>(c(0,0), c(0,1), c(0,2), c(0,3),
+                      c(1,0), c(1,1), c(1,2), c(1,3),
+                      c(2,0), c(2,1), c(2,2), c(2,3),
+                      c(3,0), c(3,1), c(3,2), c(3,3));
+#undef c
+}
+
+template<typename T>
+inline constexpr typename mat<4, 4, T>::row_type operator*
+(
+  mat<4, 4, T> const& m,
+  typename mat<4, 4, T>::col_type const& v
+)
+{
+#define c(i) (                                  \
+    m[i][0] * v[0]                              \
+  + m[i][1] * v[1]                              \
+  + m[i][2] * v[2]                              \
+  + m[i][3] * v[3] )
+
+  return typename mat<4, 4, T>::row_type(c(0), c(1), c(2), c(3));
+#undef c
+}
+
+template<typename T>
+inline constexpr mat<4, 4, T> transpose(mat<4, 4, T> const& m)
+{
+  return mat<4, 4, T>(
+    m[0][0], m[1][0], m[2][0], m[3][0],
+    m[0][1], m[1][1], m[2][1], m[3][1],
+    m[0][2], m[1][2], m[2][2], m[3][2],
+    m[0][3], m[1][3], m[2][3], m[3][3]
+  );
+}
diff --git a/math/math.hpp b/math/math.hpp
new file mode 100644
index 0000000..4e4af69
--- /dev/null
+++ b/math/math.hpp
@@ -0,0 +1,28 @@
+#pragma once
+
+template <typename T>
+constexpr T sqrt(const T n) noexcept;
+
+template <>
+constexpr float sqrt(const float n) noexcept
+{
+  return __builtin_sqrtf(n);
+}
+
+template <typename T>
+constexpr T cos(const T n) noexcept;
+
+template <>
+constexpr float cos(const float n) noexcept
+{
+  return __builtin_cosf(n);
+}
+
+template <typename T>
+constexpr T sin(const T n) noexcept;
+
+template <>
+constexpr float sin(const float n) noexcept
+{
+  return __builtin_sinf(n);
+}
diff --git a/math/vec.hpp b/math/vec.hpp
new file mode 100644
index 0000000..a63d225
--- /dev/null
+++ b/math/vec.hpp
@@ -0,0 +1,4 @@
+#pragma once
+
+template <int L, typename T>
+struct vec;
diff --git a/math/vec3.hpp b/math/vec3.hpp
new file mode 100644
index 0000000..7383cee
--- /dev/null
+++ b/math/vec3.hpp
@@ -0,0 +1,149 @@
+#pragma once
+
+#include "math.hpp"
+#include "vec.hpp"
+
+//
+// vec3
+//
+
+template <typename T>
+struct vec<3, T>
+{
+  T x, y, z;
+
+  inline constexpr vec();
+  inline constexpr vec(T scalar);
+  inline constexpr vec(T _x, T _y, T _z);
+
+  constexpr inline vec<3, T> operator-() const;
+  inline constexpr T const& operator[](int i) const;
+  inline constexpr vec<3, T>& operator=(vec<3, T> const& v);
+  inline constexpr vec<3, T>& operator+=(vec<3, T> const& v);
+  inline constexpr vec<3, T>& operator-=(vec<3, T> 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>
+constexpr inline vec<3, T> vec<3, T>::operator-() const
+{
+  return vec<3, T>(-x, -y, -z);
+}
+
+template <typename T>
+inline constexpr T const& vec<3, T>::operator[](int i) const
+{
+  switch(i)
+  {
+  default: [[fallthrough]];
+  case 0: return x;
+  case 1: return y;
+  case 2: return z;
+  }
+}
+
+template <typename T>
+inline constexpr vec<3, T>& vec<3, T>::operator=(vec<3, T> 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>
+inline constexpr vec<3, T>& vec<3, T>::operator+=(vec<3, T> const& v)
+{
+  *this = *this + vec<3, T>(v);
+  return *this;
+}
+
+template <typename T>
+inline constexpr vec<3, T>& vec<3, T>::operator-=(vec<3, T> 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 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));
+}
diff --git a/math/vec4.hpp b/math/vec4.hpp
new file mode 100644
index 0000000..2ec8b45
--- /dev/null
+++ b/math/vec4.hpp
@@ -0,0 +1,155 @@
+#pragma once
+
+#include "math.hpp"
+#include "vec.hpp"
+
+//
+// vec4
+//
+
+template <typename T>
+struct vec<4, T>
+{
+  T x, y, z, w;
+
+  inline constexpr vec();
+  inline constexpr vec(T scalar);
+  inline constexpr vec(T _x, T _y, T _z, T _w);
+
+  constexpr inline vec<4, T> operator-() const;
+  inline constexpr T const& operator[](int i) const;
+  inline constexpr vec<4, T>& operator=(vec<4, T> const& v);
+  inline constexpr vec<4, T>& operator+=(vec<4, T> const& v);
+  inline constexpr vec<4, T>& operator-=(vec<4, T> const& v);
+};
+
+template <typename T>
+inline constexpr vec<4, T>::vec()
+  : x(0), y(0), z(0), w(0)
+{}
+
+template <typename T>
+inline constexpr vec<4, T>::vec(T scalar)
+  : x(scalar), y(scalar), z(scalar), w(scalar)
+{}
+
+template <typename T>
+inline constexpr vec<4, T>::vec(T _x, T _y, T _z, T _w)
+  : x(_x), y(_y), z(_z), w(_w)
+{}
+
+template <typename T>
+constexpr inline vec<4, T> vec<4, T>::operator-() const
+{
+  return vec<4, T>(-x, -y, -z, -w);
+}
+
+template <typename T>
+inline constexpr T const& vec<4, T>::operator[](int i) const
+{
+  switch(i)
+  {
+  default: [[fallthrough]];
+  case 0: return x;
+  case 1: return y;
+  case 2: return z;
+  case 3: return w;
+  }
+}
+
+template <typename T>
+inline constexpr vec<4, T>& vec<4, T>::operator=(vec<4, T> const& v)
+{
+  this->x = static_cast<T>(v.x);
+  this->y = static_cast<T>(v.y);
+  this->z = static_cast<T>(v.z);
+  this->w = static_cast<T>(v.w);
+  return *this;
+}
+
+template <typename T>
+inline constexpr vec<4, T>& vec<4, T>::operator+=(vec<4, T> const& v)
+{
+  *this = *this + vec<4, T>(v);
+  return *this;
+}
+
+template <typename T>
+inline constexpr vec<4, T>& vec<4, T>::operator-=(vec<4, T> const& v)
+{
+  *this = *this - vec<4, T>(v);
+  return *this;
+}
+
+template <typename T>
+inline constexpr vec<4, T> operator+(vec<4, T> const& v1, vec<4, T> const& v2)
+{
+  return vec<4, T>(v1.x + v2.x,
+                   v1.y + v2.y,
+                   v1.z + v2.z,
+		   v1.w + v2.w);
+}
+
+template <typename T>
+inline constexpr vec<4, T> operator-(vec<4, T> const& v1, vec<4, T> const& v2)
+{
+  return vec<4, T>(v1.x - v2.x,
+                   v1.y - v2.y,
+                   v1.z - v2.z,
+		   v1.w - v2.w);
+}
+
+template <typename T>
+inline constexpr vec<4, T> operator*(vec<4, T> const& v1, vec<4, T> const& v2)
+{
+  return vec<4, T>(v1.x * v2.x,
+                   v1.y * v2.y,
+                   v1.z * v2.z,
+		   v1.w * v2.w);
+}
+
+template <typename T>
+inline constexpr vec<4, T> operator*(vec<4, T> const& v1, T const& scalar)
+{
+  return v1 * vec<4, T>(scalar);
+}
+
+template <typename T>
+inline constexpr vec<4, T> operator/(vec<4, T> const& v1, vec<4, T> const& v2)
+{
+  return vec<4, T>(v1.x / v2.x,
+                   v1.y / v2.y,
+                   v1.z / v2.z,
+		   v1.w / v2.w);
+}
+
+template <typename T>
+inline constexpr vec<4, T> operator/(vec<4, T> const& v1, T const& scalar)
+{
+  return v1 / vec<4, T>(scalar);
+}
+
+template <typename T>
+inline constexpr T dot(vec<4, T> const& v1, vec<4, T> const& v2)
+{
+  vec<4, T> tmp(v1 * v2);
+  return tmp.x + tmp.y + tmp.z + tmp.w;
+}
+
+template <typename T>
+inline constexpr vec<4, T> functor1(T (&func) (T const& x), vec<4, T> const& v)
+{
+  return vec<4, T>(func(v.x), func(v.y), func(v.z), func(v.w));
+}
+
+template <typename T, typename U>
+inline constexpr vec<4, U> functor1(U (&func) (T const& x), vec<4, T> const& v)
+{
+  return vec<4, U>(func(v.x), func(v.y), func(v.z), func(v.w));
+}
+
+template <typename T>
+inline constexpr T length(vec<4, T> const& v)
+{
+  return sqrt(dot(v, v));
+}