// ===== 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

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

	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
	// 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(&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)

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*
		let image = (*data).image
		mut y = (*data).ystart
		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((*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))
				image.data[y * 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 * 8ul) as thread mut*
	let threads_data = allocate(thread_count * 104ul) 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(64ul) 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()