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

  default(realprecision,2000);

  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())
      )
   );
}

solvepell(d) = if(d==1,print1("[1, 1, 0], [0, 1, -1], [2, 1, 3], [1, 2, -3], [0, 2, -4],"),if(d==2,print1("[1, 1, -1], [0, 1, -2], [2, 1, 2], [3, 2, 1], [2, 2, -4], [6, 4, 4],"),if(d==3,print1("[1, 1, -2], [0, 1, -3], [2, 1, 1], [4, 2, 4],"),if(d==4,print1("[2, 1, 0], [1, 1, -3],"),if(d==5,print1("[2, 1, -1], [1, 1, -4], [3, 1, 4], [9, 4, 1],"),if(d==6,print1("[2, 1, -2], [3, 1, 3], [5, 2, 1], [10, 4, 4],"),if(d==12,print1("[3, 1, -3], [4, 1, 4], [7, 2, 1],"),pell(d))))))))