diff --git a/examples/raytracer.png b/examples/raytracer.png
new file mode 100644
index 0000000..b5ff487
Binary files /dev/null and b/examples/raytracer.png differ
diff --git a/examples/raytracer.psl b/examples/raytracer.psl
new file mode 100644
index 0000000..5d97779
--- /dev/null
+++ b/examples/raytracer.psl
@@ -0,0 +1,829 @@
+// ===== Dynamic memory =====
+
+func allocate(size: u64) -> unit mut*:
+	foreign func malloc(size: u64) -> unit mut *
+	return malloc(size)
+
+func deallocate(ptr: unit*):
+	foreign func free(ptr: unit*)
+	free(ptr)
+
+// ===== Debug IO =====
+
+func floor(x: f32) -> i32:
+	foreign func floorf(x: f32) -> f32
+	return floorf(x) as i32
+
+// No function overloading yet...
+func print_byte(c: u8):
+	foreign func putchar(c: i32) -> i32
+	putchar(c as i32)
+
+func print_str(str: u8*):
+	// Don't use puts() because it adds a newline
+	mut p = str
+	while *p != 0ub:
+		print_byte(*p)
+		p += 1
+
+func print_i32(x: i32):
+	if x < 0:
+		print_byte('-')
+		print_i32(-x)
+		return
+	if x >= 10:
+		print_i32(x / 10)
+	print_byte('0' + (x % 10 as u8))
+
+func print_f32(x: f32):
+	if x < 0.0:
+		print_byte('-')
+		print_f32(-x)
+		return
+	// Print integer part
+	let xfloor = floor(x)
+	print_i32(xfloor)
+	print_byte('.')
+	mut y = x - (xfloor as f32)
+	// Print fractional part
+	mut i = 0
+	while i < 3:
+		y = y * 10.0
+		let yfloor = floor(y)
+		print_byte('0' + (yfloor as u8))
+		y = y - (yfloor as f32)
+		i = i + 1
+
+func print_vec3(v: vec3):
+	print_byte('(')
+	print_f32(v.x)
+	print_byte(',')
+	print_f32(v.y)
+	print_byte(',')
+	print_f32(v.z)
+	print_byte(')')
+
+func print_time(t: f32):
+	if t >= 3600.0:
+		let hours = floor(t / 3600.0)
+		print_i32(hours)
+		print_byte('h')
+		print_byte(' ')
+		let minutes = floor((t - (hours as f32) * 3600.0) / 60.0)
+		print_i32(minutes)
+		print_byte('m')
+	else if t >= 60.0:
+		let minutes = floor(t / 60.0)
+		print_i32(minutes)
+		print_byte('m')
+		print_byte(' ')
+		let seconds = t - (minutes as f32) * 60.0
+		print_f32(seconds)
+		print_byte('s')
+	else if t >= 1.0:
+		print_f32(t)
+		print_byte('s')
+	else if t >= 0.001:
+		print_f32(t * 1000.0)
+		print_byte('m')
+		print_byte('s')
+	else if t >= 0.000001:
+		print_f32(t * 1000000.0)
+		print_byte('u')
+		print_byte('s')
+	else:
+		print_f32(t * 1000000000.0)
+		print_byte('n')
+		print_byte('s')
+
+func flush():
+	foreign func fflush(stream: unit*) -> i32
+	fflush(0ul as unit*)
+
+// ===== File IO =====
+
+// Opaque empty struct for type safety
+struct file:
+
+func open(path: u8*, mode: u8) -> file*:
+	foreign func fopen(path: u8*, mode: u8*) -> file*
+	// A hack to turn a single u8 into a zero-terminated array
+	let mode_wide = mode as u32
+	return fopen(path, &mode as u8*)
+
+func close(file: file*):
+	foreign func fclose(file: file*) -> i32
+	fclose(file)
+
+func write(file: file*, data: u8*, size: u64):
+	foreign func fwrite(data: u8*, size: u64, count: u64, file: file*) -> u64
+	fwrite(data, size, 1ul, file)
+
+func write_byte(file: file*, value: u8):
+	foreign func fputc(ch: i32, file: file*) -> i32
+	fputc(value as i32, file)
+
+func write_u64(file: file*, value: u64):
+	if value >= 10ul:
+		write_u64(file, value / 10ul)
+	write_byte(file, '0' + (value % 10ul as u8))
+
+// ===== Time =====
+
+// Platform-dependent!
+struct timespec:
+	seconds: i64
+	nanoseconds: i64
+
+const TIME_UTC = 1
+
+func get_time() -> timespec:
+	foreign func timespec_get(ts: timespec mut*, base: i32) -> i32
+	mut result = timespec()
+	timespec_get(&mut result, TIME_UTC)
+	return result
+
+func time_delta(x: timespec, y: timespec) -> f32:
+	return ((x.seconds - y.seconds) as f32) + ((x.nanoseconds - y.nanoseconds) as f32) / 1000000000.0
+
+// ===== Image =====
+
+struct image:
+	width: u64
+	height: u64
+	data: u8[3] mut*
+
+func create_image(width: u64, height: u64) -> image:
+	return image(width, height, allocate(width * height * 3ul) as u8[3] mut *)
+
+// ===== PPM format helpers =====
+
+func write_ppm(path: u8*, image: image):
+	let file = open(path, 'w')
+	write_byte(file, 'P')
+	write_byte(file, '6')
+	write_byte(file, '\n')
+	write_u64(file, image.width)
+	write_byte(file, ' ')
+	write_u64(file, image.height)
+	write_byte(file, '\n')
+	write_u64(file, 255ul)
+	write_byte(file, '\n')
+	write(file, image.data as u8*, image.width * image.height * 3ul)
+	close(file)
+
+// ===== Math =====
+
+const pi = 3.141592653589793
+
+func sin(x: f32) -> f32:
+	foreign func sinf(x: f32) -> f32
+	return sinf(x)
+
+func cos(x: f32) -> f32:
+	foreign func cosf(x: f32) -> f32
+	return cosf(x)
+
+func sqrt(x: f32) -> f32:
+	foreign func sqrtf(x: f32) -> f32
+	return sqrtf(x)
+
+func tan(x: f32) -> f32:
+	foreign func tanf(x: f32) -> f32
+	return tanf(x)
+
+func atan(x: f32) -> f32:
+	foreign func atanf(x: f32) -> f32
+	return atanf(x)
+
+func pow(x: f32, y: f32) -> f32:
+	foreign func powf(x: f32, y: f32) -> f32
+	return powf(x, y)
+
+func log(x: f32) -> f32:
+	foreign func logf(x: f32) -> f32
+	return logf(x)
+
+func abs(x: f32) -> f32:
+	if x < 0.0:
+		return -x
+	return x
+
+func min(x: f32, y: f32) -> f32:
+	if x < y:
+		return x
+	return y
+
+func max(x: f32, y: f32) -> f32:
+	if x > y:
+		return x
+	return y
+
+func clamp(x: f32) -> f32:
+	return max(0.0, min(1.0, x))
+
+struct vec3:
+	x: f32
+	y: f32
+	z: f32
+
+// No operator overloading yet...
+func add(v: vec3, u: vec3) -> vec3:
+	return vec3(v.x + u.x, v.y + u.y, v.z + u.z)
+
+func sub(v: vec3, u: vec3) -> vec3:
+	return vec3(v.x - u.x, v.y - u.y, v.z - u.z)
+
+// Multiply scalar
+// No function overloading yet...
+func mults(v: vec3, s: f32) -> vec3:
+	return vec3(v.x * s, v.y * s, v.z * s)
+
+// Multiply vector (component-wise)
+func multv(v: vec3, u: vec3) -> vec3:
+	return vec3(v.x * u.x, v.y * u.y, v.z * u.z)
+
+// Divide vector (component-wise)
+func divv(v: vec3, u: vec3) -> vec3:
+	return vec3(v.x / u.x, v.y / u.y, v.z / u.z)
+
+func lerp(v: vec3, u: vec3, t: f32) -> vec3:
+	return add(v, mults(sub(u, v), t))
+
+func dot(v: vec3, u: vec3) -> f32:
+	return v.x * u.x + v.y * u.y + v.z * u.z
+
+func cross(v: vec3, u: vec3) -> vec3:
+	return vec3(v.y * u.z - v.z * u.y, v.z * u.x - v.x * u.z, v.x * u.y - v.y * u.x)
+
+func length(v: vec3) -> f32:
+	return sqrt(dot(v, v))
+
+func normalized(v: vec3) -> vec3:
+	return mults(v, 1.0 / length(v))
+
+struct ray:
+	origin: vec3
+	direction: vec3
+
+struct quaternion:
+	x: f32
+	y: f32
+	z: f32
+	w: f32
+
+const identity = quaternion(0.0, 0.0, 0.0, 1.0)
+
+func inverse(q: quaternion) -> quaternion:
+	return quaternion(-q.x, -q.y, -q.z, q.w)
+
+func rotation(axis: vec3, angle: f32) -> quaternion:
+	let s = sin(angle * 0.5)
+	let c = cos(angle * 0.5)
+	return quaternion(axis.x * s, axis.y * s, axis.z * s, c)
+
+func rotate(q: quaternion, v: vec3) -> vec3:
+	let qv = vec3(q.x, q.y, q.z)
+	let t = mults(cross(qv, v), 2.0)
+	return add(add(v, mults(t, q.w)), cross(qv, t))
+
+// ===== Random number generation =====
+
+struct rng:
+	state: u64[2]
+
+func splitmix64(x: u64) -> u64:
+	mut result = x
+	result += 11400714819323198485ul
+	result = (result ^ (result >> 30u)) * 13787848793156543929ul
+	result = (result ^ (result >> 27u)) * 10723151780598845931ul
+	return result ^ (result >> 31u)
+
+func next_u64(rng: rng mut*) -> u64:
+	func rotate_left(x: u64, k: u32) -> u64:
+		return (x << k) | (x >> (64u - k))
+	// No struct-field-by-pointer (rng->state) access yet
+	let result = rotate_left((*rng).state[0] * 5ul, 7u) * 9ul
+	(*rng).state[1] ^= (*rng).state[0]
+	(*rng).state[0] = rotate_left((*rng).state[0], 24u) ^ (*rng).state[1] ^ ((*rng).state[1] << 16u)
+	(*rng).state[1] = rotate_left((*rng).state[1], 37u)
+	return result
+
+// Uniform in [0..1)
+func next_f32(rng: rng mut*) -> f32:
+	return (next_u64(rng) as f32) / 18446744073709551615.0
+
+func next_normal(rng: rng mut*) -> f32:
+	let u = next_f32(rng)
+	let v = next_f32(rng)
+	return sqrt(- 2.0 * log(u)) * cos(2.0 * pi * v)
+
+// Uniform on unit sphere
+func next_vec3(rng: rng mut*) -> vec3:
+	let theta = 2.0 * pi * next_f32(rng)
+	let cosphi = 2.0 * next_f32(rng) - 1.0
+	let sinphi = sqrt(max(0.0, 1.0 - cosphi * cosphi))
+	return vec3(cos(theta) * sinphi, sin(theta) * sinphi, cosphi)
+
+func next_vec3_normal(rng: rng mut*) -> vec3:
+	return vec3(next_normal(rng), next_normal(rng), next_normal(rng))
+
+// ===== Camera =====
+
+// Default (non-rotated) camera uses standard OpenGL conventions:
+//    X right
+//    Y up
+//   -Z forward
+struct camera:
+	position: vec3
+	rotation: quaternion
+	fovx: f32
+	fovy: f32
+
+func compute_fovx(fovy: f32, aspect_ratio: f32) -> f32:
+	return 2.0 * atan(tan(fovy * 0.5) * aspect_ratio)
+
+// x, y are in [-1..1]
+func camera_direction(camera: camera, x: f32, y: f32) -> vec3:
+	let dir = vec3(x * tan(camera.fovx * 0.5), y * tan(camera.fovy * 0.5), -1.0)
+	return rotate(camera.rotation, normalized(dir))
+
+func camera_ray(camera: camera, x: f32, y: f32) -> ray:
+	return ray(camera.position, camera_direction(camera, x, y))
+
+// ===== Scene description =====
+
+const diffuse_tag = 1u
+const metallic_tag = 2u
+const glass_tag = 3u
+
+struct material:
+	// albedo for diffuse, reflection color for metallic
+	color: vec3
+	emission: vec3
+	type: u32
+	roughness: f32
+	ior: f32
+
+struct plane:
+	// Plane normal is defined as vec3(0,0,1) rotated by object rotation
+	// Plane origin is defined as object position
+
+struct sphere:
+	radius: f32
+
+struct box:
+	// NB: box size is 2 * extent
+	// I.e. box boundary is center +/- extent
+	extent: vec3
+
+const plane_tag = 1u
+const sphere_tag = 2u
+const box_tag = 3u
+
+// Poor man's union
+struct shape:
+	tag: u32
+	data: f32[3]
+
+struct object:
+	position: vec3
+	rotation: quaternion
+	shape: shape
+	material: material
+
+func set_plane(object: object mut*, shape: plane):
+	(*object).shape.tag = plane_tag
+
+func set_sphere(object: object mut*, shape: sphere):
+	(*object).shape.tag = sphere_tag
+	*((*object).shape.data as f32 mut* as sphere mut*) = shape
+
+func set_box(object: object mut*, shape: box):
+	(*object).shape.tag = box_tag
+	*((*object).shape.data as f32 mut* as box mut*) = shape
+
+struct object_array:
+	size: u64
+	data: object mut*
+
+struct scene:
+	background: vec3
+	objects: object_array
+
+// ===== Color utilities =====
+
+func gamma_correct(v: vec3) -> vec3:
+	const gamma = 1.0 / 2.2
+	return vec3(pow(v.x, gamma), pow(v.y, gamma), pow(v.z, gamma))
+
+func aces(x: vec3) -> vec3:
+	let a = 2.51
+	let b = 0.03
+	let c = 2.43
+	let d = 0.59
+	let e = 0.14
+	// This is why operator overloading is a must
+	let nom = multv(x, add(mults(x, a), vec3(b, b, b)))
+	let den = add(multv(x, add(mults(x, c), vec3(d, d, d))), vec3(e, e, e))
+	return divv(nom, den)
+
+func to_bytes(v: vec3) -> u8[3]:
+	func to_byte(x: f32) -> u8:
+		return clamp(x) * 255.0 as u8
+	return [to_byte(v.x), to_byte(v.y), to_byte(v.z)]
+
+// ===== Raytracing core =====
+
+struct intersection:
+	intersected: bool
+	distance: f32
+	normal: vec3
+	material: material*
+
+const max_distance = 1000000.0
+const no_intersection = intersection(false, max_distance, vec3(0.0, 0.0, 1.0), 0ul as material*)
+
+func intersect_plane(ray: ray, object: object*) -> intersection:
+	let normal = rotate((*object).rotation, vec3(0.0, 0.0, 1.0))
+	// dot(o + t * d - p, n) = 0
+	// dot(o - p, n) + t * dot(d, n) = 0
+	// t = - dot(o - p, n) / dot(d, n)
+	let t = - dot(sub(ray.origin, (*object).position), normal) / dot(ray.direction, normal)
+	return intersection(t > 0.0, t, normal, &(*object).material)
+
+func intersect_sphere(ray: ray, object: object*) -> intersection:
+	let shape = &(*object).shape.data as sphere*
+	// |o + t * d - p|^2 = r^2
+	// dot(o - p, o - p) + 2 * t * dot(o - p, d) + t^2 * dot(d, d) = r * r
+	let delta = sub(ray.origin, (*object).position)
+	// Solve quadratic
+	let a = 1.0 // assume ray.direction is normalized
+	let b = 2.0 * dot(delta, ray.direction)
+	let c = dot(delta, delta) - (*shape).radius * (*shape).radius
+	let D = b * b - 4.0 * a * c
+	if D < 0.0:
+		return no_intersection
+	let t1 = (- b - sqrt(D)) / (2.0 * a)
+	let t2 = (- b + sqrt(D)) / (2.0 * a)
+	if t2 < 0.0:
+		return no_intersection
+	mut t = t1
+	if t1 < 0.0:
+		t = t2
+	let normal = normalized(add(delta, mults(ray.direction, t)))
+	return intersection(true, t, normal, &(*object).material)
+
+func intersect_box(ray: ray, object: object*) -> intersection:
+	func swap(x: f32 mut*, y: f32 mut*):
+		let temp = *x
+		*x = *y
+		*y = temp
+
+	let shape = &(*object).shape.data as box*
+	let inverse_rotation = inverse((*object).rotation)
+	let local_delta = rotate(inverse_rotation, sub(ray.origin, (*object).position))
+	let local_direction = rotate(inverse_rotation, ray.direction)
+
+	// (o + t * d).x = +/- e.x
+	mut txmin = (- (*shape).extent.x - local_delta.x) / local_direction.x
+	mut txmax = (  (*shape).extent.x - local_delta.x) / local_direction.x
+	mut tymin = (- (*shape).extent.y - local_delta.y) / local_direction.y
+	mut tymax = (  (*shape).extent.y - local_delta.y) / local_direction.y
+	mut tzmin = (- (*shape).extent.z - local_delta.z) / local_direction.z
+	mut tzmax = (  (*shape).extent.z - local_delta.z) / local_direction.z
+
+	if txmin > txmax:
+		swap(&txmin, &txmax)
+	if tymin > tymax:
+		swap(&tymin, &tymax)
+	if tzmin > tzmax:
+		swap(&tzmin, &tzmax)
+
+	let tmin = max(txmin, max(tymin, tzmin))
+	let tmax = min(txmax, min(tymax, tzmax))
+
+	if tmin > tmax || tmax < 0.0:
+		return no_intersection
+
+	let inside = tmin < 0.0
+	mut t = tmin
+	mut tt = vec3(txmin, tymin, tzmin)
+	// No ternary if yet...
+	if inside:
+		t = tmax
+		tt = vec3(txmax, tymax, tzmax)
+
+	mut normal = vec3(0.0, 0.0, 0.0)
+
+	if t == tt.x:
+		if local_direction.x > 0.0:
+			normal = vec3(-1.0, 0.0, 0.0)
+		else:
+			normal = vec3( 1.0, 0.0, 0.0)
+	else if t == tt.y:
+		if local_direction.y > 0.0:
+			normal = vec3(0.0, -1.0, 0.0)
+		else:
+			normal = vec3(0.0,  1.0, 0.0)
+	else:
+		if local_direction.z > 0.0:
+			normal = vec3(0.0, 0.0, -1.0)
+		else:
+			normal = vec3(0.0, 0.0,  1.0)
+		
+	if inside:
+		normal = mults(normal, -1.0)
+
+	return intersection(true, tmin, rotate((*object).rotation, normal), &(*object).material)
+
+func intersect_scene(scene: scene*, ray: ray) -> intersection:
+	mut intersection = no_intersection
+	mut i = 0ul
+	while i < (*scene).objects.size:
+		mut current_intersection = no_intersection
+
+		if (*scene).objects.data[i].shape.tag == plane_tag:
+			current_intersection = intersect_plane(ray, &(*scene).objects.data[i])
+		else if (*scene).objects.data[i].shape.tag == sphere_tag:
+			current_intersection = intersect_sphere(ray, &(*scene).objects.data[i])
+		else if (*scene).objects.data[i].shape.tag == box_tag:
+			current_intersection = intersect_box(ray, &(*scene).objects.data[i])
+
+		if current_intersection.intersected && current_intersection.distance < intersection.distance:
+			intersection = current_intersection
+		
+		i += 1ul
+
+	return intersection
+
+func raytrace(scene: scene*, camera_ray: ray, rng: rng mut*) -> vec3:
+	mut current_ray = camera_ray
+	mut result = vec3(0.0, 0.0, 0.0)
+	mut factor = vec3(1.0, 1.0, 1.0)
+	let termination_probability = 0.25
+
+	mut running = true
+	while running:
+		let intersection = intersect_scene(scene, current_ray)
+
+		if intersection.intersected:
+			// Uncomment to debug normals
+			// return mults(add(intersection.normal, vec3(1.0, 1.0, 1.0)), 0.5)
+			
+			let cosine = - dot(intersection.normal, current_ray.direction)
+			let inside = cosine < 0.0
+
+			result = add(result, multv(factor, (*intersection.material).emission))
+
+			// Russian roulette ray termination
+			if next_f32(rng) < termination_probability:
+				// No break operator yet...
+				running = false
+			else:
+				mut new_direction = vec3(0.0, 0.0, 0.0)
+
+				if (*intersection.material).type == diffuse_tag:
+					// NB: albedo is assumed to be premultiplied by pi to be in [0..1] range
+					// This should also contain multiplication by cos(new_dir, normal), division by direction pdf (cos / pi)
+					// and division by pi (because of albedo normalization), but these all cancel out
+					factor = multv(factor, (*intersection.material).color)
+
+					// Cosine-weighted hemisphere direction
+					new_direction = normalized(add(next_vec3(rng), intersection.normal))
+
+				else if (*intersection.material).type == metallic_tag:
+					// This should also contain multiplication by brdf and division by direction pdf,
+					// but we'll just pretend that the random reflected ray pdf coincides with brdf and thus cancels out
+					factor = multv(factor, (*intersection.material).color)
+
+					// Compute perfect-mirror reflected direction
+					new_direction = add(current_ray.direction, mults(intersection.normal, 2.0 * cosine))
+
+					// Alter the direction based on roughness
+					new_direction = normalized(add(new_direction, mults(next_vec3_normal(rng), cosine * (*intersection.material).roughness)))
+
+				else if (*intersection.material).type == glass_tag:
+					// This should also contain multiplication by brdf and division by direction pdf,
+					// but we'll just pretend that the random refracted ray pdf coincides with brdf and thus cancels out
+					factor = multv(factor, (*intersection.material).color)
+
+					mut ior = (*intersection.material).ior
+					if inside:
+						ior = 1.0 / ior
+					
+					// Compute perfect refracted ray
+					let k = 1.0 - ior * ior * (1.0 - cosine * cosine)
+					if k >= 0.0:
+						new_direction = add(mults(current_ray.direction, ior), mults(intersection.normal, ior * abs(cosine) - sqrt(k)))
+					else:
+						// Total internal reflection
+						new_direction = add(current_ray.direction, mults(intersection.normal, 2.0 * cosine))
+
+					// Alter the direction based on roughness
+					new_direction = normalized(add(new_direction, mults(next_vec3_normal(rng), (*intersection.material).roughness)))
+
+				// Compute the new ray, and offset its origin a bit along intersection normal
+				let position = add(current_ray.origin, mults(current_ray.direction, intersection.distance))
+				mut position_shift = 0.001
+				if dot(new_direction, intersection.normal) < 0.0:
+					position_shift = -position_shift
+				current_ray = ray(add(position, mults(intersection.normal, position_shift)), new_direction)
+
+				// Account for Russian roulette
+				factor = mults(factor, 1.0 / (1.0 - termination_probability))
+		else:
+			result = add(result, multv(factor, (*scene).background))
+			running = false
+
+	return result
+
+// ===== Main =====
+
+func make_default_scene() -> scene:
+	let object_count = 8ul
+	// No sizeof yet...
+	let objects = allocate(object_count * 20ul * 4ul) as object mut*
+
+	// Filling all objects in the same function overflows the max stack size
+	// dictated by arm64 limitations (add immediate is bound by 12 bits => 4Kb)
+	// The reasonable solution: fix the compiler!
+	// The temporary solution: split into several separate functions!
+	
+	func fill_walls(objects: object mut*):
+		// Floor
+		objects[0].position = vec3(0.0, -5.0, 0.0)
+		objects[0].rotation = rotation(vec3(1.0, 0.0, 0.0), - pi * 0.5)
+		objects[0].material.color = vec3(0.6, 0.6, 0.6)
+		objects[0].material.emission = vec3(0.0, 0.0, 0.0)
+		objects[0].material.type = diffuse_tag
+		set_plane(&objects[0], plane())
+
+		// Ceiling
+		objects[1].position = vec3(0.0, 5.0, 0.0)
+		objects[1].rotation = rotation(vec3(1.0, 0.0, 0.0), pi * 0.5)
+		objects[1].material.color = vec3(0.6, 0.6, 0.6)
+		objects[1].material.emission = vec3(0.0, 0.0, 0.0)
+		objects[1].material.type = diffuse_tag
+		set_plane(&objects[1], plane())
+
+		// Left wall
+		objects[2].position = vec3(-5.0, 0.0, 0.0)
+		objects[2].rotation = rotation(vec3(0.0, 1.0, 0.0), pi * 0.5)
+		objects[2].material.color = vec3(1.0, 0.2, 0.2)
+		objects[2].material.emission = vec3(0.0, 0.0, 0.0)
+		objects[2].material.type = diffuse_tag
+		set_plane(&objects[2], plane())
+
+		// Right wall
+		objects[3].position = vec3(5.0, 0.0, 0.0)
+		objects[3].rotation = rotation(vec3(0.0, 1.0, 0.0), - pi * 0.5)
+		objects[3].material.color = vec3(0.2, 0.2, 1.0)
+		objects[3].material.emission = vec3(0.0, 0.0, 0.0)
+		objects[3].material.type = diffuse_tag
+		set_plane(&objects[3], plane())
+
+		// Back wall
+		objects[4].position = vec3(0.0, 0.0, -5.0)
+		objects[4].rotation = identity
+		objects[4].material.color = vec3(0.6, 0.6, 0.6)
+		objects[4].material.emission = vec3(0.0, 0.0, 0.0)
+		objects[4].material.type = diffuse_tag
+		set_plane(&objects[4], plane())
+
+	func fill_objects(objects: object mut*):
+		// Top lamp
+		objects[5].position = vec3(0.0, 10.0, 0.0)
+		objects[5].rotation = identity
+		objects[5].material.color = vec3(0.0, 0.0, 0.0)
+		objects[5].material.emission = vec3(20.0, 16.0, 12.0)
+		objects[5].material.type = diffuse_tag
+		set_sphere(&objects[5], sphere(5.25))
+
+		// Box
+		objects[6].position = vec3(-2.25, -2.0, -1.0)
+		objects[6].rotation = rotation(vec3(0.0, 1.0, 0.0), pi / 6.0)
+		objects[6].material.color = vec3(1.0, 1.0, 1.0)
+		objects[6].material.emission = vec3(0.0, 0.0, 0.0)
+		objects[6].material.type = glass_tag
+		objects[6].material.roughness = 0.0
+		objects[6].material.ior = 1.333
+		set_box(&objects[6], box(vec3(1.5, 3.0, 1.5)))
+
+		// Sphere
+		objects[7].position = vec3(2.5, -3.0, -1.0)
+		objects[7].rotation = rotation(vec3(0.0, 1.0, 0.0), - pi / 6.0)
+		objects[7].material.color = vec3(1.0, 1.0, 1.0)
+		objects[7].material.emission = vec3(0.0, 0.0, 0.0)
+		objects[7].material.type = metallic_tag
+		objects[7].material.roughness = 0.5
+		set_sphere(&objects[7], sphere(2.0))
+	
+	fill_walls(objects)
+	fill_objects(objects)
+
+	let background = vec3(0.0, 0.0, 0.0)
+	return scene(background, object_array(object_count, objects))
+
+const default_fovy = 2.0 * atan(0.5)
+
+func main():
+	let scene = make_default_scene()
+
+	let image = create_image(512ul, 512ul)
+	// let image = create_image(128ul, 128ul)
+
+	let aspect_ratio = (image.width as f32) / (image.height as f32)
+	let camera = camera(vec3(0.0, 0.0, 15.0), identity, default_fovy, compute_fovx(default_fovy, aspect_ratio))
+
+	// let samples_per_pixel = 4096ul
+	let samples_per_pixel = 2048ul
+	// let samples_per_pixel = 1024ul
+	// let samples_per_pixel = 512ul
+	// let samples_per_pixel = 256ul
+	// let samples_per_pixel = 16ul
+	// let samples_per_pixel = 1ul
+
+	let start_time = get_time()
+	mut last_report_time = start_time
+	mut last_report_progress = 0.0
+	mut average_speed = 0.0
+
+	func clear_line(spaces: u32):
+		mut i = 0u
+		print_byte('\r')
+		while i < spaces:
+			print_byte(' ')
+			i += 1u
+		print_byte('\r')
+
+	mut y = 0ul
+	while y < image.height:
+		mut x = 0ul
+		while x < image.width:
+			mut rng = rng([splitmix64(x | (y << 32u)), 10723151780598845931ul])
+			mut sample = 0ul
+			mut color = vec3(0.0, 0.0, 0.0)
+			while sample < samples_per_pixel:
+				let tx = - 1.0 + 2.0 * (x as f32 + next_f32(&mut rng)) / (image.width  as f32)
+				let ty =   1.0 - 2.0 * (y as f32 + next_f32(&mut rng)) / (image.height as f32)
+				let ray = camera_ray(camera, tx, ty)
+				color = add(color, raytrace(&scene, ray, &mut rng))
+				sample += 1ul
+			color = mults(color, 1.0 / (samples_per_pixel as f32))
+			image.data[y * image.width + x] = to_bytes(gamma_correct(aces(color)))
+			
+			x += 1ul
+
+			let time = get_time()
+			let delta = time_delta(time, last_report_time)
+			if delta > 0.125 || (y == 0ul && x == 1ul) || (y + 1ul == image.height && x == image.width):
+				let done = y * image.width + x
+				let total = image.width * image.height
+				let progress = (done as f32) / (total as f32)
+				let time_passed = time_delta(time, start_time)
+
+				// Running exponential-weighted average speed for
+				// accurate remaining time estimate
+				let speed_estimate = (progress - last_report_progress) / delta
+				if average_speed == 0.0:
+					average_speed = speed_estimate
+				else:
+					average_speed = average_speed + (speed_estimate - average_speed) * 0.125
+				let time_remaining = (1.0 - progress) / average_speed
+
+				let str1 = ['%', ' ', 'd', 'o', 'n', 'e', ',', ' ', '\0']
+				let str2 = [' ', 'l', 'e', 'f', 't', ' ', '\0']
+
+				clear_line(30u)
+				print_f32(100.0 * progress)
+				print_str(str1 as u8*)
+				print_time(time_remaining)
+				print_str(str2 as u8*)
+				flush()
+				last_report_time = time
+				last_report_progress = progress
+
+		y += 1ul
+
+	let total_time = time_delta(get_time(), start_time)
+	let total_samples = image.width * image.height * samples_per_pixel
+	
+	let str1 = ['R', 'e', 'n', 'd', 'e', 'r', ' ', 't', 'o', 'o', 'k', ' ', '\0']
+	let str2 = [' ', 'p', 'e', 'r', ' ', 's', 'a', 'm', 'p', 'l', 'e', '\0']
+	clear_line(30u)
+	print_str(str1 as u8*)
+	print_time(total_time)
+	print_byte('\n')
+	print_time(total_time / (total_samples as f32))
+	print_str(str2 as u8*)
+	print_byte('\n')
+
+	// No string literal support yet...
+	let path = ['o', 'u', 't', '.', 'p', 'p', 'm', '\0']
+	write_ppm(path as u8*, image)
+	deallocate(image.data as unit*)
+
+	deallocate(scene.objects.data as unit*)
+
+main()
\ No newline at end of file