#
#      commands to work with planar curves over a field
#      and divisors on them
#        (some functions only work over a finite field)
#
#
#12-11-2004
######################################################

#Read("/home/wdj/gapfiles/curves/curves.gap");

#####################################################
#
#   The data structure for a divisor on P1 
#   - a record having components:
#   * coeffs (the integer weights of the points in the support), 
#   * support (the list of points supporting the divisors), 
#   * curve
#     which has components
#     * polynomial (the function f(x,y) in F[x,y]/F[x]
#         for which f(x,y)=0 defines the curve which the divisor
#         is associated to - eg, f(x,y)=y for P^1), 
#     * ring (ie, F[x,y] - *not* as an iterated ring).
#
#   see classical_goppa_examples.gap
#
#   divisor record commands do not have type checking
#  
######################################################

MultiplicityInList:=function(L,a)
## simple list command: returns how many times a occurs in L
local mult,b;
 mult:=0;
 for b in L do
  if b=a then mult:=mult+1; fi;
 od;
 return mult;
end;

###########################################################

MatrixTransformationOnMultivariatePolynomial:=function(A,f,R)
# A is an nxn matrix with entries in a field F
# R=F[x1,...,xn]
# f is a poly in R
local m,vars,Af;
if A=() then return f; fi;
 vars:=IndeterminatesOfPolynomialRing(R);
 m:=A*vars;
 Af:=Value(f,vars,m);
return Af;
end;


##################################################

DegreeMultivariatePolynomial2:=function(f)
# f is a multivariate poly in F[x1,x2,...,xn]
# works same as one below
local partsf,monsf,varsf,deg,i,j;
 partsf:=ConstituentsPolynomial(f);
 varsf:=partsf.variables;
 monsf:=partsf.monomials;
 deg:=List([1..Length(monsf)],i->
  Sum(List([1..Length(varsf)],j->
   DegreeIndeterminate(monsf[i],varsf[j]))));
 return Maximum(deg);
end;

DegreeMultivariatePolynomial:=function(f,R)
# f is a multivariate poly in R=F[x1,x2,...,xn]
local partsf,monsf,vars,deg,i,j;
 vars:=IndeterminatesOfPolynomialRing(R);
 partsf:=ConstituentsPolynomial(f);
# varsf:=partsf.variables;
 monsf:=partsf.monomials;
 deg:=List([1..Length(monsf)],i->
  Sum(List([1..Length(vars)],j->
   DegreeIndeterminate(monsf[i],vars[j]))));
 return Maximum(deg);
end;

DegreesMultivariatePolynomial:=function(f,R)
# f is a multivariate poly in R=F[x1,x2,...,xn]
local partsf,monsf,vars,deg,i,j;
 vars:=IndeterminatesOfPolynomialRing(R);
 partsf:=ConstituentsPolynomial(f);
# varsf:=partsf.variables;
 monsf:=partsf.monomials;
 deg:=List([1..Length(monsf)],i->
  List([1..Length(vars)],j->
   [monsf[i],vars[j],DegreeIndeterminate(monsf[i],vars[j])]));
 return deg;
end;


### 
### in general, use instead  PolynomialCoefficientsPolynomial
###
CoefficientMultivariatePolynomial:=function(f,xi,power,R)
# f is a multivariate poly in R=F[x1,x2,...,xn]
# returns coeff of xi^power in f
local vars,newvars,i,f0,F;
  F:=CoefficientsRing(R); 
  f0:=f;
  vars:=IndeterminatesOfPolynomialRing(R);
  i:=Position(vars,xi);
  if Characteristic(F)=0 then
    newvars:=[];
    for i in [1..Length(vars)] do
    if vars[i]=xi then
      newvars[i]:=Zero(F);
    else
      newvars[i]:=vars[i];
    fi;
    od;
    if power=0 then 
      return Value(f,vars,newvars);
    fi;
    for i in [1..power] do
      f0:=Derivative(f0,xi);
    od;
    return Value(f0,vars,newvars)/Factorial(power);  
######## ****** works in char 0 only*******
   else
    return PolynomialCoefficientsOfPolynomial(f,i);
