Make `elliptic_j` exact when possible

[grhkm21]

When $\tau$ is imaginary quadratic, the $j$-invariant $j(\tau)$ is an algebraic integer. It would be nice if the elliptic_j method returns an element in $\overline{\mathbb{Q}}$ in this case. This probably requires implementing a new algorithm, since the current elliptic_j implementation calls pari which is purely numerical.

This is a special case of #15354.

[user202729]

Proof of concept implementation (although I don't know if the error added is sufficient to ensure answer correctness)

def symbolic_elliptic_j(t):
    """
    Implementation based on https://mathoverflow.net/a/36232 and https://mathoverflow.net/a/379727.

    EXAMPLES::

        sage: symbolic_elliptic_j(sqrt(11)*I)
        1.122662367132124?e9
        sage: symbolic_elliptic_j(QQbar.zeta(6))
        0
        sage: symbolic_elliptic_j(QQbar.zeta(6)*2)
        54000
        sage: symbolic_elliptic_j(sqrt(-3)*2)
        2.835807690422284?e9
        sage: symbolic_elliptic_j(sqrt(-3)*3)
        1.510132287045155?e14
        sage: symbolic_elliptic_j(sqrt(-3)*5)
        4.282443629598017?e23
        sage: symbolic_elliptic_j(I/2)
        287496
        sage: symbolic_elliptic_j(I*4/3)
        5138.898474224591?
        sage: symbolic_elliptic_j(sqrt(2))
        Traceback...
        AssertionError: 
        sage: symbolic_elliptic_j(2)
        Traceback...
        AssertionError: 
    """
    t = QQbar(t)
    assert t.imag()>0
    assert t.degree()==2
    f = t.minpoly()
    f *= lcm(a.denom() for a in f)  # scale by appropriate multiplier to make the polynomial in ZZ[x] and primitive (TODO might fold into dedicated method in real code)
    R.<x> = QQ[]
    l = [t1
     for f in BinaryQF_reduced_representatives(f.discriminant())
     for t1 in R(f.polynomial().subs(y=1)).roots(QQbar, multiplicities=False)
     if t1.imag()>0]  # list of t' such that j(t') is conjugate of j(t)
    F = CBF
    prec = F.prec()
    while True:
        R.<x> = F[]
        error = F(1).add_error(2^-prec)
        # We need to manually add error as above because the default constructor may decide
        # to not give any error, as in the first example below::
        #
        # sage: 1+10^-500 in CBF(CC(1))
        # False
        # sage: 1.1+10^-500 in CBF(1.1)
        # True
        # 
        # We can't just use ``add_error`` because that takes absolute error
		# while we need both absolute and relative error.
        # absolute error is for the case when result of ``elliptic_j()`` is 0 or close to it
        f = prod(x-(F(elliptic_j(r, prec+5)) * error).add_error(2^-prec) for r in l)
        # add 5 guard bits just in case
		# (TODO can put this into ``elliptic_j_ball(t, prec) -> RealBallField(prec)`` later ---
		# the only guarantee needed is that it returns higher precision result when prec increase, and the answer is correct)
        try:
            f = ZZ[x](f)
        except ValueError:
            # not accurate enough
            prec *= 2
            F = ComplexBallField(prec)
            continue
        break
    l = f.roots(QQbar, multiplicities=False)  # list of candidates of the result

    if len(l)==1:
        return l[0]

    while True:
        error = F(1).add_error(2^-prec)
        r=(F(elliptic_j(t, prec+5))*error).add_error(2^-prec) # TODO reuse ``elliptic_j_ball`` as mentioned above
        l1=[r1 for r1 in l if r1 in r]
        assert l1
        if len(l1)>1:
            # not accurate enough
            prec *= 2
            F = ComplexBallField(prec)
            continue
        return l1[0]