# Here we create a list of Nagell fundamental solutions and output a list of 
# global Stolt fundamental solutions of x^2-dy^2=n, where 4 divides n and n is nonzero.
# Also d > 0 and is nonsquare.
# We return the number of Stolt fundamental solutions.
# There are three possibilities: Let e1=x1+y1sqrt(d) be the least positive solution of
# x^2 - dy^2 = 1 and e4 = (x4+y4sqrt(d))/2 be the least positive solution of
# x^2 - dy^2 = 4. if e1=e4, then the Stolt and Nagell classes are identical.
# if e1=e4^2, each Stolt class is composed of two Nagell classes,
# if e1=e4^3, each Stolt class is composed of three Nagell classes.
# function e1versuse4(d) tells us when these occur.
# This method of finding the Stolt fundamental solutions is much quicker that the
# function stolt0 in file stolt.
# type bc stoltvialmm stoltarrange.bc squareroot patz genfacs

define stoltvialmm(d,n,printflag){
auto h,s,gpsize,i,j,x[],y[]

h=patzgen(d,n)
# if h>0, this creates arrays 
# globalfundxff[0],...,globalfundxff[h-1] and
# globalfundyff[0],...,globalfundyff[h-1].
if(h==0){
  print "no solutions\n"
  return(0)
}
s=fund4(d,0)
if(printflag){
   print "globalx=",globalx,"\n"
   print "globaly=",globaly,"\n"
}
gpsize=e1versuse4(d)
if(printflag){
   print "gpsize=",gpsize,"\n"
}
if(gpsize==1){
   if(printflag){
      print "Stolt and Nagell classes are identical\n";
   }
   for(j=0;j<h;j++){
      globalstoltvialmmx[j]=globalfundxff[j]
      globalstoltvialmmy[j]=globalfundyff[j]
      if(printflag){
         print "(globalstoltviammx[",j,"],globalstoltviammy[",j,"]=(",globalfundxff[j],",",globalfundyff[j],")\n"
      }
   }
   return(h)
}
if(printflag){
   print "Here are the Nagell fundamental solutions of x^2 - ",d,"y^2=",n,":\n"
      for(i=0;i<h;i++){
         print "(globalfundxff[",i,"],globalfundyff[",i,"])=(",globalfundxff[i],",",globalfundyff[i],")\n"
      }
}

#if(globalx/2==globalr && globaly/2==globals){
#   print "e1=e4\n"
#   return(1)
#}
#powerdd(globalx,globaly,d,2)
#print "z1=",z1,"\n"
#print "z2=",z2,"\n"
#
#if(z1/4==globalr && z2/4==globals){
#   print "e1=e4^2\n"
#   gpsize=2
#}else{
#   print "e1=e4^3\n"
#   gpsize=3
#}
#
for(i=0;i<h;i++){
  x[i]=globalfundxff[i]
  y[i]=globalfundyff[i]
}

arrange(x[],y[],h,gpsize,n/2)
for(j=0;j<h;j++){
   if(printflag){
      print "after arrange: (x[",j,"],y[",j,"]=(",x[j],",",y[j],")\n"
   }
}
if(printflag){
   print "\n"
}
j=0
for(i=0;i<h;i=i+gpsize){
 if(gpsize==2){
   if(y[i]<y[i+1]){
     index[j]=i
   }
   if(y[i]>y[i+1]){
     index[j]=i+1
   }
   if(y[i]==y[i+1]){
     if(x[i]<0){
       index[j]=i+1
     }else{
       index[j]=i
     }
   }
 }
 if(gpsize==3){
   if(y[i]<y[i+1] && y[i]<y[i+2]){
     index[j]=i
   }
   if(y[i+1]<y[i] && y[i+1]<y[i+2]){
     index[j]=i+1
   }
   if(y[i+2]<y[i] && y[i+2]<y[i+1]){
     index[j]=i+2
   }
   if(y[i]==y[i+1]){
     if(y[i]<y[i+2]){
       if(x[i]<0){
         index[j]=i+1
       }else{
         index[j]=i
       }
     }
     if(y[i]>y[i+2]){
         index[j]=i+2
     }
   }
   if(y[i+1]==y[i+2]){
     if(y[i]<y[i+1]){
         index[j]=i
     }
     if(y[i]>y[i+1]){
        if(x[i+1]>=0){
          index[j]=i+1
        }else{
          index[j]=i+2
        }
     }
   }
   if(y[i]==y[i+2]){
     if(y[i+1]<y[i]){
         index[j]=i+1
     }
     if(y[i+1]>y[i]){
        if(x[i]>=0){
          index[j]=i
        }else{
          index[j]=i+2
        }
     }
   }
 }
 globalstoltvialmmx[j]=x[index[j]]
 globalstoltvialmmy[j]=y[index[j]]
 j=j+1
}
h=h/gpsize
if(printflag){
   for(j=0;j<h;j++){
       print "(globalstoltvialmmx[",j,"],globalstoltvialmmy[",j,"]=(",globalstoltvialmmx[j],",",globalstoltvialmmy[j],")\n"
   }
   print "h=",h,"\n"
   print "gpsize=",gpsize,"\n"
}
return(h)
}

