// ===== 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(mut str: u8*) {
	// Don't use puts() because it adds a newline
	while *str != 0ub {
		print_byte(*str)
		str += 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(mut x: f32) {
	if x < 0.0 {
		print_byte('-')
		print_f32(-x)
		return
	}
	// Print integer part
	let xfloor = floor(x)
	print_i32(xfloor)
	print_byte('.')
	x -= (xfloor as f32)
	// Print fractional part
	mut i = 0
	while i < 3 {
		x = x * 10.0
		let xfloor = floor(x)
		print_byte('0' + (xfloor as u8))
		x -= (xfloor 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, 0ub]
	return fopen(path, &mode_wide 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
}

func sleep(time: timespec) {
	foreign func nanosleep(time: timespec*, rem: unit*) -> i32
	nanosleep(&time, 0ul as unit*)
}

// ===== 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 {
	return if x >= 0.0 then x else -x
}

func sqr(x: f32) -> f32 {
	return x * x
}

func min(x: f32, y: f32) -> f32 {
	return if x <= y then x else y
}

func max(x: f32, y: f32) -> f32 {
	return if x >= y then x else y
}

// I want function overloading so badly
func min_u32(x: u32, y: u32) -> u32 {
	return if x <= y then x else 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))
	}
	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
	}
	let t = if t1 < 0.0 then t2 else t1
	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 sort(x: f32 mut*, y: f32 mut*) {
		if *x > *y {
			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

	sort(&mut txmin, &mut txmax)
	sort(&mut tymin, &mut tymax)
	sort(&mut tzmin, &mut 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
	let t = if inside then tmax else tmin
	let tt = if inside then vec3(txmax, tymax, tzmax) else vec3(txmin, tymin, tzmin)

	mut normal = vec3()

	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, t, 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

	while true {
		let intersection = intersect_scene(scene, current_ray)

		if !intersection.intersected {
			result = add(result, multv(factor, scene.background))
			break
		}

		// 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 {
			break
		}

		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 {
			mut ior = intersection.material.ior
			if inside {
				ior = 1.0 / ior
			}

			// Schlick's approximation for Fresnel term
			let r0 = sqr((1.0 - ior) / (1.0 + ior))
			let reflectance = r0 + (1.0 - r0) * pow(max(0.0, 1.0 - abs(cosine)), 5.0)

			if next_f32(rng) < reflectance {
				// Compute perfect-mirror reflected direction
				new_direction = add(current_ray.direction, mults(intersection.normal, 2.0 * cosine))
			} else {
				// 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)

				// 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))
	}

	return result
}

// ===== Multithreading =====

// Opaque handle
struct thread {
	// Platform-dependent!
	handle: unit*
}

func create_thread(thread_func: unit mut* -> unit, thread_data: unit mut*) -> thread {
	foreign func pthread_create(thread: thread mut*, attr: unit*, start_routine: unit mut* -> unit, arg: unit mut*) -> i32
	mut result = thread()
	pthread_create(&mut result, 0ul as unit*, thread_func, thread_data)
	return result
}

func join_thread(thread: thread) {
	foreign func pthread_join(thread: thread, retval: unit*) -> i32
	pthread_join(thread, 0ul as unit*)
}

struct mutex {
	// Platform-dependent!
	payload: u64[8]
}

func create_mutex(mutex: mutex mut*) {
	foreign func pthread_mutex_init(mutex: mutex mut*, attr: unit*) -> i32
	pthread_mutex_init(mutex, 0ul as unit*)
}

func destroy_mutex(mutex: mutex mut*) {
	foreign func pthread_mutex_destroy(mutex: mutex mut *) -> i32
	pthread_mutex_destroy(mutex)
}

func lock_mutex(mutex: mutex mut*) {
	foreign func pthread_mutex_lock(mutex: mutex mut*) -> i32
	pthread_mutex_lock(mutex)
}

func unlock_mutex(mutex: mutex mut*) {
	foreign func pthread_mutex_unlock(mutex: mutex mut*) -> i32
	pthread_mutex_unlock(mutex)
}

// ===== Main =====

func make_default_scene() -> scene {
	let object_count = 8ul
	let objects = allocate(object_count * sizeof(object)) 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(&mut 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(&mut 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(&mut 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(&mut 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(&mut 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(&mut 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(0.7, 0.85, 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(&mut 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, 0.8, 0.6)
		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(&mut 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)

// Measured performance (apple m1 pro, 128x128, 1024 spp):
// 1 thread:  133s (~8us    per sample)
// 4 threads:  35s (~8.25us per sample)
// 8 threads:  18s (~8.5us  per sample)

func main() {
	let scene = make_default_scene()

	// let image = create_image(512ul, 512ul)
	// let image = create_image(256ul, 256ul)
	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))

	// const samples_per_pixel = 16384ul
	// const samples_per_pixel = 4096ul
	// const samples_per_pixel = 2048ul
	const samples_per_pixel = 1024ul
	// const samples_per_pixel = 512ul
	// const samples_per_pixel = 256ul
	// const samples_per_pixel = 16ul
	// const samples_per_pixel = 1ul

	struct thread_data {
		scene: scene*
		image: image
		camera: camera
		ystart: u64
		ystep: u64
		done_mutex: mutex mut*
		done: u64
	}

	func thread_func(data_raw: unit mut*) {
		let data = data_raw as thread_data mut*
		mut y = data.ystart
		while y < data.image.height {
			mut x = 0ul
			while x < data.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)) / (data.image.width  as f32)
					let ty =   1.0 - 2.0 * (y as f32 + next_f32(&mut rng)) / (data.image.height as f32)
					let ray = camera_ray((*data).camera, tx, ty)
					color = add(color, raytrace((*data).scene, ray, &mut rng))
					sample += 1ul
				}
				color = mults(color, 1.0 / (samples_per_pixel as f32))
				data.image.data[y * data.image.width + x] = to_bytes(gamma_correct(aces(color)))
				x += 1ul

				lock_mutex(data.done_mutex)
				data.done += 1ul
				unlock_mutex(data.done_mutex)
			}
			y += data.ystep
		}
	}

	let start_time = get_time()
	mut last_report_time = start_time
	mut last_report_progress = 0.0
	mut average_speed = 0.0
	mut report_count = 0u

	let thread_count = 9ul
	let threads = allocate(thread_count * sizeof(thread)) as thread mut*
	let threads_data = allocate(thread_count * sizeof(thread_data)) as thread_data mut*

	mut th = 0ul
	while th < thread_count {
		let data = threads_data + th
		data.scene = &scene
		data.image = image
		data.camera = camera
		data.ystart = th
		data.ystep = thread_count
		data.done_mutex = allocate(sizeof(mutex)) as mutex mut*
		create_mutex(data.done_mutex)
		data.done = 0ul
		threads[th] = create_thread(thread_func, data as unit mut*)
		th += 1ul
	}

	func clear_line(spaces: u32) {
		mut i = 0u
		print_byte('\r')
		while i < spaces {
			print_byte(' ')
			i += 1u
		}
		print_byte('\r')
	}

	mut done = 0ul
	let total = image.width * image.height
	while done < total {
		sleep(timespec(0l, 125000000l))
		let time = get_time()
		let delta = time_delta(time, last_report_time)
		if !(delta > 0.125 || last_report_progress == 0.0) {
			continue
		}

		done = 0ul

		mut th = 0ul
		while th < thread_count {
			lock_mutex(threads_data[th].done_mutex)
			done += threads_data[th].done
			unlock_mutex(threads_data[th].done_mutex)
			th += 1ul
		}

		if done == 0ul {
			continue
		}

		let progress = (done as f32) / (total as f32)

		// Running exponential-weighted average speed
		// for +/- accurate remaining time estimate
		let speed_estimate = (progress - last_report_progress) / delta
		average_speed = average_speed + (speed_estimate - average_speed) / (min_u32(report_count + 1u, 16u) as f32)
		let time_remaining = (1.0 - progress) / average_speed

		let str1 = ['%', ' ', 'd', 'o', 'n', 'e', ',', ' ', '\0']
		let str2 = [' ', 'p', 'a', 's', 's', 'e', 'd', ',', ' ', '\0']
		let str3 = [' ', 'l', 'e', 'f', 't', ' ', '\0']

		clear_line(50u)
		print_f32(100.0 * progress)
		print_str(str1 as u8*)
		print_time(time_delta(time, start_time))
		print_str(str2 as u8*)
		print_time(time_remaining)
		print_str(str3 as u8*)
		flush()
		last_report_time = time
		last_report_progress = progress
		report_count += 1u
	}

	// No for loops yet...
	th = 0ul
	while th < thread_count {
		join_thread(threads[th])
		th += 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(50u)
	print_str(str1 as u8*)
	print_time(total_time)
	print_byte('\n')
	print_time((thread_count as f32) * 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()