"""
Affine plane curves over a general ring. Basis computation of 
Riemann-Roch spaces of hyperelliptic curves over a finite field
for effective divisors.

AUTHORS:
    -- 2005-11-13, William Stein <wstein@ucsd.edu>
    -- 2005-11-13, David Joyner
    -- 2006-01 David Kohel <kohel@maths.usyd.edu.au>
    -- 2006-04 David Joyner (added local_coordinates_ram, fixed divisor_of_function
                             added is_ramification_point1, is_ramification_point2, 
                             riemann_roch_space, many minor bug fixes)
    -- 2007-10 David Joyner - minor modifications to make code run with sage 2.8

"""

#**********************************************************************************
#       Copyright (C) 2005 David Joyner, David Kohel, William Stein <wstein@gmail.com>
#
#  Distributed under the terms of the GNU General Public License (GPL)
#
#  The full text of the GPL is available at:
#
#                  http://www.gnu.org/licenses/
#**********************************************************************************

from sage.interfaces.all import singular
from sage.misc.sage_eval import *
from sage.misc.all import latex, generic_cmp, add

from sage.rings.all import (PolynomialRing,
                            degree_lowest_rational_function, is_PrimeField)

from sage.schemes.generic.affine_space import is_AffineSpace
from sage.schemes.generic.algebraic_scheme import AlgebraicScheme_subscheme_affine

from sage.schemes.generic.divisor import *

from sage.schemes.plane_curves.curve import Curve_generic

class AffineSpaceCurve_generic(Curve_generic, AlgebraicScheme_subscheme_affine):
    def _repr_type(self):
        return "Affine Space"
    
    def __init__(self, A, X):
        if not is_AffineSpace(A):
            raise TypeError, "A (=%s) must be an affine space"%A
        Curve_generic.__init__(self, A, X)
        d = self.dimension()
        if d != 1:
            raise ValueError, "defining equations (=%s) define a scheme of dimension %s != 1"%(X,d)

