#
# Maple programs for counting points on varieties over finite fields 
# that are defined by "highly linear" integer polynomials. Useful for
# testing Kontsevich's conjecture on the spanning tree generating
# function of a graph. 
#
# Should work with Maple V.2,3,4,5. Tested thoroughly with V.3.
#
# The fundamental data structure used here is the "Z-expression";
# i.e., a Maple expression of the form
#
#     c * q^d * Z[f_1,....,f_k],
#
# where c and d are integers, and the f_i's are multivariate integer
# polynomials. By definition, Z[f_1,....,f_k] represents the probability
# that a randomly chosen evaluation of the f_i's in F_q is zero. Thus if
# there are m variables, q^m*Z[f_1,....,f_k] is the number of 0's over F_q.
#
# A compound Z-expression is a sum of Z-expressions.
#
# For example,  q^2*Z[x^2-y^2];   Z[] ( = 1);  Z[x^2+y^2,x-y] - Z[2,x]
#
# Note: Z and q are used as global variables.
#
#########################
#
# Main procedure defined in this file:
#
# reduce  attempt to completely reduce any compound Z-expression,
#         combining randomized, greedy, and conservative methods.
#
# Examples:
#   f:=Z[x^2-y^2];  reduce(f);
#   reduce(Z[x*y+y*z+x*z]);
#
#########################
#
# FLAGS
#
# If 'savefile' is assigned a name, say foobar, then with each round
# of reduction, files named foobar1.m, foobar2.m, ... etc are saved
# with the content of the file being a single variable 'f' that is
# assigned the current compound Z-expression. This is useful when
# running a long computation that may need to be interrupted.

savefile:=NULL;        # this is NULL, so no saving will occur

# Normally reduce() runs in "quiet" mode. But there are two levels of
# verbosity available.

verbosity:=0;          # by default we don't want verbosity

# With verbosity=1, one line is printed during each round of reduction.
# The format of this diagnostic information is
#    round  nops  shape  time  space
# where
#  round = number of rounds of reductions so far,
#  nops  = number of operands in the current compound expression
#  shape = two metrics describing the "shape" (see docs below)
#  time  = number of cpu seconds used for this round
#  space = number of bytes allocated to Maple
#
# With verbosity=2, the progress of intermediate (greedy) reductions are
# reported. When attempting to reduce a compound Z-expression,
# hedge_reduce() (see below) is applied one at a time to the individual
# Z-expressions. The progress of these reductions is reported in
# the same way as above, but indented to the right by one field.
# If one of these individual reductions fails for any reason, a line
# of *'s is printed, and conservative (slow) methods take over.
# A number of Warnings about un-nice things (such as inert Z-expressions)
# are lprinted in this case as well.
#
################################
#
# Major subroutines
#
# reduce_one    given a single Z-expression, attempt an application of
#               one of the reduction relations
# hedge_reduce  given a single Z-expression, first try to reduce completely
#               If this fails, apply one round of miserly reduction.
# elim          eliminate a variable that occurs linearly in some term of
#               a Z-expression
# mysplit       simplify a Z-expression in which one of the terms has
#               a nontrivial factorization
# sqfree        the product of the distinct irreducible factors of a
#               multivariate polynomial over Z.
#
########################################################################
 
# define a space-used reporting proc. Depends on Maple version. (Yuck.)

if `+`(0)=0 then space:=proc() kernelopts(bytesalloc) end
  else space:=proc() 4*status[2] end fi;

# A proc for reporting the size and shape of intermediate results.
# The "shape" of f is a,m, where a is the average number of bytes in the
# operands of f, and m is the maximum. (Broken if f has 0 ops.)

shape:=proc(f) round(length(f)/nops(f)),max(op(map(length,[op(f)]))) end;

# A more useful version of irreduc() than the one in Maple.
# Define scalars to be "irreducible", as well as any polynomial
# that is too big for Maple to handle. On the other hand, we want
# any polynomial with content > 1 to NOT be "irreducible".

irred:=proc(f) local t;
  if degree(f)<1 then true
    elif type(f,`*`) or content(f)>1 then false
    else t:=traperror(irreduc(f));
    if t=lasterror then true else t fi
  fi
