r"""
###################################################################
#       Copyright (C) 2011 chris wuthrich
#  Distributed under the terms of the GNU General Public License ###################################################################

This file contains  contains the algorithms to compute the 
monodromy and class group pairings on elliptic
curves over number fields. This file was written to do experiments for the article

Jean Gillibert, Christian Wuthrich, The class group pairing and $p$-descent on elliptic curves, 2011.

Given an elliptic Curve $E$ over a number field $K$, there is a pairing on $E(K)$ with values in a certain group denoted by $LogPic( OK, S )$. This group contains the class group of $K$ as a subgroup and the quotient is a direct sum of copies of $\QQ/\ZZ$ one for each place in $S$, where $S$ is the set of bad places for $E$. If one of the two points has good reduction everywhere then the value of thw pairing belongs to the usual class group. The image of the pairing into the quotient is the monodromy pairing on $E$.

EXAMPLES.
Here is an elliptic curve over a number field with two interesting points in it::

    sage: E = EllipticCurve("158c1")
    sage: P = E(13,-15)
    sage: P.order()
    5
    sage: K.<t> = QuadraticField(-79)
    sage: EK = E.base_extend(K)
    sage: Q = EK.lift_x(101/9)
    sage: P_K = EK(P)

This file implements the monodromy pairing at places of bad reduciton. In our example both points have bad reduction somewhere but at different places::

    sage: monodromy_pairing(P_K,Q)
    [0, 0, 0]
    sage: monodromy_pairing(Q,Q)
    [1/4, 1/4, 0]
    sage: monodromy_pairing(P_K,P_K)
    [4/5, 4/5, 0]

There is the definition of the logarithmic class group pairing, which yields an element in the a new class defined in this file::

    sage: logarithmic_class_group_pairing(P_K,Q)
    Logarithmic ideal class  (4, 1/2*t + 1/2)
    sage: L = logarithmic_class_group_pairing(P_K+Q,Q)
    sage: L
    Logarithmic ideal class  (2, 1/2*t + 1/2) + 1/4*((2, 1/2*t + 1/2)) + 1/4*((2, 1/2*t - 1/2))
    sage: L+L
    Logarithmic ideal class  (4, 1/2*t + 1/2) + 1/2*((2, 1/2*t + 1/2)) + 1/2*((2, 1/2*t - 1/2))
    sage: L*4
    Logarithmic ideal class  (2, 1/2*t - 1/2)
    sage: L.order()
    20

Finally, in the case that two points are orthogonal under the monodromy pairing, this yields the class group pairing. These are elements of the usual class group::

    sage: class_group_pairing(P_K,Q)
    Fractional ideal class (4, 1/2*t + 1/2)
    sage: class_group_pairing(4*Q,Q)
    Trivial principal fractional ideal class
    sage: class_group_pairing(4*Q,P_K)
    Fractional ideal class (447853394671/4674, 1/738492*t + 480247863355/9348)


AUTHOR:

- chris wuthrich 2011-07


TODO:

in log_correction_term_instar(P,Q,v,n) some cases are not implemented yet.

"""

# these functions define the local contribution that corrects
# the naive pairing using the denominator ideals

def log_correction_term(P,Q, v):
    r"""
    Returns the term `-Q^*(\Phi_P)`. This term corrects the class group
    to get the logarithmic class group pairing.
    
    INPUT::
    
    - ``P`` - a point on an elliptic curve
    
    - ``Q`` - another point on the same curve
    
    - ``v`` - a prime in the base field of the curve
    
    OUTPUT::
    
    a rational number
    
    EXAMPLES::
    
        sage: E = EllipticCurve("930d2")
        sage: P = E.gens()[0]
        sage: E.tamagawa_number(5)
        14
        sage: log_correction_term(P,2*P,5)
        -10/7
        sage: log_correction_term(3*P,-8*P,5)
        -6/7
        
    """
    E = P.curve()
    K = E.base_ring()
    assert K == QQ or v.number_field() == K
    
    dav = E.local_data(v)
    if dav.has_good_reduction():
        return 0
    if dav.has_split_multiplicative_reduction():
        return log_correction_term_split_multiplicative(P,Q,v)
    if dav.has_nonsplit_multiplicative_reduction():
        return log_correction_term_nonsplit_multiplicative(P,Q,v)
   
    # else we have additive reduction
    k = dav.kodaira_symbol()._pari_code()
    if k == 2 or k == -2: #II or II*
        return 0
    if k == 3 or k == -3: #III or III*
        return log_correction_term_3(P,Q,v,k)
    if k == 4 or k == -4: #IV or IV*
        return log_correction_term_4(P,Q,v,k)
    # else type I*n
    return log_correction_term_instar(P,Q,v,-k-4)
    
