# gnubc program : patz 
# Needs squareroot 
# to run sswgeneral, one needs 
# patz squareroot2 powertest3 genfacs gcd posformrep reduceneg congruence powerdd ddzero-extra arithpartition.bc rob1 pell1

/* Typing pn(a[],n) returns p[n], where
 * p[0]=a[0],p[1]=a[0]*a[1]+1,p[i+1]=a[i+1]*p[i]+p[i-1] if i >= 1.
 * Typing qn(a[],n) returns q[n], where
 * q[0]=1,q[1]=a[1],q[i+1]=a[i+1]*q[i]+q[i-1] if i >= 1.
 * Hence pn/qn is the value of the simple continued fraction
 * [a[0];a[1],...,a[n]].
 */
define pn(a[],n){
	auto x,y,z,i
        if (n == 0)return(a[0])
        if(n==1)return(a[0]*a[1]+1)
        x=a[0];y=a[0]*a[1]+1
        for(i=2;i<=n;i++){
          z=a[i]*y+x
	  x=y
	  y=z
        }
        return(z)
}

define qn(a[],n){
        auto x,y,z,i
        if(n==0)return(1)
        if(n==1)return(a[1])
        x=1;y=a[1]
        for(i=2;i<=n;i++){
	    z=a[i]*y+x
	    x=y
	    y=z
        }
        return(z)
}

/*
 * This program finds the continued fraction expansion of a
 * quadratic irrational (u+tsqrt(d))/v, where d is not a square,
 * t,u,v integers, v nonzero.
 * After converting to (U+sqrt(D))/V, where V divides D-U^2,
 */
/* sign of an integer a */
/* sign(a)=1,-1,0, according as a>0,a<0,a=0 */

define sign(a){
	if(a>0) return(1)
	if(a<0) return(-1)
	return(0)
}

/* absolute value of an integer n */

define abs(n){
	if(n>=0) return(n)
	return(-n)
}

/*  mod(a,b)=the least non-negative remainder when an integer a is divided by a
 positive integer b */

define mod(a,b){
	auto c
	c=a%b
	if(a>=0) return(c)
	if(c==0) return(0)
        return(c+b)	
}

/* int(a,b)=integer part of a/b, a, b integers, b != 0 */

define int(a,b){
	auto c
	c=sign(b)
	a=a*c
	b=b*c
	return((a-mod(a,b))/b)
}

/*
 * This is a function for finding the period of the continued fraction
 * expansion of reduced quadratic irrational a=(u+sqrt(d))/v.
 * Here d is non-square, 1<(u+sqrt(d))/v, -1<(u-sqrt(d))/v<0.
 * The algorithm also assumes that v divides d-u*u and is based on K. Rosen,
 * Elementary Number theory and its applications, p.379-381 and Knuth's The art
 * of computer programming, Vol. 2, p. 359.  0 is returned if a is not reduced.
 */
define period(d,u,v,i){
	auto r,s,t,f,j,k,temp
	f=sqrt(d)
	if(u<=0){"a is not reduced: u<=0:";return(0)}
	if(v<=0){"a is not reduced: v<=0:";return(0)}
	if(u+f<v){"a is not reduced: a<1:";return(0)}
	if(u+v<= f){"a is not reduced: a'<=-1:";return(0)}
	if(f<u){"a is not reduced: a'>0:";return(0)}
	k=d-u*u
	if(k%v!=0){"v does not divide d - u*u:";return(0)}
	s=v
	r=u
	a_surd_length=i
/* i is created by surd(d,t,u,v) below and indexes the ith 
   convergent of (u+t*sqrt(d))/v */
	for(j=i;1;j++){
	   a=(f+u)/v
           /*  print "a[",j,"]=",a,"\n" */
	   a_surd[j]=a
	   u=a*v-u
	   temp=v
	   v=(d-u*u)/v
           u_surd[j+1]=u
           v_surd[j+1]=v
           if(u==r && v==s){
             t=j+1-i
	     a_surd_length=a_surd_length+t
             return(t)
           /* the continued fraction has period length j+1-i */
           }
	}
}

/*
 * This function uses the continued fraction algorithm expansion in K. Rosen,
 * Elementary Number theory and its applications,p.379-381 and Knuth's
 * The art of computer programming, Vol. 2, p. 359. It locates the first reduced
 * complete quotient and then uses the function period(d,u,v,i,flag)
 * to locate the period of the continued fraction expansion.
 */

define surd(d,t,u,v,print_flag){
	auto a,c,e,g,f,i,j,k,w,z,con_length
	c=d;g=u;e=v
	if(v==0){
           print "v="
           return(0)
        }
	f=sqrt(d)
	if(f*f==d){
           print "d is the square of "
           return(f)
        }
	if(t==0){
           print "t = "
           return(0)
        }
	z=sign(t)
	u=u*z
	v=v*z
	d=d*t*t
	w=abs(v)
	if((d-u*u)%w!=0){
           d=d*v*v
	   u=u*w
	   v=v*w
	}
	f=sqrt(d)
	for(j=0;1;j++){
	    u_surd[j]=u
	    v_surd[j]=v
            g=reducetest(u,v,d)
            if(g){/* start if(g) */ /* (u+sqrt(d))/v is reduced */
              i=j
              pre_period_length=i
              l=period(d,u,v,i)
              if(l%2){
                 l=period(d,u,v,j+l)
                 if(print_flag){
                    print "printing first and second period:\n"
                 }
              }
              con_length=a_surd_length-1
              for(k=0;k<=con_length;k++){
                  p_surd[k]=pn(a_surd[],k)
                  q_surd[k]=qn(a_surd[],k)
              }
              if(print_flag){
                 for(k=0;k<=con_length;k++){
                 }
                 print "\n"
                 for(k=0;k<=con_length;k++){
                     if(k==pre_period_length || l%2 && k==pre_period_length+l){
                        print "\n period:\n"
                     }
                     print "(P[",k,"]+sqrt(",d,"))/Q[",k,"]=(",u_surd[k],"+sqrt(",d,"))/",v_surd[k],", "
                     print "A[",k,"]/B[",k,"]=",p_surd[k],"/",q_surd[k],", "
                     print "a[",k,"]=",a_surd[k],"\n"
                 }
                 print "cfrac has period length ",l,"\n"
              }
              return(l)
           }/* if(g) */
           if(v>0){
              a=int(f+u,v)
           }
           if(v<0){
              a=int(f+u+1,v)
           }
           /*print "a[",j,"]=",a,"\n" */
           a_surd[j]=a
           u=a*v-u
           v=int(d-u*u,v)
        }
}

/* This finds fundamental solutions (globalfundx[h], globalfundy[h]) for the diophantine equation
 * x^2-dy^2=n, d>0, d not a perfect square n nonzero and not equal to 1. 
 * Returns the number of solution classes.
 */
define patz(d,n,print_flag){
	auto a,b,k,q0,s,t,period_length,temp,g1_surd,g2_surd,f1,f2,solution_number,pos,pos1,pos2,flagpos,g
        solution_number=0

        if(d<=1){
           print "d<=1\n"
           return(-1)
        } 
        g=sqrt(d)
        if(g^2==d){
           print "d is a perfect square\n"
           return(-1)
        }
        if(n==1){
           globalfundx[0]=1
           globalfundy[0]=0
           return(1)
        }
	q0=abs(n)
        s=sqroot(d,q0,0)
        /* s=0 means no solutions of x^2=d (mod |n|) */
        /* If s>0, we get all solutions x, 0<=x<=|n|/2 in the form
         * solution[0],...,solution[numbr-1]
         */
	if(s==0){
           if(print_flag){
	      print "x^2=",d," mod(",n,") has no solution\n"
           }
           return(solution_number)
	}
	/* Now to test the cfrac of each of w[j]=(-solution[j]+\sqrt{D})/q0
           to see if v_surd[j]=(-1)^j*sign(n) holds for some j in
           1<=j<=t+l, where l is the period of 
           w[j]=[a_surd[0],...,a_surd[t];a_surd[t+1],...,a_surd[t+l]]
           and t=pre_period_length.
         */
	for(k=0;k<numbr;k++){
            if(print_flag){
               print "processing omega[",k,"]="
               print "(",-solution[k],"+sqrt(",d,")/",q0,"\n"
            }
            period_length=surd(d,1,-solution[k],q0,print_flag)
            temp=pre_period_length+period_length
            if(period_length%2){
               temp=temp+period_length
            }
            for(a=1;a<=temp;a++){
                temp1=(-1)^a*sign(n)
                if(v_surd[a]==temp1){
                   f1=q_surd[a-1]
                   g1_surd=q0*p_surd[a-1]+solution[k]*q_surd[a-1]
                   if(print_flag){
                      print "solution: class P=",solution[k],"\n"
                      print " (x,y)=(",g1_surd,",",q_surd[a-1],")\n"
                   }
                   break
                }
            }
            if(a==temp+1){
               if(print_flag){
                  print "P=",solution[k]," yields no solution\n"
               }
               continue /* process solution[k+1] */
	    }
            if(solution[k]==0){
               if(print_flag){
                  print "Fundamental solution for x^2-",d,"y^2=",n,": ambiguous class P=0, "
               }
               pos=abs(g1_surd)
               if(print_flag){
                  print "(x,y)=(",pos,",",q_surd[a-1],")\n"
               }
               globalfundx[solution_number]=pos
               globalfundy[solution_number]=q_surd[a-1]
               solution_number=solution_number+1
            }
            if(solution[k]!=0){
               if(print_flag){
                  print "processing omega*[",k,"]=(",-solution[k],"-sqrt(",d,")/",q0,"\n"
               }
               period_length=surd(d,-1,-solution[k],q0,print_flag)
               temp=pre_period_length+period_length
               if(period_length%2){
                  temp=temp+period_length
               }
	       for(a=1;a<=temp;a++){
                   temp1=(-1)^(a+1)*sign(n)
		   if(v_surd[a]==temp1){
                      f2=q_surd[a-1]
                      g2_surd=q0*p_surd[a-1]+solution[k]*q_surd[a-1]
                      if(print_flag){
                         print "solution: class P=",solution[k]
                         print ", (x,y)=(",g2_surd,",",q_surd[a-1],")\n"
                      }
                      break
                   }
               }
               flagpos=0
               if(2*solution[k]==q0){
                  flagpos=1
                  if(print_flag){
	             print "Fundamental solution for x^2-",d,"y^2=",n,": "
                     print "ambiguous class, "
                  }
               }else{
                  if(print_flag){
	             print "Fundamental solutions for x^2-",d,"y^2=",n,": "
                     print "class "
                  }
               }
               pos1=abs(g1_surd)
               pos2=abs(g2_surd)
               if(print_flag){
                  print "P=",solution[k],", (x,y)=("
               }
               if(flagpos==1){
                  if(f1<=f2){
                     if(print_flag){
                        print pos1,",",f1
                     }
                     globalfundx[solution_number]=pos1
                     globalfundy[solution_number]=f1
                  }else{
                     globalfundx[solution_number]=pos2
                     globalfundy[solution_number]=f2
                     if(print_flag){
                        print pos2,",",f2
                     }
                  }
                  solution_number=solution_number+1
               }else{
                  if(f1<=f2){
                     if(print_flag){
                        print "+-",pos1,",",f1
                     }
                     globalfundx[solution_number]=pos1
                     globalfundy[solution_number]=f1
                     globalfundx[solution_number+1]=-pos1
                     globalfundy[solution_number+1]=f1
                  }else{
                     if(print_flag){
                        print "+-",pos2,",",f2
                     }
                     globalfundx[solution_number]=pos2
                     globalfundy[solution_number]=f2
                     globalfundx[solution_number+1]=-pos2
                     globalfundy[solution_number+1]=f2
                  }
                  solution_number=solution_number+2
               }
               if(print_flag){
                  print ")\n"
               }
            }
	}/* end of for k loop */
        return(solution_number)
}