class AffineCurve_generic(Curve_generic):
    def __init__(self, A, f):
        P = f.parent()
        if not (is_AffineSpace(A) and A.dimension != 2):
            raise TypeError, "Argument A (= %s) must be an affine plane."%A
        Curve_generic.__init__(self, A, [f])

    def _repr_type(self):
        return "Affine"

    def is_ramification_point1(self,pt):
        """
        Returns True if pt is a ramification point of self under the 
        map C --> A^1 given by (x,y) -> x.

        EXAMPLES:
            sage: F = GF(5)
            sage: A2 = AffineSpace(2, F, names = 'xy')
            sage: R = A2.coordinate_ring()
            sage: x, y = R.gens()
            sage: f = y^2 - x^9 - x
            sage: C = AffineCurve_GF(A2,f)
            sage: pt = C([0,0])
            sage: C.is_ramification_point(pt)
            True
            sage: pt = C([2,3])
            sage: C.is_ramification_point(pt)
            False

        """
        F = self.base_ring()
        f = self.defining_polynomial()
        R = f.parent()
        x,y = R.gens()
        R0 = PolynomialRing(F,4,names = [str(x),str(y),"a","t"])
        x,y,a,t = R0.gens()
        f0 = f(x,y)*t**0*a**0
        ft = f0(pt[0]+t,pt[1]+a*t,F(1),F(1))
        c = ft.coefficient(t)
        d = c.degree(a)
        if d>0:
            return False
        return True

    is_ramification_point = is_ramification_point1

    def is_ramification_point2(self,pt):
        """
        Returns True if pt is a ramification point of self under the 
        map C --> A^1 given by (x,y) -> y.

        EXAMPLES:
            sage: F = GF(5)
            sage: A2 = AffineSpace(2, F, names = 'xy')
            sage: R = A2.coordinate_ring()
            sage: x, y = R.gens()
            sage: f = y^2 - x^9 - x
            sage: C = AffineCurve_GF(A2,f)
            sage: pt = C([0,0])
            sage: C.is_ramification_point2(pt)
            False

        """
        F = self.base_ring()
        f = self.defining_polynomial()
        R = f.parent()
        x,y = R.gens()
        R0 = PolynomialRing(F,4,names = [str(x),str(y),"a","t"])
        x,y,a,t = R0.gens()
        f0 = f(x,y)*t**0*a**0
        ft = f0(pt[0]+a*t,pt[1]+t,F(1),F(1))
        c = ft.coefficient(t)
        d = c.degree(a)
        if d>0:
            return False
        return True

    def ramification_degree1(self,pt):
        """
        Returns the ramification degree d if pt is a ramification point of self under the 
        map C --> A^1 given by (x,y) -> x.

        NOTE: We assume the ramification degree is < 10000.

        EXAMPLES:
            sage: F = GF(5)
            sage: A2 = AffineSpace(2, F, names = 'xy')
            sage: R = A2.coordinate_ring()
            sage: x, y = R.gens()
            sage: f = y^2 - x^9 - x
            sage: C = AffineCurve_GF(A2,f)
            sage: pt = C([0,0])
            sage: C.is_ramification_point(pt)
            True
            self: C.ramification_degree(pt)
            2

        """
        if not(self.is_ramification_point1(pt)):
            return 0 
        F = self.base_ring()
        f = self.defining_polynomial()
        R = f.parent()
        x,y = R.gens()
        R0 = PolynomialRing(F,4,names = [str(x),str(y),"a","t"])
        x,y,a,t = R0.gens()
        f0 = f(x,y)*t**0*a**0
        for d in range(1,10000):
            ft = f0(pt[0]+t**d,pt[1]+a*t,F(1),F(1))
            c = ft.coefficient(t**d)
            if c.degree(a)>0:
                return d
        return "fail"

    ramification_degree = ramification_degree1

    def ramification_degree2(self,pt):
        """
        Returns the ramification degree d if pt is a ramification point of self under the 
        map C --> A^1 given by (x,y) -> y.

        NOTE: We assume the ramification degree is < 10000.

        EXAMPLES:
            sage: F = GF(5)
            sage: A2 = AffineSpace(2, F, names = 'xy')
            sage: R = A2.coordinate_ring()
            sage: x, y = R.gens()
            sage: f = y^2 - x^9 - x
            sage: C = AffineCurve_GF(A2,f)
            sage: pt = C([0,0])
            sage: C.is_ramification_point2(pt)
            False
            self: C.ramification_degree2(pt)
            0

        """
        if not(self.is_ramification_point2(pt)):
            return 0 
        F = self.base_ring()
        f = self.defining_polynomial()
        R = f.parent()
        x,y = R.gens()
        R0 = PolynomialRing(F,4,names = [str(x),str(y),"a","t"])
        x,y,a,t = R0.gens()
        f0 = f(x,y)*t**0*a**0
        for d in range(1,10000):
            ft = f0(pt[0]+a*t,pt[1]+t**d,F(1),F(1))
            c = ft.coefficient(t**d)
            if c.degree(a)>0:
                return d
        return "fail"


    def divisor_of_function(self, r):
        """
        Return the divisor of a function on a curve. 
        
        WARNING: The zeros and poles of r must be degree 1 (as places), smooth, and 
        unramified w.r.t. (x,y) -> x.    
        
        INPUT:
             r is a rational function on X
             
        OUTPUT:
             list -- The divisor of r represented as a list of coefficients and points.
                     (TODO: This will change to a more structural output in the future.)

        EXAMPLES:
            sage: F = GF(5)
            sage: A2 = AffineSpace(2, F, names = 'xy')
            sage: R = A2.coordinate_ring()
            sage: x, y = R.gens()
            sage: f = y^2 - x^9 - x
            sage: C = AffineCurve_GF(A2,f)
            sage: K = FractionField(R)
            sage: r = 1/x
            sage: C.divisor_of_function(r)   
            -(y, x)
            sage: r = 1/x^3
            sage: C.divisor_of_function(r)
            -3*(y, x)
        """
        F = self.base_ring()
        f = self.defining_polynomial()
        pts = self.rational_points() 
        numpts = len(pts)
        R = f.parent()
        x,y = R.gens()
        R0 = PolynomialRing(F,3,names = [str(x),str(y),"t"])
        vars0 = R0.gens()
        t = vars0[2]
        divf = []
        for pt0 in pts:
            lcs = self.local_coordinates(pt0,5)
            yt = lcs[1]
            xt = lcs[0]
            ldg = degree_lowest_rational_function(r(xt,yt),t)
            #print pt0,lcs,ldg
            if ldg[0] != 0:
                divf.append((ldg[0],self(pt0)))
        #print divf
        return self.divisor(divf)

    def local_coordinates_unram(self, pt, n):
        r"""
        Return local coordinates to precision n at the given point.

        \begin{note}
        {\bf Behaviour is flakey} - some choices of $n$ are worst that others.
        Can be very slow if n and the degree of f are large.
        \end{note}

        INPUT:
            pt -- a smooth F-rational point on X which is not a
                  point of ramification for the projection (x,y) -> x.
            n  -- the number of terms desired

        OUTPUT:
            x = x0 + t
            y = y0 + power series in t

        EXAMPLES:
            sage: F = GF(5)
            sage: A2 = AffineSpace(2, F, names='xy')
            sage: x,y = A2.gens()
            sage: f = y^2 - x^9 - x
            sage: C = AffineCurve_GF(A2,f)
            sage: pt = C([2,3])
            sage: C.local_coordinates_unram(pt, 9)
            [2 + t, 3 + 3*t^2 + t^3 + 3*t^4 + 3*t^6 + 3*t^7 + t^8 + 2*t^9 + 3*t^11 + 3*t^12]
        """
        S = singular
        f = self.defining_polynomial()
        R = f.parent()
        F = self.base_ring()
        p = F.characteristic()
        x0 = F(pt[0])
        y0 = F(pt[1])
        astr = ["a"+str(i) for i in range(1,2*n)]
        x,y = R.gens()
        R0 = PolynomialRing(F,2*n+2,names = [str(x),str(y),"t"]+astr)
        Rs = singular(R0)
        vars0 = R0.gens()
        t = vars0[2]
        yt = y0*t**0+add([vars0[i]*t**(i-2) for i in range(3,2*n+2)])
        xt = x0+t
        xt = S(xt) 
        yt = S(yt) 
        fs = str(f).replace("^","**")    ## William Stein's suggestion:
        fs = fs.replace("x","xt")        ## do multivariate polynomial
        fs = fs.replace("y","yt")        ## evaluation in singular
        ft = eval(fs)
        fstr = str(ft)
        S.eval('ring s = '+str(p)+','+str(R0.gens())+',lp;')
        S.eval('poly f = '+fstr + ';')
        c = S('coeffs(%s, t)'%fstr)
        N = int(c.size())
        b = ["%s[%s,1],"%(c.name(), i) for i in range(2,N/2-4)]
        b = ''.join(b)
        b = b[:len(b)-1] # to cut off the trailing comma
        cmd = 'ideal I = '+b
        S.eval(cmd)
        #S.eval('short = 0')    # print using *'s and ^'s.
        c = S.eval('slimgb(I)')
        d = c.split("=")
        d = d[1:]
        d[len(d)-1] += "\n"
        e = [x[:x.index("\n")] for x in d]
        vals = []
        for x in e:
            for y in vars0:
                if str(y) in x:
                    if len(x.replace(str(y),"")) != 0:
                        i = x.find("-")
                        if i>0:
                            vals.append([eval(x[1:i]),x[:i],F(eval(x[i+1:]))])
                        i = x.find("+")
                        if i>0:
                            vals.append([eval(x[1:i]),x[:i],-F(eval(x[i+1:]))])
                    else:
                        vals.append([eval(str(y)[1:]),str(y),F(0)])
        vals.sort()
        k = len(vals)
        v = [x0+t,y0+add([vals[i][2]*t**(i+1) for i in range(k)])]
        return v

    local_coordinates = local_coordinates_unram

    def local_coordinates_ram(self, pt, n):
        r"""
        Return local coordinates to precision n at the given point.
        
        \begin{note}
        {\bf Behaviour is flakey} - some choices of $n$ are worst that others.
        Can be very slow if n and the degree of f are large.
        \end{note}

        INPUT:
            pt -- a smooth F-rational point on X which is a
                  point of ramification for the projection (x,y) -> x.
            n  -- the number of terms desired

        OUTPUT:
            x = x0 + t
            y = y0 + power series in t

        EXAMPLES:
            sage: F = GF(5)
            sage: pt = (0,0)
            sage: A2 = AffineSpace(2, F, names='xy')
            sage: x,y = A2.gens()
            sage: f = y^2 - x^9 - x
            sage: C = AffineCurve_prime_GF(A2,f)
            sage: C.local_coordinates_ram(pt, 4)
            [t^2, t]

        #sage: time C.local_coordinates_ram(pt, 10)             # very time-consuming
        #CPU times: user 18.82 s, sys: 0.37 s, total: 19.19 s   # random output
        #Wall time: 4449.89                                     # random output
        #[t^2 + 4*t^18 + 4*t^19, t]
        """
        S = singular
        f = self.defining_polynomial()
        R = f.parent()
        F = self.base_ring()
        p = F.characteristic()
        x0 = F(pt[0])
        y0 = F(pt[1])
        astr = ["a"+str(i) for i in range(1,2*n)]
        x,y = R.gens()
        R0 = PolynomialRing(F,2*n+2,names = [str(x),str(y),"t"]+astr)
        Rs = singular(R0)
        vars0 = R0.gens()
        t = vars0[2]
        xt = x0*t**0+add([vars0[i]*t**(i-2) for i in range(3,2*n+2)])
        yt = y0+t
        xt = S(xt) 
        yt = S(yt) 
        fs = str(f).replace("^","**")    ## William Stein's suggestion:
        fs = fs.replace("x","xt")        ## do multivariate polynomial
        fs = fs.replace("y","yt")        ## evaluation in singular
        ft = eval(fs)                    ## (much fast than in python) 
        fstr = str(ft)
        S.eval('ring s = '+str(p)+','+str(R0.gens())+',lp;')
        S.eval('poly f = '+fstr + ';')
        S.eval('matrix C =coeffs(%s, t)'%fstr)
        c = S('coeffs(%s, t)'%fstr)
        N = int(c.size())
        cfs = [S("C[1][%s]"%i) for i in range(1,2*n)]
        b = [str(a)+"," for a in cfs]
        b = ''.join(b)
        b = b[:len(b)-1] # to cut off the trailing comma
        cmd = 'ideal I = '+b
        S.eval(cmd)
        S.eval('short=0')    # print using *'s and ^'s.
        c = S.eval('slimgb(I)')
        d = c.split("=")
        d = d[1:]
        d[len(d)-1] += "\n"
        e = [x[:x.index("\n")] for x in d]
        vals = []
        for x in e:
            for z in vars0:
                if str(z) in x:
                    if len(x.replace(str(z),"")) != 0:
                        i = x.find("-")
                        if i>0:
                            vals.append([eval(x[1:i]),x[:i],F(eval(x[i+1:]))])
                        i = x.find("+")
                        if i>0:
                            vals.append([eval(x[1:i]),x[:i],-F(eval(x[i+1:]))])
                    else:
                        vals.append([eval(str(z)[1:]),str(z),F(0)])
        vals.sort()
        k = len(vals)
        v = [x0+add([vals[i][2]*t**(i+1) for i in range(k)]),y0+t]
        return v

