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 | Group of (scene * vec * real) list

fun intersect(orig, dir, scene) =
    let fun aux((scene, center, radius), first as (l, _)) =
	    let val l' =
		    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 in
		if l' >= l then first else
		case scene of
		    Sphere => (l', unitise (orig +| l' *| dir -| center))
		  | Group scenes => foldl aux first scenes
	    end in
	aux(scene, (infinity, (0.0, 0.0, 0.0)))
    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)
	val g = dot(n, neg_light) in
	if g <= 0.0 then 0.0 else
	let val p = lambda *| dir +| delta *| n
	    val (l, _) = intersect(p, neg_light, scene) in
	    if l >= infinity then g else 0.0
	end
    end

fun create(level, r, c as (x, y, z)) =
    let val obj = (Sphere, c, 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'] in
	    (Group objs, c, 3.0 * r)
	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 ())