/* calculate_n(a,b,c,alpha,gamma,n) calulates the n on page 259, equation (2.5)
 * of the paper below.
 * It is only used for printing.
 */
define calculate_n(a,b,c,alpha,gamma,n){
auto beta,delta,q0,q00,m,tmp1,tmp2,tmp5,tmp6
        q0=abs(n)
	q00=2*q0
	beta=-gcd2(alpha,gamma)
	delta= gcd1(alpha,gamma)
   /*     tmp1=((2*a)*alpha+b*gamma)*beta
	tmp2=(b*alpha+(2*c)*gamma)*delta
    */
        tmp1=(2*a)*alpha*beta
        tmp2=2*c*gamma*delta
        tmp5=b*alpha*delta
        tmp6=b*beta*gamma
        m=tmp1+tmp2+tmp5+tmp6
	m=mod(m,q00)
        if(m>q0){
           m=m-q00
        }
        return(m)
}

/* binary0(a,b,c,n,print_flag) tests the diophantine equation
 * ax^2+bxy+cy^2=n for solubility. Here d=b^2-4ac>0 and is not a perfect
 * square. Also n is non-zero, gcd(a,b,c)=1 and gcd(a,n)=1. 
 * In the case of solubility, the solution (x,y) in each class, 
 * with least y is obtained.
 * With print_flag =1, we get verbose output.
 * See K.R. Matthews, 'The Diophantine equation ax^2+bxy+cy^2=N, 
 * D=b^2-4ac>0', J. de The'orie des Nombres de Bordeaux 14 (2002) 257-270.
 * The test for necessity is simplified by the observation that 
 * if Q_n=2(-1)^nN/|N| occurs in the cfrac of omega=(-n+sqrt(d))/q, 
 * then Q_m=2(-1)^{m+1}N/|N| occurs in the cfrac of omega*=(-n-sqrt(d))/q.
 * creates global fund_x[] and fund_y[].
 */

define binary0(a,b,c,n,print_flag){
auto no_of_solutions, d,k,m,mm,q,q0,q00,s,t,temp,x,y,ee,qdash,ndash,rdash,theta,result

        if(print_flag){
          print "Solving "
          result=printform(a,b,c,n)
        }
        if(a==n){
           fund_x[0]=1
           fund_y[0]=0
           mm=calculate_n(a,b,c,1,0,n)
           if(print_flag){
              print "class n=",mm," (mod 2)\n"
              print "solution (1,0)\n"
           }
           return(1)
        }
	d=b^2-4*a*c
	q0=abs(n)
        q00=2*q0
        q=2*a*q0
        s=quadratic(a,b,c,q0,0)
        /* s=0 means no solutions of ax^2+bx+c=0 (mod |n|) */
        /* If s>0, we get all solutions x, 0<=x<|n| in the form
         * quadratic_solution[0],...,quadratic_solution[s-1]
         */
	if(s==0){
           if(print_flag){
	      print "ax^2+bx+c=0 mod(n) has no solution\n"
           }
           return(0)
	}
	/* With m=2*a*quadratic_solution[k]+b, to test the cfrac of each of 
           w[j]=(-m+\sqrt{d})/q
           to see if v_surd[j]=(-1)^j*2**sign(n) holds for some j in
           1<=j<=t+l, where l is the period of 
           w[j]=[a_surd[0],...,a_surd[t];a_surd[t+1],...,a_surd[t+l]]
           and t=pre_period_length.
         */
        no_of_solutions=0
	for(k=0;k<s;k++){
            theta=quadratic_solution[k]
	    m=2*a*theta+b
            ee=n/q0
            qdash=ee*a*q0
            ndash=ee*m
            rdash=ee*(a*theta^2+b*theta+c)/q0
            if(print_flag){
               print "theta=",theta,": solving ",qdash,"X^2+(",ndash,")XY+(",rdash,")Y^2=1\n"
            }
            t=generalizedpell(qdash,ndash,rdash,print_flag)
            if(t){
                 x=generalizedpellfundy*theta+generalizedpellfundx*q0
                 y=generalizedpellfundy
                 fund_x[no_of_solutions]=x
                 fund_y[no_of_solutions]=y
                 if(print_flag){
                    print "least nonnegative y solution of ",qdash,"x^2+(",ndash,")xy+(",rdash,")y^2=",n," is ("  
                    print x,",",y,")\n"
                    mm=calculate_n(a,b,c,x,y,n)
                    print "solution: class n=",mm," (mod ",q00,")\n"
                 }
                 no_of_solutions=no_of_solutions+1
            }
	}
        return(no_of_solutions)
}

/* binary1(a,b,c,n,print_flag) deals with the more general case of solving
 * ax^2+bxy+cy^2=n, where gcd(a,n) may be >1. Again we assume d=b^2-4ac>0
 * and that d is not a square. Also we assume gcd(a,b,c)=1.
 * The number r of solution classes is returned.
 * verbose output if print_flag=1.
 * Printing of solution classes is suppressed if gcd(a,n)>1.
 */ 
define binary1(a,b,c,n,print_flag){
        auto aa,bb,cc,i,q00,mm,r,s,tmp1,tmp2,tmp3,tmp4,result
        if(print_flag){
          print "Solving "
          result=printform(a,b,c,n)
        }
	if(gcd(a,n)==1){
           r=binary0(a,b,c,n,print_flag)
           return(r)
	}
        /* Now we calculate aa,bb,cc, where ax^2+bxy+cy^2 is 
         * transformed into aaX^2+bbXY+ccY^2 under the transformation
         * x=alpha*X+beta*Y, y=gamma*X+delta*Y.
         */
         q00=2*abs(n)
         s=gauss(a,b,c,n)
	 aa=gauss_m
         gauss_delta=gcd1(gauss_alpha,gauss_gamma)
         gauss_beta= -gcd2(gauss_alpha,gauss_gamma)
         tmp1=a*gauss_beta*gauss_beta
         tmp2=b*gauss_beta*gauss_delta
         tmp3=c*gauss_delta*gauss_delta
         tmp4=tmp1+tmp2
         cc=tmp4+tmp3
         tmp1=2*a*gauss_alpha*gauss_beta
         tmp2=2*c*gauss_gamma*gauss_delta
         tmp3=b*gauss_alpha*gauss_delta
         tmp4=b*gauss_beta*gauss_gamma
         bb=tmp1+tmp2+tmp3+tmp4
         r=binary0(aa,bb,cc,n,print_flag)
         /* replace 0 by 1 to print fundamental solutions for reduced
          * case gcd(aa,n)=1.
          */
         if(r==0){
           return(0)
         }
         /* this creates arrays of fundamental solutions:
          * fund_x[0],...,fund_x[r-1] and
          * fund_y[0],...,fund_y[r-1].
          */
	 for(i=0;i<r;i++){
	     tmp1=gauss_alpha*fund_x[i]
	     tmp2=gauss_beta*fund_y[i]
	     tmp3=gauss_gamma*fund_x[i]
	     tmp4=gauss_delta*fund_y[i]
             fund_x[i]=tmp1+tmp2
             fund_y[i]=tmp3+tmp4
             mm=calculate_n(a,b,c,fund_x[i],fund_y[i],n)
             if(print_flag){
                print "class n=",mm," (mod ",q00,"): (x,y)=("
                print fund_x[i],",",fund_y[i]
                print ")\n"
             }
         }
         return(r)
}

/* gauss(a,b,c,n) takes a triple gcd(a,b,c)=1 and |n|>1 and produces
 * (x,y)=(gauss_alpha,gauss_gamma) and an m=gauss_m such that ax^2+bxy+cy^2=m, 
 * gcd(m,n)=1. See L.-K. Hua 'Introduction to number theory', 311-312.
 */
define gauss(a,b,c,n){
       auto absn,e,g,i,s,t,tmp1,tmp2,xx[],yy[]

       absn=abs(n)
       e=omega(absn)
       for(i=0;i<e;i++){
	   tmp1=mod(a,qglobal[i])
	   tmp2=mod(c,qglobal[i])
           s=sign(tmp1)
           t=sign(tmp2)
           if(s){
	      xx[i]=1
	      yy[i]=0
           }
           if(s==0 && t){
	      xx[i]=0
	      yy[i]=1
           }
           if(s==0 && t==0){
	      xx[i]=1
	      yy[i]=1
           }
       }
       tmp1=chineseae(xx[],qglobal[],e)
       tmp2=chineseae(yy[],qglobal[],e)
       g=gcd(tmp1,tmp2)
       gauss_alpha=int(tmp1,g)
       gauss_gamma=int(tmp2,g)
       tmp1=a*gauss_alpha*gauss_alpha
       tmp2=b*gauss_alpha*gauss_gamma
       tmp3=c*gauss_gamma*gauss_gamma
       tmp4=tmp1+tmp2
       gauss_m=tmp4+tmp3
       return
}

/* binary(a,b,c,n,print_flag) finds the primitive solution classes for 
 * ax^2+bxy+cy^2=n, d=b^2-4ac>0 and nonsquare, n nonzero, gcd(a,b,c)=1.
 * print_flag=1 gives verbose output.
 */