class AffineCurve_GF(AffineCurve_generic):
    def rational_points(self, algorithm="enum"):
        r"""
        Return sorted list of all rational points on this curve.

        Use \emph{very} naive point enumeration to find all rational points
        on this curve over a finite field.
        
        EXAMPLE:
            sage: A, (x,y) = AffineSpace(2,GF(9)).objgens()
            sage: C = AffineCurve_GF(A,x^2 + y^2 - 1)
            sage: C
            Affine Curve over Finite Field in a of size 3^2 defined by 2 + x1^2 + x0^2
            sage: C.rational_points()
            [(2*a + 2, 2*a + 2), (2*a + 2, a + 1), (a + 1, 2*a + 2),
             (a + 1, a + 1), (2, 0), (1, 0), (0, 2), (0, 1)]
        """
        f = self.defining_polynomial()
        R = f.parent()
        K = R.base_ring()
        points = []
        for x in K:
            for y in K:
                if f(x,y) == 0:
                    points.append(self((x,y)))
        points.sort()
        return points


class AffineCurve_prime_GF(AffineCurve_GF):
    # CHECK WHAT ASSUMPTIONS ARE MADE REGARDING AFFINE VS. PROJECTIVE MODELS!!!
    # THIS IS VERY DIRTY STILL -- NO DATASTRUCTURES FOR DIVISORS. 
    
    def rational_points(self, algorithm="enum"):
        r"""
        Return sorted list of all rational points on this curve.
        
        INPUT:
            algorithm -- string: 
                           'enum' -- straightforward enumeration
                           'bn' -- via Singular's Brill-Noether package.
                           'all' -- use all implemented algorithms
                                    and verify that they give the same answer,
                                    then return it

        \note{The Brill-Noether package does not always work.  When it
        fails a RuntimeError exception is raised.}
         
        EXAMPLE:
            sage: A2 = AffineSpace(2, F, names='xy')
            sage: x,y = A2.gens()
            sage: f = y^2 - x^9 - x
            sage: C = AffineCurve_GF(A2,f); C
            Affine Curve over Finite Field of size 5 defined by y^2 + 4*x + 4*x^9
            sage: C.rational_points(algorithm='bn')
            [(0, 0), (2, 2), (2, 3), (3, 1), (3, 4)]
            sage: C = AffineCurve_GF(A2,x - y + 1)
            sage: C.rational_points()
            [(0, 1), (1, 2), (2, 3), (3, 4), (4, 0)]

        Now on a different curve:

            sage: x, y = (GF(17)['x,y']).gens()
            sage: C = Curve(x^2 + y^5 + x*y - 19)
            sage: v = C.rational_points(algorithm='bn')
            sage: w = C.rational_points(algorithm='enum')
            sage: len(v)
            20
            sage: v == w
            True
            
        """
        if algorithm == "enum":

            return AffineCurve_finite_field.rational_points(self, algorithm="enum")

        elif algorithm == "bn":
            f = self.defining_polynomial()._singular_()
            singular = f.parent()
            singular.lib('brnoeth')
            try:
                X1 = f.Adj_div()
            except (TypeError, RuntimeError), s:
                raise RuntimeError, str(s) + "\n\n ** Unable to use the Brill-Noether Singular package to compute all points (see above)."
                
            X2 = singular.NSplaces(1, X1)
            X3 = singular.extcurve(1, X2)
            R = X3[1][5]
            singular.set_ring(R)

            # We use sage_flattened_str_list since iterating through
            # the entire list through the sage/singular interface directly
            # would involve hundreds of calls to singular, and timing issues with
            # the expect interface could crop up.  Also, this is vastly
            # faster (and more robust).
            v = singular('POINTS').sage_flattened_str_list()
            pnts = [self(int(v[3*i]), int(v[3*i+1])) for i in range(len(v)/3) if int(v[3*i+2])!=0]
            pnts.sort()      ##### Singular's ordering for the rational pts is random 
            return pnts      

        elif algorithm == "all":

            S_enum = self.rational_points(algorithm = "enum")
            S_bn = self.rational_points(algorithm = "bn")
            if S_enum != S_bn:
                raise RuntimeError, "Bug in rational_points -- different algorithms give different answers for curve %s!"%self
            return S_enum

        else:
            raise ValueError, "No algorithm '%s' known"%algorithm