def change_multiplicative_model(E,v): 
    r"""
    Given an elliptic curve with split multiplicative reduction at v, this gives back a Weierstrass model of the form
    `y^2 + xy = x^3 +a*x^2 +b*x+c` with a,b,c in v.
    
    INPUT:
    
    - ``E`` - an elliptic curve over a number field
    
    - ``v`` - a prime of this field
    
    OUTPUT: another elliptic curve
    
    .. note: This is nowhere near a clever or an effective implementation.
    
    EXAMPLE::
    
        sage: E = EllipticCurve("11a1")
        sage: E.has_split_multiplicative_reduction(11)
        True
        sage: change_multiplicative_model(E,11)
        Elliptic Curve defined by y^2 + x*y = x^3 - 2563/248*x^2 + 1449459/61504*x - 402788551/15252992 over Rational Field

        sage: R.<x> = QQ[]
        sage: K.<t> = NumberField(x^5-3*x+17)
        sage: E2 = EllipticCurve([1,t-1,0,3*t^2, 27+(t-1)^3])
        sage: v = K.ideal(3,t-1)
        sage: E2.has_split_multiplicative_reduction(v)
        True
        sage: E = E2.change_weierstrass_model([-7,2+t,-t,8+3*t])
        sage: E
        Elliptic Curve defined by y^2 + (2/7*t-1/7)*x*y + (-1/49*t-18/343)*y = x^3 + (-1/49*t^2+5/49*t+5/49)*x^2 + (15/2401*t^2+29/2401*t)*x + (6/117649*t^3-47/117649*t-50/117649) over Number Field in t with defining polynomial x^5 - 3*x + 17
        sage: change_multiplicative_model(E,v)
        Elliptic Curve defined by y^2 + x*y = x^3 + (8001259978752/501404359854503*t^4+24806351311776/501404359854503*t^3-4088602278702/501404359854503*t^2+1231716749135/501404359854503*t-108142895514971/501404359854503)*x^2 + (181753474741544229054610860/251406332081103939687329377009*t^4+315407562888242694425468292/251406332081103939687329377009*t^3-1466711529489151290708532812/251406332081103939687329377009*t^2-3232904831209593513239550084/251406332081103939687329377009*t+1116845256547213320346810056/251406332081103939687329377009)*x + (-6376076813528295233990069939907290808987/126056231000494521873599606475708996873321527*t^4+4546496564571637759722177401541581517287/126056231000494521873599606475708996873321527*t^3-18588947934735099934524441866003471364957/126056231000494521873599606475708996873321527*t^2-18342206588189367207812673349463704275183/126056231000494521873599606475708996873321527*t+121733795986009204084108316873244736446106/126056231000494521873599606475708996873321527) over Number Field in t with defining polynomial x^5 - 3*x + 17


        sage: E = EllipticCurve("1141a1")
        sage: E.has_split_multiplicative_reduction(7)
        True
        sage: change_multiplicative_model(E,7)
        Elliptic Curve defined by y^2 + x*y = x^3 + 77973/20752*x^2 + 24655869/6728836*x + 227822969101/139636804672 over Rational Field

        
    """
    # make it first minimal at v
    
    E = E.local_data(v).minimal_model()
    
    # then we change to a model of the form [1, 0, 0, *, *]
    
    s = (1-E.a1())/2
    r = (s**2 +E.a1() * s - E.a2())/3
    t = - (E.a1()*r + E.a3())/2
    Eq = E.change_weierstrass_model([1,r,s,t])
    # print "first step", Eq.a_invariants()
    assert Eq.a1() == 1 and Eq.a2() == 0 and Eq.a3() == 0
 
    # second we change it so that the singularity is at (0,0) with a3 = a1= 0
 
    A = Eq.a4()
    B = Eq.a6()
    th = (A-36*B)/(1-48*A)
    Eq = Eq.change_weierstrass_model([1,-2*th,0,th])
    # print "second step", Eq.a_invariants()
    assert Eq.a1() == 1 and  Eq.a3() == 0 and valuation(Eq.a6(),v) > 0 
    assert valuation(Eq.a6(),v) > 0
 
    # finally we change it so that the reduction is y^2 + xy = x^3
    
    a = Eq.a2()
    usq = 1+4*a
    # if the reduction is split multiplicaitve, then usq is a square mod v
    K = Eq.base_ring()
    if K == QQ:
        usq_red = IntegerModRing(v)(usq)
        u = usq_red.sqrt().lift()
    else:
        kv = K.ring_of_integers().residue_field(v)
        usq_red = kv(usq)
        u_red = usq_red.sqrt()
        u = kv.lift(u_red)
    s = (u-1)/2
    Eq = Eq.change_weierstrass_model([u,0,s,0])
    assert Eq.a1() == 1 and  Eq.a3() == 0 and valuation(Eq.a6(),v) > 0 and valuation(Eq.a2(),v)>0 and valuation(Eq.a4(),v)>0
    
    return Eq
 
def _shift_singularity(f,P):
    """
    Given an equation of the form
    y^2 + x*y = d* x^3 + a*x^2 + b*x + c
    with d and a having positive valuation,
    this procedure yeilds back another one such
    that (0,0) is the singular point. 
    The given point P is also moved to the new
    equation.
    
    INPUT:
    
    - ``f`` - a polynomial in two variables x and y, trusted to be of the form mentionned.
    
    - ``P`` - a solution to f
    
    OUTPUT: another polynomial and a solution to it.
    
    .. note: No checks are done.
    
    """
    R = f.parent()
    K = R.base_ring()
    x = R.gens()[0]
    y = R.gens()[1]
    d,a,_,_,b,c = (-f).coefficients()
    aa = a - 6*b*d
    bb = 12* b**2 *d - 4*a*b 
    cc = 4*a* b**2 - 8 * b**3 * d - b**2 + c
    # print "d,a,b,c", d, a, b, c
    # print "shift to ", d, aa, bb, cc
    Q =  [P[0]+2*b,P[1]-b]
    # print P, Q
    g = y**2 +x*y - d* x**3 -aa* x**2 -bb * x -cc
    # print f(P), g(Q)
    return g, Q

def _blowup(f,P,pi):
    """
    Given an equation of the form 
     y^2 + x*y = d* x^3 + a*x^2 + b*x + c
    with a singularity at (0,0),
    this procedure yeilds back the strict transform of a blowup at v
    """
    # pick a uniformizer
    R = f.parent()
    x = R.gens()[0]
    y = R.gens()[1]
    d,a,_,_,b,c = f.coefficients()
    # assert valuation(c,v) > 1 and valuation(b,v) > 0
    return y**2 + x*y + d*pi* x**3 + a* x**2 + b/pi * x + c/pi**2, [P[0]/pi,P[1]/pi]

