#
# Count points on varieties over finite fields the old-fashioned way
# (i.e., one at a time). More generally, evaluate compound Z-expressions
# at primepowers q. Use interpolation to fit the counts to polynomials,
# where possible.
#
# Should work with Maple V.3,4,5. Will not work with V.2 or earlier.
#
# 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.
#
# Examples: q^2*Z[x^2-y^2];   Z[] ( = 1);
#              Z[x^2+y^2,x-y] - Z[2,x]
#
############################################
#
# Main procedures defined in this file:
#
# evalZ    evaluate a compound Z-expression over F_q, q = any primepower.
# format   given a compound Z-expression that is alleged to be a
#          polynomial in 1/q, return a polynomial in 1/q with undetermined
#          coefficients that the Z-expression must fit.
# fitpoly  given a compound Z-expression u and a polynomial g in 1/q with
#          k undetermined coefficients that u must fit, evaluate u in F_q
#          for q = the smallest k primepowers, and find the polynomial in
#          the form of g that fits this data.
#
# Note: Z, q, t, Gf, and a1,a2,... are used as global variables.
#
# Additional procedures:
#
# ctz1        count the # of zeroes of a polynomial list over F_p, p=prime.
# ctz         count the # of zeroes of a poly. list over F_q, q=primepower.
# buildGf     build a table of data for the Galois field F_q.
# is_homog    determine whether a polynomial list is homogeneous.
# is_prim     determine whether a polynomial list is primitive.
# has_scalar  determine whether a scalar is a member of a poly. list.
#
##########################################
#
# Examples:
#
# u:=Z[x^2+x*y+y^2+x*z];
# evalZ(u,3);
# evalZ(u,2,2);
# g:=format(u);
# f:=fitpoly(u,g);
#
########################################################################
#
# 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;

# Given a list F of polynomials with integer coefficients, count
# the number of simultaneous zeroes of F over F_p, p prime.
# Fancier approaches are possible, but this is good enough if most
# of the time F is a single polynomial whose coeffs are mostly {1,-1}.

ctz1:=proc(F,p) local vars,n,ct,A,i,b;
  vars:=indets(F); n:=nops(vars);
  ct:=0; A:=table([0$n]);
  do
    b:=subs({seq(vars[i]=A[i],i=1..n)},F);
    if modp({0,op(b)},p)={0} then ct:=ct+1 fi;
    for i to n while A[i]=p-1 do A[i]:=0 od;
    if i>n then break else A[i]:=A[i]+1 fi;
  od; ct
end;

# Given a prime p and integer k>0, build a table of data about the
# field F_q, q=p^k, and return an irreducible polynomial that defines
# F_q as an extension of F_p. (Gf="Galois field").

buildGf:=proc(p,k) global Gf; local f,g,i,j,A;
  option remember;
  A:=table([0$k]);
  do
    f:=convert([seq(A[i+1]*t^i,i=0..k-1),t^k],`+`);
    if Irreduc(f) mod p then break fi;
    for i to k while A[i]=p-1 do A[i]:=0 od;
    if i>k then ERROR(`check your input`) else A[i]:=A[i]+1 fi;
  od;
  A:='A'; g:=convert([seq(A[i+1]*t^i,i=0..k-1)],`+`);
  A:=table([0$k]);
  for i from 0 do
    Gf[p,k][i]:=eval(g);
    for j to k while A[j]=p-1 do A[j]:=0 od;
    if j>k then break else A[j]:=A[j]+1 fi;
  od; f
end;

# Given a list F of polynomials with integer coefficients, count the
# number of simultaneous zeroes of F over F_q, q=p^k, p prime.

ctz:=proc(F,p,k) local f0,nz,vars,n,ct,A,f,i,b;
  f0:=buildGf(p,k); nz:=p^k-1;
  vars:=indets(F); n:=nops(vars);
  ct:=0; A:=table([0$n]);
  do ct:=ct+1;
    for f in F do
      b:=subs({seq(vars[i]=Gf[p,k][A[i]],i=1..n)},f);
      if modp(rem(b,f0,t),p)<>0 then ct:=ct-1; break fi
    od;
    for i to n while A[i]=nz do A[i]:=0 od;
    if i>n then break else A[i]:=A[i]+1 fi;
  od; ct
end;

# Given a compound Z-expression u,
#   evalZ(u,p)   evaluates u in F_p, p prime
#   evalZ(u,p,k) evaluates u in F_q, q=p^k, p prime, k>0.