define binary(a,b,c,n,print_flag){
auto s,i,g,absn,p,u,v,delta,d,result
       g=gcd3(a,b,c)
       absn=abs(n)
       if(n%absn){
          print "gcd(a,b,c) does not divide n\n"
          return(0)
       }else{
          a=a/g
          b=b/g
          c=c/g
          n=n/g
       }
       d=b^2-4*a*c
       if(print_flag){
          print "d = ",d,"\n"
       }
       if(b%2==0){
          p=b/2
          delta=p^2-a*c
          if(print_flag){
             print "delta = ",delta,"\n"
          }
       }
       s=binary1(a,b,c,n,print_flag)
       if(print_flag){
       print "The equation "
       result=printform(a,b,c,n)
       if(s==0){
             print "has no primitive solution\n"
       }else{
          if(s==1){
             print "has one primitive solution class: "
          }else{
             print "has ", s," primitive solution classes:\n";
          }
          for(i=0;i<s;i++){
              print "(",fund_x[i],",",fund_y[i],")\n"
              u=fund_x[i]
              v=fund_y[i]
          }
       }
       }
       return(s)
}
/* patz0 now simply returns the fundamental solutions and the total number of solution classes. */
define patz0(d,n,flag){
	auto a,b,k,q0,s,t,period_length,temp,g1_surd,g2_surd,f1,f2,solution_number
        solution_number=0

	q0=abs(n)
        s=sqroot(d,q0,0)
        /* s=0 means no solutions of x^2=d (mod |n|) */
        /* If s>0, we get all solutions x, 0<=x<=|n|/2 in the form
         * solution[0],...,solution[numbr-1]
         */
	if(s==0){
           return(solution_number)
	}
	/* Now to test the cfrac of each of w[j]=(-solution[j]+\sqrt{D})/q0
           to see if v_surd[j]=(-1)^j*sign(n) holds for some j in
           1<=j<=t+l, where l is the period of 
           w[j]=[a_surd[0],...,a_surd[t];a_surd[t+1],...,a_surd[t+l]]
           and t=pre_period_length.
         */
	for(k=0;k<numbr;k++){
            period_length=surd(d,1,-solution[k],q0,0)
            temp=pre_period_length+period_length
            if(period_length%2){
               temp=temp+period_length
            }
            for(a=1;a<=temp;a++){
                temp1=(-1)^a*sign(n)
                if(v_surd[a]==temp1){
                   f1=q_surd[a-1]
                   g1_surd=q0*p_surd[a-1]+solution[k]*q_surd[a-1]
                   break
                }
            }
            if(a==temp+1){
               continue /* process solution[k+1] */
	    }
            if(solution[k]==0){
              if(flag){
               print "Fundamental solution: ambiguous class P=0,"
               print "(x,y)=(",g1_surd,",",q_surd[a-1],")\n"
               patzx=g1_surd
               patzy=q_surd[a-1]
              }
               solution_number=solution_number+1
            }
            if(solution[k]!=0){
               period_length=surd(d,-1,-solution[k],q0,0)
               temp=pre_period_length+period_length
               if(period_length%2){
                  temp=temp+period_length
               }
	       for(a=1;a<=temp;a++){
                   temp1=(-1)^(a+1)*sign(n)
		   if(v_surd[a]==temp1){
                      f2=q_surd[a-1]
                      g2_surd=q0*p_surd[a-1]+solution[k]*q_surd[a-1]
                      break
                   }
               }
               if(flag){
	          print "Fundamental solution: "
               }
               if(2*solution[k]==q0){
                  if(flag){
                    print "ambiguous class, "
                    print "(x,y)=(",g1_surd,",",q_surd[a-1],")\n"
                  }
                  solution_number=solution_number+1
               }else{
                  solution_number=solution_number+2
               }
               if(flag){
                  print "P=",solution[k],", (x,y)=("
                  if(f1<=f2){
                     print g1_surd,",",f1
                     patzx=g1_surd
                     patzy=f1
                  }else{
                     print g2_surd,",",f2
                     patzx=g2_surd
                     patzy=f2
                  }
                 # print ")\n"
                  print "), "
               }
            }
	}/* end of for k loop */
        return(solution_number)
}

define selin_test(m,n,flag){
auto i,s,d,sn
for(i=m;i<=n;i++){
   s=selin(i)
   d=i^2-2
   if(flag){
     print "t=",i,":"
   }
   sn=patz0(d,s,flag)
   if(sn>4 && flag==0){
       print "t=",i,",d=",d, ",N=",s,",sn=",sn,":\n"
   }
}
}

define selin(t){
return(2*t^3-3*t^2-6*t+9)
}

define dujella_test(m,n,flag){
auto k,d,sn,temp
for(k=m;k<=n;k++){
   temp=k^2
   d=temp+1
   if(flag){
      print "k=",k,":"
   }
   sn=patz0(d,temp,flag)
   if(sn>2 && flag==0){
       print "k=",k,",d=",d, ",N=",temp,",sn=",sn,":\n"
   }
}
}

/*
define dujella_test1(m,n){
auto t,d,sn,temp
for(t=m;t<=n;t++){
   temp=f(2*t,2*t)
   d=temp^2+1
   print "t=",t,":"
   sn=patz0(d,temp^2)
   print "sn=",sn,":\n"
}
}*/

define conjecture(d,n,s){
auto p,i,q[]
   i=0
   p=s
   q[0]=p
   n[0]=n
   print "q[0]= ",s,","
   print "n[0]=",n,","
   while(n*n>d){
     i=i+1
     n=(p^2-d)/n
     n[i]=n
     print "n[",i,"]=",n,","
     p=p%n
     q[i]=p
     print "q[",i,"]=",p,","
   }
   if(n==1){
     print "reached n[",i,"]=1\n"
   }else{
     print "reached n[",i,"]=",n,"\n"
   }
}

define dujella_test2(k){
auto t,d,n,i,r
  d=k^2+1
  n=d-1
  print "k=",k,",d=k^2=",d,":"
  print "Solving congruence P^2=",1," mod(k^2)\n"
  print "1<=P<=n/2=",n/2,":\n"
  t=sqroot(1,n,0)
  t=t/2
  print "We get ",t," solutions: "
  for(i=0;i<t;i++){
    print solution[i],","
  }
  print"\n"
  for(i=0;i<t;i++){
     r=conjecture(d,n,solution[i])
  }
}

define dujella_test3(m,n){
auto k,t,cal,sn,f,s,kk,count1
   for(k=m;k<=n;k++){
       count1=0
       kk=k^2+1
       t=omega(k)
       cal=divisors(qglobal[],kglobal[],t)
       for(s=0;s<cal;s++){
         f=divisor[s]
         if(f*f<k){
           sn=patz0(kk,(k/f)^2,0) 
           if((sn>=2 && f>1)||(sn>=4 && f==1)){
              sn=patz0(kk,(k/f)^2,1) 
              print "sn:",sn,","
              print "k:",k,","
              print "f:",f,"\n"
           }
         }
       }
   }
   return
}

define dujella_test4(k){
auto d,div2k,t,s,kk,divisor2k[],divkk,ss,twok,n,a,b,g,flag,p,q

   kk=k^2+1
   twok=2*k
   t=omega(twok)
   div2k=divisors(qglobal[],kglobal[],t)
   for(s=0;s<div2k;s++){
     divisor2k[s]=divisor[s]
   }
   flag=0
   for(s=0;s<div2k;s++){
      d=divisor2k[s]
      if(d==k || d==twok){
         continue
      }
      if(k%2 && d%2){
        continue
      }
      /*print "processing divisor d=",d,":"*/
      t=omega(kk)
      divkk=divisors(qglobal[],kglobal[],t)
      n=twok/d
      /*print "2k/d=",n,"\n"*/
      for(ss=0;ss<divkk;ss++){
         b=divisor[ss]
         a=kk/b
         if(a<=2 || b<=2 || a>=b){
           continue
         }
         if(gcd(a,b)==1){
/*          g=binary1(a,0,-b,n,0)*/
            /*print "processing a=",a,",d=",d,"\n"*/
            g=rootdovera(a,k,d)
            if(g){
               print "k=",k,","
               print "d=",d,","
               print "a=",a,","
               print "b=",b,","
               print "\n"
               print "=================\n"
               flag=1
            }
         }
      }
   }
   return(flag)
}

/* for use when d=k^2+1 > 5 and a divides d. */
define rootda(d,a){
	auto h,p,q,e,f,g,x,y,k,l,m,n,u,v,r,s
        
	x=sqrt(d)
	p=0;q=a
	g=x*x
        l=0;k=1;m=1;n=0
	for(h=0;1;h++){
		y=(x+p)/q
                print y,","
                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"
                }
                l=k;m=n
                k=u;n=v
		f=p
		p=y*q-p
		e=q
		q=(d-(p*p))/q
		if(q==e){/* Q_h=Q_{h+1}, odd period 2h+1 */
                   print "Q_",h,"=Q_",h+1,"\n"
                   return(2*h+1)
		}
    }
}

# a divides dd=kk^2+1, ab=kk^2+1,2<a<b, gcd(a,b)=1.
# We get the cfrac of sqrt(dd)/a as far as the end of first period.
# dd divides 2kk, d distinct from kk and 2kk.
# returns 0 if no solution of ap^2-bq^2=pm 2k/d found, otherwise returns 1 with solutions
# arising from first period.  
define rootdovera(a,kk,dd){
	auto i,d,h,p,q,e,twokoverd,x,y,k,l,m,n,u,v,r,s,flag,b,a[],pcap[],qcap[],pu[],qv[]
        flag=0
        twokoverd=2*kk/dd
        d=kk^2+1
        b=d/a
	x=kk
	p=0;q=a
        l=0;k=1;m=1;n=0
	for(h=0;1;h++){
		y=(x+p)/q
                a[h]=y
                u=k*y+l;v=n*y+m
                       /* u/v is the i-th convergent to sqrt(d) */
                pcap[h]=p
                qcap[h]=q
                pu[h]=u
                qv[h]=v
                if(q==a && h>0){
                 if(flag){
                  for(i=0;i<=h;i++){
                    print "a[",i,"]=",a[i],",P[",i,"]=",pcap[i],",Q[",i,"]=",qcap[i],",A[",i,"]/B[",i,"]=",pu[i],"/",qv[i],"\n"
                  }
                  }
                  return(flag)
                }
                l=k;m=n
                k=u;n=v
		p=y*q-p
		e=q
		q=(d-(p*p))/q
                if(dd==1){
                  t1=abs(e-q+2*p)
                  if(t1==2*kk){
                    temp=u+l
                    temq=v+m
                     print "|Q[",h,"]-Q[",h+1,"]+2P[",h+1,"]|=",2*k,"=2k, A[",h,"]+A[",h-1,"]/B[",h,"]+B[",h-1,"]=",temp,"/",temq,"=p/q\n" 
                     print "solves ap^2-bq^2=pm2k/d=",a*(temp)^2-b*(temq)^2,"\n"
                     print "x=",dd*(a*(temp)^2+b*(temq)^2)/2,", y=",dd*(temp)*(temq),"\n"
                     flag=1
                  }
                  t2=abs(e-q-2*p)
                  if(t2==2*kk){
                     temp=u-l
                     temq=v-m
                     print "|Q[",h,"]-Q[",h+1,"]-2P[",h+1,"]|=",2*kk,"=2k, A[",h,"]-A[",h-1,"]/B[",h,"]-B[",h-1,"]=",temp,"/",temq,"=p/q\n" 
                     print "solves ap^2-bq^2=pm2k/d=",a*(u-l)^2-b*(v-m)^2,"\n"
                     print "x=",dd*(a*(u-l)^2+b*(v-m)^2)/2,", y=",dd*(u-l)*(v-m),"\n"
                     flag=1
                  }
                }
                if(q==twokoverd){
                     print "Q[",h+1,"]=",q,"=2k/d and A[",h,"]/B[",h,"]=",u,"/",v,"=p/q\n"
                     print "solves ap^2-bq^2=pm2k/d=",a*u^2-b*v^2,"\n"
                     print "x=",dd*(a*u^2+b*v^2)/2,", y=",dd*u*v,"\n"
                     flag=1
                }
       }
}