def component_of_a_point(P,v):
    r"""
    Given a point on a split multiplicative curve,
    determine the index of the component on which it lies
    
    INPUT:
    
    - ``P`` - a point on an elliptic curve over a number field
    
    - ``v`` - a place in this field
    
    OUTPUT : an integer, well defined up to sign modulo n
    
    EXAMPLES ::
    
        sage: E = EllipticCurve("930d2")
        sage: P = E.gens()[0]
        sage: component_of_a_point(P,5)
        9
        sage: component_of_a_point(2*P,5)
        4
        sage: E.tamagawa_number(5)
        14
         
        sage: R.<x> = QQ[]
        sage: K.<t> = NumberField(x^3-3*x^2+17*x-1)
        sage: E = EllipticCurve( [1,23,0,18-26*t,0])
        sage: P = E(t,1)
        sage: v = K.ideal(11,t+2)
        sage: component_of_a_point(P,v)
        2
    
    """
    E = P.curve()
    K = E.base_ring()
    n = E.tamagawa_number(v)
    if P.has_good_reduction(v):
        return ZZ(0)
    if n == 1:
        return ZZ(1)
    
    #choose a uniformiser
    if K == QQ:
        pi = v
    else:
        pi = [pi for pi in v.gens() if pi.valuation(v) == 1][0]
    
    # now n > 2 and we need to blow up
    
    # change to a model of the form [1,*,0,*,*] with reduction y^2 + xy = x^3
    Eq = change_multiplicative_model(E,v)
    # move P and Q over to the new model
    iso = E.isomorphism_to(Eq)
    P = iso(P)
    # print "point ",P," on Weierstrass equation changed to ", Eq.a_invariants() 
    assert valuation(P[0],v) > 0 and valuation(P[1],v) > 0
 
    # initiate
    R = PolynomialRing(K, 2, ["x", "y"])
    x, y = R.gens()
    a = Eq.a2()
    b = Eq.a4()
    c = Eq.a6()
    d = 0
    f = y**2 + x*y - x**3 - a * x**2 - b*x - c
    assert valuation(a,v) >0 and valuation(b,v) > 0 and valuation(c,v) > 1
    # print valuation(P[0],v), valuation(P[1],v), valuation(P[0]+P[1],v)
    assert f(P[0],P[1]) == 0
        
    
    while True: 
        f, P = _blowup(f,P, pi)
        #print "point P is %s and strict transform is %s"%(P,f) 
        #print valuation(P[0],v), valuation(P[1],v), valuation(P[0]+P[1],v)
        #print f.coefficient({x:0,y:0})
        #print type(f.coefficient({x:0,y:0}))
        #print f.coefficient({x:0,y:0}) in K
        assert f(P[0],P[1]) == 0
        d += 1
        _,a,_,_,b,c = (-f).coefficients()
        assert valuation(a,v) > 0
        if valuation(b**2 - c,v) == 0:
            # the blowup is smooth
            assert n == 2*d
            return ZZ(d)    
        f, P = _shift_singularity(f,P)
        # print "point P is %s and shifted eq is %s"%(P,f) 
        # print valuation(P[0],v), valuation(P[1],v), valuation(P[0]+P[1],v)
        assert f(P[0],P[1]) == 0
        
        assert valuation(c,v) >0 
        
        if valuation(b**2-c,v) == 1:
            # the blowup consits of two line, but we won't need a further blowup
            assert n == 2*d + 1
            if valuation(P[1],v) > 0:
                return d + 1
            elif valuation(P[0]+P[1],v) >0:
                return d
            else:
                raise RunTimeError
        
        # check if the point is no longer singular
        if valuation(P[0],v) <1 or valuation(P[1],v) < 1:
            if valuation(P[1],v) > 0:
                return n-d
            elif valuation(P[0]+P[1],v) >0:
                return d
            else:
                raise RunTimeError
    
        assert n >= 2*(d+1) # breaks to make sure the loop if finite if there are bugs.
        
        #do the next blowup
                
    

def log_correction_term_split_multiplicative(P,Q,v):
    r"""    
    Returns the term -Q*(Phi_P). This term corrects the class group
    to get the logarithmic class group pairing in the case of a
    place of split multiplicative reduction.
    
    INPUT::
    
    - ``P`` - a point on an elliptic curve
    
    - ``Q`` - another point on the same curve
    
    - ``v`` - a prime in the base field of the curve
    
    OUTPUT::
    
    a rational number
    
    EXAMPLES::
    
        sage: E = EllipticCurve("930d2")
        sage: P = E.gens()[0]
        sage: log_correction_term_split_multiplicative(P,-P,5)
        -25/14

    """
    E = P.curve()
    if P.has_good_reduction(v) or Q.has_good_reduction(v):
        return 0
    # number of components. = tamagawa number
    n = -4 + E.local_data(v).kodaira_symbol()._pari_code()
    if n == 2:
        return -ZZ(1)/2
    
    j = component_of_a_point(P,v)
    k = component_of_a_point(Q,v)
    if j <= k :
        return ZZ(j)*(k-n)/n
    else:
        return ZZ(k)*(j-n)/n
        
