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

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

  v=[];
  r=[];
  s=[];

  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,r=concat(r,[2*x,2*y]));
         if(z==1,s=concat(s,[2*x,2*y]);
      if(!setsearch(Set(v),-4) && setsearch(Set(v),-1),print1([r[1],r[2],-4], ", "));
      if(!setsearch(Set(v),4),print1([s[1],s[2],4], ", "));
      return())
      )
   );
}