define dujella_cfrac_test(m,n){
auto k,t
   for(k=m;k<=n;k++){
      t=dujella_test4(k)
      if(t==0){
        continue
      }
   }
}

/* This lists the primitive fundamental solutions of x^2 - dy^2=n. */
define patztest(d,n){
auto h,g
   g=patz(d,n,0)
   print "Fundamental solutions:\n"
   for(h=0;h<g;h++){
     print "("
     print globalfundx[h],",",globalfundy[h]
     print ")\n"
   }
}

/* This creates the fundamental solutions of x^2 - dy^2=n as a global array. */
define patzgen(d,n){
auto absn,t,g,countpatzgen,temp,h,f
   absn=abs(n)
   if(absn!=1){
      t=omega(absn)
      g=divisors(qglobal[],kglobal[],t)
   }else{
      g=1
      divisor[0]=1
   }
   g=divisors(qglobal[],kglobal[],t)
   countpatzgen=0
   for(f=0;f<g;f++){
      temp=divisor[f]^2
      if(n%temp==0){
         t=patz(d,n/temp,0)
         for(h=0;h<t;h++){
            globalfundxff[countpatzgen+h]=divisor[f]*globalfundx[h]
            globalfundyff[countpatzgen+h]=divisor[f]*globalfundy[h]
         }
         countpatzgen=countpatzgen+t
      }
   }
   return(countpatzgen)
}

/* This lists the fundamental solutions of x^2 - dy^2=n. */
define patzgentest(d,n){
auto h,g
   g=patzgen(d,n)
   print g," Fundamental solutions (x,y):\n"
   for(h=0;h<g;h++){
     print "("
     print globalfundxff[h],",",globalfundyff[h]
     print ")\n"
   }
   return(g)
}


/* This finds all solutions (primitive and imprimitive) of ax^2+bxy+cy^2=n as a global array. */
/* b^2-4ac>0, a!=0,n!=0, gcd(a,b,c)=1. */
define binarygen(a,b,c,n,print_flag){
auto absn,t,g,countbinarygen,temp,h,f,result
   g=gcd3(a,b,c)
   if(g>1){
      print "gcd(a,b,c)>1\n"
      return(-1)
   }
   if(b^2-4*a*c<=0){
     print "b^2-4ac<=0\n"
     return(-1)
   }
   absn=abs(n)
   if(absn!=1){
      t=omega(absn)
      g=divisors(qglobal[],kglobal[],t)
   }else{
      g=1
      divisor[0]=1
   }
   countbinarygen=0
   for(f=0;f<g;f++){
      temp=divisor[f]^2
      if(n%temp==0){
/*       print "divisor[",f,"]=",divisor[f],"\n"*/
         temp1=n/temp
         if(print_flag){
            print "solving "
            result=printform(a,b,c,temp1)
            print " for primitive classes\n"
         }
         t=binary(a,b,c,temp1,print_flag)
         if(t>0){
            for(h=0;h<t;h++){
               globalbinarygenx[countbinarygen+h]=divisor[f]*fund_x[h]
               globalbinarygeny[countbinarygen+h]=divisor[f]*fund_y[h]
            }
            countbinarygen=countbinarygen+t
         }
      }
   }
   return(countbinarygen)
}

/* This lists the solution classes of ax^2+bxy+cy^2=n, b^2-4ac> 0, a!=0, n!=0. */
define binarygenlist(a,b,c,n,print_flag){
auto h,g,d,bmod2,bigdelta,phi,psi,null
# global globalx,globaly,globalr,globals

   g=binarygen(a,b,c,n,print_flag)
   if(g==-1){
      return(-1)
   }
   if(g==0){
      if(print_flag){
         print "there are no solutions for ",a,"x^2+(",b,")xy+(",c,")y^2=",n,":\n"
      }
      return(0)
   }
   d=b^2-4*a*c
   bmod2=b%2
   if(bmod2==0){
      bigdelta=d/4
      print "(b^2-4ac)/4=",bigdelta,"\n"
   }else{
      print "b^2-4ac=",d,"\n"
   }
   if(bmod2){
      null=fund4(d,0)
      phi=globalx;psi=globaly
      print "least positive solution (",phi,",",psi,") of  Pell's equation\n"
      print "phi^2 - ",d,"psi^2=4\n";
   }else{
      null=fund1(bigdelta)
      phi=globalr;psi=globals
      print "least positive solution (",phi,",",psi,") of  Pell's equation\n"
      print "phi^2 - ",bigdelta,"psi^2=1\n";
   }
   if(print_flag){
      if(g==1){
         print "one solution class for ",a,"x^2+(",b,")xy+(",c,")y^2=",n,":\n"
      }else{
         print g," solution classes for ",a,"x^2+(",b,")xy+(",c,")y^2=",n,"\n"
      }
   }
   for(h=0;h<g;h++){
     print "class ",h+1,": (",globalbinarygenx[h],",",globalbinarygeny[h],")\n"
   }
   return(g)
}


/* Here a is nonzero, d=b^2-4ac>0 and is nonsquare. 
 * We solve ax^2+bxy+cy^2=1. 
 */
define generalizedpell1(a,b,c,print_flag){
auto d,i,k,q,period_length,temp,g1_surd,g2_surd,f1,f2,u,v,x,y,ss,delta,p,returnflag,t,result

        if(print_flag){
           print "Solving "
           t=printform(a,b,c,1)
        }
        returnflag=0
        if(a==1){
           generalizedpellfundx=1
           generalizedpellfundy=0
	   print "solution generalized pell equation "
           result=printform(a,b,c,1)
           if(b%2==0){
              p=b/2
              delta=p^2-a*c
              print "delta = ",delta,"\n"
              print "x = Phi - (",p,")Psi\n"
              print "y = Psi,\n"
              print "where Phi^2 - ",delta,"Psi^2 = 1\n"
            }else{
              print "x = (Phi - (",b,")Psi)/2\n"
              print "y = Psi,\n"
              print "where Phi^2 - ",d,"Psi^2 = 4\n"
            }
            return(1)
        }
	d=b^2-4*a*c
        q=2*a
	/* to test the cfrac of w=(-b+\sqrt(d))/q
           to see if v_surd[j]=2(-1)^j holds for some j in
           1<=j<=t+l, where l is the period of 
           w=[a_surd[0],...,a_surd[t];a_surd[t+1],...,a_surd[t+l]]
           and t=pre_period_length.
	   also to test the cfrac of w*=(-b-\sqrt(d))/q
           to see if v_surd[j]=2(-1)^(j+1) holds for some j in
           1<=j<=t+l, where l is the period of 
           w*=[a_surd[0],...,a_surd[t];a_surd[t+1],...,a_surd[t+l]]
           and t=pre_period_length.
         */
        if(print_flag){
           print "processing omega=(",-b,"+sqrt(",d,")/",q,"\n"
        }
        period_length=surd(d,1,-b,q,print_flag)
        temp=pre_period_length+period_length
        if(period_length%2){
           temp=temp+period_length
        }
        for(i=1;i<=temp;i++){
            temp1=(-1)^i*2
            if(v_surd[i]==temp1){
               returnflag=1
               solution_flag=1
               f1=q_surd[i-1]
               g1_surd=p_surd[i-1]
               if(print_flag){
                  print "solution (x,y)=(",g1_surd,",",f1,")\n"
                  print "of generalized Pell equation "
                  result=printform(a,b,c,1)
               }
               break
            }
        }
        if(print_flag){
           print "processing omega*=(",-b,"-sqrt(",d,")/",q,"\n"
        }
        period_length=surd(d,-1,-b,q,print_flag)
        temp=pre_period_length+period_length
        if(period_length%2){
	   temp=temp+period_length
	}
	for(i=1;i<=temp;i++){
            temp1=(-1)^(i+1)*2
	    if(v_surd[i]==temp1){
	       f2=q_surd[i-1]
	       g2_surd=p_surd[i-1]
               if(print_flag){
	          print "solution (x,y)=(",g2_surd,",",f2,")\n"
                  print "of generalized Pell equation "
                  result=printform(a,b,c,1)
               }
	       break
            }
	}
        if(d!=5 || (d==5 && a>0)){
           if(returnflag==0){
              return(0)
           }
        }
        if(d!=5 || (d==5 && a>0)){
           if(f1<=f2){
              u=g1_surd
	      v=f1
	   }else{
              u=g2_surd
              v=f2
	   }
	   generalizedpellfundx=u
	   generalizedpellfundy=v
           if(print_flag){
	      print "solution with least positive y (", u,",",v,") "
              print "of generalized Pell equation "
              result=printform(a,b,c,1)
           }
        }
        if(d==5 && a<0){
           ss=pre_period_length
	   if(ss>1){
	      tmp1=p_surd[ss-1]-p_surd[ss-2]
	      tmp2=q_surd[ss-1]-q_surd[ss-2]
	   }else{
	      tmp1=p_surd[0]-1
	      tmp2=q_surd[0]
	   }
           returnflag=1
	   generalizedpellfundx=tmp1
	   generalizedpellfundy=tmp2
           u=tmp1
           v=tmp2
           if(print_flag){
	      print "exceptional solution (", tmp1,",",tmp2,") with least positive y\n"
              print "of generalized Pell equation "
              result=printform(a,b,c,1)
           }
	}
        if(print_flag){
           print "general solution of generalized pell equation:\n"
           if(b%2==0){
              p=b/2
              delta=p^2-a*c
              print "delta = ",delta,"\n"
              print "x=",u,"Phi - (",c*v+p*u,")Psi\n"
              print "y=",v,"Phi + (",p*v+a*u,")Psi\n"
              print "Phi^2 - ",delta,"Psi^2 = 1\n"
           }else{
              print "d = ",d,"\n"
              print "x=",u,"Phi/2 - (",2*c*v+b*u,")Psi/2\n"
              print "y=",v,"Phi/2 + (",b*v+2*a*u,")Psi/2\n"
              print "Phi^2 - ",d,"Psi^2 = 4\n"
           }
        }
         return(returnflag)
}
/* Here a is nonzero, d=b^2-4ac>0 and is nonsquare. 
 * We solve ax^2+bxy+cy^2=1. 
 */