def log_correction_term_nonsplit_multiplicative(P,Q,v):
    r"""
    Returns the term -Q*(Phi_P). This term corrects the class group
    to get the logarithmic class group pairing in the case of a
    place of non-split multiplicative reduction.
    
    INPUT::
    
    - ``P`` - a point on an elliptic curve
    
    - ``Q`` - another point on the same curve
    
    - ``v`` - a prime in the base field of the curve
    
    OUTPUT::
    
    a rational number
    
    EXAMPLES::
            
        sage: E = EllipticCurve("1058c2")
        sage: E.local_data(2).kodaira_symbol()
        I6
        sage: E.local_data(2).has_nonsplit_multiplicative_reduction()
        True
        sage: P = E.gens()[0]
        sage: Q = E.gens()[1]
        sage: log_correction_term_nonsplit_multiplicative(P,Q,2)
        -3/2
        sage: log_correction_term_nonsplit_multiplicative(P,P,2)
        -3/2
        sage: log_correction_term_nonsplit_multiplicative(P,Q,2)
        -3/2
        sage: log_correction_term_nonsplit_multiplicative(P,2*P,2)
        0

    """
    if P.has_good_reduction(v) or Q.has_good_reduction(v):
        return 0
    #otherwise they are on the same component
    E = P.curve()
    n = -4 + E.local_data(v).kodaira_symbol()._pari_code()
    assert n % 2 == 0
    return -ZZ(n)/4  
   
def log_correction_term_3(P,Q,v,k):
    r"""
    Returns the term -Q*(Phi_P). This term corrects the class group
    to get the logarithmic class group pairing in the case of a
    place of additive reduction of type III and III*.
    
    INPUT::
    
    - ``P`` - a point on an elliptic curve
    
    - ``Q`` - another point on the same curve
    
    - ``v`` - a prime in the base field of the curve
    
    - ``k`` - either -3 or 3 meaning III or III*
    
    OUTPUT::
    
    a rational number
    
    EXAMPLES::
    
        sage: E = EllipticCurve("4312l2")
        sage: P = E.gens()[0]
        sage: Q = E.gens()[1]
        sage: [(dav.prime().gens()[0], dav.kodaira_symbol()) for dav in E.local_data()]
        [(2, III*), (7, III), (11, I1)]
        sage: log_correction_term_3(P,Q,2,-3)
        -3/2
        sage: log_correction_term_3(P,Q,3,3)
        0

    """
    # the type is III or III*, so there are two connected components
    if P.has_good_reduction(v) or Q.has_good_reduction(v):
        return 0
    #otherwise they are on the same component
    if k == 3: # type III
        return -ZZ(1)/2
    else:      # type III*
        return -ZZ(3)/2
        
def log_correction_term_4(P,Q,v,k):
    r"""
    Returns the term -Q*(Phi_P). This term corrects the class group
    to get the logarithmic class group pairing in the case of a
    place of additive reduction of type IV and IV*.
    
    INPUT::
    
    - ``P`` - a point on an elliptic curve
    
    - ``Q`` - another point on the same curve
    
    - ``v`` - a prime in the base field of the curve  
    
    - ``k`` - either -4 or 4 meaning IV or IV*

    OUTPUT::
    
    a rational number
    
    EXAMPLES::
    
        sage: E = EllipticCurve("916c1")
        sage: E.local_data(2).kodaira_symbol()
        IV
        sage: P = E.gens()[0]
        sage: Q = E.gens()[1]
        sage: log_correction_term_4(P,Q,2,4)
        0
        sage: P.has_good_reduction(2)
        True
        sage: Q.has_good_reduction(2)
        False
        sage: log_correction_term_4(Q,Q,2,4)
        -2/3
        sage: log_correction_term_4(Q,2*Q,2,4)
        -1/3
        
        sage: E  = EllipticCurve("1431a1")
        sage: E.local_data(3).kodaira_symbol()
        IV*
        sage: P = E.gens()[0]
        sage: Q = E.gens()[1]
        sage: log_correction_term_4(P,Q,3,-4)
        -4/3
        sage: log_correction_term_4(2*P,Q,3,-4)
        -2/3
        sage: log_correction_term_4(3*P,Q,3,-4)
        0

    """
    if P.has_good_reduction(v) or Q.has_good_reduction(v):
        return 0
    #otherwise they are either on the same component 
    if (P-Q).has_good_reduction(v):
        if k == 4: #type IV
            return -ZZ(2)/3
        else: #type IV*
            return -ZZ(4)/3
    # or not on the same component
    if k == 4: # type IV
        return -ZZ(1)/3
    else:
        return -ZZ(2)/3
        
    
def log_correction_term_instar(P,Q,v,n):
    r"""    
    Returns the term -Q*(Phi_P). This term corrects the class group
    to get the logarithmic class group pairing in the case of a
    place of additive reduction of type I*_n.
    
    INPUT::
    
    - ``P`` - a point on an elliptic curve
    
    - ``Q`` - another point on the same curve
    
    - ``v`` - a prime in the base field of the curve
    
    - ``n`` - an integer such that the type has to be I*_n
    
    OUTPUT::
    
    a rational number
    
    EXAMPLES::
    
        sage: E = EllipticCurve("944e1")
        sage: E.local_data(2).kodaira_symbol()
        I2*
        sage: E.local_data(2).tamagawa_number()
        4
        sage: P = E.gens()[0]
        sage: Q = E.gens()[1]
        sage: log_correction_term_instar(P,Q,2,2)
        -1/2
        sage: log_correction_term_instar(2*P,Q,2,2)
        0
        sage: log_correction_term_instar(P,P,2,2)
        Traceback (most recent call last): 
        ...  
        NotImplementedError: The awkward case is not yet fully implemented.

        sage: E = EllipticCurve("1088j2")
        sage: E.local_data(2).tamagawa_number()
        4
        sage: E.local_data(2).kodaira_symbol()
        I5*
        sage: P = E.gens()[0]
        sage: Q = E.gens()[1]
        sage: log_correction_term_instar(P,Q,2,5)
        0
        sage: log_correction_term_instar(Q,2*Q,2,5)
        -1/2
        sage: log_correction_term_instar(Q,3*Q,2,5)
        -1/2
        sage: log_correction_term_instar(Q,4*Q,2,5)
        0

    """
    if P.has_good_reduction(v) or Q.has_good_reduction(v):
        return 0
    if not (P-Q).has_good_reduction(v):
        return -ZZ(1)/2
    # or not on the same component
    # either they are both on Gamma_1 or not
    if P.curve().local_data(v).tamagawa_number() == 2:
        return -ZZ(1)
    if n == 0:
        return -ZZ(1)
    else:
        raise NotImplementedError, "The awkward case is not yet fully implemented."
    # we have to move the equation such that a_i have valuation
    # [>= 1, 1, >=2, >=3, >=4]
    # then the valuation of the x-coordinate should tell us if P and Q are both on Gamma_1 or on one of the others.
    # if it is Gamma_1 then we return -1, else it is -1-n/4.
    
