function d(P0,P1)
   a := (P0[1]-P1[1])^2 + (P0[2]-P1[2])^2;
   return Sqrt(RealField()!a);
end function;

intrinsic GraphComponent(C::Crv, P0::Pt, t::FldPrElt, max::RngIntElt) -> SeqEnum
{}
   AssertAttribute(FldPr, "OutputPrecision", 3);
   ans := [P0];
   P := P0;
   A := AmbientSpace(C);
   F := DefiningEquation(C);
   R<x,y> := Parent(F);
   R<Z> := PolynomialRing(ComplexField());

   Fx := Derivative(F,x);
   Fy := Derivative(F,y);
   last_dir_x := 0;
   last_dir_y := 0;
   while #ans lt max and (#ans lt 5 or d(P0,P) gt 1.1*t) do
      print "#ans = ",#ans;
      FxP := Evaluate(Fx,Eltseq(P));
      FyP := Evaluate(Fy,Eltseq(P));
//      print "FxP =", FxP, "   FyP =", FyP;
      if FyP eq 0 then
         if last_dir_x gt 0 then
            Pnew := A![P[1],P[2]+t];
         elif last_dir_x lt 0 then
            Pnew := A![P[1],P[2]-t];
         else
            error "Confused by vertical tangent!";
         end if;
      else
         m := -FxP/FyP;   // the slope
         dir_x := t/Sqrt(1+m^2);
         dir_y := m*dir_x;
//         print "m = ", m, "  dir = ", [dir_x, dir_y];
//         print "x^2 + y^2 = ", dir_x^2 + dir_y^2;
         if d([dir_x,dir_y],[last_dir_x,last_dir_y]) gt t/2 then
            dir_x := -dir_x;
            dir_y := -dir_y;
         end if;
         if Abs(dir_x) ge Abs(dir_y) then
            x_new := P[1]+dir_x;
            g := Evaluate(F,[x_new,Z]);
            y_new := Sort([<Abs(RealField()!(Real(r[1])-(P[2]+dir_y))),Real(r[1])> :
                                              r in Roots(g)])[1][2];
         else 
            y_new := P[2]+dir_y;
            g := Evaluate(F,[Z,y_new]);
            x_new := Sort([<Abs(RealField()!(Real(r[1])-(P[1]+dir_x))),Real(r[1])> :
                                              r in Roots(g)])[1][2];
         end if;   
         Pnew := A![x_new, y_new];
      end if;
      last_dir_x := Pnew[1]-P[1];
      last_dir_y := Pnew[2]-P[2];
      P := Pnew;
      Append(~ans,P);
      print P;
   end while;
   return ans cat [P0];
end intrinsic;

intrinsic pscon(g::SeqEnum) -> MonStgElt
{}
   ans := "";
   for P in g do
      ans := ans*Sprintf("(%o,%o)\n",P[1],P[2]);
   end for;
   return ans;
end intrinsic;
