def CRT_nf(reslist, Ilist, check=True):
    r"""
    Solve Chinese remainder problem over a number field.

    INPUT:

    - ``reslist`` -- a list of residues, i.e. integral number field elements

    - ``Ilist`` -- a list of integral ideas, assumed pairsise coprime

    - ``check`` (boolean, default True) -- if True, result is checked

    OUTPUT:

    An integral element x such that x-reslist[i] is in Ilist[i] for all i.
    """
    n = len(reslist)
    if n == 0:
        # we have no parent field so this is all we can do
        return 0
    if n == 1:
        return reslist[0]
    if n > 2:
        # use induction / recursion
        x = CRT_nf([reslist[0], CRT_nf(reslist[1:], Ilist[1:])],
                   [Ilist[0], prod(Ilist[1:])])
        if check:
            check_CRT_nf(reslist, Ilist, x)
        return x
    # now n=2
    r = Ilist[0].element_1_mod(Ilist[1])
    x = ((1 - r) * reslist[0] + r * reslist[1]).mod(prod(Ilist))
    if check:
        check_CRT_nf(reslist, Ilist, x)
    return x

    def dbd(self, d):
        """
        Return matrix of <d>.

        INPUT:

            - `d` -- integer

        OUTPUT:

            - a matrix modulo 2

        EXAMPLES::

            sage: from mdsage import *
            sage: C = KamiennyCriterion(29)
            sage: C.dbd(2)
            22 x 22 dense matrix over Finite Field of size 2 (use the '.str()' method to see the entries)
            sage: C.dbd(2)^14==1
            True
            sage: C.dbd(2)[0]
            (0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 0, 0, 0, 1, 1, 0, 1, 0, 1, 1, 1)
        """
        d=ZZ(d)
        if self.verbose: tm = cputime(); mem = get_memory_usage(); print "dbd start"
        try: return self._dbd[d % self.p]
        except AttributeError: pass
        # Find generators of the integers modulo p:
        gens = Integers(self.p).unit_gens()
        orders = [g.multiplicative_order() for g in gens]
        # Compute corresponding <z> operator on integral cuspidal modular symbols
        
        X = [self.M.diamond_bracket_operator(z).matrix() for z in gens]
        if self.verbose: print "time and mem", cputime(tm), get_memory_usage(mem), "create d"
        X = [x.restrict(self.S_integral, check=False) for x in X]
        if self.verbose: print "time and mem", cputime(tm), get_memory_usage(mem), "restrict d"
        
        X = [matrix_modp(x) for x in X]
        if self.verbose: print "time and mem", cputime(tm), get_memory_usage(mem), "mod d"
        # Take combinations to make list self._dbd of all dbd's such that
        # self._dbd[d] = <d>
        from itertools import product
        v = [None] * self.p
        for ei in product(*[range(i) for i in orders]):
            di = prod(g**e for e,g in zip(ei,gens)).lift()
            m = prod(g**e for e,g in zip(ei,X))
            v[di] = m
        if self.verbose: print "time and mem", cputime(tm), get_memory_usage(mem), "mul"

        assert v.count(None) == (self.p-euler_phi(self.p))
        self._dbd = v
        if self.verbose: print "time and mem", cputime(tm), get_memory_usage(mem), "bdb finnished"
        return v[d % self.p]