def e_zero(E,v):
    r"""
    The term e_v(O) at the place v.
    
    INPUT:
   
    - ``E`` - an elliptic curve over a number field
    
    - ``v`` - a prime of this number field
    
    OUTPUT: an integer.

    EXAMPLES::
    
        sage: E = EllipticCurve([0,7])
        sage: K.<t> = NumberField(x^6-7)
        sage: EK = E.base_extend(K)
        sage: v = K.ideal(t)
        sage: v.is_prime()
        True
        sage: e_zero(EK,v)
        1

    """
    # uu is the quotient of the Neron differential at pp divided by
    # the differential associated to this particular equation E
    K = E.base_ring()
    if K.degree() == 1:
        return 0
    dav = E.local_data(v)
    uu = E.isomorphism_to(dav.minimal_model()).u
    hv = valuation(uu,v)
    return hv
    
###############################

_fp = lambda x: x-floor(x)

def monodromy_pairing(P,Q,v=None):
    r"""
    Given two points on an elliptic curve over a number field, 
    this returns the monodromy pairing between them. Optionally,
    the prime v in the field can be given to specify the fibre
    in which we compute; otherwise, the vector of all values at
    bad places in given.
    
    The result is given back as a rational number between 0 and 1. 
    It should be considered as an element of QQ/ZZ.
    
    INPUT:
    
    - ``P`` - a point on an ellliptic curve over a number field
    
    - ``Q`` - another point on the same curve
    
    - ``v`` - (optional) a prime in the field
    
    OUTPUT: a rational number or a list of rational numbers.
    
    EXAMPLES::
    
        sage: E = EllipticCurve("930d2")
        sage: P = E.gens()[0]
        sage: monodromy_pairing(P,P)
        [0, 0, 11/14, 0]
        sage: monodromy_pairing(P,2*P,5)
        4/7
        sage: monodromy_pairing(P,7*P,5)
        1/2
        sage: monodromy_pairing(P,14*P,5)
        0
        sage: (14*P).has_good_reduction()
        True
        
        sage: R.<x> = QQ[]
        sage: K.<t> = NumberField(x^3-3*x^2+17*x-1)
        sage: E = EllipticCurve( [1,23,0,18-26*t,0])
        sage: P = E(t,1)
        sage: v = K.ideal(11,t+2)
        sage: monodromy_pairing(P,P)
        [0, 0, 0, 1/3, 0, 0, 0]
       
    """
    E = P.curve()
    if Q.curve() != E:
        raise ValueError, "The points must lie on the same curve"
    if v == None:
        vs = [dav.prime() for dav in E.local_data() if dav.has_bad_reduction()]
        if E.base_ring() == QQ:
            vs = [v.gens()[0] for v in vs]
        return [monodromy_pairing(P,Q,v) for v in vs]
    else:
        return _fp( log_correction_term(P,Q,v) )

###

  
###############################
# the class to hold our pairing
###############################

