ClosestPair := module()

local
    ModuleApply := proc(L::list,$)
    local Lx, Ly, out;
        Ly := sort(L, 'key'=(i->i[2]), 'output'='permutation');
        Lx := sort(L, 'key'=(i->i[1]), 'output'='permutation');
        out := Recurse(L, Lx, Ly, 1, numelems(L));
        return sqrt(out[1]), out[2];
    end proc; # ModuleApply

local
    BruteForce := proc(L, Lx, r1:=1, r2:=numelems(L), $)
    local d, p, n, i, j;
        d := infinity;
        for i from r1 to r2-1 do
            for j from i+1 to r2 do
                n := dist( L[Lx[i]],  L[Lx[j]] );
                if n < d then
                    d := n;
                    p := [ L[Lx[i]], L[Lx[j]] ];
                end if;
            end do; # j
        end do; # i
        return (d, p);
    end proc; # BruteForce

local dist := (p, q)->(( (p[1]-q[1])^2+(p[2]-q[2])^2 ));

local Recurse := proc(L, Lx, Ly, r1, r2)
    local m, xm, rDist, rPair, lDist, lPair, minDist, minPair, S, i, j, Lyr, Lyl;

    if r2-r1 <= 3 then
        return BruteForce(L, Lx, r1, r2);
    end if;

    m := ceil((r2-r1)/2)+r1;
    xm := (L[Lx[m]][1] + L[Lx[m-1]][1])/2;

    (Lyr, Lyl) := selectremove( i->L[i][1] < xm, Ly);

    (rDist, rPair) := thisproc(L, Lx, Lyr, r1, m-1);
    (lDist, lPair) := thisproc(L, Lx, Lyl, m, r2);

    if rDist < lDist then
        minDist := rDist;
        minPair := rPair;
    else
        minDist := lDist;
        minPair := lPair;
    end if;

    S := [ seq( `if`(abs(xm - L[i][1])^2< minDist, L[i], NULL ), i in Ly ) ];

    for i from 1 to nops(S)-1 do
        for j from i+1 to nops(S) do
            if abs( S[i][2] - S[j][2] )^2 >= minDist then
                break;
            elif dist(S[i], S[j]) < minDist then
                minDist := dist(S[i], S[j]);
                minPair := [S[i], S[j]];
            end if;
        end do;
    end do;

    return (minDist, minPair);

    end proc; #Recurse

end module; #ClosestPair