#from sage.schemes.hyperelliptic_curves.hyperelliptic_generic import *
from sage.schemes.plane_curves.projective_curve import ProjectiveCurve_generic
from sage.rings.all import (factor, FractionField, PolynomialRing, degree_lowest_rational_function)
from sage.misc.all import add, mul
from sage.schemes.generic.divisor import *
from sage.schemes.plane_curves.constructor import *
from sage.sets.set import *
import copy

def riemann_roch_space(C, D, E):
        """
        INPUT:
            D, E -- effective divisors on a hyperelliptic curve C
                    defined over a finite field, D<E
                    The support of E must 
                     (a) consist only of affine, unramified points
                     (b) have distinct x-coordinates (ie, 
                         (x0,y0) and (x0,-y0) cannot both 
                         belong to supp(E)).

        OUTPUT:
           A sub-basis of L(E) no element of which is
           contained in $L(D)$.  Moreover, every element of
           this basis is of the form a(x)+y*b(x), for some
           rational functions a(x), b(x).

        Using \code{self.riemann_roch_space(D,E)} and
        \code{self.riemann_roch_space([0,...,0],D)}, one can compute
        bases $B$ of $L(D)$ and $C$ of $L(E)$ such that $B\subset C$.

        This is an implementation of the pivot method, described in
        Joyner-Ksir, {\em Bases for Riemann-Roch spaces of hyperelliptic
        curves: Effective case}. This method is basically a version
        of Coates' algorithm in this special case.

        Let $X$ denote a hyperelliptic curve over an algebraically
        closed field $F$ and let $D<E$ denote effective divisors on
        $X$.  See the docstring of the function \code{rr_space1_3} for
        more details of the algorithm.
        
        EXAMPLES:
            sage: F = GF(17)
            sage: A2 = AffineSpace(2, F, names='xy')
	    sage: x,y = A2.gens()
            sage: f = x**9 + x
            sage: C = AffineCurve_GF(A2,y**2-f)
            sage: pts = C.rational_points()
            sage: D = C.divisor([(0,pts[1]),(0,pts[6])])
            sage: E = C.divisor([(3,pts[1]),(4,pts[6])])
            sage: b = riemann_roch_space(C,D,E); b
            [(16 + y + 14*x + 8*x^2 + 6*x^3 + x^4 + 16*x^5 + x^6)/(16 + 2*x + 6*x^2 + 15*x^3 + 15*x^4 + 13*x^6 + x^7),
             (12 + y + 6*x + 3*x^4 + 7*x^5)/(4 + 8*x + 8*x^2 + 6*x^3 + 15*x^4 + 9*x^5 + x^6),
             (12 + 16*y + 2*x + 12*x^2 + 16*x^3 + 15*x^4)/(16 + 11*x + 8*x^2 + 5*x^3 + 12*x^4 + 16*x^5)]
            sage: C.divisor_of_function(b[0])
            (y, 12 + x) + (15 + x, 2 + y) - 4*(3 + y, 4 + x) - 3*(16 + x, 6 + y) + (1 + x, 7 + y)
            sage: C.divisor_of_function(b[1])
            -3*(3 + y, 4 + x) - 3*(16 + x, 6 + y)
            sage: C.divisor_of_function(b[2])
            -2*(3 + y, 4 + x) - 3*(16 + x, 6 + y)

        REFERENCES:

        """
        Esupp = E.support()
        Ecoeffs = [E.coeff(p) for p in Esupp]
        Dsupp = D.support()
        Dcoeffs = [D.coeff(p) for p in Dsupp]
        E0 = [[E.coeff(p),p] for p in Esupp]
        D0 = [[D.coeff(p),p] for p in Dsupp]
        f = C.defining_polynomial()
        R = f.parent()
        if f.coefficient(R.gen(1)) != 0:
            raise NotImplementedError, "The defining equation (=%s) must be of the form y^2 - h(x) = 0"%f
        F = C.base_ring()
        p = F.characteristic()
        x = R.gens()[0]
        deg = f.degree()
        genus = C.geometric_genus()
        dim = sum(Ecoeffs)+1-genus        ##  by the RR thrm *if* div(E)>= 2g-1
        suppDE = []
        suppDEidx = []
        pts = C.rational_points() 
        for i in range(len(E0)):
            if Ecoeffs[i]>Dcoeffs[i]:
                for j in range(Ecoeffs[i]-Dcoeffs[i]):
                    suppDE.append([pts[i],i])
                suppDEidx.append(i)
        basis = []
        Enew = E
        while (dim>1 and len(suppDE)>0):
            b = rr_space1_3(f,Esupp,Ecoeffs,Dsupp,Dcoeffs)
            basis.append(b[0])
            PP = suppDE.pop()
            P = PP[0]
            idxP = PP[1]
            Ecoeffs[idxP] = Ecoeffs[idxP]-1
            Enew = Enew - C.divisor([(1,P)])
            dim = sum(Ecoeffs)+1-genus
        return basis