######## ****** works in char p as well *******
  fi;
end;


HomogenizeMultivariatePolynomial:=function(f,R1,R2)
#
# eg, f is in R1=F[x,y]
# result is homogeneous in R2=F[x,y,z]
# GAP does not allow R2=F[X,Y,Z]
#
 local Hf,vars1,vars2,numvars,HfoverZd;
 vars2:=IndeterminatesOfPolynomialRing(R2);
 vars1:=IndeterminatesOfPolynomialRing(R1);
 if not(CoefficientsRing(R1)=CoefficientsRing(R2)) then
   Error("\n\n Base fields of rings must be same.\n\n");
 fi; 
 numvars:=Length(vars1);
 HfoverZd:=Value(f,vars1,List([1..numvars],i->vars2[i]/vars2[numvars+1]));
 return NumeratorOfRationalFunction(HfoverZd);
end;

HomogenizeCoordinates:=function(Pt,R1,R2)
#
# eg, Pt is in Spec(R1), R1=F[x,y]
# result is homogeneous in Spec(R2), R2=F[x,y,z]
#
 local x,i,j,F,vars1,vars2,numvars,Hpt;
 F:=CoefficientsRing(R1);
 vars2:=IndeterminatesOfPolynomialRing(R2);
 vars1:=IndeterminatesOfPolynomialRing(R1);
 numvars:=Length(vars1);
 if numvars<>Length(Pt) then
   Error("\n\n Point must be in affine space of first ring.\n\n");
 fi; 
 Hpt:=[];
 for x in vars2 do
   i:=Position(vars2,x);
   if x in vars1 then
     j:=Position(vars1,x);
     Hpt[i]:=Pt[j];
    else
     Hpt[i]:=One(F);
   fi;
 od;
 return Hpt;
end;

DehomogenizeMultivariatePolynomial:=function(h,R2,R1)
#
# eg, h is in R2=F[x,y,z]
# result is homogeneous in R1=F[x,y]
#
 local f,vars1,vars2,numvars,vars,F;
 F:=CoefficientsRing(R1);
 vars2:=IndeterminatesOfPolynomialRing(R2);
 vars1:=IndeterminatesOfPolynomialRing(R1);
#if not(IsSubset(Set(vars2),Set(vars1))) then
#  Error("\n\n Second ring must be subring of first.\n\n");
#fi; 
 numvars:=Length(vars2);
 vars:=ShallowCopy(vars2); 
 vars[numvars]:=One(F);
 return Value(h,vars2,vars);
end;

ProjectivePoint2PointTransformation:=function(Pt1,Pt2)
# Input: Pt1, Pt2 are points in P^2
# Output: a 3x3 invertible matrix A, APt1=Pt2
local A,F,A1,A2,infty;
 F:=DefaultField(Pt1[1]);
 A:=NullMat(3,3,F);
 if (Pt1[1]=Zero(F) and Pt1[2]=Zero(F) and Pt1[3]=Zero(F)) then
    Error("\n\n First point is not in P^2\n\n");
 fi;
 if (Pt2[1]=Zero(F) and Pt2[2]=Zero(F) and Pt2[3]=Zero(F)) then
    Error("\n\n Second point is not in P^2\n\n");
 fi;
 if (Pt1[1]=Zero(F) and Pt1[2]=Zero(F) and Pt1[3]<>Zero(F)) then
    if Pt2[3]<>Zero(F) then
        A[1][1]:=One(F);
        A[2][2]:=One(F);
        A[1][3]:=Pt2[1]/Pt1[3];
        A[2][3]:=Pt2[2]/Pt1[3];
        A[3][3]:=Pt2[3]/Pt1[3];
        return A;
    fi;
    if Pt2[2]<>Zero(F) then
        A[1][2]:=One(F);
        A[3][1]:=One(F);
        A[1][3]:=Pt2[1]/Pt1[3];
        A[2][3]:=Pt2[2]/Pt1[3];
        A[3][3]:=Pt2[3]/Pt1[3];
        return A;
    fi;
    if Pt2[1]<>Zero(F) then
        A[2][2]:=One(F);
        A[3][1]:=One(F);
        A[1][3]:=Pt2[1]/Pt1[3];
        A[2][3]:=Pt2[2]/Pt1[3];
        A[3][3]:=Pt2[3]/Pt1[3];
        return A;
    fi;
 fi;
 infty:=[Zero(F),Zero(F),One(F)];
 A1:=ProjectivePoint2PointTransformation(infty,Pt1);
 A2:=ProjectivePoint2PointTransformation(infty,Pt2);
 A:=A2*A1^(-1);
 return A;