# This function returns 1 if e1=e4, 2 if e1=e4^2 and 3 if e1=e4^3. 
# This is from Stolt's paper on the Diophantine equatioins u^2-Dv^2=\pm4Nm
# Arkiv for Mathematik, Band 2, Nr 1, 1951, 1-23.
define e1versuse4(d){
auto s,y4
  if(d%8==1 || d%4 == 2 || d%4==3){
     return(1)
  }
  s=fund4(d,0)
  y4=globaly
  if(d%8==5){
    if(y4%2==0){
      return(1)
    }else{
      return(3)
    }
  }
  if(d%4==0){
    if(y4%2==0){
      return(1)
    }else{
      return(2)
    }
  }
}


# We use completion of the square to solve ax^2 + bxy + cy^2 = n, d=b^2-4ac>0 and squarefree.
# We reduce the equation to X^2 -dY^2=4an, where X=2ax+y and Y=y and then solve this equation
# to get the Stolt fundamental solutions.
define binaryviasfs(a,b,c,n,printflag){
auto d,fouran,alpha,beta,twoa,t,count,x,y,u,temp,temp1,temp2,g,h,absa,j
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)
}
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
twoa=2*a
d=b^2-4*a*c
g=stoltvialmm(d,fouran,printflag)
if(printflag){
   print g," is the number of Stolt fundamental classes of x^2 - ",d,"y^2=",fouran,"\n"
}
absa=abs(a)
   count=0
   for(h=0;h<g;h++){
      alpha=globalstoltvialmmx[h]
      beta=globalstoltvialmmy[h]
      if((alpha-b*beta)%twoa==0){
         #print "alpha-b*beta is divisible by 2a\n"
         y=globalbinaryviastolty[count]=beta
         x=(alpha-b*beta)/twoa
         temp1=-(a*x+b*beta)
         if(temp1%absa==0){
           u=temp1/a
           temp2=equivfunc(u,beta,alpha,beta,n)
           if(temp2==1){
              if(u>x){
                 x=u
              }
              if(printflag){
                 print "interchanging u and x\n"
              }
           }
         }
         globalbinaryviastoltx[count]=x
         if(printflag){
            print "Stolt fundamental solution[",count,"]=(",x,",",y,")\n"

            print "(",x,",",beta,")\n"
            print"\tx=("
            temp=-(b*x+2*c*y)
            t=printauplusbv(x,temp)
            print")/2\n"
            print"\ty=("
            temp=2*a*x+b*y
            t=printauplusbv(y,temp)
            print")/2\n"
         }
         count=count+1
      }
   }
   if(printflag){
      if(count){
         print "Here u^2-",d,"v^2=4\n"
      }
   }
   if(printflag){
      for(j=0;j<count;j++){
         print "(globalbinaryviastoltx[",j,"],globalbinaryviastolty[",j,"]=(",globalbinaryviastoltx[j],",",globalbinaryviastolty[j],")\n"
      }
      print "the number of solution families of \n"
      print a,"x^2+",b,"xy+",c,"y^2=",n," is ",count,"\n"
   }
   return(count)
}