end;

# Return the square-free part of f; i.e., the product of the distinct
# irreducible factors of f over Z.

sqfree:=proc(f) local g;
  if not type(f,polynom(integer)) then ERROR(`bad input`,f) fi;
  g:=sqrfree(f);
  g:=abs(g[1])*convert(map(u->op(1,u),op(2,g)),`*`);
  if type(g,polynom(integer)) then g else ERROR(`bad output`,g) fi;
end;

# Simplify a Z-expression in which one of the terms has a nontrivial
# factorization. Parameters:
#   F = a list of multivariate polynomials over Z,
#   j = an index such that F[j] factors nontrivially,
#   d = the coefficient of the Z-expression where F came from
# Note that Maple sometimes returns factors with rational, non-integer
# coefficients. Doh! Must work with primitive parts.

mysplit:=proc(F,j,d) local g,h,i,c;
  if type(F[j],`*`) then
    g:=F[j]; c:=abs(content(g)); h:=c*primpart(op(1,g));
    g:=convert([seq(primpart(op(i,g)),i=2..nops(g))],`*`);
  else
    g:=factors(F[j]); c:=abs(content(F[j]));
    g:=map(u->primpart(u[1]),g[2]);
    if nops(g)>1 then
      h:=c*g[1]; g:=convert([seq(g[i],i=2..nops(g))],`*`);
    elif nops(g)=1 and c>1 then
      h:=c; g:=g[1];
    else g:=[op(g),1];
      RETURN(d*Z[op(subsop(j=c*g[1],F))]) 
    fi;
  fi; 
  d*Z[op(subsop(j=g,F))]+d*Z[op(subsop(j=h,F))]
    -d*Z[op(subsop(j=(g,h),F))]
end;

# Given that we haven't found any opportunity to apply mysplit or elim,
# clean up our Z-expression and put it into canonical form:
#   each term must have positive leading coeff, duplicates must be
#   removed, and the scalar (if any) must be prime.
# With an optional third argument (a flag), reject with an ERROR any
# irreducible Z-expression other than Z[] or Z[m] for some integer m.
# (Used in conservative mode for deciding complete reducibility.)

endplay:=proc(F,d) local i,G,m,g;
  G:=NULL;
  for i to nops(F) do
    g:=F[i]; if lcoeff(g)<0 then g:=-g fi;
    if not member(g,[G]) then G:=G,g fi
  od; G:=[G];
  if nops(G)=nops(F) and not type(G,list(integer)) then
    if nargs>2 then ERROR(`too ugly`)
      elif verbosity>1 then lprint(`WARNING: irreducible: `,G) fi
  fi;
  if type(G[1],'integer') then
    m:=numtheory[factorset](G[1]);
    convert([seq(d*Z[op(subsop(1=i,G))],i=m)],`+`);  
  else
    d*Z[op(G)]
  fi;
end;

# The preferred ordering of members of a polynomial list. Use with sort().

pref:=proc(u,v)
  if degree(u)=degree(v) then evalb(length(u)>=length(v))
    else evalb(degree(u)<=degree(v)) fi
end;

# Eliminate a variable that occurs linearly in some polynomial. Parameters:
#  args[4] = (unnamed) a list of polynomials
#     y    = the variable
#     i    = the index of args[4] where y occurs linearly
#     d    = the coeff. of the Z-expression where the list came from
# With an optional 5th argument (a flag), trigger an error whenever
# the objects get too big (>10^5 bytes).

elim:=proc(y,i,d) local F,g,G,u;
  F:=args[4]; g:=collect(F[i],y);
  g:=[coeff(g,y),subs(y=0,g)];
  F:=subsop(i=NULL,F);
  G:=[seq(collect(numer(subs(y=y/q,u)),[q,y]),u=F)];
  G:=map(sqfree,subs(y=-g[2],q=g[1],G));
  g:=map(sqfree,g);
  if nargs>4 and length(g[1])+length(G)>100000 then
    ERROR(`object too messy`) fi;
  d*Z[op(g),op(F)]+(d/q)*Z[op(G)]-(d/q)*Z[g[1],op(G)];
end;