define generalizedpell(a,b,c,print_flag){
auto d,i,k,q,period_length,temp,g1_surd,g2_surd,f1,f2,u,v,x,y,ss,delta,p,returnflag,t,result

        if(print_flag){
           print "Solving "
           t=printform(a,b,c,1)
        }
        returnflag=0
        if(a==1){
           generalizedpellfundx=1
           generalizedpellfundy=0
	   print "solution generalized pell equation "
           result=printform(a,b,c,1)
           print "with smallest nonnegative y: (1,0)\n"
           return(1)
        }
	d=b^2-4*a*c
        q=2*a
	/* to test the cfrac of w=(-b+\sqrt(d))/q
           to see if v_surd[j]=2(-1)^j holds for some j in
           1<=j<=t+l, where l is the period of 
           w=[a_surd[0],...,a_surd[t];a_surd[t+1],...,a_surd[t+l]]
           and t=pre_period_length.
	   also to test the cfrac of w*=(-b-\sqrt(d))/q
           to see if v_surd[j]=2(-1)^(j+1) holds for some j in
           1<=j<=t+l, where l is the period of 
           w*=[a_surd[0],...,a_surd[t];a_surd[t+1],...,a_surd[t+l]]
           and t=pre_period_length.
         */
        if(print_flag){
           print "processing omega=(",-b,"+sqrt(",d,")/",q,"\n"
        }
        period_length=surd(d,1,-b,q,print_flag)
        temp=pre_period_length+period_length
        if(period_length%2){
           temp=temp+period_length
        }
        for(i=1;i<=temp;i++){
            temp1=(-1)^i*2
            if(v_surd[i]==temp1){
               returnflag=1
               solution_flag=1
               f1=q_surd[i-1]
               g1_surd=p_surd[i-1]
               if(print_flag){
                  print "solution (x,y)=(",g1_surd,",",f1,")\n"
                  print "of generalized Pell equation "
                  result=printform(a,b,c,1)
               }
               break
            }
        }
        if(print_flag){
           print "processing omega*=(",-b,"-sqrt(",d,")/",q,"\n"
        }
        period_length=surd(d,-1,-b,q,print_flag)
        temp=pre_period_length+period_length
        if(period_length%2){
	   temp=temp+period_length
	}
	for(i=1;i<=temp;i++){
            temp1=(-1)^(i+1)*2
	    if(v_surd[i]==temp1){
	       f2=q_surd[i-1]
	       g2_surd=p_surd[i-1]
               if(print_flag){
	          print "solution (x,y)=(",g2_surd,",",f2,")\n"
                  print "of generalized Pell equation "
                  result=printform(a,b,c,1)
               }
	       break
            }
	}
        if(d!=5 || (d==5 && a>0)){
           if(returnflag==0){
              return(0)
           }
        }
        if(d!=5 || (d==5 && a>0)){
           if(f1<=f2){
              u=g1_surd
	      v=f1
	   }else{
              u=g2_surd
              v=f2
	   }
	   generalizedpellfundx=u
	   generalizedpellfundy=v
           if(print_flag){
	      print "solution with least positive y (", u,",",v,") "
              print "of generalized Pell equation "
              result=printform(a,b,c,1)
           }
        }
        if(d==5 && a<0){
           ss=pre_period_length
	   if(ss>1){
	      tmp1=p_surd[ss-1]-p_surd[ss-2]
	      tmp2=q_surd[ss-1]-q_surd[ss-2]
	   }else{
	      tmp1=p_surd[0]-1
	      tmp2=q_surd[0]
	   }
           returnflag=1
	   generalizedpellfundx=tmp1
	   generalizedpellfundy=tmp2
           u=tmp1
           v=tmp2
           if(print_flag){
	      print "exceptional solution (", tmp1,",",tmp2,") with least positive y\n"
              print "of generalized Pell equation "
              result=printform(a,b,c,1)
           }
	}
        if(print_flag){
           print "general solution of generalized pell equation:\n"
           if(b%2==0){
              p=b/2
              delta=p^2-a*c
              print "delta = ",delta,"\n"
              print "x=",u,"Phi - (",c*v+p*u,")Psi\n"
              print "y=",v,"Phi + (",p*v+a*u,")Psi\n"
              print "Phi^2 - ",delta,"Psi^2 = 1\n"
           }else{
              print "d = ",d,"\n"
              print "x=",u,"Phi/2 - (",2*c*v+b*u,")Psi/2\n"
              print "y=",v,"Phi/2 + (",b*v+2*a*u,")Psi/2\n"
              print "Phi^2 - ",d,"Psi^2 = 4\n"
           }
        }
         return(returnflag)
}

define printform(a,b,c,n){
print a,"x^2+(",b,")xy+(",c,")y^2=",n,"\n" 

} 

