val delta = Math.sqrt (Real.nextAfter(1.0, 2.0) - 1.0)
val infinity = Real.posInf
type vec = real * real * real
infix 7 *| fun s *| (x, y, z) : vec = (s*x, s*y, s*z)
infix 6 +| fun (x1, y1, z1) +| (x2, y2, z2) : vec = (x1+x2, y1+y2, z1+z2)
infix 6 -| fun (x1, y1, z1) -| (x2, y2, z2) : vec = (x1-x2, y1-y2, z1-z2)
fun dot((x1, y1, z1), (x2, y2, z2)) : real = x1*x2 + y1*y2 + z1*z2
fun unitise r = (1.0 / Real.Math.sqrt (dot(r, r))) *| r
datatype scene = Sphere of vec * real | Group of vec * real * scene list

fun length r = Math.sqrt(dot(r, r))

fun ray_sphere(orig, dir, center, radius) =
    let val v = center -| orig val b = dot(v, dir)
	val disc = b * b - dot(v, v) + radius * radius in
	if disc < 0.0 then infinity else
	let val disc = Real.Math.sqrt disc
	    val t2 = b + disc in
	    if t2 < 0.0 then infinity else
	    let val t1 = b - disc in
		if t1 > 0.0 then t1 else t2 end
	end
    end

fun intersect(orig, dir, scene) =
    let fun aux(Sphere (center, radius), first as (l, _)) =
	    let val l' = ray_sphere(orig, dir, center, radius) in
		if l' >= l then first else
		(l', unitise (orig +| l' *| dir -| center))
	    end
	  | aux(Group (center, radius, scenes), first as (l, _)) =
	    if ray_sphere(orig, dir, center, radius) >= l then first else
	    foldl aux first scenes
    in aux(scene, (infinity, (0.0, 0.0, 0.0))) end

fun sintersect(orig, dir, scene) =
    let fun aux (Sphere (center, radius)) =
	    ray_sphere(orig, dir, center, radius) < infinity
	  | aux (Group (center, radius, scenes)) =
	    ray_sphere(orig, dir, center, radius) < infinity andalso
	    List.exists aux scenes
    in aux scene end

val neg_light = unitise (1.0, 3.0, ~2.0)

fun ray_trace(dir, scene) =
    let val (lambda, n) = intersect((0.0, 0.0, 0.0), dir, scene) in
	if lambda >= infinity then 0.0 else
	let val g = dot(n, neg_light) in
	    if g <= 0.0 then 0.0 else
	    let val p = lambda *| dir +| delta *| n in
		if sintersect(p, neg_light, scene) then 0.0 else g
	    end
	end
    end

fun bound (Sphere (c', r'), (c, r)) = (c, Real.max(r, (length (c -| c') + r')))
  | bound (Group (_, _, l), (c, r)) = foldl bound (c, r) l

fun create(level, r, (x, y, z)) =
    let val obj = Sphere ((x, y, z), r) in
	if level = 1 then obj else
	let val r' = 3.0 * r / Real.Math.sqrt 12.0
	    fun aux x' z' = create(level-1, 0.5 * r, (x-x', y+r', z+z'))
	    val objs = [obj, aux (~r') (~r'), aux r' (~r'),
			aux (~r') r', aux r' r']
	    val (c, r) = foldl bound ((x, y+r, z), 0.0) objs in
	    Group (c, r, objs)
	end
    end
fun loop(f, accu, (x, y, n)) =
    if y=n then accu else
    if x=n then loop(f, accu, (0, y+1, n)) else
    loop(f, (f(accu, real x, real y)), (x+1, y, n))
fun eye_ray(scene, n, ss, x, y) (g, dx, dy) =
    let val dir = unitise(x + dx/ss, y + dy/ss, n) in
	g + ray_trace(dir, scene)
    end
fun pixel(scene, n, ss) ((), x, y) =
    let val (x, y) = (x - n / 2.0, (n - 1.0) / 2.0 - y)
	val g = loop((eye_ray(scene, n, real ss, x, y)), 0.0, (0, 0, ss)) in
	print (String.str(Char.chr(Real.trunc (255.0 * g / real (ss*ss)))))
    end
fun aux(level, n) =
    let val ss = 4 and scene = create(level, 1.0, (0.0, ~1.0, 4.0)) in
	(fn s => print ("P5\n"^s^" "^s^"\n255\n")) (Int.toString n);
	loop((pixel(scene, real n, ss)), (), (0, 0, n));
	OS.Process.success
    end
fun main (_, [level, n]) =
    aux (getOpt(Int.fromString level, 9), getOpt(Int.fromString n, 512))
  | main (_, _) = aux (9, 512)
val _ = main ("", CommandLine.arguments ())