Also, the C++ code includes some optimisations not in the C code because it
already exceeded the 100 LOC limit of the shootout. Specifically,
specialised calls for shadow rays. This shaves off another 30% or so.

(* The Great Computer Language Shootout
   http://shootout.alioth.debian.org/
   Contributed by Jon Harrop, 2005
   Compile: ocamlopt -inline 100 ray.ml -o ray *)

(* This implementation differs from the original in several ways:

   Uses an implicit scene, generated as a ray is traced rather than being
   precalculated and stored explicitly in a tree.

   Specialized shadow-ray intersection functions. *)

let delta = sqrt epsilon_float and pi = 4. *. atan 1.

(* 3D vector and associated functions *)
type vec = {x:float; y:float; z:float}
let ( *| ) s r = {x = s *. r.x; y = s *. r.y; z = s *. r.z}
let ( +| ) a b = {x = a.x +. b.x; y = a.y +. b.y; z = a.z +. b.z}
let ( -| ) a b = {x = a.x -. b.x; y = a.y -. b.y; z = a.z -. b.z}
let dot a b = a.x *. b.x +. a.y *. b.y +. a.z *. b.z
let unitise r = (1. /. sqrt (dot r r)) *| r

(* A semi-infinite ray starting at "orig" and with direction "dir". *)
type ray = { orig: vec; dir: vec }

(* Calculate the parametric intersection of the given ray with the given
sphere. *)
let ray_sphere orig dir center radius =
  let v = center -| orig in
  let b = dot v dir in
  let disc = b *. b -. dot v v +. radius *. radius in
  if disc < 0. then infinity else
    let disc = sqrt disc in
    (fun t2 -> if t2 < 0. then infinity else
       ((fun t1 -> if t1 > 0. then t1 else t2) (b -. disc))) (b +. disc)

(* Calculate whether or not the given ray intersects the given sphere. *)
let ray_sphere' orig dir center radius =
  let v = center -| orig in
  let b = dot v dir in
  let disc = b *. b -. dot v v +. radius *. radius in
  if disc < 0. then false else b +. sqrt disc >= 0.

(* Ratio of the radii of one level of spheres to the next. *)
let s = 6. /. sqrt 12.

