71 lines
1.9 KiB
Plaintext
71 lines
1.9 KiB
Plaintext
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
|