/* Nagell's criterion for fundamental solution of u^2-Dv^2=n, n>0 */
define nagell(d,n){
auto bound,y,t,t1,x,e,f,start,nn,flag,x1,y1,tt,flag1,g
if(d<1){
       print "d=",d," is < 1\n"
       return(0)
}
x=sqrt(d)
g=x*x
if(d==g){
       print "d=",d," is a perfect square\n"
       return(0)
}
nn=abs(n)
tt=fund(d)
if(n==1){
   print "(1, 0) is a fundamental solution in the sense of Nagell\n"
   print" and defines an ambiguous class.\n"
   print "(",globalr,",", globals,") is the solution of Pell's equation x^2 - ",d,"y^2 = 1, with least positive y.\n"
   return;
}

x1=globalr
y1=globals
e=y1^2*nn
if(n>0){
   f=2*(x1+1)
}else{
   f=2*(x1-1)
}
t=e/f
bound=sqrt(t)
print "bound=",bound,"\n"
if(n>0){
  start=0
}else{
  start=1
}
flag1=0
for(y=start;y<=bound;y++){
   t1=d*y^2+n
   if(t1<0){
     continue
   }else{
       x=sqrt(t1)
   if(x^2==t1){
     flag1=1
     print "(",x,",",y,") is a fundamental solution\n"
     if(y==0){
        print" and defines an ambiguous class\n"
        flag1=0
     }
     if(y>0 && ((y^2*f==y1^2*nn) || x==0)){
        print" and defines an ambiguous class\n"
        flag1=0
     }
     if(flag1){
        print "(",-x,",",y,") is a fundamental solution\n"
        flag1=0
     }
   }
}
}
return
}

/* this finds the fundamental solution globalr+sqrt(d)globals  of Pell's equation x^2-dy^2=1
   using the midperiod approach. */
define fund(d){
	auto h,p,q,e,f,g,x,y,k,l,m,n,u,v,r,s,oldl,oldm,oldu,oldv,oldn
	x=sqrt(d)
	p=0;q=1
	g=x*x
	/*if(d==g+1){*//* period length = 1 */
        /*return(1)
	}*/
        l=0;k=1;m=1;n=0
	for(h=0;1;h++){
		y=(x+p)/q
                u=k*y+l;v=n*y+m
                       /* u/v is the i-th convergent to sqrt(d) */
 /*             if(h>=0){
                   print "A[",h,"]/B[",h,"]=",u,"/",v,"\n"
                }*/
                oldl=l;oldm=m;oldn=n
                oldu=u;oldv=v;
                l=k;m=n
                k=u;n=v
                
		f=p
		p=y*q-p
		e=q
		q=(d-(p*p))/q
		if(p==f){/* P_h=P_{h+1}, even period 2h */
                /* print "P_",h,"=P_",h+1,"\n"*/
                   globalr=oldn*(oldu+oldl)+(-1)^h
                   globals=oldn*(oldv+oldm)
                /* print "globalr=",globalr,",globals=",globals,"\n"*/
                   return(2*h)
		}
		if(q==e){/* Q_h=Q_{h+1}, odd period 2h+1 */
 /*                print "Q_",h,"=Q_",h+1,"\n"*/
                   r=k*n+l*m
                   s=n^2+m^2
                   globalr=r*r+d*s*s
                   globals=2*r*s
                /* print "globalr=",globalr,",globals=",globals,"\n" */
                   return(2*h+1)
		}
    }
}

/* absolute value of an integer n */
define abs(n){
	if(n>=0) return(n)
	return(-n)
}
define bcmul3(a,b,c){
auto temp
temp=a*b
temp=temp*c
return(temp)
}

/* This is a version of nagell(d,n) without printing. 
 * The H fundamental solutions, if any, are (nagellx[h],nagelly[h]) for h=0,,,,H-1. */
