function ClosePoints_mzip % MATLAB zip file, a self-extracting MATLAB archive. % Usage: Run this file to recreate the original directory. fname = mfilename; fin = fopen([fname '.m'],'r'); dname = fname(1:find(fname=='_',1,'last')-1); mkdir(dname); mkdir([dname filesep 'lib']) addpath(dname) addpath([dname filesep 'lib']) L = fgetl(fin); while length(L) < 2 || ~isequal(L(1:2),'%%') L = fgetl(fin); end while ~isequal(L,'%% EOF') F = [dname filesep L(4:end)]; disp(F) fout = fopen(F,'w'); L = fgetl(fin); while length(L) < 2 || ~isequal(L(1:2),'%%') fprintf(fout,'%s\n',L); L = fgetl(fin); end fclose(fout); end fclose(fin); end %% Center.m function [d,kout,jout] = Center(z,din) % Center(z,d) is used by DivCon to examine the % strip of half-width d about the center point. n = length(z); m = floor(n/2); xh = real(z(m)); d = din; [~,p] = sort(imag(z)); z = z(p); s = []; ks = []; ns = 0; for k = 1:n if abs(real(z(k)) - xh) < d ns = ns+1; s(ns,1) = z(k); ks(ns,1) = k; end end kout = NaN; jout = NaN; for k = 1:ns for j = k+1:ns if (imag(s(j)) - imag(s(k))) < d && abs(s(k) - s(j)) < d d = abs(s(k) - s(j)); kout = p(ks(k)); jout = p(ks(j)); end end end end %% DivCon.m function [d,kout,jout] = DivCon(z,~) % [d,k,j] = DivCon(z) finds indices k and j so that % z(k) and z(j) minimizes d = abs(z(k)-z(j)). % % DivCon employs a recursive divide and conquer algorithm with % complexity O(n*log(n)) where n = length(z). % [d,k,j] = DivCon(z,true) is a recursive call with ascending real(z). % % Center(z,d) examines the central strip of half-width d. % Pairs(z) computes the same thing as DivCon, but has complexity O(n^2). n = length(z); if n <= 3 [d,kout,jout] = Pairs(z); else if nargin < 2 [~,p] = sort(real(z)); z = z(p); else p = (1:n)'; end m = floor(n/2); % Left half [dl,kl,jl] = DivCon(z(1:m),true); % Right half [dr,kr,jr] = DivCon(z(m+1:end),true); if dl < dr d = dl; k = kl; j = jl; else d = dr; k = kr+m; j = jr+m; end % Center strip [ds,ks,js] = Center(z,d); if ds < d d = ds; k = ks; j = js; end kout = min(p(k),p(j)); jout = max(p(k),p(j)); end end %% Pairs.m function [d,kout,jout] = Pairs(z) % [d,k,j] = Pairs(z) finds indices k and j so that % z(k) and z(j) minimizes d = abs(z(k)-z(j)). % % Pairs employs a brute force algorithm that examines all possible % pairs. This has complexity O(n^2) where n = length(z). % DivCon(z) computes the same thing, but has complexity O(n*log(n)) % and consequently is much faster. n = length(z); d = Inf; for k = 1:n for j = k+1:n if abs(z(k) - z(j)) < d d = abs(z(k) - z(j)); kout = k; jout = j; end end end end %% xDivCon.m function xDivCon(z,seed) if nargin == 0 disp('------') for n = 4:20 disp(n) for seed = 0:20 rng(seed) z = rand(n,1) + rand(n,1)*1i; xDivCon(z,seed) end end return end [db,kb,jb] = Pairs(z); [d,k,j] = DivCon(z); if ~(d == db && k == kb && j == jb) n = length(z); disp([n seed kb jb k j]) fprintf('Pairs: %8.4f %6.4f+%6.4fi %6.4f %6.4fi\n', ... db,real(z(kb)),imag(z(kb)),real(z(jb)),imag(z(jb))) fprintf('DivCon: %8.4f %6.4f+%6.4fi %6.4f %6.4fi\n', ... d,real(z(k)),imag(z(k)),real(z(j)),imag(z(j))) end end %% EOF