#define ssw2(a,b,c,d,e,f){
#auto gcount,h,dd,u,v,k,delta0,delta,gamma,zeta,epsilon,alpha,beta,psi,phi,p,aa,m,temp,temp1,x1,y1,null,twodelta0
#   if(a==0){
#      print "a=0\n"
#      return(-1)
#   }
#   gcount=0
#   dd=b^2-4*a*c
#   if(dd<0){
#      print "D<0\n"
#      return(-1)
#   }
#   if(dd==0){
#      print "D=0\n"
#      return(-1)
#   }
#   print "D=",dd,"\n"
#   if(b%2){
#      u=b*e-2*c*d
#      v=b*d-2*a*e
#   }else{
#      u=(b/2)*e-c*d
#      v=(b/2)*d-a*e
#   }
#   alpha=-u
#   beta=-v
#   if(b%2==0){
#      p=b/2
#      delta0=p^2-a*c
#      twodelta0=2*delta0
#      print "delta0 = ",delta0,"\n"
#      null=fund1(delta0)
#      phi=globalr
#      psi=globals
#      print "phi^2 - ",delta0,"psi^2 = 1\n"
#   }
#  if(b%2){
#     null=fund4(dd,0)
#     phi=globalx
#     psi=globaly
#     print "phi^2 - ",dd,"psi^2 = 4\n"
#  }
#   print "phi=",phi,", psi=",psi,"\n"
#   print "alpha=",alpha,", beta=",beta,"\n"
#  if(b%2==0){
#     k=-(a*u^2+b*u*v+c*v^2-d*twodelta0*u-e*twodelta0*v+f*twodelta0^2)
#     kk=-delta0*(a*e^2-b*e*d+c*d^2+4*f*delta0)
#   }else{
#     k=-dd*(a*e^2-b*e*d+c*d^2+f*dd)
#   }
#   print "k=",k,"\n"
#   print "kk=",kk,"\n"
#   if(k==0){
#      print "k=0\n"
#      r=-u
#      s=-v
#      if(r%dd==0 && s%dd==0){
#         print "(x,y)=(",r/dd,",",s/dd,")\n"
#         gcount=1
#         return(1)
#      }else{
#         print "no integer solution\n"
#         return(0)
#      }
#  }
#  if(b%2){
#     print dd,"x=X+(",alpha,")," 
#     print dd,"y=Y+(",beta,")\n" 
#     print "\n"
#  }else{
#     print twodelta0,"x=X+(",alpha,")," 
#     print twodelta0,"y=Y+(",beta,")\n" 
#     print "\n"
#  }
#
# 
#  print "First find the solution classes of ",a,"X^2+",b,"XY+",c,"Y^2=",k,"\n"
#  g=binarygenlist(a,b,c,k,0)
#  for(h=0;h<g;h++){
#      /*print "h=",h,"\n"*/
#      gamma=globalbinarygenx[h]
#      epsilon=globalbinarygeny[h]
#      if(b%2==0){
#         delta=-(c*epsilon+p*gamma)
#         zeta=p*epsilon+a*gamma
#         a1=2*delta*globals
#         a2=2*zeta*globals
#
#         b1=-(alpha+gamma)*globalr
#         b2=-(beta+epsilon)*globalr
#print "1:a1=",a1,",a2=",a2,"\n"
#print "1:b1=",b1,",b2=",b2,"\n"
#         r=legendresimulcong(a1,b1,a2,b2,delta0,delta0)
#print "1:r=",r,"\n"
#         if(r){
#              print twodelta0,"x=(",gamma,")F+(",delta,")G+(",alpha,")\n"
#              print twodelta0,"y=(",epsilon,")F+(",zeta,")G+(",beta,")\n"
#              aa=globalchineselegsol
#              m=globalchineselegmod
#              null=powerdd(phi,psi,delta0,2*aa)
#              x1=gamma*z1+delta*z2
#              y1=epsilon*z1+zeta*z2
#              print "(x,y)=(",(x1+alpha)/delta0,",",(y1+beta)/twodelta0,")\n"
#              print "F+Gsqrt(",delta0,")=(",phi,"+(",psi,")sqrt(",delta0,"))^"
#              null=printlegendre(aa,m,0)
#              gcount=gcount+1
#              print "\n"
#         }
#
#         b1=-(alpha+gamma*globalr+delta*globals)
#         b2=-(beta+epsilon*globalr+zeta*globals)
#print "2:a1=",a1,",a2=",a2,"\n"
#print "2:b1=",b1,",b2=",b2,"\n"
#         r=legendresimulcong(a1,b1,a2,b2,delta0,delta0)
#print "2:r=",r,"\n"
#         if(r){
#              print twodelta0,"x=(",gamma,")F+(",delta,")G+(",alpha,")\n"
#              print twodelta0,"y=(",epsilon,")F+(",zeta,")G+(",beta,")\n"
#              aa=globalchineselegsol
#              m=globalchineselegmod
#              null=powerdd(phi,psi,delta0,2*aa+1)
#              x1=gamma*z1+delta*z2
#              y1=epsilon*z1+zeta*z2
#              print "(x,y)=(",(x1+alpha)/twodelta0,",",(y1+beta)/twodelta0,")\n"
#              print "F+Gsqrt(",delta0,")=(",phi,"+(",psi,")sqrt(",delta0,"))^"
#              null=printlegendre(aa,m,1)
#              gcount=gcount+1
#              print "\n"
#         }
#         b1=(alpha-gamma)*globalr
#         b2=(beta-epsilon)*globalr
#print "3:a1=",a1,",a2=",a2,"\n"
#print "3:b1=",b1,",b2=",b2,"\n"
#         r=legendresimulcong(a1,b1,a2,b2,delta0,delta0)
#print "3:r=",r,"\n"
#         if(r){
#              print twodelta0,"x=(",gamma,")F+(",delta,")G+(",alpha,")\n"
#              print twodelta0,"y=(",epsilon,")F+(",zeta,")G+(",beta,")\n"
#              aa=globalchineselegsol
#              m=globalchineselegmod
#              null=powerdd(phi,psi,delta0,2*aa)
#print "3:aa=",aa,",m=",m,"\n"
#print "3:z1=",z1,",z2=",z2,"\n"
#print "3:gamma=",gamma,",epsilon=",epsilon,"\n"
#print "3:delta=",delta,",zeta=",zeta,"\n"
#              x1=gamma*(-z1)+delta*(-z2)
#              y1=epsilon*(-z1)+zeta*(-z2)
#print "3:x1=",x1,",y1=",y1,"\n"
#print "3:alpha=",alpha,",beta=",beta,"\n"
#              print "(x,y)=(",(x1+alpha)/twodelta0,",",(y1+beta)/twodelta0,")\n"
#              print "F+Gsqrt(",delta0,")=-(",phi,"+(",psi,")sqrt(",delta0,"))^"
#              null=printlegendre(aa,m,0)
#              gcount=gcount+1
#              print "\n"
#         }
#
#         b1=alpha-gamma*globalr-delta*globals
#         b2=beta-epsilon*globalr-zeta*globals
#print "4:a1=",a1,",a2=",a2,"\n"
#print "4:b1=",b1,",b2=",b2,"\n"
#         r=legendresimulcong(a1,b1,a2,b2,twodelta0,twodelta0)
#print "4:r=",r,"\n"
#         if(r){
#              print twodelta0,"x=(",gamma,")F+(",delta,")G+(",alpha,")\n"
#              print twodelta0,"y=(",epsilon,")F+(",zeta,")G+(",beta,")\n"
#              aa=globalchineselegsol
#              m=globalchineselegmod
#              null=powerdd(phi,psi,delta0,2*aa+1)
#              x1=gamma*(-z1)+delta*(-z2)
#              y1=epsilon*(-z1)+zeta*(-z2)
#              print "(x,y)=(",(x1+alpha)/twodelta0,",",(y1+beta)/twodelta0,")\n"
#              print "F+Gsqrt(",delta0,")=-(",phi,"+(",psi,")sqrt(",delta0,"))^"
#              null=printlegendre(aa,m,1)
#              gcount=gcount+1
#              print "\n"
#         }
#      }else{
#         delta=-(2*c*epsilon+b*gamma)
#         zeta=b*epsilon+2*a*gamma
#         a1=2*delta*globaly
#         a2=2*zeta*globaly
#         b1=-(alpha+gamma)*globalx
#         b2=-(beta+epsilon)*globalx
#         r=legendresimulcong(a1,b1,a2,b2,dd,dd)
#         if(r){
#              print dd,"X=(",gamma,")F+(",delta,")G+(",alpha,")\n"
#              print dd,"Y=(",epsilon,")F+(",zeta,")G+(",beta,")\n"
#              aa=globalchineselegsol
#              m=globalchineselegmod
#              temp=2*aa
#              temp1=2^temp
#              null=powerdd(phi,psi,dd,temp)
#              x1=(gamma*z1+delta*z2)/temp1
#              y1=(epsilon*z1+zeta*z2)/temp1
#              print "(x,y)=(",(x1+alpha)/dd,",",(y1+beta)/dd,")\n"
#              print "F+Gsqrt(",dd,")=(",phi,"/2+(",psi,"/2)sqrt(",dd,"))^"
#              null=printlegendre(aa,m,0)
#              gcount=gcount+1
#              print "\n"
#         }
#         b1=-(2*alpha+gamma*globalx+delta*globaly)
#         b2=-(2*beta+epsilon*globalx+zeta*globaly)
#         r=legendresimulcong(a1,b1,a2,b2,dd,dd)
#         if(r){
#              print dd,"x=(",gamma,")F+(",delta,")G+(",alpha,")\n"
#              print dd,"y=(",epsilon,")F+(",zeta,")G+(",beta,")\n"
#              aa=globalchineselegsol
#              m=globalchineselegmod
#              temp=2*aa+1
#              temp1=2^temp
#              null=powerdd(phi,psi,dd,temp)
#              x1=(gamma*z1+delta*z2)/temp1
#              y1=(epsilon*z1+zeta*z2)/temp1
#              print "(x,y)=(",(x1+alpha)/dd,",",(y1+beta)/dd,")\n"
#              print "F+Gsqrt(",dd,")=(",phi,"/2+(",psi,"/2)sqrt(",dd,"))^"
#              null=printlegendre(aa,m,1)
#              gcount=gcount+1
#              print "\n"
#         }
#         b1=(alpha-gamma)*globalx
#         b2=(beta-epsilon)*globalx
#         r=legendresimulcong(a1,b1,a2,b2,dd,dd)
#         if(r){
#              print dd,"x=(",gamma,")F+(",delta,")G+(",alpha,")\n"
#              print dd,"y=(",epsilon,")F+(",zeta,")G+(",beta,")\n"
#              aa=globalchineselegsol
#              m=globalchineselegmod
#              temp=2*aa
#              temp1=2^temp
#              null=powerdd(phi,psi,dd,temp)
#              x1=(gamma*(-z1)+delta*(-z2))/temp1
#              y1=(epsilon*(-z1)+zeta*(-z2))/temp1
#              print "(x,y)=(",(x1+alpha)/dd,",",(y1+beta)/dd,")\n"
#              print "F+Gsqrt(",dd,")=-(",phi,"/2+(",psi,"/2)sqrt(",dd,"))^"
#              null=printlegendre(aa,m,0)
#              gcount=gcount+1
#              print "\n"
#         }
#         b1=2*alpha-gamma*globalx-delta*globaly
#         b2=2*beta-epsilon*globalx-zeta*globaly
#         r=legendresimulcong(a1,b1,a2,b2,dd,dd)
#         if(r){
#              print dd,"x=(",gamma,")F+(",delta,")G+(",alpha,")\n"
#              print dd,"y=(",epsilon,")F+(",zeta,")G+(",beta,")\n"
#              aa=globalchineselegsol
#              m=globalchineselegmod
#              temp=2*aa+1
#              temp1=2^temp
#              null=powerdd(phi,psi,dd,temp)
#              x1=(gamma*(-z1)+delta*(-z2))/temp1
#              y1=(epsilon*(-z1)+zeta*(-z2))/temp1
#              print "(x,y)=(",(x1+alpha)/dd,",",(y1+beta)/dd,")\n"
#              print "F+Gsqrt(",dd,")=-(",phi,"/2+(",psi,"/2)sqrt(",dd,"))^"
#              null=printlegendre(aa,m,1)
#              gcount=gcount+1
#              print "\n"
#         }
#      }
#  }
#  print "There are ",gcount," solution classes\n"
#  return(gcount)
#}

/* this finds the fundamental solution globalr+sqrt(d)globals  of Pell's equation x^2-dy^2=1
   using the midperiod approach. */
define fund1(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)
		}
    }
}
/* This finds the fundamental solution (globalx,globaly) of Pell's equation x^2-dy^2=4,
 * where d>0 and is nonsquare. 
 * We use the fact that a positive solution with gcd(x,y)=1 must be a convergent of
 * the continued fraction expansion of sqrt(d). If there is no such convergent, we use th
 * familiar midpoint method of finding the least solution of Pell's equations X^2-dY^2=1
 * and then (globalx,globaly)=(2X,2Y).
 */
define fund4(d,printflag){
	auto h,p,q,e,f,g,x,y,k,l,m,n,u,v,r,s,oldl,oldm,oldu,oldv,t,flag
        flag=0
        if(d==2){
          globalx=6
          globaly=4
          flag=1
        }
        if(d==3){
          globalx=4
          globaly=2
          flag=1
        }
        if(d==5){
          globalx=3
          globaly=1
          flag=1
        }
        if(d==6){
          globalx=10
          globaly=4
          flag=1
        }
        if(d==7){
          globalx=16
          globaly=6
          flag=1
        }
        if(d==8){
          globalx=6
          globaly=2
          flag=1
        }
        if(d==10){
          globalx=38
          globaly=12
          flag=1
        }
        if(d==11){
          globalx=20
          globaly=6
          flag=1
        }
        if(d==12){
          globalx=4
          globaly=1
          flag=1
        }
        if(d==13){
          globalx=11
          globaly=3
          flag=1
        }
        if(d==14){
          globalx=30
          globaly=8
          flag=1
        }
        if(d==15){
          globalx=8
          globaly=2
          flag=1
        }
        x=sqrt(d)
        if(d==x^2+1 && d >17){
          globalx=4*x^2+2
          globaly=4*x
          flag=1
        }
        if(printflag && flag){
           print "globalx=",globalx,",globaly=",globaly,"\n"
        }
        if(flag){
           return(d)
        }
	x=sqrt(d)
	p=0;q=1
	g=x*x
        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
/*print "n=",n,",y=",y,",m=",m,",v=",v,"\n"*/
                if(printflag){
                      print "A[",h,"]/B[",h,"]=",u,"/",v,"\n"
                }
                oldl=l;oldm=m
                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
                t=(-1)^(h+1)
                if(q==4 && t==1){
                   globalx=k
                   globaly=n
                   if(printflag){
                      print "globalx=",globalx,",globaly=",globaly,"\n"
                   }
                   return(4)
                }
                if(q==4 && t!=1){
                   globalx=(k^2+d*n^2)/2
                   globaly=k*n
                   if(printflag){
                      print "globalx=",globalx,",globaly=",globaly,"\n"
                   }
                   return(-4)
                }
		if(p==f){/* P_h=P_{h+1}, even period 2h */
                   if(printflag){
                      print "P_",h,"=P_",h+1,"\n"
                   }
                   globalx=2*(m*(oldu+oldl)+(-1)^h)
                   globaly=2*m*(oldv+oldm)
                   if(printflag){
                      print "globalx=",globalx,",globaly=",globaly,"\n"
                   }
                   return(1)
		}
		if(q==e){/* Q_h=Q_{h+1}, odd period 2h+1 */
                   if(printflag){
                      print "Q_",h,"=Q_",h+1,"\n"
                   }
                   r=k*n+l*m
                   s=n^2+m^2
                   globalx=2*(r*r+d*s*s)
                   globaly=4*r*s
                   if(printflag){
                      print "globalx=",globalx,",globaly=",globaly,"\n"
                   }
                   return(-1)
		}
       }
}

