with javascript_semantics
function bruteForceClosestPair(sequence s)
    atom {x1,y1} = s[1], {x2,y2} = s[2],
         dx = x1-x2, dy = y1-y2,
         mind = dx*dx+dy*dy
    sequence minp = s[1..2]
    for i,si in s to length(s)-1 do
        {x1,y1} = si
        for sj in s from i+1 do
            {x2,y2} = sj
            dx = x1-x2
            dx = dx*dx
            if dx<mind then
                dy = y1-y2
                dx += dy*dy
                if dx<mind then
                    {mind,minp} = {dx,{si,sj}}
                end if
            end if
        end for
    end for
    minp = sort(minp) -- (see note below)
    return {sqrt(mind),minp}
end function

sequence testset = sq_rnd(repeat({1,1},10000))
atom t0 = time()
atom {d, {{x1,y1},{x2,y2}}} = bruteForceClosestPair(testset)
-- (Sorting the final point pair makes brute/dc more likely to tally. Note however
--  when >1 equidistant pairs exist, brute and dc may well return different pairs;
--  it is only a problem if they decide to return different minimum distances.)
printf(1,"Closest pair: {%f,%f} {%f,%f}, distance=%f (%3.2fs)\n",{x1,y2,x2,y2,d,time()-t0})

t0 = time()
constant X = 1, Y = 2
sequence xP = sort(deep_copy(testset)),
         yP = sort_columns(deep_copy(testset),{Y})

function distsq(sequence p1,p2)
    atom {x1,y1} = p1, {x2,y2} = p2
    x1 -= x2
    y1 -= y2
    return x1*x1 + y1*y1
end function

function closestPair(sequence xP, yP)
-- where xP is P(1) .. P(N) sorted by x coordinate, and
--       yP is P(1) .. P(N) sorted by y coordinate (ascending order)
    integer N = length(xP),
         midN = floor(N/2)
    assert(length(yP)=N)
    if N<=3 then
        return bruteForceClosestPair(xP)
    end if
    sequence xL = xP[1..midN],
             xR = xP[midN+1..N],
             yL = {},
             yR = {}
    atom xm = xP[midN][X]
    for yi in yP do
        if yi[X]<=xm then
            yL = append(yL,yi)
        else
            yR = append(yR,yi)
        end if
    end for
    {atom dmin, sequence pairMin} = min(closestPair(xL, yL),
                                        closestPair(xR, yR))
    sequence yS = {}
    for yi in yP do
        if abs(xm-yi[X])<dmin then
            yS = append(yS,yi)
        end if
    end for
    atom closest = dmin*dmin
    for i,yi in yS to length(yS)-1 do
        for yk in yS from i+1 do
            if yk[Y]-yi[Y]>=dmin then exit end if
            d = distsq(yk,yi)
            if d<closest then
                {closest, pairMin} = {d, {yk, yi}}
            end if
        end for
    end for
    pairMin = sort(pairMin) -- (see note above)
    return {sqrt(closest), pairMin}
end function

{d,{{x1,y1},{x2,y2}}} = closestPair(xP,yP)
printf(1,"Closest pair: {%f,%f} {%f,%f}, distance=%f (%3.2fs)\n",{x1,y2,x2,y2,d,time()-t0})