# Find a solution of a quadratic systems of equations using weak
# Grobner-type reductions, given that the variety is 0-dimensional,
# and that all solutions have coordinates of the form sqrt(rationals).
# If there is a solution, return a linear system for it in triangular form

quadsolve:=proc(eq) local red,u,sa,v,new;
  red:=[`quadsolve/reduce`(eq)];
  if nops(red)=0 then RETURN() fi;
  for u in red[1] do
    sa:=traperror(factors(u)[2]);
    if sa=lasterror then next    #this trap is for a bug in Maple9
      elif sa[1][2]=2 or nops(sa)>1 then
        sa:=map(x->op(1,x),sa)
      elif nops(indets(sa[1]))=1 then
        v:=indets(sa[1])[1];
        sa:=map((x,y)->y-x,[solve(sa[1][1])],v);
      else next
    fi;
    for v in sa do
      new:=procname({v,op(red[1])});
      if nops([new])>0 then RETURN({op(new),op(red[2])}) fi
    od;
    RETURN()
  od;
  if nops(red[1])>0 then
    if nargs>1 then assign(args[2],red) fi;
    ERROR(`full reduction failed`)
  else red[2] fi;
end;

# given a set of quadratic expressions, eliminate the first
# variable in each linear expression, and reduce a quadratic
# if it has a leading term that matches that of another quadratic,
# and repeat until no more reductions of these types are possible

`quadsolve/reduce`:=proc(eq) local foo,u,i,j,v,lm,co,lin,new,notdone;
  foo:=[op(eq)]; lin:=NULL; notdone:=true;
  while notdone do
    foo:=`quadsolve/linreduce`(foo,'new');
    if foo=NULL then RETURN() else lin:=lin,op(new) fi;
    foo:=sort(foo,(x,y)->evalb(nops(indets(x))<=nops(indets(y))));
    lm:=[]; notdone:=false;
    for i to nops(foo) do
      u:=foo[i]; v:=`quadsolve/leadm`(u,'co');
      while member(v,lm,'j') and v<>1 do
        if not type(co^2,'rational') then ERROR(`bad residuals`) fi;
        u:=simplify(expand(u/co-foo[j]),'sqrt');
        v:=`quadsolve/leadm`(u,'co'); notdone:=true;
      od;
      if not type(co^2,'rational') then ERROR(`bad residuals`)
        elif co<>1 and v<>1 then u:=simplify(expand(u/co),'sqrt') fi;
      foo:=subsop(i=u,foo); lm:=[op(lm),v];
    od;
  od; foo,{lin};
end:

# Given a set of polynomials defining an ideal, eliminate the first
# variable in each linear polynomial, and repeat until there are no
# more linear polynomials available. Assign to args[2] the set of linear
# equations used to do the eliminations. Return NULL if the system has
# a (linear) inconsistency.
# It would be nice to use Maple's solve, but they haven't even bothered
# to correctly implement linear algebra over quadratic number fields. ;-(

`quadsolve/linreduce`:=proc(eq) local foo,lin,new,v,co;
  foo:=eq; lin:=NULL;
  do
    foo:=subs(0=NULL,simplify(expand(foo),'sqrt'));
    new:=map(proc(x) if degree(x)<2 then x fi end,[op(foo)]);
    if nops(new)=0 then assign(args[2],{lin}); RETURN(foo) fi;
    new:=sort(new,(x,y)->evalb(nops(indets(x))<=nops(indets(y))));
    v:=`quadsolve/firstvar`(new[1]);
    if v=1 then RETURN() else co:=coeff(new[1],v) fi;
    if not type(co^2,'rational') then ERROR(`bad pivots`) fi;
    new:=simplify(expand(v-new[1]/co),'sqrt');
    lin:=lin,v=new; foo:=subs(v=new,foo);
  od;
  ERROR(`this cannot happen`);
end:

# Given a triangular set of linear equations of the form 'var=other stuff',
# solve the system by substitution, and then apply sqrt simplifications.

tri_solve:=proc(eq) local res,new,sys;
  res:=NULL; sys:=eq;
  do
    new:=map(x->if nops(indets(x))=1 then x fi,sys);
    sys:=map(x->if nops(indets(x))>0 then x fi,subs(new,sys));
    res:=res,op(new);
    if nops(sys)=0 then break elif nops(new)=0 then
      ERROR(`not triangular`) fi;
  od;
  simplify({res},'sqrt');
end;

# extract the "first" variable that appears in an expression

`quadsolve/firstvar`:=proc(eq) local vars,u,v;
  vars:=indets(eq);
  if nops(vars)=0 then v:=1 else v:=vars[1] fi;
  for u in vars do
    if op(1,u)<op(1,v) then v:=u
      elif op(1,u)=op(1,v) and op(2,u)<op(2,v) then v:=u fi
  od; v
end:

# extract the leading monomial (and coefficient) of a polynomial

`quadsolve/leadm`:=proc(u) local a,d,res,v;
  a:=u; d:=degree(a); res:=1;
  if type(a,`+`) then a:=map(proc(x,y) if degree(x)=y
    then x else 0 fi end,a,d) fi;
  while degree(a)>0 do
    v:=`quadsolve/firstvar`(a); d:=degree(a,v);
    res:=res*v^d; a:=coeff(a,v,d); 
  od;
  if nargs>1 then assign(args[2],a) fi; res
end: