.NET 4's improved tail call elimination

Existing versions of the .NET Framework sometimes unnecessarily fail to eliminate tail calls. Grant Richins of Microsoft recently published a blog post explaining how this has been addressed in .NET 4. The improvements make it entirely feasible to write idiomatic pure and impure functional F# programs relying upon tail call elimination without having to worry about stack overflows.

Such changes in the very fabric of .NET also demonstrate just how committed Microsoft are to their new F# programming language.

New VS2010 release date

Microsoft have updated the product page for Visual Studio 2010 (which includes F#!) with a new release date of 12th April 2010 under the "Buy" section.

24 days and counting!!!

Writing an article editor with Windows Presentation Foundation

The F#.NET Journal just published an article about the design and implementation of a complete GUI application:

"This article describes the article editing tool we are creating to help automate the process of generating journal content. The editor provides a simple text editor where the author writes a form of markdown. The markdown is interactively parsed and pretty printed in a separate pane and the editor allows the result to be exported to HTML for inclusion in the journal. The application is based entirely upon Windows Presentation Foundation..."

To read this article and more, subscribe to The F#.NET Journal today!

F# vs Unmanaged C++ for parallel numerics

We obtained a surprising performance result when comparing optimized parallel ray tracers written in F# and C++ recently. The following two programs render the same highly complex scenes containing over a million objects. Surprisingly, the 136-line managed F# program runs slightly faster at 17s than the 168-line unmanaged C++ which takes 18s.

Here is the F# program:

let t = System.Diagnostics.Stopwatch.StartNew()

let delta = sqrt epsilon_float

[<Struct>]
type vec =
  val x : float
  val y : float
  val z : float

  new(x, y, z) = {x=x; y=y; z=z}

  static member ( + ) (a : vec, b : vec) = new vec(a.x+b.x, a.y+b.y, a.z+b.z)
  static member ( - ) (a : vec, b : vec) = new vec(a.x-b.x, a.y-b.y, a.z-b.z)
  static member ( * ) (a, b : vec) = new vec(a*b.x, a*b.y, a*b.z)

let vec x y z = new vec(x, y, z)

let dot (a : vec) (b : vec) = a.x * b.x + a.y * b.y + a.z * b.z
let length r = sqrt(dot r r)
let unitise r = 1. / length r * r

type scene = Scene of vec * float * scene_type
and scene_type =
  | Sphere
  | Group of scene * scene * scene * scene * scene

let infinity = System.Double.PositiveInfinity

let inline ray_sphere (d : vec) (v : vec) r =
  let vx = v.x
  let vy = v.y
  let vz = v.z
  let disc = vx * vx + vy * vy + vz * vz - r * r
  if disc < 0. then infinity else
    let b = vx * d.x + vy * d.y + vz * d.z
    let b2 = b * b
    if b2 < disc then infinity else
      let disc = sqrt(b2 - disc)
      let t1 = b - disc in
      if t1 > 0. then t1 else b + disc

let ray_sphere' (o : vec) (d : vec) (c : vec) r =
  let vx = c.x - o.x
  let vy = c.y - o.y
  let vz = c.z - o.z
  let vv = vx * vx + vy * vy + vz * vz
  let b = vx * d.x + vy * d.y + vz * d.z
  let disc = b * b - vv + r * r
  disc >= 0. && b + sqrt disc >= 0.

type hit = {mutable l: float; mutable nx: float; mutable ny: float; mutable nz: float}

let intersect_sphere (dir: vec) hit l' (center: vec) =
    let x = l' * dir.x - center.x
    let y = l' * dir.y - center.y
    let z = l' * dir.z - center.z
    let il = 1. / sqrt(x * x + y * y + z * z)
    hit.l <- l'
    hit.nx <- il * x
    hit.ny <- il * y
    hit.nz <- il * z

let rec intersect dir hit = function
  | Scene(center, radius,typ) ->
      let l' = ray_sphere dir center radius
      if l' < hit.l then
          match typ with
          | Sphere ->
              intersect_sphere dir hit l' center
          | Group(a, b, c, d, e) ->
              intersect dir hit a
              intersect dir hit b
              intersect dir hit c
              intersect dir hit d
              intersect dir hit e

let rec intersect' orig dir = function
  | Scene(center, radius, typ) ->
      ray_sphere' orig dir center radius &&
      match typ with
      | Sphere -> true
      | Group (a, b, c, d, e) ->
          intersect' orig dir a ||
          intersect' orig dir b ||
          intersect' orig dir c ||
          intersect' orig dir d ||
          intersect' orig dir e

let neg_light = unitise(vec 1. 3. -2.)

let rec ray_trace dir scene =
  let hit = {l=infinity; nx=0.; ny=0.; nz=0.}
  intersect dir hit scene;
  if hit.l = infinity then 0. else
    let n = vec hit.nx hit.ny hit.nz in
    let g = dot n neg_light in
    if g < 0. then 0. else
      if intersect' (hit.l * dir + delta * n) neg_light scene then 0. else g

let fold5 f x a b c d e = f (f (f (f (f x a) b) c) d) e

let rec create level c r =
  let obj = Scene(c, r,Sphere)
  if level = 1 then obj else
    let a = 3. * r / sqrt 12.
    let rec bound (c, r) = function
      | Scene(c', r',Sphere) -> c, max r (length (c - c') + r')
      | Scene(_, _, Group(v, w, x, y, z)) -> fold5 bound (c, r) v w x y z
    let aux x' z' = create (level - 1) (c + vec x' a z') (0.5 * r)
    let w, x, y, z = aux a (-a), aux (-a) (-a), aux a a, aux (-a) a
    let c, r = fold5 bound (c + vec 0. r 0., 0.) obj w x y z in
    Scene(c, r, Group(obj, w, x, y, z))

let level, n = 11, 2048

let scene = create level (vec 0. -1. 4.) 1.
let ss = 4

do
  use ch = System.IO.File.Create("C:\Users\Jon\Documents\Visual Studio 10\Projects\F#\RayTracer\image.pgm", n * n + 32)
  sprintf "P5\n%d %d\n255\n" n n
  |> String.iter (ch.WriteByte << byte)
  let image = Array.create (n*n) 0uy
  System.Threading.Tasks.Parallel.For(0, n, fun y ->
    for x = 0 to n - 1 do
      let mutable g = 0.
      for dx = 0 to ss - 1 do
        for dy = 0 to ss - 1 do
          let aux x d = float x - 0.5 * float n + float d / float ss
          let dir = unitise(vec (aux x dx) (aux y dy) (float n))
          g <- g + ray_trace dir scene
      image.[x+y*n] <- 0.5 + 255. * g / float (ss*ss) |> byte)
  |> ignore
  ch.Write(image, 0, n*n)
  printf "Took %gs" t.Elapsed.TotalSeconds

And here is the C++ program:

#include "stdafx.h"
#include <vector>
#include <iostream>
#include <fstream>
#include <limits>
#include <cmath>

using namespace std;

numeric_limits<double> real;

double delta = sqrt(real.epsilon()), infinity = real.infinity();

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

Vec operator+(const Vec &a, const Vec &b)
{ return Vec(a.x+b.x, a.y+b.y, a.z+b.z); }

Vec operator-(const Vec &a, const Vec &b)
{ return Vec(a.x-b.x, a.y-b.y, a.z-b.z); }

Vec operator*(double a, const Vec &b) { return Vec(a*b.x, a*b.y, a*b.z); }

double dot(const Vec &a, const Vec &b) { return a.x*b.x + a.y*b.y + a.z*b.z; }

double length(const Vec &a) { return sqrt(dot(a, a)); }

Vec unitise(const Vec &a) { return (1 / sqrt(dot(a, a))) * a; }

struct Hit {
 double lambda;
  Vec normal;
  Hit() : lambda(infinity), normal(Vec(0, 0, 0)) {}
};

struct Sphere;

struct Scene {
  virtual ~Scene() {};
  virtual void intersect(Hit &, const Vec &) const = 0;
  virtual bool intersect(const Vec &, const Vec &) const = 0;
  virtual Sphere bound(Sphere b) const = 0;
};

struct Sphere : public Scene {
  Vec center;
  double radius;
  Sphere(Vec c, double r) : center(c), radius(r) {}
  ~Sphere() {}
  double ray_sphere(const Vec &dir) const {
    double b = center.x*dir.x + center.y*dir.y + center.z*dir.z;
    double disc = b*b - dot(center, center) + radius * radius;
    if (disc > 0) {
      double d = sqrt(disc), t2 = b + d;
      if (t2 > 0) {
  double t1 = b - d;
  return (t1 > 0 ? t1 : t2);
      }
    }
    return infinity;
  }
  bool sray_sphere(const Vec &orig, const Vec &dir) const {
    Vec v = center - orig;
    double b = dot(v, dir), disc = b*b - dot(v, v) + radius * radius;
    return (disc < 0 ? false : b + sqrt(disc) >= 0);
  }
  void intersect(Hit &hit, const Vec &dir) const {
    double lambda = ray_sphere(dir);
    if (lambda < hit.lambda) {
      hit.lambda = lambda;
      double
  nx = lambda*dir.x - center.x,
  ny = lambda*dir.y - center.y,
  nz = lambda*dir.z - center.z;
      double il = 1/sqrt(nx*nx + ny*ny + nz*nz);
      hit.normal.x = il*nx;
      hit.normal.y = il*ny;
      hit.normal.z = il*nz;
    }
  }
  bool intersect(const Vec &orig, const Vec &dir) const
  { return sray_sphere(orig, dir); }
  Sphere bound(Sphere b) const {
    return Sphere(b.center, max(b.radius, length(center - b.center) + radius));
  }
};

typedef vector<Scene *> Scenes;

struct Group : public Scene {
  Sphere b;
  Scenes child;
  Group(Sphere b, Scenes c) : b(b), child(c) {}
  ~Group() {
    for (Scenes::const_iterator it=child.begin(); it!=child.end(); ++it)
      delete *it;
  }
  void intersect(Hit &hit, const Vec &dir) const {
    if (b.ray_sphere(dir) < hit.lambda)
      for (Scenes::const_iterator it=child.begin(); it!=child.end(); ++it)
 (*it)->intersect(hit, dir);
  }
  bool intersect(const Vec &orig, const Vec &dir) const {
    if (!b.sray_sphere(orig, dir)) return false;
    for (Scenes::const_iterator it=child.begin(); it!=child.end(); ++it)
      if ((*it)->intersect(orig, dir)) return true;
    return false;
  }
  Sphere bound(Sphere b) const {
    Sphere b2 = b;
    for (Scenes::const_iterator it=child.begin(); it!=child.end(); ++it)
      b2 = (*it)->bound(b2);
    return b2;
  }
};

double ray_trace(const Vec &neg_light, const Vec &dir, const Scene &s) {
  Hit hit;
  s.intersect(hit, dir);
  if (hit.lambda == infinity) return 0;
  double g = dot(hit.normal, neg_light);
  if (g < 0) return 0.;
  Vec p = hit.lambda*dir + delta*hit.normal;
  return (s.intersect(p, neg_light) ? 0 : g);
}

Scene *create(int level, Vec c, double r) {
  Scene *s = new Sphere(c, r);
  if (level == 1) return s;
  Scenes child;
  child.reserve(5);
  child.push_back(s);
  double rn = 3*r/sqrt(12.);
  for (int dz=-1; dz<=1; dz+=2)
    for (int dx=-1; dx<=1; dx+=2)
      child.push_back(create(level-1, c + rn*Vec(dx, 1, dz), r/2));
  Sphere b2(c + Vec(0, r, 0), r);
  for (Scenes::const_iterator it=child.begin(); it!=child.end(); ++it)
    b2 = (*it)->bound(b2);
  return new Group(b2, child);
}

int _tmain(int argc, char* argv[]) {
  int level = (argc==3 ? atoi(argv[1]) : 11),
    n = (argc==3 ? atoi(argv[2]) : 2048), ss = 4;
  Vec neg_light = unitise(Vec(1, 3, -2));
  Scene *s(create(level, Vec(0, -1, 4), 1));
  ofstream out;
  std::vector<char> image(n*n);
  #pragma omp parallel for
  for (int y=0; y<n; ++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) {
    Vec dir(unitise(Vec(x+dx*1./ss-n/2., y+dy*1./ss-n/2., n)));
    g += ray_trace(neg_light, dir, *s);
  }
      image[x+(n-1-y)*n] = .5 + 255. * g / (ss*ss);
    }
  out.open("C:\\Users\\Jon\\Documents\\image.pgm", std::ios_base::binary);
  out << "P5\n" << n << " " << n << "\n255\n";
  out.write(&image[0], n*n);
  return 0;
}

An introduction to asynchronous workflows

The F#.NET Journal just published an article about concurrent programming:

"Asynchronous workflows were introduced into the F# programming language in 2007 and have recently matured pending the imminent release of the first officially-supported F# in Visual Studio 2010. This article describes the basic concepts that underpin this language feature and walks through some simple examples that demonstrate both basic and advanced use of this language feature..."

To read this article and more, subscribe to The F#.NET Journal today!