/* gnubc program : patz */
/* Needs squareroot */

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

define gcd3(a,b,c){
auto t
t=gcd(a,b)
t=gcd(t,c)
return(t)
}

/* posrep(a,b,c,n,printflag) 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. Here a>0 and n>0.
 * The number 2*r of solutions is returned.
 */ 
define posrep(a,b,c,n,printflag){
        auto aa,bb,cc,i,r,s,tmp1,tmp2,tmp3,tmp4,discriminant
       if(gcd3(a,b,c)>1){
           print "gcd(a,b,c)>1\n"
           return(-1)
       }
        if(a<=0){
           print "a<=0\n"
           return(-1)
        }
        discriminant=b^2-4*a*c
        if(printflag){
           print "discriminant=",discriminant,"\n"
        }
        if(discriminant>=0){
           print "discriminant >=0\n"
           return(-1)
        }
	if(gcd(a,n)==1){
           r=posrep1(a,b,c,n,printflag)
           for(i=0;i<r;i++){
	       print "solution[",2*i,"]: (t,u)=(",posrep_solutionx0[i],",",posrep_solutiony0[i],")\n"
	       print "solution[",2*i+1,"]: (t,u)=(",-posrep_solutionx0[i],",",-posrep_solutiony0[i],")\n"
           }
           return(2*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=posrep1(aa,bb,cc,n,printflag)
         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
	     print "solution[",2*i,"]: (t,u)=(",posrep_solutionx[i],",",posrep_solutiony[i],")\n"
	     print "solution[",2*i+1,"]: (t,u)=(",-posrep_solutionx[i],",",-posrep_solutiony[i],")\n"
         }
         return(2*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=chinesea(xx[],qglobal[],e)
       tmp2=chinesea(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 solving 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. In the case of solubility,
 * we get one, two or three non-negative solutions, according as D < -4, D=-4 or D= -3.
 * We first deals with exceptional cases: otherwise we express -q/2p as a continued fraction
 *  using Euclid's algorithm. 
 * We only consider convergents agloba[i]/bgloba[i] where bglobal[i]<= sqrt(4p/-D). 
 * We test f[i]=p*aglobal[i[^2+q*aglobal[i]*bglobal[i[+c*bglobal[i]^2 in qformvalues()
 * to see if it hits the value 1.  If not, there is no solution (u,y) for this value of theta 
 * and we proceed to the next value of theta.  If we do hit a value 1, and D < -4, we proceed 
 * to the next value of theta; if D=-4, we know we will have f[i+1]=1 and then we
 * exit the while loop;  finally, if D=-3, we know that f[i+1]=1 and f[i+2]=1 and after considering  * convergents aglobal[i]/bglobal[i], aglobal[i+1]/bglobal[i+1],aglobal[i+2]/bglobal[i+2], we exite  * the while loop. At each successful exit, we calculate and record the solutions (t,u) of the 
 * original equation.
 * printflag=1 prints out the intermediate results.
 */

define posrep1(b,c,d,a,printflag){
auto theta,p,q,r,q0,k,tt,u,y,solnumber,t,s,dd,n,number,aover2
    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){
      print "b.theta^2+c.theta+d=0 (mod a) has no solution\n"
      return(0)
   }
   dd=c^2-4*b*d
   solnumber=0
   aover2=a/2
   for(k=0;k<s;k++){
      theta=quadratic_solution[k]
      if(theta>aover2){
         theta=theta-a
      }
      if(printflag){
         print "theta[",k,"]=",theta,"\n"
      }
      p=(b*theta^2+c*theta+d)/a
      q=-2*b*theta-c
      r=a*b
      print "(p,q,r)=(",p,",",q,",",r,")\n"
      if(dd<-4 && p==1){
         if(printflag){
            print "exceptional case D < -4 and P = 1\n"
         }
         u=1
         y=0
         tt=theta
         if(printflag){
            print "solution0[",solnumber,"]: (t,u)=(",tt,",",u,")\n"
         }
          posrep_solutionx0[solnumber]=tt
          posrep_solutiony0[solnumber]=u
          solnumber=solnumber+1
          continue
      }
      if(dd==-4){
         if(p==1){
            if(printflag){
               print "exceptional case D = -4 and P = 1\n"
            }
            u=1
            y=0
            tt=theta
            if(printflag){
               print "solution0[",solnumber,"]: (t,u)=(",tt,",",u,")\n"
            }
            posrep_solutionx0[solnumber]=tt
            posrep_solutiony0[solnumber]=u
            solnumber=solnumber+1
            u=-q/2
            y=1
            tt=theta*u-a
            if(printflag){
               print "solution0[",solnumber,"]: (t,u)=(",tt,",",u,")\n"
            }
            posrep_solutionx0[solnumber]=tt
            posrep_solutiony0[solnumber]=u
            solnumber=solnumber+1
            continue
        }
        if(p==2){
           if(printflag){
              print "exceptional case D = -4 and P = 2\n"
           }
           n=q/2
           u=(-n+1)/2
           y=1
           tt=theta*u-a
           if(printflag){
              print "solution0[",solnumber,"]: (t,u)=(",tt,",",u,")\n"
           }
           posrep_solutionx0[solnumber]=tt
           posrep_solutiony0[solnumber]=u
           solnumber=solnumber+1
           u=(-n-1)/2
           y=1
           tt=theta*u-a
           if(printflag){
              print "solution0[",solnumber,"]: (t,u)=(",tt,",",u,")\n"
           }
           posrep_solutionx0[solnumber]=tt
           posrep_solutiony0[solnumber]=u
           solnumber=solnumber+1
           continue
        }
      }/*end of dd=-4 section*/
      if(dd==-3){
        n=(q-1)/2
        if(p==1){
           if(printflag){
              print "exceptional case D = -3 and P=1\n"
            }
            u=1
            y=0
            tt=theta*u
            if(printflag){
               print "(u,y)=(1,0)\n"
               print "solution0[",solnumber,"]: (t,u)=(",tt,",",u,")\n"
            } 
            posrep_solutionx0[solnumber]=tt
            posrep_solutiony0[solnumber]=u
            solnumber=solnumber+1
            u=-n
            y=1
            tt=theta*u-a
            if(printflag){
               print "(u,y)=(",u,",",y,")\n"
               print "solution0[",solnumber,"]: (t,u)=(",tt,",",u,")\n"
            }
            posrep_solutionx0[solnumber]=tt
            posrep_solutiony0[solnumber]=u
            solnumber=solnumber+1
            u=-n-1
            y=1
            tt=theta*u-a
            if(printflag){
               print "(u,y)=(",u,",",y,")\n"
               print "solution0[",solnumber,"]: (t,u)=(",tt,",",u,")\n"
            }
            posrep_solutionx0[solnumber]=tt
            posrep_solutiony0[solnumber]=u
            solnumber=solnumber+1
            continue
        }
        if(p==3){
           if(printflag){
              print "exceptional case D = -3 and P=3\n"
            }
            u=(-n+1)/3
            y=1
            tt=theta*u-a
            if(printflag){
               print "(u,y)=(",u,",",y,")\n"
               print "solution0[",solnumber,"]: (t,u)=(",tt,",",u,")\n"
            }
            posrep_solutionx0[solnumber]=tt
            posrep_solutiony0[solnumber]=u
            solnumber=solnumber+1
            u=-q/3
            y=2
            tt=theta*u-a*y
            if(printflag){
               print "(u,y)=(",u,",",y,")\n"
               print "solution0[",solnumber,"]: (t,u)=(",tt,",",u,")\n"
            }
            posrep_solutionx0[solnumber]=tt
            posrep_solutiony0[solnumber]=u
            solnumber=solnumber+1
            u=-(n+2)/3
            y=1
            tt=theta*u-a
            if(printflag){
               print "(u,y)=(",u,",",y,")\n"
               print "solution0[",solnumber,"]: (t,u)=(",tt,",",u,")\n"
            }
            posrep_solutionx0[solnumber]=tt
            posrep_solutiony0[solnumber]=u
            solnumber=solnumber+1
            continue
        }
      }/*end of dd=-3 section*/
      number=qformvalues(p,q,r,a,theta,solnumber,printflag)
      solnumber=solnumber+number
  }
  print "number of solutions is ", 2*solnumber,"\n"
  return(solnumber)
}

/* this examines the continued fraction for -q/2p, where convergents B[i]=by 
 * do not exceed Lagrange's upper bound sqrt(4p/-disc). It is only used for non-exceptional 
 * solutions. See the algorithm at Section 7 of the author's paper.
 */
define qformvalues(p,q,r,a,theta,solnumber,printflag){
auto i,m,n,t,num,den,x,y,d,fourp,disc,aa,bb,ax,ay,bx,by,az,bz,qq,u,rr,fourpoverdisc
disc=4*p*r-q^2
d=-disc
fourp=4*p
m=-q
n=2*p
if(printflag){
   print "-q/2p=",m,"/",n,"\n"
}
i=0
aa=m;bb=n
ax=0;bx=1;ay=1;by=0
fourpoverdisc=fourp/disc
while(by^2<=fourpoverdisc){
   if(i){
     aa=bb
     bb=rr
   }
   rr=mod(aa,bb)
   qq=(aa-rr)/bb
   az=qq*ay+ax
   bz=qq*by+bx
   u=az;y=bz
   ax=ay;ay=az
   bx=by;by=bz
   f[i]=p*u^2+q*u*y+r*y^2
   if(printflag){
      print "A[",i,"]/B[",i,"]=",u,"/",y,": "
      print "f[",i,"]=",f[i],"\n"
   }
   if(f[i]==1){
     posrep_solutionx0[solnumber]=u*theta-a*y
     posrep_solutiony0[solnumber]=u
     if(printflag){
        print "(u,y)=(",u,",",y,")\n"
     }
     print "solution0[",solnumber,"]: (t,u)=(",posrep_solutionx0[solnumber],",",u,")\n"
     solnumber=solnumber+1
     if(d<-4){
        return(1)
     }
     if(d==-4 || d==-3){
        aa=bb
        bb=rr
        rr=aa%bb
        qq=(aa-rr)/bb
        az=qq*ay+ax
        bz=qq*by+bx
        u=az;y=bz
        ax=ay;ay=az
        bx=by;by=bz
        if(printflag){
           i=i+1
           print "A[",i,"]/B[",i,"]=",u,"/",y,": "
           print "(u,y)=(",u,",",y,")\n"
        }
        posrep_solutionx0[solnumber]=u*theta-a*y
        posrep_solutiony0[solnumber]=u
        print "solution0[",solnumber,"]: (t,u)=(",posrep_solutionx0[solnumber],",",u,")\n"
        solnumber=solnumber+1
        if(d==-4){
           return(2)
        }
        if(d==-3){
           aa=bb
           bb=rr
           rr=aa%bb
           qq=(aa-rr)/bb
           az=qq*ay+ax
           bz=qq*by+bx
           u=az;y=bz
           if(printflag){
              i=i+1
              print "A[",i,"]/B[",i,"]=",u,"/",y,": "
              print "(u,y)=(",u,",",y,")\n"
           }
           posrep_solutionx0[solnumber]=u*theta-a*y
           posrep_solutiony0[solnumber]=u
           print "solution0[",solnumber,"]: (t,u)=(",posrep_solutionx0[solnumber],",",u,")\n"
           solnumber=solnumber+1
           return(3)
        }
     }
   }
   i=i+1
}
if(printflag){
   print "no solution found for theta=",theta,"\n"
}
return(0)
}