# Attempt to apply one round of reduction to a single Z-expression,
# choosing a variable at random to eliminate. If this is not possible,
# try a "mysplit". 
# With a second argument (a flag), try to find the "best" reduction
# available, measuring success by the length of the generated expression.
# (Can be expensive!)

reduce_one:=proc(a) local d,F,y,i,u,vars,b,b0,cost,st;
  F:=op(indets(a,'indexed'));
  if nops([F])<>1 then ERROR(`bad input`,a) fi;
  d:=coeff(a,F); F:=subs(0=NULL,[op(F)]);
  if nops(F)=0 then RETURN(d*Z[]) fi;
  if member(1,F) or member(-1,F) then RETURN(0) fi;
  F:=sort(F,pref); u:=0;
  for i to nops(F) while degree(F[i])=0 do u:=igcd(u,F[i]) od;
  if abs(u)=1 then RETURN(0) elif u=0 then u:=NULL
    else u:=convert(numtheory[factorset](u),`*`) fi;
  F:=[u,op(i..nops(F),F)];
  if nargs>1 then # be conservative
    cost:=infinity; b:=NULL; st:=time();
    for i to nops(F) do
      if irred(F[i]) then next fi;
      b:=traperror(mysplit(F,i,d));
      if b=lasterror then 
        if verbosity>0 then lprint(`WARNING (mysplit): `,b) fi
        elif length(b)<cost then b0:=b; cost:=length(b) fi;
    od;
    if cost=infinity then
      vars:=[seq(op(map(proc(v,f,ii) if degree(f,v)=1 then [v,ii] fi end,
        indets(F[i]),F[i],i)),i=1..nops(F))];
      for y in vars do
        b:=traperror(elim(op(y),d,F));
        if b=lasterror then 
          if verbosity>0 then lprint(`WARNING (elim): `,b) fi
          elif length(b)<cost then b0:=b; cost:=length(b) fi;
      od;
    fi;
    if cost<infinity then 
      if verbosity>1 then lprint(` `,cost,time()-st,space()) fi; b0
      elif b=NULL then endplay(F,d) 
      else ERROR(failed)
    fi;
  else # be profligate
    for i to nops(F) do
      y:=map(proc(v,f) if degree(f,v)=1 then v fi end,indets(F),F[i]);
      if nops(y)>0 then
        y:=y[rand(1..nops(y))()]; RETURN(elim(y,i,d,F,0)) fi
    od;
    for i to nops(F) do
      if not irred(F[i]) then RETURN(mysplit(F,i,d)) fi;
    od;
    endplay(F,d,0)  # disallow ugly results
  fi
end;

# First try to completely reduce a single Z-expression cheaply, setting
# error traps for ugliness, largeness (individual summands that use >10^5
# bytes), or exceeding Maple's capacity. If it fails at any point along
# the way, restart and apply just one round of conservative reduction.
# Takes only one argument: a single Z-expression.

hedge_reduce:=proc() local i,a,a0,st;
  a:=args[1]; a0:=Z[];
  for i while a<>a0 and a<>0 do
    a0:=a; st:=time();
    if type(a,`+`) then a:=traperror(map(reduce_one,a))
      else a:=traperror(reduce_one(a)) fi;
    if a=lasterror or shape(a)[2]>100000 then 
      if verbosity>1 then lprint(` `,i,`*`,`*`,time()-st,space()) fi;
      RETURN(reduce_one(args[1],0))
    elif verbosity>1 then
      lprint(` `,i,nops(a),shape(a),time()-st,space()) fi;
  od; a
end;  

# The main procedure. Reduces any compound Z-expression to inert form,
# working hard trying to get a complete reduction.

reduce:=proc() local i,f,f0,st;
  interface(quiet=true);
  f:=args[1]; f0:=Z[];
  for i while f0<>f do
    st:=time(); f0:=f;
    if type(f,`+`) then f:=map(hedge_reduce,f) else f:=hedge_reduce(f) fi;
    if verbosity>0 then lprint(i,nops(f),shape(f),time()-st,space()) fi;
    if savefile<>NULL then save f, cat(savefile,i,`.m`) fi;
  od;
  interface(quiet=false); f
end;