def upper_bound_index_cusps_in_JG_torsion(G,d, bound = 60):
    """
    INPUT:
        
    - G - a congruence subgroup
    - d - integer, the size of the rational cuspidal subgroup
    - bound (optional, default = 60) - the bound for the primes p up to which to use
      the hecke matrix `T_p - <p> - p` for bounding the torsion subgroup
    
    OUTPUT:
        
    - an integer `i` such that `(\#J_G(\QQ)_{tors})/d` is a divisor of `i`.
    
    EXAMPLES::

        sage: from mdsage import *
        sage: d = rational_cuspidal_classgroup(Gamma1(23)).cardinality()
        sage: upper_bound_index_cusps_in_JG_torsion(Gamma1(23),d)
        1

    """
    N = G.level()
    M=ModularSymbols(G);
    Sint=cuspidal_integral_structure(M)
    kill_mat=(M.star_involution().matrix().restrict(Sint)-1)
    kill=kill_mat.transpose().change_ring(ZZ).row_module()
    for p in prime_range(3,bound):
        if not N % p ==0:
            kill+=kill_torsion_coprime_to_q(p,M).restrict(Sint).change_ring(ZZ).transpose().row_module()
        if kill.matrix().is_square() and kill.matrix().determinant()==d:
            #print p
            break
    kill_mat=kill.matrix().transpose()
    #print N,"index of torsion in stuff killed",kill.matrix().determinant()/d
    if kill.matrix().determinant()==d:
        return 1
        
    pm=integral_period_mapping(M)
    period_images1=[sum([M.coordinate_vector(M([c,infinity])) for c in cusps])*pm for cusps in galois_orbits(G)]
    
    m=(Matrix(period_images1)*kill_mat).stack(kill_mat)
    diag=m.change_ring(ZZ).echelon_form().diagonal()
    #print diag,prod(diag)
    assert prod(diag)==kill.matrix().determinant()/d
    
    period_images2=[M.coordinate_vector(M([c,infinity]))*pm for c in G.cusps() if c != Cusp(oo)]
    m=(Matrix(period_images2)*kill_mat).stack(kill_mat)
    m,denom=m._clear_denom()
    diag=(m.change_ring(ZZ).echelon_form()/denom).diagonal()
    #print diag
    #print prod(i.numerator() for i in diag),"if this is 1 then :)"
    return prod(i.numerator() for i in diag)

    def cardinality(self):
        r"""
        Return the cardinality of self

        EXAMPLES::

            sage: S = Subsets([1,1,2,3],submultiset=True)
            sage: S.cardinality()
            12
            sage: len(S.list())
            12

            sage: S = Subsets([1,1,2,2,3],submultiset=True)
            sage: S.cardinality()
            18
            sage: len(S.list())
            18

            sage: S = Subsets([1,1,1,2,2,3],submultiset=True)
            sage: S.cardinality()
            24
            sage: len(S.list())
            24
        """
        from sage.all import prod
        return Integer(prod(k+1 for k in self._d.values()))

def get_T2(K, S, unit_first=True, verbose=False):
   u = -1 if K==QQ else  K(K.unit_group().torsion_generator())
   from KSp import IdealGenerator
   Sx = [IdealGenerator(P) for P in S]
   if unit_first:
      Sx = [u] + Sx
   else:
      Sx = Sx + [u]
   r = len(Sx)
   r2 = r*(r-1)//2
   N = prod(S,1)
   T2 = {}
   primes = primes_iter(K,1)
   p = primes.next()
   while len(T2)<r+r2:
      p = primes.next()
      while p.divides(N):
         p = primes.next()
      # Compute the values alpha_p(Delta) for Delta in Sx:
      ro = alphalist(p,Sx)
      # Compute the set I(P) of i for which alpha_p(Delta_i)=1:
      ij = Set([i for i,vi in enumerate(ro) if vi])
      if verbose:
         print("P={}, I(P)={}".format(p,ij))
      # Keep I(p) and p if #I(p)=1 or 2 and it's a new [email protected]
      if len(ij) in [1,2] and not ij in T2:
         T2[ij] = p
   return T2

def NonCubicSet(K,S, verbose=False):
    u = -1 if K==QQ else  K(K.unit_group().torsion_generator())
    from KSp import IdealGenerator
    Sx = [u] + [IdealGenerator(P) for P in S]
    r = len(Sx)
    d123 = r + binomial(r,2) + binomial(r,3)
    vecP = vec123(K,Sx)
    A = Matrix(GF(2),0,d123)
    N = prod(S,1)

    primes = primes_iter(K,None)
    T = []
    while A.rank() < d123:
        p = primes.next()
        while p.divides(N):
            p = primes.next()
        v = vecP(p)
        if verbose:
            print("v={}".format(v))
        A1 = A.stack(vector(v))
        if A1.rank() > A.rank():
            A = A1
            T.append(p)
            if verbose:
                print("new A={} with {} rows and {} cols".format(A,A.nrows(),A.ncols()))
                print("T increases to {}".format(T))
    return T

