I recently tried my hand at writing a ray tracer (something I haven't done
for 15 years!) and, in particular, I to write the shortest possible ray
tracer which was still comprehensible and interesting. The initial result
was a 222-line OCaml program which incrementally renders a sphere-flake via
OpenGL:

  http://www.ffconsultancy.com/free/ray_tracer/index.html

I then cut this program down to a 66-line OCaml program and ported it into a
97-line C++ program. These programs are compared on this page:

  http://www.ffconsultancy.com/free/ray_tracer/comparison.html

The cut-down versions don't do reflections, color or such a pretty
sphere-flake and output a greyscale PGM but they do still use hierarchical
spherical bounding volumes to accelerate rendering enough that they can
render scenes containing millions of spheres.

Somewhere on an opteron 146 running debian64:
# g++ -g -O3 ray.cpp -o ray
# time ./ray >tt

real    0m13.959s
user    0m13.946s
sys     0m0.009s

# objdump -x ray
[2mn later]
# diff ray_orig.cpp ray.cpp
12,16c12,18
< Vec operator+(Vec a, Vec b) { return Vec(a.x + b.x, a.y + b.y, a.z +
b.z); }
< Vec operator-(Vec a, Vec b) { return Vec(a.x - b.x, a.y - b.y, a.z -
b.z); }
< Vec operator*(double a, Vec b) { return Vec(a * b.x, a * b.y, a *
b.z); }
< double dot(Vec a, Vec b) { return a.x*b.x + a.y*b.y + a.z*b.z; }
< Vec unitise(Vec a) { return (1 / sqrt(dot(a, a))) * a; }
---

real    0m10.493s
user    0m10.483s
sys     0m0.005s

Note also that he mentions on his web page that he's doing this for "the
computer language shootout benchmarks", http://shootout.alioth.debian.org/.
That page clearly points out all the caveats around benchmarks like these
(to the point that they can be effectively meaningless in  the real world.)

Specifically:

In fact if you do the same for other hotpath functions as in:
19c21
<   virtual void intersect(double &, Vec &, Ray) = 0;
---

You'll notice i haven't even fixed the gratitious use of virtual
functions, and other details that force standard compliant c++
compilers to pessimize a lot when facing such... hmm.. source.

I'm sure OCaml is a wonderful & expressive language, but crafting such
a fishy c++ equivalent to make it shine (because you have an agenda) is
a despicable practice.

But i guess the main point of the whole operation wasn't fair
benchmarking, education, or fun but selling books: "This program has
demonstrated some of the ways that OCaml improves upon C++. For a
thorough introduction to OCaml, read our book "Objective CAML for
Scientists".

Or put another way, that's a sad PR stunt.

Moreover, your choices of initial optimisation are decidedly suboptimal. I'd
have gone for:

1. Terminating the intersection of shadow rays when the first intersection
is found (instead of finding the closest intersection).

2. Use an implicit scene, to avoid storing it explicitly.

3. Use single-precision storage (or no storage at all).

Both of these optimisations will give much bigger performance improvements
when averaged over different platforms and architectures.

As it happens, I had already implemented all of your optimisations and all
of these optimisations to both implementations before I made my original
post. None of them add anything new. Sometimes C++ is faster, sometimes
OCaml is faster. OCaml is always much more succinct. You can keep doing
this until you are blue in the face but I don't think you'll learn much of
interest.

As OCaml is for symbolic use, it performs no inlining by default (unlike
g++) so it is common to ask OCaml to do some inlining on numerical code.
Thus, specifying "-inline 100" actually makes for a fairer comparison.

In fact, I'd already implemented your optimisations (and many more effective
optimisations) and had to throw them out because of this limit. Note that
the OCaml has 30 more lines with which to optimise. My optimised <100 LOC
OCaml program is much faster than my optimised <100 LOC C++ on all
platforms.

I chose an inheritance hierarchy because this is the closest C++ equivalent
to a variant type which I see in typical C++ programs.

Another optimisation that I made was to replace inheritance with a single
Scene struct which represented a group of child scenes or a sphere when it
had no children. I threw this out as well because it is not a fair
equivalent to OCaml's variant type. For example, you could not specify a
colour for each sphere without also specifying colours for all bounding
volumes.

Also, it generates (numerically) different results, which makes it trickier
to validate against the OCaml.

  The keyword 'inline' is not a command for the compiler that says it must
inline the function or else it will be deleted immediately. And the other
way around, the lack of 'inline' keyword by no means stops the compiler
from inlining the function if it can and it estimates it would be useful.

  In *theory* 'inline' is a hint to the compiler. In practice it has nothing
to do with inlining. Most of the modern compilers will simply ignore that
keyword when they are estimating whether a function should be inlined or
not. In this context 'inline' is a no-op. It's like the 'register'
keyword which is nowadays obsolete and a no-op.

  What 'inline' *does* have some effect on is the linker. If a function
