{{{id=1|
def can_height(P, F, N=10, prec=53, debug=false):
    r""""
        Evaluates the canonical height of the projective point with respect to the map F
        using a telescoping sum of Green's functions, as in Silverman ADS ex. 5.29
    """
    from sage.rings.real_mpfr import is_RealField
    K = P.domain().base()
    Kdeg = K.absolute_degree()
    if debug: print 'Base field K is ', K
    if debug: print 'Degree of base field K is ', Kdeg

    # TODO: Check K is a number field

    R = RealField(prec=prec)

    # TODO: Compute N based on prec, rather than as a defaulting parameter

    P.normalize_coordinates()
    if debug: print 'We have normalized the coordinates of P to ', P
    d = F.degree()
    if debug: print 'deg(F) = ', d
    D=F.codomain().ambient_space().dimension_relative()
    if debug: print 'Ambient dimension D = ', D
    localht = R(0)

    # Archimedean local heights
    emb = K.places(prec=prec)
    Ptup = list(P)
    Ftup = list(F)
    for v in emb:
        # :: WARNING: If places is fed the default Sage precision of 53 bits,
        # it uses Real or Complex Double Field in place of RealField(prec) or ComplexField(prec)
        # the function is_RealField does not identify RDF as real, so we test for that ourselves.    
        if is_RealField(v.codomain()) or v.codomain() is sage.rings.real_double.RDF:
            dv = 1
        else:
            dv = 2
        if debug: print 'Current v is ', v
        x = Ptup
        if debug: print 'Current x is ', x
        if debug: print 'dv = ', dv, ' and Kdeg = ', Kdeg
        if debug: print 'For this v, dv/Kdeg = ', R(dv/Kdeg)
        for i in range(N+1):
            if debug: print 'Starting i = ', i
            Pv = [ v(t).abs() for t in x ]
            if debug: print 'Pv = ', Pv
            m = -1
            j = -1
            #compute the maximum absolute value of entries of a, and where it occurs
            for n in range(D+1):
                if Pv[n] > m:
                    j = n
                    m = Pv[n]
            if debug: print 'max coord abs val is now m = ', m
            # add to sum for the Green's function
            localht += R(dv/Kdeg)*((1/R(d))**R(i))*(R(m).log())
            if debug: print 'localht is now = ', localht
            if i < N:
                if debug: print 'About to scale x = ', x
                x = [ t/x[j] for t in x ]
                if debug: print 'Scaled x is = ', x
                x = [ Ff(x) for Ff in Ftup ]
                if debug: print 'F(x) is now ', x
    bad_primes = []
    for b in Ptup:
        bad_primes += b.denominator_ideal().prime_factors()
    if D==1:
        bad_primes += K(f.resultant()).support()
    else:
        from sage.rings.polynomial.macaulay_resultant import macaulay_resultant
        bad_primes += K(macaulay_resultant(list(F))).support()
    bad_primes = list(Set(bad_primes))
    for p in bad_primes:
        if debug: print 'Starting bad prime p = ', p
        x = Ptup
        if debug: print 'Current x is ', x
        for i in range(N+1):
            if debug: print 'Starting i = ', i
            Pv = [ R(t.abs_non_arch(p)) for t in x ]
            if debug: print 'Pv = ', Pv
            m = -1
            j = -1
            #compute the maximum absolute value of entries of a, and where it occurs
            for n in range(D+1):
                if Pv[n] > m:
                    j = n
                    m = Pv[n]
            if debug: print 'max coord abs val is now m = ', m
            # add to sum for the Green's function
            localht += R(p.residue_class_degree()*p.absolute_ramification_index()/Kdeg)*((1/R(d))**R(i))*(R(m).log())
            if debug: print 'localht is now = ', localht
            if i < N:
                if debug: print 'About to scale x = ', x
                x = [ t/x[j] for t in x ]
                if debug: print 'Scaled x is = ', x
                x = [ Ff(x) for Ff in Ftup ]
                if debug: print 'F(x) is now ', x
    return localht
///
}}}

{{{id=2|
z = var('z')
K.<a> = NumberField(z^3-z-3); K
PS.<x,y> = ProjectiveSpace(K,1)
H = Hom(PS,PS)
f = H([5*x^2+y^2,y^2])
print f
Q = PS(3,2)
Q
///
Scheme endomorphism of Projective Space of dimension 1 over Number Field in a with defining polynomial z^3 - z - 3
  Defn: Defined on coordinates by sending (x : y) to
        (5*x^2 + y^2 : y^2)
(3/2 : 1)
}}}

{{{id=18|
%time
can_height(Q,f,N=11)
///
2.75017626386820
CPU time: 0.21 s,  Wall time: 0.21 s
}}}

{{{id=59|
K.<a> = NumberField(z^3-z-3); K
P2K.<X,Y,Z> = ProjectiveSpace(K,2)
H2K = Hom(P2K,P2K)
f2 = H2K([5*X^2+Y^2+Z^2,3*Y^2,Z^2])
Q2 = P2K(a^2,2,3); Q2
///
(1/3*a^2 : 2/3 : 1)
}}}

<p>Differs by a bounded amount depending on the map $f$ from the global height:</p>

{{{id=58|
%time
can_height(Q2,f2,N=11)
///
1.76676716451166
CPU time: 2.99 s,  Wall time: 2.99 s
}}}

{{{id=64|
%time
Q2.global_height()
///
1.09861228866811
CPU time: 0.12 s,  Wall time: 0.12 s
}}}

{{{id=40|
%time
can_height(Q,f,N=12)
///
2.75056919304604
CPU time: 0.44 s,  Wall time: 0.44 s
}}}

{{{id=50|
PQ.<z,w> = ProjectiveSpace(QQ,1)
HQ = Hom(PQ,PQ)
gq = HQ([5*z^2+w^2,w^2])
Qq = PQ([3,2])
///
}}}

{{{id=44|
%time
Qq.canonical_height(gq, N=12)
///
2.7505691930460382934277268982
CPU time: 0.02 s,  Wall time: 0.02 s
}}}

{{{id=43|
%time
can_height(Q,f,N=15)
///
}}}

{{{id=6|
can_height(Q,f)-Qq.canonical_height(gq)
///
-4.44089209850063e-16
}}}

<p>Let's do something interesting now, where the point comes from our number field. Let's recall what K is:</p>

{{{id=21|
K
///
Number Field in a with defining polynomial z^3 - z - 3
}}}

{{{id=52|
Q = PS(a+1,2)
///
}}}

{{{id=53|
%time
can_height(Q,f,N=11)
///
1.68337344250905
CPU time: 1.49 s,  Wall time: 1.49 s
}}}

<p>Differs by a bounded amount from the global height, of course:</p>

{{{id=54|
Q.global_height()
///
0.789669763582290
}}}

{{{id=57|

///
}}}