end;

##############################################################

AffineCurve:=function(poly,ring)
## crv has a polynomial component and a ring component
## this does some type checking...
local crv;
 crv:=rec();
 if IsPolynomialRing(ring) then crv.ring:=ring; fi; 
 if not(IsPolynomialRing(ring)) then 
   Error("\n 4th argument must be a polynomial ring (eg, F[x,y])\n");
 fi; 
 if poly in ring then crv.polynomial:=poly; fi; 
 if not(poly in ring) then 
   Error("\n 3rd argument must be a function in the polynomial ring (eg, y in F[x,y] for P^1)\n");
 fi; 
 return crv;
end;

OnCurve:=function(Pts,crv)
# Pts is a list of points: returns true if they are all on crv
#                       returns false otherwise
local p,f,R,F,vars,val,values;
 f:=crv.polynomial;
 R:=crv.ring;
 F:=CoefficientsRing(R);
 vars:=IndeterminatesOfPolynomialRing(R);
 values:=List(Pts,p->Value(f,vars,p));
 if f in vars then   ###   P^1 case 
   for p in Pts do
     if not(p in F) then return false; fi;
   od;
   return true;
 fi;
 for p in Pts do
   for val in values do
    if val<>Zero(F) then return false; fi;
   od;
 od;
 return true;
end;


CurveFiberFF1:=function(alpha,crv)
#
# computes the fiber over x=alpha of rational points
# on a curve over a FF. Eg, if crv is x^q+x-y^(q+1)=0 on GF(q^2):
# the fiber is the set of beta such that
# beta^q + beta = alpha^(q+1)
#
local F,f,R,vars,beta,fiber;
  fiber:=[];
  f:=crv.polynomial;
  R:=crv.ring;
  F:=CoefficientsRing(R);
  vars:=IndeterminatesOfPolynomialRing(R);
  for beta in Elements(F) do
    if Value(f,vars,[alpha,beta])=Zero(F) then
       fiber:=Concatenation(fiber,[beta]);
    fi;
  od;
  return Set(fiber);
end;

CurveFiberFF2:=function(beta,crv)
#
# computes the fiber over y=alpha of rational points
# on a curve over a FF. Eg, if crv is x^q+x-y^(q+1)=0 on GF(q^2):
# the fiber is the set of alpha such that
# beta^q + beta = alpha^(q+1)
#
local F,f,R,vars,alpha,fiber;
  fiber:=[];
  f:=crv.polynomial;
  R:=crv.ring;
  F:=CoefficientsRing(R);
  vars:=IndeterminatesOfPolynomialRing(R);
  for alpha in Elements(F) do
    if Value(f,vars,[alpha,beta])=Zero(F) then
       fiber:=Concatenation(fiber,[alpha]);
    fi;
  od;
  return Set(fiber);
end;

############### stopped here ........

PointsOnAffineCurveFF:=function(crv,extF)
## returns list of points X(E), where
## X is a curve defined over F and 
## exF/F is an extension of finite fields
 local R,F,vars,f,x0,y0,solns;
 solns:=[];
 f:=crv.polynomial;
 R:=crv.ring;
 F:=CoefficientsRing(R);
 vars:=IndeterminatesOfPolynomialRing(R);
 for x0 in Elements(extF) do
  for y0 in Elements(extF) do
    if Value(f,vars,[x0,y0])=Zero(extF) then 
     solns:=Concatenation(solns,[[x0,y0]]);
    fi;
   od;
 od;
 return solns;
