\\
\\ solve Pell's Equation |x^2-Dy^2|=1
\\     D is square-free
\\

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

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

  bb=contfrac(sqrt(d));

  len=length(bb);
  for(i=1,len,
      b=vecextract(bb,2^i-1);
\\      print(b);
      m=contfracpnqn(b);  
\\      print("m=",m);
      x=m[1,1];
      y=m[2,1];
      z=x^2-d*y^2;
\\      print("[x,y]=",[x,y,z]);
      if(z==1,
\\       print("i=",i);
\\       print("x=",x);
\\       print("y=",y);
\\       print("x^2-",d,"*y^2=",z);
         return([x,y,z])
      );
      if(z==-1,
         x2=x*x+d*y*y;
         y2=2*x*y;
         return([x2,y2,1])
      )
  );
  return([])
}

\\
\\ find (x,y) \in N
\\ x^2-Dy^2=-1 or 1
\\
pell0(d,pr)=
{
  local(bb,b,len,m,x,y,z);

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

  bb=contfrac(sqrt(d));

  len=length(bb);
  for(i=1,len,
      b=vecextract(bb,2^i-1);
\\      print(b);
      m=contfracpnqn(b);  
\\      print("m=",m);
      x=m[1,1];
      y=m[2,1];
      z=x^2-d*y^2;
\\      print("[x,y]=",[x,y,z]);
      if(abs(z)==1,
\\       print("i=",i);
\\       print("x=",x);
\\       print("y=",y);
\\       print("x^2-",d,"*y^2=",z);
         return([x,y,z])
      )
  );
  return([])
}

findpell(n1,n2,pr)=
{
  local(d2,pp,pre);

  for(d=n1,n2,
      d2=sqrtint(d);
      if(d2*d2<d,
         pre=pr;
         pp=pell0(d,pre);
         while(pp==[],
               pre=2*pre;
               pp=pell0(d,pre)
         );
         print("d=",d,"; ",pp)
     )
  )
}

findpellbig(n1,n2,pr)=
{
  local(d2);

  for(d=n1,n2,
      d2=sqrtint(d);
      if(d2*d2<d && pell0(d,pr)==[],
        print([d])
      )
  )
}

\\
\\ find (x,y) \in N
\\  x^4-Dy^2=1
\\     1 < n1 <= D <= n2
\\
findpellx2(n1,n2,pr)=
{
  local(d2,pp,pre,x,y,z,x1,y1,x2,y2);

  for(d=n1,n2,
      d2=sqrtint(d);
      if(d2*d2<d,
         pre=pr;
         pp=pell(d,pre);
         while(pp==[],
               pre=2*pre;
               pp=pell(d,pre)
         );
         if(pp!=[],
            x=pp[1];
            y=pp[2];
            z=pp[3];
            if(z==-1,
               x2=x*x+d*y*y;
               y2=2*x*y;
               x=x2;
               y=y2
            );
\\ x+y*sqrt(d)=e_d
            x1=sqrtint(x);
            if(x1*x1==x,
               print("x^4-",d,"y^2=1; ",[x1,y])
            );
\\ x+y*sqrt(d)=(e_d)^2
            x2=x*x+d*y*y;
            y2=2*x*y;
            x=x2;
            y=y2;
            x1=sqrtint(x);
            if(x1*x1==x,
               print("x^4-",d,"y^2=1; ",[x1,y])
            )
         )
      )
  )
}

\\
\\ find (x,y) \in N
\\  x^2-Dy^4=1
\\     1 < n1 <= D <= n2
\\
findpelly2(n1,n2,pr)=
{
  local(d2,pp,pre,x,y,z,x1,y1,x2,y2);

  for(d=n1,n2,
      d2=sqrtint(d);
      if(d2*d2<d,
         pre=pr;
         pp=pell(d,pre);
         while(pp==[],
               pre=2*pre;
               pp=pell(d,pre)
         );
         if(pp!=[],
\\ fundamental solution e_d
            x=pp[1];
            y=pp[2];
            z=pp[3];
\\            print(pp);
            if(z==-1,
               x2=x*x+d*y*y;
               y2=2*x*y;
               x=x2;
               y=y2
            );
\\ x+y*sqrt(d)=e_d
\\            print([x,y]);
            y1=sqrtint(y);
            if(y1*y1==y,
               print("x^2-",d,"y^4=1; ",[x,y1])
            );
\\ x+y*sqrt(d)=(e_d)^4
            x2=x*x+d*y*y;
            y2=2*x*y;
            x=x2*x2+d*y2*y2;
            y=2*x2*y2;
            y1=sqrtint(y);
            if(y1*y1==y,
               print("x^2-",d,"y^4=1; ",[x,y1])
            )
         )
      )
  )
}

\\
\\ find (x,y) \in N
\\  x^4-Dy^4=1
\\     1 < n1 <= D <= n2
\\
findpellx2y2(n1,n2,pr)=
{
  local(d2,pp,pre,x,y,z,x1,y1,x2,y2);

  for(d=n1,n2,
      d2=sqrtint(d);
      if(d2*d2<d,
         pre=pr;
         pp=pell(d,pre);
         while(pp==[],
               pre=2*pre;
               pp=pell(d,pre)
         );
         if(pp!=[],
\\ fundamental solution e_d
            x=pp[1];
            y=pp[2];
            z=pp[3];
\\            print(pp);
            if(z==-1,
               x2=x*x+d*y*y;
               y2=2*x*y;
               x=x2;
               y=y2
            );
\\ x+y*sqrt(d)=e_d
\\            print([x,y]);
            x1=sqrtint(x);
            y1=sqrtint(y);
            if(x1*x1==x && y1*y1==y,
               print("x^4-",d,"y^4=1; ",[x1,y1])
            )
         )
      )
  )
}

\\
\\ findthue(x,y) \in N
\\  x^4-Dy^4=1 or -1
\\   D: positive integer, D != m^4.
\\    1 < n1<= D <= n2 .
\\    
findthue(d)=
{

  local(th,d2,pp);

  d2=sqrtint(sqrtint(d));
  if(d2^4<d && matsize(factor(x^4-d))[1]==1,
     th=thueinit(x^4-d);
     pp=thue(th,1);
     if(pp!=[] && pp!=[[1,0],[-1,0]],
        print("x^4-",d,"y^4=1; ",pp)
\\      ,print("d=",d)
     );
     pp=thue(th,-1);
     if(pp!=[],
        print("x^4-",d,"y^4=-1; ",pp)
\\      ,print("d=",d)
     )
  )
}