function [px py]=interpfrac(x,y,d,npoints)
  n=length(x);
  b=x(n)-x(1);

  for i=2:n
    aw(i-1)=(x(i)-x(i-1))/b;
    ew(i-1)=(x(n)*x(i-1)-x(1)*x(i))/b;
    cw(i-1)=(y(i)-y(i-1)-d(i)*(y(n)-y(1)))/b;
    fw(i-1)=(x(n)*y(i-1)-x(1)*y(i)-d(i)*(x(n)*y(1)-x(1)*y(n)))/b;
  endfor

  oldx=0;
  oldy=0;
  px=oldx;
  py=oldy;
  for j=1:npoints
    k=floor((n-1)*rand)+1;
    newx=aw(k)*oldx+ew(k);
    newy=cw(k)*oldx+d(k)*oldy+fw(k);
    oldx=newx;
    oldy=newy;
    px=[px; newx];
    py=[py; newy];
  endfor
endfunction