def global_minimality_class(E):
    r"""
    Returns the ideal class representing the obstruction to this
    elliptic curve having a global minimal model.

    INPUT:

    - ``E`` -- an elliptic curve over a number field

    OUTPUT:

    An ideal class of the base number field, which is trivial if and
    only if E has a global minimal model, and which can be used to
    find global and semi-global minimal models.
    """
    K = E.base_field()
    Cl = K.class_group()
    if K.class_number() == 1:
        return Cl(1)
    D = E.discriminant()
    dat = E.local_data()
    primes = [d.prime() for d in dat]
    vals = [d.discriminant_valuation() for d in dat]
    I = prod([P ** ((D.valuation(P) - v) // 12) for P, v in zip(primes, vals)],
             E.base_field().ideal(1))
    return Cl(I)

def get_T0_mod3(K,S, flist=None, verbose=False):
   # Compute all absolutely irreducible quartics unramified outside S
   # if not supplied:
   if flist == None:
      from S4 import abs_irred_extensions
      flist = abs_irred_extensions(K,S)
   if verbose:
      print("quartics: {}".format(flist))

# Append a poly with lam3=0
   x = polygen(K)
   flist0 = flist + [x**3]
   n = len(flist)

# Starting with no primes, compute the lambda matrix
   plist = []
   vlist = [lamvec(f,plist,lam3) for f in flist0]
   ij = equal_vecs(vlist)
   if verbose:
      print("With plist={}, vlist={}, ij={}".format(plist,vlist,ij));
   N = prod(S,1) * prod([f.discriminant() for f in flist])
# As long as the vectors in vlist are not all distinct, find two
# indices i,j for which they are the same and find a new prime which
# distinguishes these two, add it to the list and update.  The primes
# must not be in S or divide the discriminant of any of the quartics.
   while ij:
      i,j = ij
      if j==n:
         p = get_p_1(K,flist[i],N, lam3)
      else:
         p = get_p_2(K,flist[i],flist[j],N, lam3)
      plist = plist + [p]
      if verbose:
         print("plist = {}".format(plist))
      vlist = [lamvec(f,plist,lam3) for f in flist0]
      ij = equal_vecs(vlist)
      if verbose:
         print("With plist={}, vlist={}, ij={}".format(plist,vlist,ij));
   vlist = vlist[:-1]

   # Sort the primes into order and recompute the vectors:
   plist.sort()
   vlist = [lamvec(f,plist,lam3) for f in flist]
   return flist, plist, vlist

def an_dict_from_ap(ap, N, B):
    r"""
    Give a dict ``ap`` of the `a_p`, for primes with norm up to `B`,
    return the dictionary giving all `a_I` for all ideals `I` up to
    norm `B`.

    NOTE: This code is specific to Dirichlet series of elliptic
    curves.

    INPUT:
        - ``ap`` -- dictionary of ap, as output, e.g., by the ap_dict function
        - `N` -- ideal; conductor of the elliptic curve
        - `B` -- positive integer
    
    OUTPUT: dictionary mapping reduced rep of primes (N(p),p.reduced_gens()) of ideals to integers

    NOTE: This should be really, really slow.  It's really just a toy
    reference implementation.
    """
    from sage.all import prod   # used below
 
    F = N.number_field()
    A = F.ideals_of_bdd_norm(B)

    an = dict(ap)
    
    for n in sorted(A.keys()):
        X = A[n]
        for I in X:
            if an.has_key(reduced_rep(I)):
                # prime case, already done
                pass
            else:
                # composite case
                fac = I.factor()
                if len(fac) == 0:
                    # unit ideal
                    an[reduced_rep(I)] = ZZ(1)
                elif len(fac) > 1:
                    # not a prime power, so just multiply together
                    # already known Dirichlet coefficients, for
                    # prime power divisors (which are all known).
                    an[reduced_rep(I)] = prod(an[reduced_rep(p**e)] for p, e in fac)
                else:
                    p, e = fac[0]
                    # a prime power
                    if p.divides(N):
                        # prime divides level
                        an[reduced_rep(I)] = an[reduced_rep(p)]**e
                    else:
                        # prime doesn't divide conductor: a_{p^e} = a_p*a_{p^(e-1)} - Norm(p)*a_{p^(e-2)}
                        assert e >= 2
                        an[reduced_rep(I)] = (an[reduced_rep(p)] * an[reduced_rep(p**(e-1))]
                                                - p.norm()*an[reduced_rep(p**(e-2))])
                        
    return an

def factorization(original_poly):
    poly = ZZT(original_poly)
    assert poly[0] == 1
    if poly == 1:
        return [1]
    try:
        facts = poly.factor()
    except NotImplementedError:
        try:
            # try to sort out the memory leak
            PolynomialRing(ZZ, 'T', implementation='NTL')([1] + [0]*(len(original_poly) - 2) + [1])._pari_with_name().factor()
            facts = PolynomialRing(ZZ, 'T', implementation='NTL')(original_poly).factor()
        except NotImplementedError:
            raise
    # if the factor is -1+T^2, replace it by 1-T^2
    # this should happen an even number of times, mod powers
    out = [[-g if g[0] == -1 else g, e] for g, e in facts]
    assert prod( g**e for g, e in out ) == poly, "%s != %s" % (prod( [g**e] for g, e in out ), poly)
    return [[g.list(), e] for g,e  in out]

def euler_p_factor(root_list, PREC):
    ''' computes the coefficients of the pth Euler factor expanded as a geometric series
      ax^n is the Dirichlet series coefficient p^(-ns)
    '''
    PREC = floor(PREC)
    # return satake_list
    R = LaurentSeriesRing(CF, 'x')
    x = R.gens()[0]
    ep = prod([1 / (1 - a * x) for a in root_list])
    return ep + O(x ** (PREC + 1))

def ppower_norm_ideal_from_label(K,lab):
    r""" return the ideal of prime-power norm from its label.
    """
    n, i = [int(c) for c in lab.split(".")]
    p, f = ZZ(n).factor()[0]
    make_keys(K,p)
    PP = K.primes_dict[p]
    ff = [P.residue_class_degree() for P in PP]
    vec = exp_vec_wt(f,ff)[i-1]
    return prod([P**v for P,v in zip(PP,vec)])

def level_attributes(level):
    # returns level_radical, level_primes, level_is_prime, level_is_prime_power, level_is_squarefree, level_is_square
    fact = Integer(level).factor()
    level_primes = [elt[0] for elt in fact]
    level_radical = prod(level_primes)
    level_is_prime_power = len(fact) == 1
    level_is_prime = level_is_prime_power and level_radical == level
    level_is_square = all( elt[1] % 2 == 0 for elt in fact)
    level_is_squarefree = all( elt[1] == 1 for elt in fact)
    return [level_radical, level_primes, level_is_prime, level_is_prime_power, level_is_squarefree, level_is_square]

def get_coeffs_p_over_nf(curve, prime_number, accuracy=20 , conductor=None):
    """                                                                         
    Computes the inverse of product of L_prime on all primes above prime_number,
    then returns power series of L_p up to desired accuracy.                    
    But will not return power series if need not do so (depends on accuracy).   
    """
    if conductor is None:
        conductor = curve.conductor()
    primes = curve.base_field().prime_factors(prime_number)
    series_p = [get_factor_over_nf(curve, prime_id, prime_number, conductor, accuracy) for prime_id in primes]
    return ( prod(series_p).O(accuracy) )**(-1)

def ideals_of_norm(K,n):
    r""" Return a list of all ideals of norm n (sorted).  Cached.
    """
    if not hasattr(K,'ideal_norm_dict'):
        K.ideal_norm_dict = {}
    if not n in K.ideal_norm_dict:
        if n==1:
            K.ideal_norm_dict[n] = [K.ideal(1)]
        else:
            K.ideal_norm_dict[n] = [prod(Q) for Q in cartesian_product_iterator([ppower_norm_ideals(K,p,e) for p,e in n.factor()])]
    return K.ideal_norm_dict[n]

    def alpha_alt(roots, p):
        """ Computes invariant alpha in the 'ALT' case
        """
        from sage.all import prod, Permutation
        tmp = prod(roots[i] - roots[j] for i in range(len(roots)) for j in range(i))
        try:
            frob_perm = Permutation(alternating_group(frobenius_permutation(roots, p)))
        except TypeError:
            raise TypeError("The Frobenius element does not generate an element of the alternating group")

        sign = -(-1) ** are_conjugate_in_alternating(frob_perm, data_perm)
        return tmp * sign

def X_to_kappa(f,kappa):
    """
    f -- A polynomial in X. The constant term is ignored.
    kappa -- A function so that kappa(a) is kappa_a.
    
    Returns::
        A polynomial in kappas, converting ...
    """
    psi_list = []
    for expon,coef in f.dict().items():
        #print expon[0], coef
        if expon[0] !=0:
            psi_list += [expon[0]]*coef
        
    #print psi_list
    result = 0 
    for p in SetPartitions(range(len(psi_list))):   
        #print p  
        #print [ kappa( sum((psi_list[i] for i in s)) ) for s in p]
        result +=   prod((factorial(len(s) - 1) for s in p)) * prod(( kappa( sum((psi_list[i] for i in s)) ) for s in p ))
    return result

    def test_spanning_trees(self):
        r"""
        Test coset representatives obtained from spanning trees for even
        subgroup (Kulkarni's method with generators ``S2``, ``S3`` and Verrill's
        method with generators ``L``, ``S2``).

        EXAMPLES::

            sage: from sage.modular.arithgroup.tests import Test
            sage: Test().test_spanning_trees() #random
        """
        from sage.all import prod
        from all import SL2Z
        from arithgroup_perm import S2m,S3m,Lm

        G = random_even_arithgroup(self.index)

        m = {'l':Lm, 's':S2m}
        tree,reps,wreps,gens = G._spanning_tree_verrill()
        assert reps[0] == SL2Z([1,0,0,1])
        assert wreps[0] == ''
        for i in xrange(1,self.index):
            assert prod(m[letter] for letter in wreps[i]) == reps[i]
        tree,reps,wreps,gens = G._spanning_tree_verrill(on_right=False)
        assert reps[0] == SL2Z([1,0,0,1])
        assert wreps[0] == ''
        for i in xrange(1,self.index):
            assert prod(m[letter] for letter in wreps[i]) == reps[i]

        m = {'s2':S2m, 's3':S3m}
        tree,reps,wreps,gens = G._spanning_tree_kulkarni()
        assert reps[0] == SL2Z([1,0,0,1])
        assert wreps[0] == []
        for i in xrange(1,self.index):
            assert prod(m[letter] for letter in wreps[i]) == reps[i]
        tree,reps,wreps,gens = G._spanning_tree_kulkarni(on_right=False)
        assert reps[0] == SL2Z([1,0,0,1])
        assert wreps[0] == []
        for i in xrange(1,self.index):
            assert prod(m[letter] for letter in wreps[i]) == reps[i]

    def _prod_f_A_G(Gs, kappa, psi, psi2):
        """
        Used in the computation of product. Do not call directly.
        """
        result = 1
        G = Gs.G.strataG
        A = Gs.A.strataG
        for v in range(1,G.num_vertices()+1):
            #print "Gs", Gs
            #print "G,g",v
            result *= prod(( sum((kappa.get((j,w),0) for w in Gs.alpha_inv[v]))**f for j,f in G.kappa_on_v(v))) * prod(( Gs.psi_no_loop(edge,v,ex,psi) for edge, ex in G.psi_no_loop_on_v(v) )) * prod(( Gs.psi_loop(edge,v, ex1, ex2, psi,psi2) for edge, ex1, ex2 in G.psi_loop_on_v(v) ))

        return result                    

def ppower_norm_ideals(K,p,f):
    r""" Return a sorted list of ideals of K of norm p^f with p prime
    """
    make_keys(K,p)
    if not hasattr(K,'ppower_dict'):
        K.ppower_dict = {}
    if not (p,f) in K.ppower_dict:
        PP = K.primes_dict[p]
        # These vectors are sorted, first by unweighted weight (sum of
        # values) then lexicographically with the reverse ordering on Z:
        vv = exp_vec_wt(f,[P.residue_class_degree() for P in PP])
        Qs = [prod([P**v for P,v in zip(PP,v)]) for v in vv]
        K.ppower_dict[(p,f)] = Qs
    return K.ppower_dict[(p,f)]