def rr_space1_3(f,Esupp,Ecoeffs,Dsupp,Dcoeffs):
    r"""
    Main subprocedure of riemann_roch_space.
    
    The hyperelliptic curve $C$ is defined by $f(x,y) = y^2-h(x) = 0$,
    h(x) is a polynomial over a finite field F.  $D \leq E$ are
    effective divisors supported on unramified F-rational (degree 1) places on
    $X$. The divisors $D$ and $E$ are represented by

    Let $E = \sum_i e_i P_i, D = \sum_i d_i P_i$
        $L = \{ i \ |\ e_i > d_i \} = \{ i_0 < i_1 < ... < i_k \}$

    \begin{enumerate}
    \item Pick $j = i_0$, pivot pt = $P_j^*$, 
       seed fcn = $f_0(x,y) = (y+y_j)/\prod_i (x-x_i)$,
       pivot functions = $f_{j,a}(x,y) = (x - x_j)^{-a}$

    \item Compute local coords x = x(t), y = y(t) about pivot point 

    \item Compute $c_1,...,c_{e_j}$ such that
          $f^1_0(x,y)= f_0(x,y)-c_{e_j}f_{j,e_j}(x,y)-...-c_1f_{j,1}(x,y)$
       has not pole at the pivot pt $P_j^*$
    
    \item Replace j by next element of L, replace $f_0(x,y)$ by 
       $f^1_0(x,y)$, and redefine the pivot fcns.
       Repeat steps 2, 3. 

    \item Repeat step 4 until elements of L are exhausted and all poles at
       $P_{i_0}^*, ..., P_{i_k}^*$ are  removed. Add the resulting fcn 
       $f^k_0(x,y)$ to the basis.

    \item Replace E by E-P, where P in supp(E)-supp(D) is arbitrary. 

    \item Repeat 1-6 until E = D.
    \end{enumerate}

    Steps 1--3 are implemented in this function. The remaining steps
    in the riemann_roch_space function above.
    
    EXAMPLES:
        sage: F = GF(5)
        sage: A2 = AffineSpace(2, F, names='xy')
        sage: x,y = A2.gens()
        sage: f = x**9 + x
        sage: C = AffineCurve_GF(A2,y^2-f)
        sage: pts = C.rational_points()
        sage: D = C.divisor([(0,pts[1]),(0,pts[3])])
        sage: E = C.divisor([(4,pts[1]),(2,pts[3])])
        sage: b = riemann_roch_space(C,D,E); b
        [(y + 3*x + 2*x^2 + x^3 + x^4 + 2*x^5)/(4 + x + 4*x^2 + 2*x^3 + x^4 + x^5 + x^6),
         (4 + 4*y + x + 2*x^4)/(3 + 3*x + 4*x^2 + 2*x^3 + x^4 + 4*x^5)]
        sage: C.divisor_of_function(b[1])
        -4*(3 + y, 3 + x) - (2 + x, 4 + y)
        sage: C.divisor_of_function(b[0])
        (y, x) - 4*(3 + y, 3 + x) - 2*(2 + x, 4 + y)

    """
    Esuppx = [Esupp[i][0] for i in range(len(Esupp))]
    E0 = [[Ecoeffs[i],Esupp[i]] for i in range(len(Esupp))]
    D0 = [[Dcoeffs[i],Dsupp[i]] for i in range(len(Dsupp))]
    F = f.base_ring()
    C = Curve(f)
    pts = C.rational_points() 
    numpts = len(pts)
    R = f.parent()
    FR = FractionField(R)
    x,y = R.gens()
    genus = C.geometric_genus()
    dim = sum(Ecoeffs)+1-genus   ## RR thrm *if* div(E)>= 2g-1
    suppE = []
    suppEidx = []
    for p in pts:     #### finds the pts in Esupp with coeff>0
        if list(p) in [list(pt[1]) for pt in E0]:
            i = Esupp.index(p)
            if Ecoeffs[i]>0:
                suppE.append(pts[i])
                suppEidx.append(i)
    for i in range(len(suppEidx)):  ########### start here with new for loop
        ptindx = suppEidx[i]
        pt0 = [Esupp[ptindx][0],Esupp[ptindx][1]]
        pivot_pt = [Esupp[ptindx][0],-Esupp[ptindx][1]]
        lcs = C.local_coordinates(pivot_pt,5)
        yt = lcs[1]
        xt = lcs[0]
        if i==0:
            denom = [(x-Esuppx[int(k)])**Ecoeffs[int(k)] for k in range(len(Esupp))]
            f0 = (y-pivot_pt[1])/prod(denom)
            ############## this is the WRONG pivot function if Esupp has more than one
            ############## point in an iota orbit (ie, a (x0,y0) and a (xo,-y0). 
            ############## In that case, group the pts in Esupp into iota
            ############## orbits and only take those which have maximum exponent
            glist = [f0]
        if i>0:
            f0 = add(glist)
        f1 = [1/(x-Esupp[ptindx][0])**j for j in range(1,Ecoeffs[ptindx]+1)]
        f1.reverse()
        R0 = PolynomialRing(F,3,names = [str(x),str(y),"t"])
        vars0 = R0.gens()
        t = vars0[2]
        g = f0
        coeff0 = degree_lowest_rational_function(g(xt,yt),t)[1]
        for j in range(Ecoeffs[ptindx]):                        ## for j: subtracts the pivot function from
            h = f1[j]                                           ## the seed function to reduce the 
            ldg = degree_lowest_rational_function(g(xt,yt),t)   ## order of the pole at pivot_pt
            d0 = ldg[0]                                         
            ldh = degree_lowest_rational_function(h(xt,yt),t)
            d1 = ldh[0]
            if d1==0:                                          ## this case seems not to be used ?
                return g,glist
            if d1==d0: 
                coeff1 = R(ldh[1])
                fcn_matching_pole = h*R(ldg[1])/coeff1
                g = g - fcn_matching_pole
                glist.append(-fcn_matching_pole)
    return g,glist