end;

OnePointRiemannRochSpaceAtInfinity:=function(m,crv)
## 
## crv is a curve record representing a plane curve
##  (poly for curve must be monomial in y)
## m is an integer > 0
## Computes a basis for the RR space
## L(m*infinity) using Proposition VI.4.1 in Stichtenoth
## (see also Proposition III.10.5)
##
 local F,vars,pt,i,j,G,C,degx,degy,basisLD,f,xx,yy,R;
 R:=crv.ring;
 F:=CoefficientsRing(R);
 vars:=IndeterminatesOfPolynomialRing(R);
 xx:=vars[1];  yy:=vars[2]; 
 f:=crv.polynomial;
 degy:=DegreeIndeterminate(f,yy);
 degx:=DegreeIndeterminate(f,xx);
####### in Prop VI.4.1 degy and degx are rel prime
####### and degy is a power of p, where F=GF(q) has char p
 basisLD:=[];
 for i in [0..(degx-1)] do
  for j in [0..m] do
   if degx*j+degy*i<m+1 then
    basisLD:=Concatenation([xx^i*yy^j],basisLD);
   fi;
  od;
 od;
return Set(basisLD);
end;


DivisorOnAffineCurve:=function(cdiv,sdiv,crv)
# cdiv is the list of coeffs
# sdiv is the list of points on crv in the support
local div,F,vars,R;
 R:=crv.ring;
 F:=CoefficientsRing(R);
 vars:=IndeterminatesOfPolynomialRing(R);
 div:=rec();
 if (IsList(cdiv) and cdiv[1] in Integers) then div.coeffs:=cdiv; fi;
 if (not(IsList(cdiv)) or not(cdiv[1] in Integers)) then 
    Error("\n 1st argument is not a list of integers\n"); 
 fi;
 if ((crv.polynomial in vars) and IsList(sdiv) and sdiv[1] in F) then ### this is for P^1
   div.support:=sdiv; 
 fi;
 if (not(crv.polynomial in vars) and IsList(sdiv) and sdiv[1] in F) then ### not P^1
   div.support:=sdiv; 
 fi;
 if (not(IsList(sdiv)) or not(sdiv[1] in F)) then 
    Error("\n 2nd argument is not a list of points\n"); 
 fi;
 if Length(sdiv)<>Length(cdiv) then 
    Error("\n 1st and 2nd arguments must have same length\n"); 
 fi;
 div.curve:=crv;
 return div;
end;

DivisorAddition:=function(D1,D2)
local c1,c2,supp1,supp2,pos1,pos2,sumc,sums,pt;
 if not(D1.curve.ring=D2.curve.ring) then
      Error("\n 1st and 2nd divisor must have the same curve\n");
 fi; 
 if not(D1.curve.polynomial=D2.curve.polynomial) then
      Error("\n 1st and 2nd divisor must have the same curve\n");
 fi; 
 c1:=D1.coeffs;
 c2:=D2.coeffs;
 supp1:=D1.support;
 supp2:=D2.support;
 sumc:=[];
 sums:=[];
 for pt in Union(supp1,supp2) do
  if (pt in supp1) and (pt in supp2) then
    pos1:=PositionSublist(supp1,[pt]);
    pos2:=PositionSublist(supp2,[pt]);
    sumc:=Concatenation(sumc,[c1[pos1]+c2[pos2]]);
    sums:=Concatenation(sums,[pt]);
  fi;
  if (pt in supp1) and not(pt in supp2) then
    pos1:=PositionSublist(supp1,[pt]);
    sumc:=Concatenation(sumc,[c1[pos1]]);
    sums:=Concatenation(sums,[pt]);
  fi;
  if (pt in supp2) and not(pt in supp1) then
    pos2:=PositionSublist(supp2,[pt]);
    sumc:=Concatenation(sumc,[c2[pos2]]);
    sums:=Concatenation(sums,[pt]);
  fi;
 od;
 return rec(coeffs:=sumc,support:=sums,curve:=D1.curve);
end;


DivisorDegree:=function(div)
local c;
 c:=div.coeffs;
 return Sum(c);
