pell(d,pr)=
{
  local(bb,b,len,m,x,y,z,v);

  if(pr>28,
     default(realprecision,pr)
  );

  v=[];

  bb=contfrac(sqrt(d));

  len=length(bb);
  for(i=1,len,
      b=vecextract(bb,2^i-1);
      m=contfracpnqn(b);  
      x=m[1,1];
      y=m[2,1];
      z=x^2-d*y^2;
      if(abs(z)<=4 && !setsearch(Set(v),z),
         print1([x,y,z], ", ");
         v=concat(v,z);
         if(z==1,return())
      )
  );
}