/* gnubc program : posformrep */
/* Needs squareroot and reduceneg */

/* 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)
}



/* posrep(a,b,c,n) deals with the more general case of finding the primitive soluitons of
 * 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. Here a>0 and n>0.
 * The number r of primitive solutions is returned.
 */ 
define posrep(a,b,c,n,printflag){
        auto aa,bb,cc,i,r,s,tmp1,tmp2,tmp3,tmp4
	if(gcd(a,n)==1){
           r=posrep0(a,b,c,n,printflag)
           for(i=0;i<r;i++){
               if(printflag){
	          print "solution[",i,"]: (x,y)=(",posrep_solutionx0[i],",",posrep_solutiony0[i],")\n"
               }
               posrep_solutionx[i]=posrep_solutionx0[i]
               posrep_solutiony[i]=posrep_solutiony0[i]
           }
           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.
         */
         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=posrep0(aa,bb,cc,n,printflag)
         if(printflag){
            print "(global_a,global_b,global_c)=(",global_a,",",global_b,",",global_c,")\n"
         }
         if(r==0){
           return(0)
         }
         /* this creates arrays of solutions:
          * posrep_solutionx0[0],...,posrep_solutionx0[r-1] and
          * posrep_solutiony0[0],...,posrep_solutiony0[r-1].
          */
	 for(i=0;i<r;i++){
	     tmp1=gauss_alpha*posrep_solutionx0[i]
	     tmp2=gauss_beta*posrep_solutiony0[i]
	     tmp3=gauss_gamma*posrep_solutionx0[i]
	     tmp4=gauss_delta*posrep_solutiony0[i]
             posrep_solutionx[i]=tmp1+tmp2
             posrep_solutiony[i]=tmp3+tmp4
             if(printflag){
	        print "solution[",i,"]: (x,y)=(",posrep_solutionx[i],",",posrep_solutiony[i],")\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
}

/* We are finding the primitive solutions of the quadratic diophantine equation b.t^2_c.tu+d.u^2=a.
/* Here gcd(b,c,d)=1=gcd(b,a) and c^2-4bc<0. Also b>0,d>0,a>0.
 * We first solve the congruence b.theta^2+c.theta+d=0 (mod a).
 * We let t=theta.u-ay to get p.u^2+quy+ry^2=1, where p=(b.theta^2+c.theta+d)/a,
 * q=-(2b.theta+c) and r=ab. We find the reduced form (global_a,global_b,global_c)
 * which is equivalent to (p,q,r). If global_a>1, we continue the loop. 
 * If global_a=1,  denote the unimodular transformation employed by
 * u=a11.X+a12.Y, y=a21.X+a22.Y. we get two solutions (u,y)=(a11,a21) and (-a11,-a21). 
 * If global_a=global_c, we get two more solutions (u,y)=(a12,a22) and (-a12,-a22).
 * Then we get corresponding solutions (t,u), where t=theta.u-ay.
 * The number of primitive solutions is returned.
 * printflag=1 prints out the intermediate results.
 * See http://www.numbertheory.org/pdfs/posrep.pdf.
 */

define posrep0(b,c,d,a,printflag){
auto theta,p,q,r,q0,k,tt,u,y,solnumber,t,s,denomp

    s=quadratic(b,c,d,a,0)
        /* s=0 means no solutions of bx^2+cx+d=0 (mod a) */
        /* If s>0, we get all solutions x, 0<=x<a as
         * quadratic_solution[0],...,quadratic_solution[s-1]
         */
   if(s==0){
      if(printflag){
         print b,"theta^2+",c,"theta+",d,"=0 (mod ",a,") has no solution\n"
      }
      return(0)
   }
   solnumber=0
   for(k=0;k<s;k++){
      theta=quadratic_solution[k]
      if(printflag){
         print "theta[",k,"]=",theta,"\n"
      }
denomp=b*theta^2+c*theta+d
if(denomp%a){
   print denomp, " is not divisible by ",a,"\n"
   return(0)
}
      p=(b*theta^2+c*theta+d)/a
      q=-2*b*theta-c
      r=a*b
      t=reduce1(p,q,r) /* produces  global variables (global_a, global_b,global_c) and
                          unimodular transforming matrix 
                          [globala11,globala12,globala21,globala22]. */
      if(printflag){
         print "(a11,a12,a21,a22)=(",globala11,",",globala12,",",globala21,",",globala22,")\n"
         print "(global_a,global_b,global_c)=(",global_a,",",global_b,",",global_c,")\n"
      }
      if(global_a>1){
        continue
      }else{
           u=globala11
           y=globala21
           tt=theta*u-a*y
           if(printflag){
              print "solution0[",solnumber,"]: (t,u)=(",tt,",",u,")\n"
           }
           posrep_solutionx0[solnumber]=tt
           posrep_solutiony0[solnumber]=u
           solnumber=solnumber+1
           if(printflag){
              print "solution0[",solnumber,"]: (t,u)=(",-tt,",",-u,")\n"
           }
           posrep_solutionx0[solnumber]=-tt
           posrep_solutiony0[solnumber]=-u
           solnumber=solnumber+1
           if(global_a==global_c){
              u=globala12
              y=globala22
              tt=theta*u-a*y
              if(printflag){
                 print "solution0[",solnumber,"]: (t,u)=(",tt,",",u,")\n"
              }
              posrep_solutionx0[solnumber]=tt
              posrep_solutiony0[solnumber]=u
              solnumber=solnumber+1
              if(printflag){
                 print "solution0[",solnumber,"]: (t,u)=(",-tt,",",-u,")\n"
              }
              posrep_solutionx0[solnumber]=-tt
              posrep_solutiony0[solnumber]=-u
              solnumber=solnumber+1
              if(global_b==1){
                 u=globala11-globala12
                 y=globala21-globala22
                 tt=theta*u-a*y
                 if(printflag){
                    print "solution0[",solnumber,"]: (t,u)=(",tt,",",u,")\n"
                 }
                 posrep_solutionx0[solnumber]=tt
                 posrep_solutiony0[solnumber]=u
                 solnumber=solnumber+1
                 if(printflag){
                    print "solution0[",solnumber,"]: (t,u)=(",-tt,",",-u,")\n"
                 }
                 posrep_solutionx0[solnumber]=-tt
                 posrep_solutiony0[solnumber]=-u
                 solnumber=solnumber+1
              }
           }
      }
   }
   return(solnumber)
}

/* 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. */
define posrepgen(a,b,c,n){
auto absn,t,g,countpatzgen,temp,h,f
   if(a<0){
      a=-a;b=-b;c=-c;n=-n;
   }
   if(n<0){
     print "no solutions\n"
     return(0)
   }
   if(n==0){
     print "one solution (x,y)=(0,0)\n"
     globalposrepx[0]=0
     globalposrepy[0]=0
     return(1)
   }
   absn=abs(n)
   if(absn!=1){
      t=omega(absn)
      g=divisors(qglobal[],kglobal[],t)
   }else{
      g=1
      divisor[0]=1
   }
   countpatzgen=0
   for(f=0;f<g;f++){
      temp=divisor[f]^2
      if(n%temp==0){
        #print "divisor[",f,"]=",divisor[f],"\n"
        #print "solving ", a,"x^2+(",b,"xy+(",c,")y^2=",n/temp,"\n"
         t=posrep(a,b,c,n/temp,0)
         if(t>0){
         for(h=0;h<t;h++){
            globalposrepx[countpatzgen+h]=divisor[f]*posrep_solutionx[h]
            globalposrepy[countpatzgen+h]=divisor[f]*posrep_solutiony[h]
         }
         countpatzgen=countpatzgen+t
         }
      }
   }
   return(countpatzgen)
}

/* This lists the solutions of ax^2+bxy+cy^2=n, b^2-4ac< 0, a>0, n>0. */
define posrepgenlist(a,b,c,n){
auto h,g
   g=posrepgen(a,b,c,n)
   for(h=0;h<g;h++){
     print "("
     print globalposrepx[h],",",globalposrepy[h]
     print ")\n"
   }
   return(g)
}

/* This solves the diophantine equation ax^2+bxy+cy^2+dx+ey+f=0, where a>0, b^2-4ac<0. 
 * We use the transformation of Hua's book, p.278.
 */
define ssw1(a,b,c,d,e,f){
auto dd,gcount,k,u,v,h,r,s,g,kk

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"
u=b*e-2*c*d
v=b*d-2*a*e
k=-(a*u^2+b*u*v+c*v^2-d*dd*u-e*dd*v+f*dd^2)
kk=-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)
}
print dd,"x=X-(",u,")," 
print dd,"y=Y-(",v,")\n" 
print "Find the solutions of ",a,"X^2+",b,"XY+",c,"Y^2=",k,"\n"
g=posrepgen(a,b,c,k)
print "There are ",g," solutions (X,Y)\n"
for(h=0;h<g;h++){
    r=globalposrepx[h]-u
    s=globalposrepy[h]-v
    if(r%dd==0 && s%dd==0){
       print "(x,y)=(",r/dd,",",s/dd,")\n"
       gcount=gcount+1
    }
}
print "number of solutions (x,y) found is "
return(gcount)
}