end;

DivisorIsEffective:=function(div)
local c,a;
 c:=div.coeffs;
 for a in c do
   if a<0 then return false; fi;
 od;
 return true;
end;

DivisorNegate:=function(div)
local c,s;
 c:=div.coeffs;
 s:=div.support;
 return rec(coeffs:=(-1)*c,support:=s,curve:=div.curve);
end;

DivisorIsZero:=function(div)
local c,a;
 c:=div.coeffs;
 for a in c do
   if a<>0 then return false; fi;
 od;
 return true;
end;

DivisorEqual:=function(div1,div2)
local div;
 div:=DivisorAddition(div1,DivisorNegate(div2));
 return DivisorIsZero(div);
end;

DivisorGCD:=function(D1,D2)
local c1,c2,supp1,supp2,pos1,pos2,gcdcoeffs,gcdsupp,pt;
 if not(D1.curve.ring=D2.curve.ring) then
      Error("\n 1st and 2nd divisor must have the same curve\n");
 fi; 
 if not(D1.curve.polynomial=D2.curve.polynomial) then
      Error("\n 1st and 2nd divisor must have the same curve\n");
 fi; 
 c1:=D1.coeffs;
 c2:=D2.coeffs;
 supp1:=D1.support;
 supp2:=D2.support;
 gcdcoeffs:=[];
 gcdsupp:=[];
 for pt in Union(supp1,supp2) do
  if (pt in supp1) and (pt in supp2) then
    pos1:=PositionSublist(supp1,[pt]);
    pos2:=PositionSublist(supp2,[pt]);
    gcdcoeffs:=Concatenation(gcdcoeffs,[Minimum(c1[pos1],c2[pos2])]);
    gcdsupp:=Concatenation(gcdsupp,[pt]);
  fi;
  if (pt in supp1) and not(pt in supp2) then
    pos1:=PositionSublist(supp1,[pt]);
    gcdcoeffs:=Concatenation(gcdcoeffs,[Minimum(c1[pos1],0)]);
    gcdsupp:=Concatenation(gcdsupp,[pt]);
  fi;
  if (pt in supp2) and not(pt in supp1) then
    pos2:=PositionSublist(supp2,[pt]);
    gcdcoeffs:=Concatenation(gcdcoeffs,[Minimum(0,c2[pos2])]);
    gcdsupp:=Concatenation(gcdsupp,[pt]);
  fi;
 od;
 return rec(coeffs:=gcdcoeffs,support:=gcdsupp,curve:=D1.curve);
end;



###########################################################
#
#            P^1 only stuff
#
#  ..... bases of L(D) on P^1 ...
#  if D is effective then the basis is easy...
#


RiemannRochSpaceBasisFunctionsP1:=function(P,k,R2)
# input:
# R2 is a polynomial ring in x,y
# P is an element of the base field
# k is an integer
# output: 1/(x-P)^k
local x,vars;
  vars:=IndeterminatesOfPolynomialRing(R2);
  x:=vars[1];
  return 1/(x-P)^k;
end;

RiemannRochSpaceBasisEffectiveP1:=function(div)
# R is poly ring
# div is divisor record over base field
local F,n,basis,pt,cdiv,sdiv,i,j,k,pos,R;
  R:=div.curve.ring;
  F:=CoefficientsRing(R);
  if not(DivisorIsEffective(div)) then 
    Error("\n divisor must be effective \n");
  fi;
  basis:=[RiemannRochSpaceBasisFunctionsP1(Zero(F),0,R)];
  cdiv:=div.coeffs;
  sdiv:=div.support;
  n:=Length(cdiv);
  for k in [1..n] do
   for i in [1..cdiv[k]] do
    basis:=Concatenation(basis,
            [RiemannRochSpaceBasisFunctionsP1(sdiv[k],i,R)]);
   od;
  od;
  return basis;
end;


