MATLAB colon operator source code, obtained from http://www.mathworks.com/support/solutions/attachment.html?resid=1-4FLIOY&solution=1-4FLI96 as of 2013-11-09.

colonop.m

function v = colonop(a,d,b)
% COLONOP  Demonstrate how the built-in a:d:b is constructed.
%
%   v = colonop(a,b) constructs v = a:1:b.
%   v = colonop(a,d,b) constructs v = a:d:b.
%
%   v = a:d:b is not constructed using repeated addition.  If the
%   textual representation of d in the source code cannot be
%   exactly represented in binary floating point, then repeated
%   addition will appear to have accumlated roundoff error.  In
%   some cases, d may be so small that the floating point number
%   nearest a+d is actually a.  Here are two imporant examples.
%
%   v = 1-eps : eps/4 : 1+eps is the nine floating point numbers
%   closest to v = 1 + (-4:1:4)*eps/4.  Since the spacing of the
%   floating point numbers between 1-eps and 1 is eps/2 and the
%   spacing between 1 and 1+eps is eps,
%   v = [1-eps 1-eps 1-eps/2 1 1 1 1 1+eps 1+eps].
%
%   Even though 0.01 is not exactly represented in binary,
%   v = -1 : 0.01 : 1 consists of 201 floating points numbers
%   centered symmetrically about zero.
%
%   Ideally, in exact arithmetic, for b > a and d > 0,
%   v = a:d:b should be the vector of length n+1 generated by
%   v = a + (0:n)*d where n = floor((b-a)/d).
%   In floating point arithmetic, the delicate computatations
%   are the value of n, the value of the right hand end point,
%   c = a+n*d, and symmetry about the mid-point.

if nargin < 3
    b = d;
    d = 1;
end

tol = 2.0*eps*max(abs(a),abs(b));
sig = sign(d);

% Exceptional cases.

if ~isfinite(a) || ~isfinite(d) || ~isfinite(b)
   v = NaN;
   return
elseif d == 0 || a < b && d < 0 || b < a && d > 0
   % Result is empty.
   v = zeros(1,0);
   return
end

% n = number of intervals = length(v) - 1.

if a == floor(a) && d == 1
   % Consecutive integers.
   n = floor(b) - a;
elseif a == floor(a) && d == floor(d)
   % Integers with spacing > 1.
   q = floor(a/d);
   r = a - q*d;
   n = floor((b-r)/d) - q;
else
   % General case.
   n =  round((b-a)/d);
   if sig*(a+n*d - b) > tol
      n = n - 1;
   end
end

% c = right hand end point.

c = a + n*d;
if sig*(c-b) > -tol
   c = b;
end

% v should be symmetric about the mid-point.

v = zeros(1,n+1);
k = 0:floor(n/2);
v(1+k) = a + k*d;
v(n+1-k) = c - k*d;
if mod(n,2) == 0
   v(n/2+1) = (a+c)/2;
end