define nagell0(d,n){
auto bound,y,t,t1,x,e,f,start,nn,flag,x1,y1,tt,flag1,count
count=0
x=sqrt(d)
nn=abs(n)
tt=fund(d)
if(n==1){
   nagellx[0]=1
   nagelly[0]=0
   count=1
   return(count);
}

x1=globalr
y1=globals
e=y1^2*nn
if(n>0){
   f=2*(x1+1)
}else{
   f=2*(x1-1)
}
t=e/f
bound=sqrt(t)
/*print "bound=",bound,"\n"*/
if(n>0){
  start=0
}else{
  start=1
}
flag1=0
for(y=start;y<=bound;y++){
   t1=d*y^2+n
   if(t1<0){
     continue
   }else{
       x=sqrt(t1)
   if(x^2==t1){
     flag1=1
     nagellx[count]=x
     nagelly[count]=y
     count=count+1
     if(y==0){
        flag1=0
     }
     if(y>0 && ((y^2*f==y1^2*nn) || x==0)){
        flag1=0
     }
     if(flag1){
        nagellx[count]=-x
        nagelly[count]=y
        count=count+1
        flag1=0
     }
   }
}
}
return(count)
}

define test(d,n){
auto count,h
count=nagell0(d,n)
for(h=0;h<count;h++){
   print "(x,y)=(",nagellx[h],",",nagelly[h],")\n"
}
print "The number of fundamental solutions is "
return(count)
}

/* This uses the rather slow nagell0(d,n) to solve ax^2+bxy+cy^2+dx+ey+f=0,
 * where a is nonzero and b^2-4ac>0 is nonsauare.
 */
define ssw(a,b,c,d,e,f){
auto dd,ee,ff,n,g,h,t,u,twoa,count,x,y,xx,yy,temp,flag

count=0
dd=b^2-4*a*c
if(dd<=0){
   print "D<=0\n"
   return
}
ee=b*d-2*a*e
ff=d^2-4*a*f
n=ee^2-dd*ff
print "E=",ee,"\n"
print "F=",ff,"\n"

print "D=",dd,"\n"
print "N=",n,"\n"
if(n==0){
   print "n=0\n"
   return
}
twoa=2*a
print "X=",dd,"y+",ee,",Y=",twoa,"x+",b,"y+",d,"\n"
print "Find the fundamental solutions of X^2-",dd,"Y^2=",n,"\n"
g=nagell0(dd,n)
print "There are ",g," fundamental solutions (FundX,FundY)\n"
t=globalr
u=globals
print "(t,u)=(",t,",",u,")\n"
for(h=0;h<g;h++){
   flag=0
   x=nagellx[h]
   y=nagelly[h]
   if((x-ee)%dd && (t*x-ee)%dd){
      continue
   }
   if((x-b*y)%twoa){
      continue
   }
   if(((d*dd-b*ee-dd*y+b*x)/twoa)%dd==0 && (x-ee)%dd==0){
      yy=(x-ee)/dd
      xx=(y-b*yy-d)/twoa
   print "X=",x,", Y=",y,": "
      print "type 1: (x,y)=(",xx,",",yy,")\n"
      count=count+1
      flag=1
   } 
   temp=(d*dd-b*ee-t*(dd*y-b*x))/twoa
   if(temp%dd==0 && (t*x-ee)%dd==0){
      if(flag){
         print "type 1 & type 2\n"
         continue
      }
      yy=(x*t+y*u*dd-ee)/dd
      xx=(x*u+y*t-b*yy-d)/twoa
   print "X=",x,", Y=",y,": "
      print "type 2: (x,y)=(",xx,",",yy,")\n"
      count=count+1
   }
   continue
}
print "number of solutions (x,y) found is "
return(count)
}

/* This uses the fast patzgen(d,n) to solve ax^2+bxy+cy^2+dx+ey+f=0,
 * where a is nonzero and b^2-4ac>0 is nonsquare.
 */
