% This code draws a 2-D grid walk with N steps 
% and then searches for short paths that connect 

function paths5(N)
tic;
% set number of steps in the random walk
%N=1000000;
t=randi(4,1,N);
u=round(cos(.5*pi*t));
v=round(sin(.5*pi*t));
% define the random walk
x=cumsum(u)-u(1);
y=cumsum(v)-v(1);

% find min and max value that path reaches
minx=min(x);
maxx=max(x);
miny=min(y);
maxy=max(y);
% compute size of smallest grid containing path
xrange=maxx-minx;
yrange=maxy-miny;
% define an array big enough to hold path
% plus a empty `frame' around path
visited=zeros(xrange+4,yrange+4);

% redefine x and y vectors to lie in grid
% but not touch outer most sites
newx=x-minx+2;
newy=y-miny+2;

for k=1:N
     visited(newx(k),newy(k))=1;
end % for k

% define distance to starting point
% if distance is already define do nothing,
% otherwise distance is one more that
% the minimum of four surrounding sites
% maxpath keeps track of how far from
% origin we can get in k steps

grid=zeros(xrange+4,yrange+4);
checked=zeros(xrange+4,yrange+4);
upper=1000000;
maxpath=0;
grid=grid+upper;
x1=newx(1);
y1=newy(1);
grid(x1,y1)=0;
checked(x1,y1)=1;
% listx and listy contain coordinates of set of points
% distance k from starting point
clear listx;
clear listy;
clear newlistx;
clear newlisty;
newlistx(1)=x1;
newlisty(1)=y1;
listlength=1;
flag=1;
dist=0;
while flag==1
   dist=dist+1;
   counter=0;
   listlength=length(newlistx);
   listx=newlistx;
   listy=newlisty;
   clear newlistx;
   clear newlisty;
   for k=1:listlength
       x1=listx(k);
       y1=listy(k);
       % check each neighbor if in trace and not checked
       nextx=x1+1;
       nexty=y1;
       if (checked(nextx,nexty)==0)&&(visited(nextx,nexty)==1)
         % add to newlist, record distance and check
         counter=counter+1;
         newlistx(counter)=nextx;
         newlisty(counter)=nexty;
         grid(nextx,nexty)=dist;
         checked(nextx,nexty)=1;
       end % if
       nextx=x1-1; nexty=y1;
       if (checked(nextx,nexty)==0)&&(visited(nextx,nexty)==1)
         % add to newlist, record distance and check
         counter=counter+1;
         newlistx(counter)=nextx;
         newlisty(counter)=nexty;
         grid(nextx,nexty)=dist;
         checked(nextx,nexty)=1;
       end % if
       nextx=x1; nexty=y1+1;
       if (checked(nextx,nexty)==0)&&(visited(nextx,nexty)==1)
         % add to newlist, record distance and check
         counter=counter+1;
         newlistx(counter)=nextx;
         newlisty(counter)=nexty;
         grid(nextx,nexty)=dist;
         checked(nextx,nexty)=1;
       end % if
       nextx=x1; nexty=y1-1;
       if (checked(nextx,nexty)==0)&&(visited(nextx,nexty)==1)
         % add to newlist, record distance and check
         counter=counter+1;
         newlistx(counter)=nextx;
         newlisty(counter)=nexty;
         grid(nextx,nexty)=dist;
         checked(nextx,nexty)=1;
       end % if
   end % for k
   if counter==0
      flag=0;
   end % if
end % while


maxpath=0;
for k=1:N
    if grid(newx(k),newy(k)) > maxpath
       maxpath=grid(newx(k),newy(k));
       maxx=newx(k);
       maxy=newy(k);
    end  % if grid
end % for k

% find path between origin and maximal distance path

   clear savex;
   clear savey;
   oldx=maxx;
   oldy=maxy;
   savex(1)=oldx;
   savey(1)=oldy;
   for k=1:maxpath
      dist=maxpath-k;
      if   grid(oldx-1,oldy)==dist
           nextx=oldx-1;
           nexty=oldy;
      end % if
      if   grid(oldx+1,oldy)==dist
           nextx=oldx+1;
           nexty=oldy;
      end % if
      if   grid(oldx,oldy-1)==dist
           nextx=oldx;
           nexty=oldy-1;
      end % if
      if   grid(oldx,oldy+1)==dist
           nextx=oldx;
           nexty=oldy+1;
      end % if
      savex(k+1)=nextx;
      savey(k+1)=nexty;
      oldx=nextx;
      oldy=nexty;
   end % for k
   figure;
   axis equal;
   hold on;
   set(gca,'visible','off');
   plot(newx,newy);

   figure;
   axis equal;
   hold on;
   set(gca,'visible','off');
   plot(newx,newy);
   plot(savex,savey,'k','LineWidth',3);
   toc;
   length(savex)
   log(length(savex)/
return