class LogPic(SageObject):
    """
    The class which is a logarithmic class group
    """
    
    def __init__(self,F,S):
        r"""
        The logarithmic class group of a number field with
        respect to a finite set of places. The group is additive.
        
        INPUT:
        
        - ``F`` - a number field
        
        - ``S`` - a list of primes of F
        
        OUTPUT: A LogPic class.
        
        EXAMPLES::
        
            sage: K.<t> = NumberField(x^5-9*x^4+22*x-71)
            sage: E = EllipticCurve("30a1").base_extend(K)
            sage: S = [dav.prime() for dav in E.local_data()]
            sage: A = LogPic(K,S)
            sage: A
            Logarithmic class group of the Number Field in t with defining polynomial x^5 - 9*x^4 + 22*x - 71 with respect to the finite set [(3, t - 1), (2, t^2 - 7*t - 3), (5, t^2 - 6*t - 4), (2, t^3 - 8*t^2 - 7*t - 3), (3, t^3 - 7*t^2 - 15*t - 5), (5, t^3 - 8*t^2 - 9*t - 1)]

            sage: c = K.class_group().gens()[0]
            sage: di = dict([[v,1/v.smallest_integer()] for v in S])            
            sage: L = A(c,di)
            sage: L
            Logarithmic ideal class  (2, t^2 - 7*t - 3) + 1/3*((3, t - 1)) + 1/5*((5, t^3 - 8*t^2 - 9*t - 1)) + 1/3*((3, t^3 - 7*t^2 - 15*t - 5)) + 1/2*((2, t^2 - 7*t - 3)) + 1/5*((5, t^2 - 6*t - 4)) + 1/2*((2, t^3 - 8*t^2 - 7*t - 3))
            sage: L+L
            Logarithmic ideal class  (4, t^2 - 7*t - 3) + 2/3*((3, t - 1)) + 2/5*((5, t^3 - 8*t^2 - 9*t - 1)) + 2/3*((3, t^3 - 7*t^2 - 15*t - 5)) + 2/5*((5, t^2 - 6*t - 4))
            sage: L*15
            Logarithmic ideal class  (1) + 1/2*((2, t^2 - 7*t - 3)) + 1/2*((2, t^3 - 8*t^2 - 7*t - 3))
            sage: L*30
            Logarithmic ideal class  (1)
            
            sage: L == L*31
            True
            sage: L == L*16
            False
            sage: L*5 == 0
            False
            sage: L*30 == 0
            True

            sage: -L
            Logarithmic ideal class  (4, t^2 - 7*t - 3) + 2/3*((3, t - 1)) + 4/5*((5, t^3 - 8*t^2 - 9*t - 1)) + 2/3*((3, t^3 - 7*t^2 - 15*t - 5)) + 1/2*((2, t^2 - 7*t - 3)) + 4/5*((5, t^2 - 6*t - 4)) + 1/2*((2, t^3 - 8*t^2 - 7*t - 3))
            sage: M = L*7
            sage: M - L
            Logarithmic ideal class  (1) + 1/5*((5, t^3 - 8*t^2 - 9*t - 1)) + 1/5*((5, t^2 - 6*t - 4))

        """
        self._F = F
        self._S = S
        self._clgp = F.class_group()
        
    def _repr_(self):
        Srep = "[" 
        te = 0
        for v in self._S:
            if te == 0:
                te = 1
            else:
                Srep += ", "
            Srep +=  v._repr_short() 
        Srep += "]"
        return "Logarithmic class group of the %s with respect to the finite set %s"%(self._F,Srep)
        
    #def __contains__(self,el):   
        
    def _element_class(self):
        return LogPicElement
    
    def __call__(self,c,li=None):
        """
        Given an element of the class group and (optionnally) a dictionnary of 
        {place in S:rational number} this coerces this into the
        logarithmic class group.
        """
        if c == 0:
            return self.zero_element()
        if c not in self._clgp:
            c = self._clgp(c)
            # raise ValueError, "Not in the class group of the corresponding field"
        return LogPicElement(c,li,parent=self)
        
    def zero_element(self):
        return self.__call__(self._clgp(self._F.ideal(1)), dict([[v,0] for v in self._S]))

    def number_field(self):
        return self._F
    
    def logarithmic_places(self):
        return self._S
        
class LogPicElement(SageObject):
    """
    The class of elements in a logarithmic class group.
    """
    
    def _reduce(self):
        r"""
        This shifts the integer part of the S-part into the class group.
        The resulting representation is unique.
        """
        S = self._S
        di = self._di
        c = self._c
        A = c.parent()
        newdi = []
        for v in S:
            nv = di[v]
            mv = ZZ(floor(nv))
            if mv != 0 and not v.is_principal():
                vc = A(v)
                # print vc, type(vc), mv, type(mv), c, type(c), vc.order()
                # bug in clgroups is solved by # 11234
                vcm = vc**mv
                c = c * vcm
            newdi.append([v,_fp(nv)])
                
        self._di = dict(newdi)
        self._c = c
      
    def __init__(self, c, di, parent=None):
        r"""      
        Given an element of the class group and (optionnally) a dictionnary of 
        {place in S:rational number} this coerces this into the
        logarithmic class group (written additively).
 
        INPUT:
        
        - ``c`` - an ideal class in an number field
        
        - ``di`` - a dictionary of places:rational numbers
        
        EXAMPLES::
            
            sage: K.<t> = QuadraticField(-5)
            sage: c = K.class_group().gens()[0]
            sage: S = union(K.primes_above(3),K.primes_above(11))
            sage: di = dict([[v,1/v.residue_class_degree()] for v in S])
            sage: L = LogPic(K,S)(c,di)
            sage: L
            Logarithmic ideal class  (2, t + 1) + 1/2*((11))
            sage: L*2
            Logarithmic ideal class  (1)
            sage: L+L
            Logarithmic ideal class  (1)
            sage: -L
            Logarithmic ideal class  (2, t + 1) + 1/2*((11))
            sage: di = dict([[v,1/v.residue_class_degree()/5] for v in S])
            sage: M = LogPic(K,S)(c,di)
            sage: M
            Logarithmic ideal class  (2, t + 1) + 1/5*((3, t + 1)) + 1/5*((3, t + 2)) + 1/10*((11))
            sage: M+L
            Logarithmic ideal class  (1) + 1/5*((3, t + 1)) + 3/5*((11)) + 1/5*((3, t + 2))
            sage: M-L
            Logarithmic ideal class  (1) + 1/5*((3, t + 1)) + 3/5*((11)) + 1/5*((3, t + 2))
            sage: (M-L)*5
            Logarithmic ideal class  (1)
            
            sage: L.order()
            2
            sage: M.order()
            10
            
            sage: di = dict([[v,v.residue_class_degree()/7] for v in S])
            sage: L2 = LogPic(K,S)(c,di)
            sage: L2.order()
            14
            sage: L2
            Logarithmic ideal class  (2, t + 1) + 1/7*((3, t + 1)) + 1/7*((3, t + 2)) + 2/7*((11))
            sage: L2*7
            Logarithmic ideal class  (2, t + 1)
            sage: (L2*7).is_in_class_group()
            True
            sage: (L2*7).class_group_part()
            Fractional ideal class (2, t + 1)

        """
        I = c.ideal()
        F = I.number_field()
        S = di.keys()
        if parent is None:
            parent = LogPic(F,S)
        self._S = S
        self._F = F
        self._di = di
        self._c = c
        self._parent = parent
        self._reduce()
    
    def parent(self):
        return self._parent
        
    def number_field(self):
        return self._F
    
    def logarithmic_places(self):
        return self._S
    
    def logarithmic_part(self):
        self._reduce()
        return self._di
        
    def class_group_part(self):
        r"""
        This is only well defined if all places in S are principal. 
        But our representation of it is unique.
        """
        self._reduce()
        return self._c
        
    def _repr_(self):
        self._reduce()
        if self._c.is_one():
            s = "(1)"
        else:
            s = self._c.ideal()._repr_short()
        for v in self._S:
            if  self._di[v] != 0:
                s += " + " + self._di[v].str() + "*(" + v._repr_short()+")"
        return "Logarithmic ideal class  %s"%s
        
    def __eq__(self,other):
        self._reduce()
        if other == 0:
            other = self.parent().zero_element()
        else:
            other._reduce()
        return self._F == other._F and self._S == other._S and self._c == other._c and self._di == other._di 
        
    def __ne__(self, other):
        return not self.__eq__(other)

    def __nonzero__(self):
        return not self.__eq__(self.parent().zero_element())
        
    def is_zero(self):
        return self.__eq__(self.parent().zero_element())
        
    def __add__(self,other):
        self._reduce()
        other._reduce()
        c = self._c * other._c
        di = dict([ [v,self._di[v] + other._di[v]] for v in self._S])
        return self.parent()(c,di)
    
    def __sub__(self,other):
        self._reduce()
        other._reduce()
        c = self._c / other._c
        di = dict([ [v,self._di[v] - other._di[v]] for v in self._S])
        return self.parent()(c,di)
  
    
    def __neg__(self):
        self._reduce()
        c = (self._c)**(-1)
        di = dict([ [v,-self._di[v]] for v in self._S])
        return self.parent()(c,di)
    
    
    def _rmul_(self,scalar):
        self._reduce()
        c = (self._c)**(scalar)
        di = dict([ [v,self._di[v]*scalar] for v in self._S])
        return self.parent()(c,di)
  
    def _lmul_(self,scalar):
        return self._rmul_(scalar)
        
    def __mul__(self,scalar):
        return self._rmul_(scalar)
        
    def order(self):
        r"""
        Returns the additive order of this logarithmic class group element
        """
        m = LCM( [self._di[v].denominator() for v in self._S])
        L = self*m
        k = L._c.order()
        return m*k
        
    def is_in_class_group(self):
        r"""
        Tests weather the element belongs to the usual class group of the field.
        """
        self._reduce()
        return all( self._di[v] == 0 for v in self._S)