define ssw0(a,b,c,d,e,f){
auto dd,ee,ff,n,g,h,t,u,twoa,gcount,x,y,xx1,yy1,xx2,yy2,temp,flag,r,flag1

if(a==0){
   print "a==0\n"
   return
}
gcount=0
dd=b^2-4*a*c
if(dd<=0){
   print "Here D<=0\n"
   return(-1)
}
ee=b*d-2*a*e
ff=d^2-4*a*f
n=ee^2-dd*ff
print "E=",ee,"\n"
print "F=",ff,"\n"

print "D=",dd,"\n"
g=sqrt(dd)
if(g*g==dd){
    print "D=",g,"^2\n"
    return(-1)
}
print "N=",n,"\n"
twoa=2*a
if(n==0){
   if(ee%dd){
      print "no solution\n"
      return(0)
   }else{
      y=(-ee)/dd
      if((b*y+d)%twoa){
          print "no solution\n"
          return(0)
      }else{
          x=-(b*y+d)/twoa
          print "(x,y)=(",x,",",y,") is the only solution\n"
          return(1)
      }
   }
}
print "X=",dd,"y+",ee,",Y=",twoa,"x+",b,"y+",d,"\n"
print "Find the fundamental solutions of X^2-",dd,"Y^2=",n,"\n"
g=patzgen(dd,n)
print "There are ",g," fundamental solutions (FundX,FundY)\n"
r=fund(dd)
t=globalr
u=globals
print "(t,u)=(",t,",",u,")\n"
for(h=0;h<g;h++){
   flag=0
   flag1=0
   x=globalfundxff[h]
   y=globalfundyff[h]
   if((x-ee)%dd && (t*x-ee)%dd){
      continue
   }
   if((x-b*y)%twoa){
      continue
   }
   if(((d*dd-b*ee-dd*y+b*x)/twoa)%dd==0 && (x-ee)%dd==0){
      yy1=(x-ee)/dd
      xx1=(y-b*yy1-d)/twoa
      print "FundX=",x,", FundY=",y,": "
      /*print "type 1: (x,y)=(",xx1,",",yy1,")\n"
      print "X+Y\\sqrt(",dd,")=+-(",x,"+",y,"\\sqrt(",dd,")(",t,"+",u,"\\sqrt(",dd,")^2n\n"*/
      gcount=gcount+1
      flag=1
   } 
   temp=(d*dd-b*ee-t*(dd*y-b*x))/twoa
   if(temp%dd==0 && (t*x-ee)%dd==0){
      yy2=(x*t+y*u*dd-ee)/dd
      xx2=(x*u+y*t-b*yy2-d)/twoa
      print "FundX=",x,", FundY=",y,": "
      /*print "type 2: (x,y)=(",xx2,",",yy2,")\n"
      print "X+Y\\sqrt(",dd,")=+-(",x,"+",y,"\\sqrt(",dd,")(",t,"+",u,"\\sqrt(",dd,")^(2n+1)\n"*/
      gcount=gcount+1
      flag1=1
   }
   if(flag && flag1){
      print "X+Y\\sqrt(",dd,")=+-(",x,"+",y,"\\sqrt(",dd,")(",t,"+",u,"\\sqrt(",dd,")^n, "
      print "(x,y)=(",xx1,",",yy1,")\n"
   }
   if(flag && flag1==0){
      print "X+Y\\sqrt(",dd,")=+-(",x,"+",y,"\\sqrt(",dd,")(",t,"+",u,"\\sqrt(",dd,")^2n, "
      print "(x,y)=(",xx1,",",yy1,")\n"
   }
   if(flag==0 && flag1){
      print "X+Y\\sqrt(",dd,")=+-(",x,"+",y,"\\sqrt(",dd,")(",t,"+",u,"\\sqrt(",dd,")^(2n+1), "
      print "(x,y)=(",xx2,",",yy2,")\n"
   }
}
print "number of solutions (x,y) found is "
return(gcount)
}