[Haskell-cafe] Implementing Mathematica
Jon Harrop
jon at ffconsultancy.com
Wed May 30 18:56:30 EDT 2007
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
More information about the Haskell-Cafe
mailing list