• A Mathematica Optimization Problem
  http://xahlee.org/UnixResource_dir/writ/Mathematica_optimization.html

delta = Sqrt[$MachineEpsilon];
RaySphere[o_, d_, c_, r_] := 
      Block[{v, b, disc, t1, t2}, 
        v = c - o; 
b = v . d; 
disc = 
          Sqrt[b^2 - v . v + r^2]; 
t2 = 
          b + disc; 
If[
          Im[disc] ≠ 0 || t2 ≤ 0, ∞, 
          t1 = b - disc; 
If[t1 > 0, t1, 
            t2]]]

    
Intersect[o_, d_][{lambda_, n_}, Sphere[c_, r_]] := 
      Block[{lambda2 = RaySphere[o, d, c, r]}, 
        If[lambda2 ≥ lambda, {lambda, n}, {lambda2, 
            Normalize[o + lambda2 d - c]}]]

    
Intersect[o_, d_][{lambda_, n_}, Bound[c_, r_, s_]] := 
      Block[{lambda2 = RaySphere[o, d, c, r]}, 
        If[lambda2 ≥ lambda, {lambda, n}, 
          Fold[Intersect[o, d], {lambda, n}, s]]]

    
neglight = N@Normalize[{1, 3, -2}];
nohit = {∞, {0, 0, 0}};

    
RayTrace[o_, d_, scene_] := 
      Block[{lambda, n, g, 
          p}, {lambda, n} = Intersect[o, d][nohit, 
            scene]; 
If[lambda == ∞, 0, 
          g = n . neglight; 
If[g ≤ 0, 
            0, {lambda, 
                n} = Intersect[o + lambda d + delta n, neglight][nohit, 
                scene]; 
If[lambda < ∞, 0, 
              g]]]]

    
Create[level_, c_, r_] := 
      Block[{obj = Sphere[c, r]}, 
        If[level == 1, obj, 
          Block[{a = 3*r/Sqrt[12], Aux}, 
            Aux[x1_, z1_] := 
              Create[level - 1, c + {x1, a, z1}, 
                0.5 r]; 
Bound[c, 
              3 r, {obj, Aux[-a, -a], Aux[a, -a], 
                Aux[-a, a], Aux[a, a]}]]]]

    
scene = Create[1, {0, -1, 4}, 1];
Main[level_, n_, ss_] := 
      Block[{scene = Create[level, {0, -1, 4}, 1]}, 
        Table[Sum[
              RayTrace[{0, 0, 0}, 
                N@Normalize[{(x + s/ss/ss)/n - 
                        1/2, (y + Mod[s, ss]/ss)/n - 1/2, 1}], 
                scene], {s, 0, ss^2 - 1}]/ss^2, {y, 0, n - 1}, {x, 0, 
            n - 1}]]

    
AbsoluteTiming[
      Export["image.pgm", 
        Graphics@Raster@Main[9, 512, 4]]]