RiemannRochSpaceBasisP1:=function(div)
# R is poly ring
# div is a non-zero divisor record over 
#   base field having degree 0
local R,vars,x,div0,deg,f,F,basis,pt,cdiv,sdiv,i,j,k,pos;
  R:=div.curve.ring;
  F:=CoefficientsRing(R);
  if DivisorIsZero(div) then 
    return [One(F)];
  fi;
  deg:=DivisorDegree(div);
  if deg<0 then 
    return [Zero(F)];
  fi;
  if DivisorIsEffective(div) then 
    return RiemannRochSpaceBasisEffectiveP1(div);
  fi;
  vars:=IndeterminatesOfPolynomialRing(R);
  x:=vars[1];
  div0:=Immutable(div); ### unnecessary...
  cdiv:=div.coeffs;
  sdiv:=div.support;
  k:=Length(cdiv);  
  cdiv[k]:=cdiv[k]-deg;   ## pick the last point in div to subtract away
  f:=One(F);
  for i in [1..Length(cdiv)] do
    f:=f*(x-sdiv[i])^(-cdiv[i]);
  od;
  basis:=Concatenation([One(F)],List([0..deg],i->f*(x-sdiv[k])^(-i)));
  cdiv[k]:=cdiv[k]+deg; ## restores divisor to original
 return basis;
end;


DivisorOfRationalFunction:=function(f,R)
## R = F[x,y], f is a ratnl fcn of x
local crv,vars,y,n1,n2,suppdiv,coeffdiv,i,divf,rootsd,rootsn,den,num;
  vars:=IndeterminatesOfPolynomialRing(R); y:=vars[1];
  num:=NumeratorOfRationalFunction(f); 
  rootsn:=RootsOfUPol(num);
  den:=DenominatorOfRationalFunction(f); 
  rootsd:=RootsOfUPol(den);
  n1:=Length(Set(rootsn)); n2:=Length(Set(rootsd));
  coeffdiv:=Concatenation(List([1..n1],
    i->MultiplicityInList(rootsn, Set(rootsn)[i])),
       List([1..n2],i->-MultiplicityInList(rootsd, Set(rootsd)[i])));
  suppdiv:=Concatenation(Set(rootsn),Set(rootsd));
  crv:=AffineCurve(y,R);
  divf:=rec(coeffs:=coeffdiv,support:=suppdiv,curve:=crv);
  return divf;
end;


###############################################
#
#  Group action on RR space and associate AG code
#   for the curve P^1
#
###################################################3


MoebiusTransformation:=function(A,R)
# A is a 2x2 matrix with entries in a field F
# R=F[x]
local var,f,x,a,b,c,d,F;
 var:=IndeterminatesOfPolynomialRing(R); 
 F:=CoefficientsRing(R);
 x:=var[1];
 a:=A[1][1];
 b:=A[1][2];
 c:=A[2][1];
 d:=A[2][2];
 if c=Zero(F) and d=Zero(F) then return "infinity"; fi;
 f:=(a*x+b)/(c*x+d);
 return f;
end;


ActionMoebiusTransformationOnFunction:=function(A,f,R2)
# A is a 2x2 matrix with entries in a field F
# R=F[x,y]
# f is a rational function in F(x)
local m,numf,var,p,denf,F,R1;
if A=() then return f; fi;
 F:=CoefficientsRing(R2);
 var:=IndeterminatesOfPolynomialRing(R2);
 R1:= PolynomialRing(F,[var[1]]);
 var:=IndeterminatesOfPolynomialRing(R1);
 m:=MoebiusTransformation(A,R1);
 denf:=DenominatorOfRationalFunction(f);
 numf:=NumeratorOfRationalFunction(f);
# return Value(numf,var,[m])/Value(denf,var,[m]);
 return Value(f,var,[m]);
end;