def divisor_of_function(C,r):
        """
        Return the divisor of a function on a curve.
        
        INPUT:
             r is a rational function on X
             
        OUTPUT:
             list -- The divisor of r represented as a list of coefficients and points.
                     (TODO: This will change to a more structural output in the future.)

        WARNING: Produced correct output if the support of the divisor consists only of
        *unramified* points.

        EXAMPLES:
            sage: F = GF(5)
            sage: A2 = AffineSpace(2, F, names = 'xy')
            sage: R = A2.coordinate_ring()
            sage: x, y = R.gens()
            sage: f = y^2 - x^9 - x
            sage: C = AffineCurve_GF(A2,f)
            sage: K = FractionField(R)
            sage: r = 1/x
            sage: divisor_of_function(C,r) 
                  [[-1, (0, 0)]]
            sage: r = 1/x^3
            sage: divisor_of_function(C,r) 
                  [[-3, (0, 0)]]
        """
        self = C
        F = self.base_ring()
        f = self.defining_polynomial()
        pts = self.rational_points() 
        numpts = len(pts)
        R = f.parent()
        x,y = R.gens()
        R0 = PolynomialRing(F,3,names = [str(x),str(y),"t"])
        vars0 = R0.gens()
        t = vars0[2]
        divf = []
        for pt0 in pts:
                #if pt0[1] != F(0):
                lcs = self.local_coordinates(pt0,5)
                yt = lcs[1]
                xt = lcs[0]
                ldg = degree_lowest_rational_function(r(xt,yt),t)
                if ldg[0] != 0:
                    divf.append([ldg[0],pt0])
        return divf