define printlegendre(aa,m,parity){
if(parity==0){
   if(aa==0){
      print 2*m,"k\n"
   }else{
      print "(",2*m,"k+",2*aa,")\n"
   }
}else{
   if(aa==0){
      print "(",2*m,"k+1)\n"
   }else{
      print "(",2*m,"k+",2*aa+1,")\n"
   }
}
}

# This solves the diophantine equation ax^2+bxy+cy^2+dx+ey+f=0, where  (a,b,c)!=(0,0,0).
# It uses a method of John Robertson, communicated on March 25, 2015 to deal with the 
# general hyperbolic case.
# needs squareroot2 powertest3 genfacs gcd posformrep reduceneg congruence powerdd
define sswgeneral(a,b,c,d,e,f){
auto gg1,gg2,gg,gcount,h,dd,k,delta,gamma,zeta,epsilon,alpha,beta,psi,phi,x1,y1,null,twoa,absdd,twoadd,bigdelta,g,sol[]

  dd=b^2-4*a*c
  print "D=",dd,"\n"
  twoa=2*a
  if(dd>0){
     g=squaretest(dd)
     if(g<0 && a^2==1 && b==0 && d==0 && e==0 && f!=0){
        if(a==-1){
           a=1;c=-c;f=-f
        }
        g=patzgentest(-c,-f)
        t=fund1(-c)
        print "general solutions +-(x+ysqrt(",-c,"))(",globalr,"+",globals,"sqrt(",-c,")^n\n"
        return(g)
     }
     if(g<0 && d==0 && e==0 && f!=0){
        g=binarygen(a,b,c,-f,0)
        for(h=0;h<g;h++){
          gamma=globalbinarygenx[h]
          epsilon=globalbinarygeny[h]
          if(b%2){
             delta=-(2*c*epsilon+b*gamma)
             zeta=b*epsilon+2*a*gamma
             print "x=(",gamma,")F+(",delta,")G\n"
             print "y=(",epsilon,")F+(",zeta,")G\n"
          }else{
             delta=-(c*epsilon+(b/2)*gamma)
             zeta=(b/2)*epsilon+a*gamma
             print "x=(",gamma,")F+(",delta,")G\n"
             print "y=(",epsilon,")F+(",zeta,")G\n"
          }
          print "\n"
      }
      if(b%2){
         null=fund4(dd,0)
         phi=globalx
         psi=globaly
         print "phi^2 - ",dd,"psi^2 = 4\n"
         print "phi=",phi,", psi=",psi,"\n"
         print "F+Gsqrt(",dd,")=+-(",phi,"/2+(",psi,"/2)sqrt(",dd,"))^m\n"
      }else{
         bigdelta=dd/4
         null=fund1(bigdelta)
         phi=globalr
         psi=globals
         print "phi^2 - ",dd,"psi^2 = 1\n"
         print "phi=",phi,", psi=",psi,"\n"
         print "F+Gsqrt(",bigdelta,")=+-(",phi,"+",psi,"sqrt(",dd,"))^m\n"
      }
      return(g)
     }
  }
  if(dd<0 && d==0 && e==0){
     g=posrepgenlist(a,b,c,-f)
     return(g)
  }
  if(dd==0){
     t=ddzero(a,b,c,d,e,f,0)
     return
  }
  if(dd>0){
    g=squaretest(dd)
  }
  gg=gcd3(a,b,c)
  gcount=0
  alpha=-b*e+2*c*d
  beta=-b*d+2*a*e
  print "alpha=",alpha,", beta=",beta,"\n"
  print dd,"x=X+(",alpha,")," 
  print dd,"y=Y+(",beta,")\n" 
  print "\n"
  k=-dd*(a*e^2-b*e*d+c*d^2+f*dd)
  print "solving (",a,")X^2+(",b,")XY+(",c,")Y^2=",k,"\n"
   if(k%gg){
      print "no solution\n"
      return(0)
   }
   if(dd>0){
      g=squaretest(dd)
   }
   absdd=abs(dd)
   if(dd<0 || (dd>0 && g<0)){
      if(k==0){
         print "k=0\n"
         if(alpha%absdd==0 && beta%absdd==0){
            print "(x,y)=(",alpha/dd,",",beta/dd,")\n"
            return(1)
         }else{
            print "no integer solution\n"
            return(0)
         }
     }
  }
  
  if(dd<0 && k){# ellipse
      a=a/gg;b=b/gg;c=c/gg;k=k/gg
      gcount=ellipse(a,b,c,k,dd,alpha,beta)
      return(gcount)
  }
  if(dd>0 && g<0 && k!=0){# hyperbola
     g=bigu(a,b,c,d,e,f,0)
     return
  }
  gcount=ddsquare(a,b,c,k,g,twoa,dd,alpha,beta)
  return(gcount)
}

define hyperbola(a,b,c,k,alpha,beta,dd){# gcd(a,b,c) may be > 1, dd=b^2-4ac.
auto gg,h,t,g,aa,bb,cc,kk,bmod2,delta,zeta, deltover2,zetaover2,gamma,epsilon

  gg=gcd3(a,b,c)
  aa=a/gg;bb=b/gg;cc=c/gg;kk=k/gg
  g=binarygenlist(aa,bb,cc,kk,0) # this function needs gcd(aa,bb,cc)=1
  for(h=0;h<g;h++){
      /*print "h=",h,"\n"*/
      gamma=globalbinarygenx[h]
      epsilon=globalbinarygeny[h]
      t=functionu(a,b,c,k,alpha,beta,dd,gamma,epsilon,0)
      #print "(gamma,epsilon)=(",gamma,",",epsilon,")"
      if(t!=0){
        print " gives the solution family:\n"
        delta=-(2*cc*epsilon+bb*gamma)
        zeta=bb*epsilon+2*aa*gamma
        bmod2=bb%2
        if(bmod2==0){
          deltaover2=delta/2
          zetaover2=zeta/2
        }
        if(bmod2){
           print dd,"x=(",gamma,")F+(",delta,")G+(",alpha,")\n"
           print dd,"y=(",epsilon,")F+(",zeta,")G+(",beta,")\n"
        }else{
           print dd,"x=(",gamma,")F+(",deltaover2,")G+(",alpha,")\n"
           print dd,"y=(",epsilon,")F+(",zetaover2,")G+(",beta,")\n"
        }
         gcount=gcount+1
         print "\n"
     }
  }
  if(gcount>1 || gcount==0){
     print "There are ",gcount," solution families\n"
  }else{
     print "There is 1 solution family\n"
  }
  return(gcount)
}


#dd>0 && g>0
define ddsquare(a,b,c,k,g,twoa,dd,alpha,beta){
auto g1,g2,fourak,h,tcount,gg,t,g1timesg2,flag,xhplusalpha,yhplusalpha
#global globallinearnx,globallinearny

      fourak=4*a*k
print "k=",k,"\n"
      if(a!=0){
        g1=gcd(twoa,b+g)
        g2=gcd(twoa,b-g)
        g1timesg2=g1*g2
        if(fourak%g1timesg2){# g1g2 does not divide 4ak
           print "no solution\n"
           return(0)
        }else{
            print "i.e., solving (",twoa/g1,"X+(",(b+g)/g1,")Y)(",twoa/g2,"X+(",(b-g)/g2,")Y)=",fourak/g1timesg2,"\n"
        }
        if(k){# g1g2 divides 4ak and k is nonzero
           print "This gives finitely many solutions on factorising the rhs\n"
           t=lineareqnn(twoa/g1,(b+g)/g1,twoa/g2,(b-g)/g2,0,0,fourak/g1timesg2)
           tcount=0
           for(h=0;h<t;h++){
             xhplusalpha=globallinearnx[h]+alpha
             yhplusbeta=globallinearny[h]+beta
             if(xhplusalpha%dd==0 && yhplusbeta%dd==0){
                print "solution[",tcount,"]: (",xhplusalpha/dd,",",yhplusbeta/dd,")\n"
                tcount=tcount+1
             }
           }
           print tcount," solutions\n"
           return(tcount)
        }else{#k=0
          #print "This gives infinitely many solutions corresponding to\n"
          #print twoa,"X+(",b+g,")Y=0 and ",twoa,"X+(",b-g,")Y)=0\n"
           t=twoa*alpha+(b+g)*beta
           if(t%(g1*dd)){
              print "no solution\n"
              return(0)
           }else{
              temp=t/(g1*dd)
              print "infinitely many solutions:", twoa/g1,"x+(",(b+g)/g1,")y=",temp,"\n"
           }
           t=twoa*alpha+(b-g)*beta
           if(t%(g2*dd)){
              print "no solution\n"
              return(0)
           }else{
              temp=t/(g2*dd)
              print "infinitely many solutions: ",twoa/g2,"x+(",(b-g)/g2,")y=",temp,"\n"
              return(-2)
           }
        }#end of k=0
      }else{#a=0 and we are solving Y(bX+Y)=k
        print "Y(",b,"X+(",c,")Y)=",k,"\n"
        gg=gcd(b,c)
        if(k%gg){# gcd(b,c) does not divide k
           print "no solution\n"
           return(0)
        }else{
           if(gg>1){
               b=b/gg;c=c/gg;k=k/gg
               print "i.e., Y(",b,"X+(",c,")Y=",k,"\n"
           }
        }
        if(k){# k is nonzero and now gcd(b,c)=1
           t=lineareqnn(0,1,b,c,0,0,k)
           tcount=0
           for(h=0;h<t;h++){
               xhplusalpha=globallinearnx[h]+alpha
               yhplusbeta=globallinearny[h]+beta
               if(xhplusalpha%dd==0 && yhplusbeta%dd==0){
                  print "solution[",tcount,"]: (",xhplusalpha/dd,",",yhplusbeta/dd,")\n"
                  tcount=tcount+1
               }
           }
           print tcount," solutions\n"
           return(tcount)
        }else{#k=0 and we are solving "Y(",b,"X+(",c,")Y=0\n"
           print "i.e. line of solutions , ",dd,"y=",beta,"\n"
           if(beta%dd){# d does not divide beta
              flag=0
           }else{# D divides beta
              print "line of solutions , ",dd,"y=",beta/dd,"\n"
              flag=1
           }
           t=b*alpha+c*beta
           if(t%dd){# D does not divide b.alpha+c.beta
              return(-2)
           }else{ # D divides b.alpha+c.beta
              if(flag){# D diides beta and we have solutions y=beta/D
                  print "and the\n"
              }
              print "line of solutions ",b,"x+(",c,")y=",t/dd,"\n"
              return(-2)
           }
        }
      }#end of a=0
}
# This function is based on AA's book, page 84.
define aa1(a,b,c){
   auto d,g,alpha,beta,count,x
   if(a<=0){
      print "a<=0\n"
      return(-1)
   }
   if(b<=0){
      print "b<=0\n"
      return(-1)
   }
   d=a*b
   x=sqrt(d)
   if(x^2==d){
     print "ab is a perfect square\n"
     return(-1)
   }
   if(c==0){
      print "rc==0\n"
      return(-1)
   }
   count=0
   g=patzgen(d,a*c)
   print g," Fundamental solutions (x,y):\n"
   for(h=0;h<g;h++){
      alpha=globalfundxff[h]
      beta=globalfundyff[h]
      if(alpha%a==0){
        print "(",alpha,",",beta,")\n"
        print"\tx=", (alpha/a),"u+",b*beta,"v\n"
        print"\ty=", beta,"u+",alpha,"v\n"
        count=count+1
      }
   }
   print"\n"
   print "Here u^2-",d,"v^2=1\n"
   print "the number of solution families of \n"
   print a,"x^2-",b,"y^2=",c," is "
   return(count)
}

