[Haskell-cafe] Implementing Mathematica

[Jon Harrop]

On Wednesday 30 May 2007 22:15:55 Andrew Coppin wrote:
> Jon Harrop wrote:
> > I wrote a toy Mathematica implementation in OCaml while I
> > waited to be viva'd for my PhD. It garnered so much interest that Wolfram
> > Research bought it from me for £4,500 and gave me several free copies of
> > Mathematica.
>
> Are you serious?! o_O

Yes.

> > 1. Numerical libraries: you should be able to reuse existing libraries
> > like GMP, BLAS, LAPACK, FFTW and so on. These are often much faster than
> > Mathematica's. For example, FFTW was about 4x faster than Mathematica's
> > FFT the last time I benchmarked it. However, they do not support interval
> > arithmetic.
>
> Now this is interesting. The claim is that Mathematica is the fastest,
> most powerful software on planet earth, second to none.

If you write a simple, numerically-intensive program that runs in the 
Mathematica rewriter then its performance is about 100-1,000x slower than 
that of a native-code compiled language like Haskell. Mathematica is often 
30x slower than interpreted OCaml bytecode.

Take this ray tracer, for example:

scene = {Sphere[{0., 0., 4.}, 1.], Sphere[{-1., 1., 4.}, 1.], 
         Sphere[{-1., -1., 4.}, 1.], Sphere[{1., 1., 4.}, 1.], 
         Sphere[{1., -1., 4.}, 1.]};

\[Delta] = Sqrt[$MachineEpsilon];

Unitise[p_] := p/Sqrt[p.p]

RaySphere[o_, d_, c_, r_] :=
  Block[{v = c - o, b = v.d, disc = b^2 - v.v + r^2},
    If[disc <= 0., \[Infinity],
      disc = Sqrt[disc];
      Block[{t2 = b + disc},
        If[t2 <= 0., \[Infinity],
          Block[{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, Unitise[o + lambda2 d - c]}]
    ]
Intersect[o_, d_][hit_, list_List] := Fold[Intersect[o, d], hit, list]

nohit = {\[Infinity], {0., 0., 0.}};

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

Timing[image = 
      Table[Table[
          RayTrace[{0., 0., -2.5}, Unitise[{x, y, 128.}], scene], {y, -64, 
            64}], {x, -64, 64}];]

This program takes 4.8s to run here. I bet if we translate it into Haskell it 
will run much faster than that.

As a guide, this Haskell ray tracer is much more advanced and it can render a 
bigger (100x100) image in only 0.2s:

  http://www.nobugs.org/developer/htrace/

Incidentally, when I try to recompile with optimizations turned on, GHC 
refuses to work:

$ ghc htrace.hs -o htrace
$ ghc -O2 htrace.hs -o htrace
compilation IS NOT required

I must delete the target or edit the source to get it to recompile. I assume 
this is a known bug?

> Actually it 
> turns out that at least for factoring moderately big integers, pari/gp
> seems to be a fair bit faster (50% or so). I have no idea about the rest.

Right, but that is just calling an internal function that is written in C. 
Provided you only ever call a few library functions, performance will be 
excellent in Mathematica. But when you cannot phrase your program in terms of 
the built-in library functions, performance is terrible and this is when 
everyone reaches for a more efficient tool.

To me, performance is way down on the list of things that make Mathematica 
great. Integrated graphics is probably the main benefit. I mostly use MMA as 
a glorified graph plotter.

> Note that (as I understand it) GHC implements Haskell's Integer type
> using the GMP. And for some reason or other, they want to remove this
> feature...

Arbitrary precision integers are quite a performance burden and they are 
rarely used. I would not expect a language that is trying to be efficient to 
impose arbitrary precision integers (or floats).

> Erm... have you seen Mathematica 6?

Yes.

> That's OpenGL accelerated too.

Yes. Similarly, making graphics fast takes a lot more than just "using 
OpenGL".

> I've just been playing with it in fact - it's pretty fast as far as I can
> tell= 

It still invokes a completely generic term rewriter for everything it 
evaluates. You can really feel this when you play with some of their 
interactive demos. Even the simple ones are notably sluggish. Try translating 
the Tiger demo from our site into Mathematica, for example.

Many of their demos will be trivial to write in Haskell but performance would 
be a lot better. I'd like to write some graphics demos in Haskell using 
OpenGL...

-- 
Dr Jon D Harrop, Flying Frog Consultancy Ltd.
OCaml for Scientists
http://www.ffconsultancy.com/products/ocaml_for_scientists/?e