(* Find the first intersection point of the given ray with the scene. *)
let intersect level orig dir =
  let rec of_scene center radius lambda normal level =
    if level = 1 then
      let lambda' = ray_sphere orig dir center radius in
      if lambda' >= lambda then lambda, normal else
        lambda', unitise (orig +| lambda' *| dir -| center)
    else
      if ray_sphere orig dir center (3. *. radius) >= lambda
      then lambda, normal else
        let accu = of_scene center radius lambda normal 1 in
        let r = 0.5 *. radius and l = level - 1 in let r' = s *. r in
        let aux dx dz (lambda, normal) =
          of_scene (center +| {x=dx; y=r'; z=dz}) r lambda normal l in
        let mr' = -.r' in
        aux r' mr' (aux r' r' (aux mr' r' (aux mr' mr' accu))) in
  of_scene {x=0.; y= -1.; z=0.} 1. infinity {x=0.; y=0.; z=0.} level

(* Find if the given ray intersects the scene. *)
let intersect' level orig dir =
  let rec of_scene center radius level =
    if level = 1 then ray_sphere' orig dir center radius else
      (* Exploit short-circuit evaluation of boolean comparisons to
terminate
         this function early. *)
      ray_sphere' orig dir center (3. *. radius) &&
        (of_scene center radius 1 ||
           let r = 0.5 *. radius and l = level - 1 in let r' = s *. r in
           of_scene (center +| {x= -.r'; y=r'; z= -.r'}) r l ||
             of_scene (center +| {x= r'; y=r'; z= -.r'}) r l ||
             of_scene (center +| {x= -.r'; y=r'; z= r'}) r l ||
             of_scene (center +| {x= r'; y=r'; z= r'}) r l) in
  of_scene {x=0.; y= -1.; z=0.} 1. level

(* Trace a single ray by casting it into the scene and, if it intersects
   anything, casting a second ray toward the light to determine occlusion.
*)
let rec ray_trace l light orig dir =
  let lambda, n = intersect l orig dir in
  if lambda = infinity then 0. else
    let g = -. dot n light in
    (* If we are on the shadowed side of a sphere then don't bother casting
a
       shadow ray as we know it will intersect the same sphere. *)
    if g <= 0. then 0. else
      let p = orig +| lambda *| dir +| delta *| n in
      if intersect' l p ({x=0.; y=0.; z=0.} -| light) then 0. else g

(* Ray trace the scene at the given resolution. *)
let () =
  (* Resolution *)
  let n = match Sys.argv with [| _; l |] -> int_of_string l | _ -> 256 in
  (* Light direction *)
  let light = unitise {x= -1.; y= -3.; z=2.} in
  (* Number of levels of spheres, and oversampling. *)
  let level = 6 and ss = 4 in

  Printf.printf "P5\n%d %d\n255\n" n n;
  for y = n - 1 downto 0 do
    for x = 0 to n - 1 do
      (* Average each pixel over ss*ss separate rays. *)
      let g = ref 0. in
      for dx = 0 to ss - 1 do
        for dy = 0 to ss - 1 do
          (* Calculate the origin and direction of this ray. *)
          let orig = {x=0.; y=0.; z= -4.} in
          let dir = unitise {x = float (x - n / 2) +. float dx /. float ss;
                             y = float (y - n / 2) +. float dy /. float ss;
                             z = float n} in
          g := !g +. ray_trace level light orig dir
        done
      done;
      let g = int_of_float (0.5 +. 255. *. !g /. float (ss*ss)) in
      Printf.printf "%c" (char_of_int g)
    done
  done

Here's the C++:

// The Great Computer Language Shootout
// http://shootout.alioth.debian.org/
// Contributed by Jon Harrop, 2005
// Compile: g++ -Wall -O3 -ffast-math ray.cpp -o ray

#include <vector>
#include <iostream>
#include <limits>
#include <cmath>

using namespace std;

numeric_limits<double> real;
double delta = sqrt(real.epsilon()), infinity = real.infinity(), pi = M_PI;

// 3D vector
struct Vec {
  double x, y, z;
  Vec(double x2, double y2, double z2) : x(x2), y(y2), z(z2) {}

// Semi-infinite ray
struct Ray { Vec orig, dir; Ray(Vec o, Vec d) : orig(o), dir(d) {} };

// Scene tree
// In this implementation, a node in the scene tree is represented by a
single
// struct which is either a group of scene trees with a spherical bound or,
// implicitly, a single sphere if the group has no children.
// This is not equivalent to the variant type used to represent a node in
the
// original OCaml implementation because this representation cannot
associate
// data with only leaf nodes (such as color, reflectivity etc.) but requires
// significantly less C++ code and gives room to implement a specialized
// shadow-ray intersection algorithm.
struct Scene {
  vector<Scene> child; // Child nodes in the scene tree
  Vec center; // Center of the sphere or spherical bound
  double radius; // Radius of the sphere or spherical bound

  Scene(Vec c, double r) : child(), center(c), radius(r) {}

  std::cout << "P5\n" << n << " " << n << "\n255\n";
  for (int y=n-1; y>=0; --y) {
    for (int x=0; x<n; ++x) {
      double g=0;
      for (int dx=0; dx<ss; ++dx)
        for (int dy=0; dy<ss; ++dy) {
          // We use "dx*1." instead of "double(dx)" to save space
          Vec d(x+dx*1./ss-n/2., y+dy*1./ss-n/2., n);
          g += ray_trace(1, light, Ray(Vec(0, 0, -4), unitise(d)), s);
        }
      std::cout << char(.5 + 255. * g / (ss*ss));
    }
  }

  // In this implementation the Scene destructor recursively deallocates the
  // entire scene for us, obviating the need for manually-defined
destructors.

  return 0;