ActionMoebiusTransformationOnDivisorP1:=function(A,div)
# A is a 2x2 matrix with entries in a field F
# div is a divisor over F
local f,sdiv,Adiv,var,p,denf,F,R,xx,R1;
if A=() then return div; fi;
 R:=div.curve.ring;
 F:=CoefficientsRing(R);
 var:=IndeterminatesOfPolynomialRing(R);
 xx:=X(F,var);
 R1:= PolynomialRing(F,[xx]);
 var:=IndeterminatesOfPolynomialRing(R1);
 Adiv:=ShallowCopy(div);
 sdiv:=div.support;
 f:=MoebiusTransformation(A,R1);
 denf:=DenominatorOfRationalFunction(f);
 for p in sdiv do
  if Value(denf,var,[p])=Zero(F) then
   Print("\n f.l.t. = ",f,",   point = ",p,"\n\n");
   Error("\n Sorry, action on this divisor is undefined\n\n");
  fi;
 od;
 Adiv.support:=List(sdiv,p->Value(f,var,[p]));
 return Adiv;
end;

ActionMoebiusTransformationOnDivisorDefinedP1:=function(A,div)
# A is a 2x2 matrix with entries in a field F
# div is a divisor over F
local f,sdiv,Adiv,var,p,denf,F,R,R1,xx;
if A=() then return div; fi;
 R:=div.curve.ring;
 F:=CoefficientsRing(R);
 var:=IndeterminatesOfPolynomialRing(R);
 xx:=X(F,var); #### this be called more than once:-)
 R1:= PolynomialRing(F,[xx]);
 var:=IndeterminatesOfPolynomialRing(R1);
 Adiv:=ShallowCopy(div);
 sdiv:=div.support;
 f:=MoebiusTransformation(A,R1);
 denf:=DenominatorOfRationalFunction(f);
 for p in sdiv do
  if Value(denf,var,[p])=Zero(F) then
   return false;
  fi;
 od;
 return true;
end;


DivisorAutomorphismGroup:=function(div)
### only works for divisors on P1 *** very slow ***
local R,F,A,autgp,sdiv,G,Adiv,eG;
 sdiv:=div.support;
 autgp:=[];
 R:=div.curve.ring;
 F:=CoefficientsRing(R);
 G:=GL(2,F); 
 eG:=Elements(G);
 for A in eG do
#  f:=MoebiusTransformation(A,R);
   if ActionMoebiusTransformationOnDivisorDefinedP1(A,div) then
      Adiv:= ActionMoebiusTransformationOnDivisorP1(A,div);
      if DivisorEqual(div,Adiv) then 
        autgp:=Concatenation(autgp,[A]);
#        eG:=Difference(eG,Elements(Group(autgp)));
#  leaving the above in slows it down!
      fi;
   fi;
 od;
 if autgp<>[] then return Group(autgp); fi;
 return Group(());
end;

MatrixRepresentationOnRiemannRochSpaceP1:=function(g,div)
#
# input: g in G subgp Aut(D) subgp Aut(X)
#       D a divisor on a curve X
# output: a dxd matrix, where d = dim L(D),
#          representing the action of g on L(D).
# Note: g sends L(D) to r*L(D), where
#       r is a polynomial of degree 1 depending on 
#       g and D
#
local i,n,R,F,f,B,gB,num,gBgood,basisLD,LD,coeffs_g,xx,R1,var;
 R:=div.curve.ring;
 var:=IndeterminatesOfPolynomialRing(R);
 F:=CoefficientsRing(R);
 B:=RiemannRochSpaceBasisP1(div);
 n:=Length(B);
 LD:=VectorSpace(F,B);
 basisLD:=Basis(LD,B);
 xx:=X(F,var);
 R1:= PolynomialRing(F,[xx]);
 gB:=List(B,f->ActionMoebiusTransformationOnFunction(g,f,R));
# this ring R for gB must be same ring as for B
 if not(DivisorIsEffective(div)) then
  coeffs_g:=[];
  for i in [1..n] do
    coeffs_g[i]:=Coefficients( basisLD, gB[i] ); 
  od;
 fi;
 if DivisorIsEffective(div) then  # Coefficients can't handle a 
                                  # constant function so doing this by hand
  coeffs_g:=[List(basisLD,f->Zero(F))];
  coeffs_g[1][1]:=One(F);
  for i in [2..n] do
    coeffs_g[i]:=Coefficients( basisLD, gB[i] ); 
  od;
 fi;
 return coeffs_g;
end;