###############################
### The pairing
###############################

def denominator_ideal(P):
    r"""
    Let `P` be a point on an elliptic curve `E` defined over a number field `K`.
    Returns the integral ideal `e` in `K` such that the `x`-coordinate of `E`
    generates the factional ideal `a\cdot e^{-2}` for some integral ideal `a`
    coprime to `e`. The `y`-coordinates give `b\cdot e^{-3}`.
    
    If P is O, then the definition of this ideal is slightly different.
    
    INPUT: a point P on an elliptic curve over a number field
    
    OUTPUT: An ideal of the number field
    
    EXAMPLES::
    
        sage: K.<t> = NumberField(x^3-3*x^2+17*x-1)
        sage: E = EllipticCurve( [1,23,0,18-26*t,0])
        sage: P = E(t,1)
        sage: denominator_ideal(P)
        Fractional ideal (1)
        sage: Q = 7*P
        sage: denominator_ideal(Q)
        Fractional ideal (9918992371934701017391856165112865945001108889, 3*t - 2803221975214349290753197246757953520429824)
        sage: denominator_ideal(E(0))
        Fractional ideal (1)
        
        sage: E = EllipticCurve([0,7])
        sage: K.<t> = NumberField(x^6-7)
        sage: EK = E.base_extend(K)
        sage: denominator_ideal(EK(0))
        Fractional ideal (7, t)

    """
    E = P.curve()
    K = E.base_ring()
    OK = K.ring_of_integers()
    
    if P == E(0):
        S = [dav.prime() for dav in E.local_data()]
        I = prod( v **  e_zero(E,v) for v in S)
        return I
    # else:
        
    # wrong method
    # I = K.fractional_ideal(P[0])
    # ff = I.factor()
    # e = OK.ideal(1)
    # for fa in ff:
    #     if fa[1] < 0:
    #        assert fa[1] % 2 == 0, "oooops"
    #        e *= fa[0]^(-fa[1]//2)
    
    # correct method
    xx = P[0]
    eesq = xx.denominator()
    aa = eesq*P[0]
    esq = OK.ideal(eesq)
    a = OK.ideal(aa)
    # print "xx = ",xx
    # print "e_squared = ", esq.factor()
    # print "a =         ", a.factor()
    esq = esq/(esq+a)
    
    yy = P[1]
    eecu = yy.denominator()
    bb = eecu*P[1]
    ecu = OK.ideal(eecu)
    b = OK.ideal(bb)
    # print "yy = ", yy
    # print "e_cubed   = ", ecu.factor()
    # print "b         = ", b.factor()
    ecu = ecu/(ecu+b)

    # print "e_cubed   = ", ecu.factor()
    # print "e_squared = ", esq.factor()
    ee = ecu/esq
    assert ee.is_integral()
    return ee
 