evalZ:=proc(u,p) local alg,k,vars,v,g;
  interface(quiet=true);
  if nargs=2 or args[3]=1 then alg:=ctz1; k:=1
    else alg:=ctz; k:=args[3] fi;
  vars:=indets(u,indexed);
  g:=subs({seq(v=v/q^nops(indets([op(v)])),v=vars)},u);
  g:=subs({q=p^k,seq(v=alg([op(v)],p,k),v=vars)},g);
  interface(quiet=false); g
end;

# Return true if the polynomial list F is homogeneous.
# Could fail if the terms are not collected or expanded.

is_homog:=proc(F) local f;
  for f in F do
    if ldegree(f)<>degree(f) then RETURN(false) fi
  od; true
end;

# return true if the polynomial list is primitive

is_prim:=proc(F) evalb(igcd(op(map(content,F)))=1) end;

# Return true if there is a scalar appearing among the members of a
# polynomial list. Could fail if the terms are not expanded/collected.

has_scalar:=proc(F) local f;
  for f in F do if type(f,'integer') then RETURN(true) fi od; false
end;

# Given a compund Z-expression u that is alleged to be a polynomial
# in 1/q, return a polynomial in 1/q with undetermined coefficients
# that the Z-expression must fit, using
# 1. The scalar trick: 
#      If f_i is a nonzero scalar then Z[f_1,...,f_k] can be deleted
#      from u, since f_i<>0 in sufficiently large characteristics.
# 2. The homogeneous trick:
#      If every Z-expression (after #1) has homogeneous components,
#      then the value of u at q=1 can be obtained by substituting
#      Z[f_1,...]=1 for all constituent Z-expressions (and q=1) in u.
# 3. The top/bottom trick:
#      The top of c * q^-d * Z[f_1,...,f_k] is d+1 if [f_1,...,f_k] is
#      primitive. Otherwise, it is d. The bottom is d+m, where m is the
#      number of variables that appear. For a compound Z-expression,
#      the top and bottom are the min and max of the tops and bottoms
#      of the constituents. The only undetermined coefficients we need
#      to worry about are those of q^-i, i=top..bottom. 
# Uses a1,a2,... as names for the undetermined coefficients.

format:=proc(u) local b,vars,homog,low,high,v,F,c,m,d,res,i,eq;
  b:=subs(Z[]=1,u); vars:=indets(b,indexed);
  homog:=true; low:=NULL; high:=NULL;
  for v in vars do
    F:=subs(0=NULL,map(expand,[op(v)]));
    if has_scalar(F) then b:=subs(v=0,b); next fi;
    if homog and not is_homog(F) then homog:=false fi;
    c:=coeff(b,v); m:=nops(indets(F)); d:=-degree(c);
    if is_prim(F) then d:=d+1 fi;
    low:=min(low,d); high:=max(high,-ldegree(c)+m);
  od;
  res:=subs({seq(v=0,v=vars)},b);
  if low=NULL then RETURN(res) fi;
  if res=0 then d:=high else d:=max(high,-ldegree(res,q)) fi;
  res:=[seq(coeff(res,q,-i),i=0..d)];
  res:=subsop(seq(i+1=cat('a',i),i=low..high),res);
  res:=convert([seq(res[i+1]/q^i,i=0..d)],`+`);
  if homog then
    eq:=subs({q=1,seq(v=1,v=vars)},b-res);
    res:=subs(solve(eq,{cat('a',high)}),res)
  fi; res;
end;

# Given a compound Z-expression u and a polynomial g in 1/q with k
# unknown coefficients, evaluate u for q = the k smallest primepowers,
# and find the unique poly in the form of g that fits this data.
#
# Intermediate results are lprinted along the way in the format:
#  q   value_at_q   cpu_seconds_used_to_find_value   bytes_allocated

fitpoly:=proc(u,g) local eqns,k,i,p,st,v,primepowers;
  primepowers:=[[2,1],[3,1],[2,2],[5,1],[7,1],[2,3],[3,2],
    [11,1],[13,1],[2,4],[17,1],[19,1],[23,1],[5,2],[3,3]];
  eqns:=NULL; k:=nops(indets(g) minus {q});
  if k>nops(primepowers) then ERROR(`too big`) fi;
  for i to k do
    p:=primepowers[i]; st:=time():
    v:=evalZ(u,op(p));
    lprint(p[1]^p[2],v,time()-st,space());
    eqns:=eqns,subs(q=p[1]^p[2],g)=v
  od; lprint();
  subs(solve({eqns}),g)
end;