with an identical declaration is implemented in more than one place
then in normal circumstances the linker will give an error (because it
can't know which one of them it should use). However, if the functions
were declared 'inline', then the compiler will tell the linker that
all of those functions are actually one and the same, and then the linker
will just use one of them.

  And by the way, if you are afraid of "litter" in C++ code, then you
should switch to another programming language.
  That "litter", as you call it, is not useless junk as you seem to imply.
It has a reason and a significance.

  If your opinion of "clearer" and "shorter" is that the & symbol should
not be used, then you might want to either check your priorities or
change your programming language (because C++ is not for you).

ray.cpp: In member function `virtual void Sphere::intersect(double&, Vec&,
Ray)
   ':
ray.cpp:39: error: no match for 'operator+' in 'ray.Ray::orig +
   operator*(double, Vec)(((&ray) + 24))'
ray.cpp:12: error: candidates are: Vec operator+(Vec&, Vec&)
/usr/include/c++/3.3/bits/stl_bvector.h:261: error:
   std::_Bit_const_iterator std::operator+(int, const
   std::_Bit_const_iterator&)
/usr/include/c++/3.3/bits/stl_bvector.h:202: error:
   std::_Bit_iterator std::operator+(int, const std::_Bit_iterator&)
ray.cpp: In function `double ray_trace(double, Vec, Ray, Scene*)':
ray.cpp:62: error: no match for 'operator+' in 'ray.Ray::orig +
   operator*(double, Vec)(((&ray) + 24))'
ray.cpp:12: error: candidates are: Vec operator+(Vec&, Vec&)
/usr/include/c++/3.3/bits/stl_bvector.h:261: error:
   std::_Bit_const_iterator std::operator+(int, const
   std::_Bit_const_iterator&)
/usr/include/c++/3.3/bits/stl_bvector.h:202: error:
   std::_Bit_iterator std::operator+(int, const std::_Bit_iterator&)

You must either use "const ... &" instead, or make sure that all values are
created explicitly as local variables in the caller (rather than being
inline subexpressions in the function call).

It is precisely this obfuscation (and the fact that these are tiny
structures, and that this is not necessarily an optimisation, and that it
could be done by the compiler) that made me choose pass by value.

  You very seldom need non-const references.

So you're saying intersections aren't part of the hotpath?
Well you maybe right, after all it's only called a few million times.

And i suppose that the fact that this emasculated C++ version performed
slower than its OCaml counterpart was merely a side effect?

What's next? You're going to try to sell me a bridge?
Oops. A book first, then a bridge.

As there are currently only two values being returned by the intersect
routine, I could have used STL pair<> instead. I might investigate that (I
think it will be simpler but slower). I have seen code which uses
pair<pair<...>, ...> to return three values but I don't much like it.

To add a new kind of object to the C++, you create a new class derived from
Scene which implements its own "intersect" member function.

To add a new kind of object to the OCaml, you supplement the variant type
"scene" with another constructor and supplement the pattern match in the
"intersect" function with a case to handle your new object.

Similarly for new types of bounding volume.

Perhaps it would be interesting to add a new kind of object to both
implementations, to see how this is done in each language?

As I've said, if you want to optimise my ray tracer (or any other program)
then I would recommend algorithmic optimisations over low-level ones like
SIMD, cache coherency etc. There are plenty more algorithmic optimisations
to be done (although I've already done the main one - hierarchical bounding
volumes).

In this case, it looks as though you'd want to write specialised C code for
certain platforms. This could be called from C++ or OCaml and it would
speed both programs up by roughly the same amount. This might be a good
idea though, as OCaml is clearly preferable for the higher-level, not-
performance-critical bulk of the code. So you could get rid of C++
altogether. :-)

  When there's a logical reason to return many values at once, returning
a class or struct as the return value of the function may in fact be a
good idea. If you call the function like this:

ReturnValues values = theFunction(whatever);

where ReturnValues is a class or struct, C++ compilers (or at least gcc;
probably also most of the other compilers) will actually make that very
efficiently (they will perform the so-called return value optimization).

  The return value optimization means that, in a code like the above,
the function being called actually builds its return value directly
into 'values' and will not create any temporaries nor will it call any
copy constructors. Perhaps a bit surprisingly this works even if the
function in question is not inline and its implementation is in a
completely separate compilation unit (ie. object file) which is just
linked later to the binary (and in fact it works even if the function
in question is in a dynamically loaded library).

  It might even be that the return value optimization technique produces
slightly faster code than the function returning its values by writing
to references given to it as parameters.

  However, if you can design the program so that the function to be
called can be resolved at compilation time (instead of at runtime as
happens with virtual functions or switch-case blocks) then you might
get a slight speedup, specially if the function is short and called
millions of times.

  Have you ever heard of abstraction?

real    0m5.862s
user    0m5.800s
sys     0m0.006s
# cmp pix.orig.ppm pix.ppm; echo $?
0
# composite -compose difference pix.orig.ppm pix.ppm diff.ppm
# display diff.ppm
[contemplate a black picture]

So, on my box, the original version took ~13s then ~10s after some
patching, ~8s after some more patching and now we're at ~6s with a
const correct non virtual rule compliant version.

Before you start nitpicking please consider that:
a) this version complies with the shootout rules AFAIK.
b) it outputs the exact same picture.
c) i had to raise the epsilon a bit to avoid acnea, it has to do with
the way i compute normals but i'm not going to fix that.
d) i removed the virtualness in the most straightforward way and it's
certainly not aesthetically pleasing
e) i'm not going to spend ages on this, hence c) & d).
f) it's algorithmically equivalent to your version (ie primary and
shadow rays aren't special cased etc...)

For example, the following OCaml function returns the given number as well
as its square and cube:

let f n = n, n*n, n*n*n;;

The C++ equivalent of this using a struct is:

struct res { int n, n2, n3; };

res f(int n) {
  res ans;
  ans.n = n;
  ans.n2 = n*n;
  ans.n3 = n*n*n;
  return ans;

pair<int, pair<int, int> > f(int n) {
  return pair<int, pair<int, int> >(n, n*n, n*n*n);

I would be very interested to see a highly optimised version of my ray
tracer. I think a lot of people could learn from a set of web pages which
show how to optimise the programs. If you're willing to help and I can find
the time then I'd like to do this sometime.

  a = _mm_load_ps(inp_sse1);
  b = _mm_load_ps(inp_sse2);
  c = _mm_add_ps(a, b);
  _mm_store_ps(out_sse, c);