def logarithmic_class_group_pairing(P,Q):
    r"""
    Computes the pairing on the Mordell-Weil group with values in the
    logarithmic class group. Given two points `P` and `Q` on an elliptic
    curve over a number field, this returns an element of the 
    logarithmic class grop of this field with respect to the finite
    set of bad primes of the curve.
    
    INPUT:
    
    - ``P`` - a point on an elliptic curve over a number field
    
    - ``Q`` - another point on the same curve
    
    OUTPUT: an element of the class LogPic
    
    EXAMPLES::
        
        sage: E = EllipticCurve('11a1')
        sage: K.<t> = NumberField(x^2+47)
        sage: EK = E.base_extend(K)
        sage: P = EK(-2, -1/2*t - 1/2)   # the two generator of the free part of E(K)
        sage: Q = EK(-1, -1/2*t - 1/2)
        sage: T = EK(5,5)                # a generator of the torsion part of order 5
        sage: logarithmic_class_group_pairing(P,Q)
        Logarithmic ideal class  (1)                                                           
        sage: logarithmic_class_group_pairing(P,T)                                            
        Logarithmic ideal class  (7, 1/2*t + 3/2)                                              
        sage: logarithmic_class_group_pairing(Q,T)
        Logarithmic ideal class  (6, 1/2*t + 1/2)
       
    second example::
        
        sage: E = EllipticCurve('14a1')
        sage: K.<t> = NumberField(x^2+23)
        sage: EK = E.base_extend(K)
        sage: P = EK(-7/8*t - 5/8, -7/16*t - 101/16)
        sage: Q = EK(1/2*t + 5/2, -3/2*t - 11/2)
        sage: logarithmic_class_group_pairing(P,Q)
        Logarithmic ideal class  (1)
        sage: logarithmic_class_group_pairing(Q,Q)
        Logarithmic ideal class  (2, 1/2*t - 1/2) + 1/2*((2, 1/2*t - 1/2))
        sage: logarithmic_class_group_pairing(P,P)
        Logarithmic ideal class  (18, 1/2*t + 25/2) + 1/3*((7))

    yet another example::
    
        sage: E = EllipticCurve("1141a1")
        sage: P = E.gens()[0]
        sage: Q = E.gens()[1]
        sage: K.<t> = NumberField(x^3-7)
        sage: EK = E.base_extend(K)
        sage: P = EK(P)
        sage: Q = EK(Q)
        sage: logarithmic_class_group_pairing(P,Q)
        Logarithmic ideal class  (1)
        sage: logarithmic_class_group_pairing(P,P)
        Logarithmic ideal class  (1)
        sage: logarithmic_class_group_pairing(Q,Q)
        Logarithmic ideal class  (1) + 2/5*((t))
        sage: logarithmic_class_group_pairing(Q,2*Q)-logarithmic_class_group_pairing(Q,Q)*2
        Logarithmic ideal class  (1)
    
    """
    E = P.curve()
    if E != Q.curve():
        raise ValueError, "The points must lie on the same curve"
        
    eP = denominator_ideal(P)
    eQ = denominator_ideal(Q)
    eP_plus_Q = denominator_ideal(P+Q)
    ezero = denominator_ideal(E(0))
    K = E.base_field()
    CK = K.class_group()
    c = CK(eP_plus_Q/eP/eQ*ezero)
    
    S  = [dav.prime() for dav in E.local_data() if dav.has_bad_reduction()]
    if E.base_ring() == QQ:
        vs = [v.gens()[0] for v in vs]
        
    di = dict([[v,log_correction_term(P,Q, v)] for v in S])
    Gr = LogPic(K,S)
    return Gr(c,di)
    
    
def class_group_pairing(P,Q):
    r"""
    Given two points `P` , `Q` on an elliptic curve `E` over a number field `K`,
    this returns the value of the class group pairing of `P` and `Q` in the class group of `K`. This is only well-defined if the two points are orthogonal under the monodromy pairing.
    
    INPUT:
        
    - ``P`` - a point on an elliptic curve over a number field
    
    - ``Q`` - another point on the same curve
    
    OUTPUT: an element of the class group of the field

    EXAMPLES::
        
        sage: E = EllipticCurve('11a1')
        sage: K.<t> = NumberField(x^2+47)
        sage: EK = E.base_extend(K)
        sage: P = EK(-2, -1/2*t - 1/2)   # the two generator of the free part of E(K)
        sage: Q = EK(-1, -1/2*t - 1/2)
        sage: T = EK(5,5)                # a generator of the torsion part of order 5
        sage: class_group_pairing(P,Q)
        Trivial principal fractional ideal class
        sage: class_group_pairing(P,T)                                            
        Fractional ideal class (7, 1/2*t + 3/2)
        sage: class_group_pairing(Q,T)
        Fractional ideal class (6, 1/2*t + 1/2)
      
    second example::
        
        sage: E = EllipticCurve('14a1')
        sage: K.<t> = NumberField(x^2+23)
        sage: EK = E.base_extend(K)
        sage: P = EK(-7/8*t - 5/8, -7/16*t - 101/16)
        sage: Q = EK(1/2*t + 5/2, -3/2*t - 11/2)
        sage: class_group_pairing(P,Q)
        Trivial principal fractional ideal class
        sage: class_group_pairing(Q,Q)
        Traceback (most recent call last):
        ...
        ValueError: The points must be orthogonal under the monodromy pairing. Use the logarithmic pairing instead.
 
    """    
    if P.curve() != Q.curve():
        raise ValueError, "The points must lie on the same curve"

    L = logarithmic_class_group_pairing(P,Q)
    if not L.is_in_class_group():
        raise ValueError, "The points must be orthogonal under the monodromy pairing. Use the logarithmic pairing instead."
    
    return L.class_group_part()