# This function is based on AA's book, page 94-95.  We use Stolt equivalence instead of
# Nagell equivalence.
define aa2(a,b,c,n){
auto d,fouran,x1,y1,r,an,alpha,beta,twoa,bound,foura,x,count
if(a==0){
   print "a=0\n"
   return(-1)
}
if(b==0){
   print "b=0\n"
   return(-1)
}
if(n==0){
   print "n=0\n"
   return(-1)
}
if(a<0){
   b=-b;c=-c;n=-n
}
d=b^2-4*a*c
if(d<=0){
   print "d<=0\n"
   return(-1)
}
x=sqrt(d)
if(x^2==d){
   print "d is a perfect square\n"
   return(-1)
}
fouran=4*a*n
foura=4*a
twoa=2*a
d=b^2-4*a*c
g=stolt0(d,fouran,0)
   for(h=0;h<g;h++){
      alpha=globalstoltfundx[h]
      beta=globalstoltfundy[h]
      if((alpha-b*beta)%twoa==0){
         print "(",alpha,",",beta,")\n"
         if(d%2){
            print"\tx=("
            x=printauplusbv((alpha-b*beta)/twoa,(d*beta-b*alpha)/twoa)
            print")/2\n"
            print"\ty=("
            x= printauplusbv(beta,alpha)
            print")/2\n"
         }else{
            print"\tx="
            x=printauuplusbv((alpha-b*beta)/twoa,(d*beta-b*alpha)/foura)
            print "\n"
            print"\ty="
            x=printauuplusbv(beta,alpha/2)
            print "\n"
      
         }
         count=count+1
      }
   }
   if(count){
      if(d%2){
         print "Here u^2-",d,"v^2=4\n"
      }else{
         print "Here U^2-",d/4,"v^2=1\n"
      }
   }
   print "the number of solution families of \n"
   print a,"x^2+",b,"xy+",c,"y^2=",n," is "
   return(count)
}

define printauplusbv(a,b){
   if(a>1){
     print a,"u"
   }
   if(a==1){
     print "u"
   }
   if(a==-1){
     print "-u"
   }
   if(a<-1){
     print a,"u"
   }
   if(b>1){
     print "+",b,"v"
   }
   if(b==1){
     print "+v"
   }
   if(b==-1){
     print "-v"
   }
   if(b<-1){
     print b,"v"
   }
}
define printauuplusbv(a,b){
   if(a>1){
     print a,"U"
   }
   if(a==1){
     print "U"
   }
   if(a==-1){
     print "-U"
   }
   if(a<-1){
     print a,"U"
   }
   if(b>1){
     print "+",b,"v"
   }
   if(b==1){
     print "+v"
   }
   if(b==-1){
     print "-v"
   }
   if(b<-1){
     print b,"v"
   }
}

define lineareqn(a,b,c,d,e,f){
auto delta1,delta2,absdelta,delta
 delta1=e*d-b*f
 delta2=a*f-e*c
 delta=a*d-b*c
 absdelta=abs(delta)
 if(delta1%absdelta==0 && delta2%absdelta==0){
    globallinearx=delta1/delta
    globallineary=delta2/delta
    return(1)
 }else{
    return(0)
 }
}

define lineareqnn(a,b,c,d,r,s,n){
auto g,t,tcount,f,twog
   tcount=0
   absn=abs(n)
   if(absn!=1){
      t=omega(absn)
      g=divisors(qglobal[],kglobal[],t)
   }else{
      g=1
      divisor[0]=1
   }
   for(f=0;f<g;f++){
      divisor[g+f]=-divisor[f]
   }
   twog=2*g
   for(f=0;f<twog;f++){
       g=lineareqn(a,b,c,d,divisor[f]-r,n/divisor[f]-s)
       if(g){
         globallinearnx[tcount]=globallinearx
         globallinearny[tcount]=globallineary
          tcount=tcount+1
       }
   }
   return(tcount)
}

# Here dd=b^2-4ac=0. We use completion of the square as in L.K. Hua, p.278.
define ddzero(a,b,c,d,e,f,flag){
auto g,aa,bb,cc,roota,rootc,t,u,v,absu,s,r,j,ksol,twoa,x,y,tcount

   if(b==0){# then a=0 or c=0, but not both.
      if(a==0){
         print "We are interchanging x and y to solve\n"
         print "(",c,")x^2+(",e,")x+(",d,")y+(",f,")=0\n"
         a=c
         temp=d
         d=e
         e=temp
      }
      if(flag){
         print "solving (",a,")x^2+(",d,")x+(",e,")y+(",f,")=0\n"
      }
   }
   twoa=2*a
   u=2*(b*d-2*a*e)
   v=d^2-4*a*f
   # solving (t+d)^2=uy+v
   if(u==0){  
      #print "we have (t+d)^2=v\n"
      if(v==0){
         if(b==0){
            if(d%twoa==0){
               print "we get a line of solutions x=",-d/twoa,", y arbitrary\n"
               return(-2)
            }else{
               print "no solutions\n"
               return(0)
            }
         }else{
            t=gcd(twoa,b)
            if(d%t==0){
               print "we get a line of solutions ",twoa/t,"x+(",b/t,")y=",-d/t,"\n"
               return(-2)
            }else{
               print "no solutions\n"
               return(0)
            }
         }
      }else{#v!=0
        if(v<0){
           print "no solutions\n"
           return(0)
        }else{#v=d^2-4af>0
           g=sqrt(v)
           if(g^2==v){
              if(b==0){
                 if((g-d)%twoa==0){
                    print "solution x=",(g-d)/twoa," y, arbitrary\n"
                 }
                 if((-g-d)%twoa==0){
                    print "solution x=",-(g+d)/twoa,", y arbitrary\n"
                 }
                 return(-2)
              }else{#b!=0
                    t=gcd(twoa,b)
                    tcount=0
                    if((-d+g)%t==0){
                         print "we get a line of solutions ",twoa/t,"x+(",b/t,")y=",(-d+g)/t,"\n"
                         tcount=1
                    }
                    if((-d-g)%t==0){
                       print "and a line of solutions ",twoa/t,"x+(",b/t,")y=",(-d-g)/t,"\n"
                       tcount=tcount+1
                    }
                    if(tcount>0){
                       return(-2)
                    }else{
                       print "no solutions\n"
                       return(0)
                    }
              }
           }else{
              print "no solutions\n"
              return(0)
           }
        }#end of v>0
      }#end of v!=0
      return
   }#end of if u==0
   # we can now proceed with the knowledge that u is nonzero below.
   # solving (t+d)^2=uy+v, u nonzero
absu=abs(u)
if(flag){
print "absu=",absu,"\n"
print "u=",u,"\n"
print "v=",v,"\n"
print "solve T^2=",v," (mod ",absu,")\n"
}
g=squareroot0(v,absu)
if(g==0){
   print "no solutions\n"
   return(0)
}
tcount=0
for(h=0;h<g;h++){
if(flag){
    print "solution[",h,"]=",solution[h],"\n"
    print "t = ",solution[h]-d,"+(",u,")K, K arbitrary\n"
}
    s=2*solution[h]
    r=(solution[h]^2-v)/u
if(flag){
    print "r=",r,"\n"
    print "s=",s,"\n"
    print "y=",r,"+(",s,")K+(",u,")K^2\n"
}
    j=u-b*s
    ksol=d+b*r-solution[h]
if(flag){
    print "j=",j,"\n"
    print "ksol=",ksol,"\n"
    print "Now solving ",j,"K=",ksol,"(mod ",twoa,")\n"
}
    l=legendrecong(j,ksol,twoa)
    if(l==0){
       if(flag){
          print "l=0\n"
       }
       continue
    }
    x=globalcongsol
    y=globalcongmod
    if(flag){
        print "checking to see if 2a divides jy:\n"
        print "jy%  = ",(j*y)%twoa,"\n" 
        print "K = ",x," + ",y,"w, where w is arbitrary\n"
    }
    temp=2*b*e-4*c*d
    print "SOL[",tcount,"]=",solution[h],"\n"
    sol[tcount]=solution[h]
    print "x=",(x*j-ksol)/twoa+temp*x^2,"+(",y*j/twoa+2*temp*x*y,")w+(",temp*y^2,")w^2\n"
    print "y=",r+s*x+u*x^2,"+(",s*y+2*u*x*y,")w+(",u*y^2,")w^2\n"
    print "\n"
    tcount=tcount+1
}#end of for(h)
if(tcount){
   if(tcount>1){
      print tcount," families of solutions\n"
   }else{
      print "one family of solutions\n"
   }
   print "here is the initial solution array:"
       print_array(sol[],tcount)
   improvement(u,v,sol[],tcount,a,b,d,e)
}else{
   print "no solutions\n"
}
}#end of define ddzero

define reducetest(u,v,d){
auto f
  f=sqrt(d)
  if(v>0 && u>0 && u<=f && f<u+v && v-u<=f){
    return(1)
  }else{
    return(0)
  }
}

# This finds all solutions of x^2 = d (mod n) in the range -n/2 < x <= n/2. 
# outputs global solution[h]. used to be in file juk. Added to patz 1st December 2019.

define squareroot0(d,n){
auto flag1,flag2,g,h,solcount
    
    flag1=0
    flag2=0
    g=sqroot(d,n,0)
    if(g==0){
       return(0)
    }
    solcount=numbr
    if(2*solution[numbr-1]==n){
      flag1=1
    }
    if(solution[0]==0){
      flag2=1
    }
    for(h=0;h<numbr;h++){
       if(flag1==0 && flag2==0){
          solution[numbr+h]=-solution[h]
          solcount=solcount+1
       }
    }
   #for(h=0;h<solcount;h++){
   #    print "solution[",h,"]=",solution[h],"\n"
   #}
   #print "\n"
    return(solcount)
}