function siam2(T,ndig)

if nargin<2
 ndig = 100;
end

mp.Digits(ndig);
trem = mp(T);


p = [mp('0.5'); mp('1/10')];
pth = p;

v = [mp(1); mp(0)];
s = 0;
s = mp(s);
x0 = mp([0;0]);

options.Display = 'off';
t = mp(0);

while trem>0
   m = round(p+mp(2)*v/mp(3));
   s = fzero(@(t) dot(p+t*v-m,p+t*v-m)-mp(1)/mp(9),mp('0'),options);
   if s<0
     s = 1e10;
   end
   if s<trem
     disp('Reflection:');
     q = p+s*v;
     fprintf('Actual point: (%16.14f,%16.14f)\n',q(1),q(2));
     d = q-m;
     H = mmat(d(1),d(2));
     p = p+s*v;
     v = H*v;
   else
     s = min(trem,mp(2)/mp(3));
     p = p+s*v;
   end
   trem = trem-s;
   pth = [pth p];
end
fprintf('Final point: (%16.14f,%16.14f)\n',p(1),p(2));

xc = -1:1;
yc = -1:2;
clf
hold on
axis equal
dis = linspace(0,2*pi,200);
onec = 1/3*[cos(dis);sin(dis)];
for i=1:length(xc)
 for j=1:length(yc)
   plot(xc(i)+onec(1,:),yc(j)+onec(2,:),'LineWidth',2);
 end
end
plot(pth(1,:),pth(2,:),'-r','LineWidth',2);
gph = [];
for i = 1:size(pth,2)
  ap = pth(:,i);
  if abs(norm(ap-round(ap))-1/3)<1e-6
    plot(ap(1),ap(2),'ok');
  end
end

    plot(p(1),p(2),'.','markersize',20);

hold off
xs = xlim;
ys = ylim;
w = 0.01;
axis([xs(1)-w xs(2)+w ys(1)-w ys(2)+w]);

axis off
fprintf(['Result: %.',num2str(ndig),'f\n'],norm(pth(:,end)));

function H = mmat(a,b)

H = 9*[b^2-a^2 -2*a*b; -2*a*b a^2-b^2];
