def _is_a_cover(mt0, mt1):
    r"""
    Define the cover relations.

    Return ``True`` if and only if the second argument is a cover of
    the first one.

    EXAMPLES::

        sage: import sage.combinat.alternating_sign_matrix as asm
        sage: asm._is_a_cover([[1,2,3],[1,2],[1]], [[1,2,3],[1,3],[1]])
        True
        sage: asm._is_a_cover([[1,2,3],[1,3],[2]], [[1,2,3],[1,2],[1]])
        False
    """
    diffs = 0
    for (a,b) in itertools.izip(flatten(mt0), flatten(mt1)):
        if a != b:
            if a+1 == b:
                diffs += 1
            else:
                return False
        if diffs > 1:
            return False
    return diffs == 1

    def module_generator(self, shape):
        """
        This yields the module generator (or highest weight element) of a classical
        crystal of given shape. The module generator is the unique tableau with equal
        shape and content.

        EXAMPLES::

            sage: T = crystals.Tableaux(['D',3], shape = [1,1])
            sage: T.module_generator([1,1])
            [[1], [2]]

            sage: T = crystals.Tableaux(['D',4],shape=[2,2,2,-2])
            sage: T.module_generator(tuple([2,2,2,-2]))
            [[1, 1], [2, 2], [3, 3], [-4, -4]]
            sage: T.cardinality()
            294
            sage: T = crystals.Tableaux(['D',4],shape=[2,2,2,2])
            sage: T.module_generator(tuple([2,2,2,2]))
            [[1, 1], [2, 2], [3, 3], [4, 4]]
            sage: T.cardinality()
            294
        """
        type = self.cartan_type()
        if type[0] == 'D' and len(shape) == type[1] and shape[type[1]-1] < 0:
            invert = True
            shape = shape[:-1] + (-shape[type[1]-1],)
        else:
            invert = False
        p = Partition(shape).conjugate()
        # The column canonical tableau, read by columns
        module_generator = flatten([[val-i for i in range(val)] for val in p])
        if invert:
            module_generator = [(-x if x == type[1] else x) for x in module_generator]
        return self(list=[self.letters(x) for x in module_generator])

def _prepare_coefficient_by_restriction(precision, weight_parity, relation_precision, S) :
    r"""
    Provide input data to ``_coefficient_by_restriction__with_restriction_matrix``.

    INPUT:
    
    - ``precision`` -- A filter for Jacobi forms of arbitrary index.
    
    - ``weight_parity`` -- An integer.
    
    - ``relation_precision`` -- A filter for Jacobi forms.

    - ``S`` -- A list of vectors.
    """
    L = precision.jacobi_index()

    rand = Random()
    max_S_length = max([L(s) for s in S])
    relation_S_pre = flatten(L.short_vector_list_up_to_length(max_S_length + 1, True)[1:])
    relation_S = list()
    for _ in range(4 * L.det()) :
        s = rand.choice(relation_S_pre)
        if s not in relation_S :
            relation_S.append(s)

    (global_restriction_matrix__big, row_groups, row_labels, column_labels) = _global_restriction_matrix(precision, S, weight_parity)
    (global_relation_matrix, column_labels_relations) = _global_relation_matrix(relation_precision, relation_S, weight_parity )
    global_restriction_matrix__big.change_ring(QQ)
    global_relation_matrix.change_ring(QQ)

    return ( global_restriction_matrix__big, row_groups, row_labels, column_labels,
             global_relation_matrix, column_labels_relations )

    def YoungsLattice(n):
        """
        Return Young's Lattice up to rank `n`.

        In other words, the poset of partitions
        of size less than or equal to `n` ordered by inclusion.

        INPUT:

        - ``n`` -- a positive integer

        EXAMPLES::

            sage: P = Posets.YoungsLattice(3); P
            Finite meet-semilattice containing 7 elements
            sage: P.cover_relations()
            [[[], [1]],
             [[1], [1, 1]],
             [[1], [2]],
             [[1, 1], [1, 1, 1]],
             [[1, 1], [2, 1]],
             [[2], [2, 1]],
             [[2], [3]]]
        """
        from sage.combinat.partition import Partitions, Partition
        from sage.misc.flatten import flatten
        partitions = flatten([list(Partitions(i)) for i in range(n + 1)])
        return JoinSemilattice((partitions, Partition.contains)).dual()

     def bruhat_interval(self, x, y):
         """
         Returns the list of t such that x <= t <= y.
 
         EXAMPLES::
 
             sage: W = WeylGroup("A3", prefix="s")
             sage: [s1,s2,s3]=W.simple_reflections()
             sage: W.bruhat_interval(s2,s1*s3*s2*s1*s3)
             [s1*s2*s3*s2*s1, s2*s3*s2*s1, s3*s1*s2*s1, s1*s2*s3*s1, s1*s2*s3*s2, s3*s2*s1, s2*s3*s1, s2*s3*s2, s1*s2*s1, s3*s1*s2, s1*s2*s3, s2*s1, s3*s2, s2*s3, s1*s2, s2]
             sage: W = WeylGroup(['A',2,1], prefix="s")
             sage: [s0,s1,s2]=W.simple_reflections()
             sage: W.bruhat_interval(1,s0*s1*s2)
             [s0*s1*s2, s1*s2, s0*s2, s0*s1, s2, s1, s0, 1]
         """
         if x == 1:
             x = self.one()
         if y == 1:
             y = self.one()
         if x == y:
             return [x]
         ret = []
         if not x.bruhat_le(y):
             return ret
         ret.append([y])
         while ret[-1] != []:
             nextlayer = []
             for z in ret[-1]:
                 for t in z.bruhat_lower_covers():
                     if t not in nextlayer:
                         if x.bruhat_le(t):
                             nextlayer.append(t)
             ret.append(nextlayer)
         return flatten(ret)

    def module_generator(self, shape):
        """
        This yields the module generator (or highest weight element) of a classical
        crystal of given shape. The module generator is the unique tableau with equal
        shape and content.

        EXAMPLE::

            sage: T = CrystalOfTableaux(['D',3], shape = [1,1])
            sage: T.module_generator([1,1])
            [[1], [2]]
        """
        type = self.cartan_type()
        if type[0] == 'D' and len(shape) == type[1] and shape[type[1]-1] < 0:
            invert = True
            shape = shape[:-1]+(-shape[type[1]-1],)
        else:
            invert = False
        p = Partition(shape).conjugate()
        # The column canonical tableau, read by columns
        module_generator = flatten([[p[j]-i for i in range(p[j])] for j in range(len(p))])
        if invert:
            for i in range(type[1]):
                if module_generator[i] == type[1]:
                    module_generator[i] = -type[1]
        return self(list=[self.letters(x) for x in module_generator])

    def polish_notation(self):
        r"""
        Convert the calling boolean formula into polish notation.

        OUTPUT:

        A string representation of the formula in polish notation.

        EXAMPLES:

        This example illustrates converting a formula to polish notation::

            sage: import sage.logic.propcalc as propcalc
            sage: f = propcalc.formula("~~a|(c->b)")
            sage: f.polish_notation()
            '|~~a->cb'

            sage: g = propcalc.formula("(a|~b)->c")
            sage: g.polish_notation()
            '->|a~bc'

        AUTHORS:

        - Paul Scurek (2013-08-03)
        """
        return ''.join(flatten(logicparser.polish_parse(repr(self))))

    def minimal_composition_filter(self, ls, rs) :
        if len(ls) == 0 or len(rs) == 0 :
            return SiegelModularFormGnFilter_diagonal_lll(self.__n, 0, self.__reduced)

        maxd = flatten( map(lambda (ml, mr): [ml[i,i] + mr[i,i] for i in xrange(self.__n)],
                                itertools.product(ls, rs) ) )

        return SiegelModularFormGnFilter_diagonal_lll(self.__n, maxd + 1, self.__reduced)

    def plot_cluster_fan_stereographically(self, northsign=1, north=None, right=None, colors=None):
        from sage.plot.graphics import Graphics
        from sage.plot.point import point
        from sage.misc.flatten import flatten
        from sage.plot.line import line
        from sage.misc.functional import norm

        if self.rk !=3:
            raise ValueError("Can only stereographically project fans in 3d.")
        if not self.is_finite() and self._depth == infinity:
            raise ValueError("For infinite algebras you must specify the depth.")

        if north == None:
            if self.is_affine():
                north = vector(self.delta())
            else:
                north = vector( (-1,-1,-1) )
        if right == None:
            if self.is_affine():
                right = vector(self.gamma())
            else:
                right = vector( (1,0,0) )
        if colors == None:
            colors = dict([(0,'red'),(1,'green'),(2,'blue'),(3,'cyan'),(4,'yellow')])
        G = Graphics()

        roots = list(self.g_vectors())
        compatible = []
        while roots:
            x = roots.pop()
            for y in roots:
                if self.compatibility_degree(x,y) == 0:
                    compatible.append((x,y))
        for (u,v) in compatible:
            G += _stereo_arc(vector(u),vector(v),vector(u+v),north=northsign*north,right=right,thickness=0.5,color='black')

        for i in range(3):
            orbit = self.ith_orbit(i)
            for j in orbit:
                G += point(_stereo_coordinates(vector(orbit[j]),north=northsign*north,right=right),color=colors[i],zorder=len(G))

        if self.is_affine():
            tube_vectors = map(vector,flatten(self.affine_tubes()))
            for v in tube_vectors:
                G += point(_stereo_coordinates(v,north=northsign*north,right=right),color=colors[3],zorder=len(G))
            if north != vector(self.delta()):
                G += _stereo_arc(tube_vectors[0],tube_vectors[1],vector(self.delta()),north=northsign*north,right=right,thickness=2,color=colors[4],zorder=0)
            else:
                # FIXME: refactor this before publishing
                tube_projections = [
                        _stereo_coordinates(v,north=northsign*north,right=right)
                        for v in tube_vectors ]
                t=min((G.get_minmax_data()['xmax'],G.get_minmax_data()['ymax']))
                G += line([tube_projections[0],tube_projections[0]+t*(_normalize(tube_projections[0]-tube_projections[1]))],thickness=2,color=colors[4],zorder=0)
                G += line([tube_projections[1],tube_projections[1]+t*(_normalize(tube_projections[1]-tube_projections[0]))],thickness=2,color=colors[4],zorder=0)
        G.set_aspect_ratio(1)
        G._show_axes = False
        return G

 def contained_partitions(l):
     """
     Nested function returning those partitions contained in
     the partition `l`
     """
     if l == Partition([]):
         return l
     return flatten([l, [contained_partitions(m)
                         for m in lower_covers(l)]])

    def _call_(self, x):
        r"""
        Return the image of ``x`` in the tableau model of `B(\infty)`.

        EXAMPLES::

            sage: T = crystals.infinity.Tableaux(['A',3])
            sage: RC = crystals.infinity.RiggedConfigurations(['A',3])
            sage: phi = T.coerce_map_from(RC)
            sage: x = RC.an_element().f_string([2,2,1,1,3,2,1,2,1,3])
            sage: y = phi(x); y.pp()
              1  1  1  1  1  2  2  3  4
              2  2  3  4
              3
            sage: (~phi)(y) == x
            True
        """
        lam = [sum(nu)+1 for nu in x]
        ct = self.domain().cartan_type()
        I = ct.index_set()
        if ct.type() == 'D':
            lam[-2] = max(lam[-2], lam[-1])
            lam.pop()
            l = sum([ [[r+1,1]]*v for r,v in enumerate(lam[:-1]) ], [])
            n = len(I)
            l = l + sum([ [[n,1], [n-1,1]] for k in range(lam[-1])], [])
        else:
            if ct.type() == 'B':
                lam[-1] *= 2
            l = sum([ [[r,1]]*lam[i] for i,r in enumerate(I) ], [])

        RC = RiggedConfigurations(ct.affine(), reversed(l))
        elt = RC(x)
        if ct.type() == 'A':
            bij = RCToKRTBijectionTypeA(elt)
        elif ct.type() == 'B':
            bij = RCToMLTBijectionTypeB(elt)
        elif ct.type() == 'C':
            bij = RCToKRTBijectionTypeC(elt)
        elif ct.type() == 'D':
            bij = RCToMLTBijectionTypeD(elt)
        else:
            raise NotImplementedError("bijection of type {} not yet implemented".format(ct))
        y = bij.run()

        # Now make the result marginally large
        y = [list(c) for c in y]
        cur = []
        L = CrystalOfLetters(ct)
        for i in I:
            cur.insert(0, L(i))
            c = y.count(cur)
            while c > 1:
                y.remove(cur)
                c -= 1
        return self.codomain()(*flatten(y))

    def loops_iterator(self, other=None):
        r"""
        INPUT:

             - ``other`` -- a perfect matching of the same set of ``self``.
               (if the second argument is empty, the method :meth:`an_element` is
               called on the parent of the first)

        OUTPUT:

            If we draw the two perfect matchings simultaneously as edges of a
            graph, the graph obtained is a union of cycles of even lengths.
            The function returns an iterator for these cycles (each cycle is
            given as a list).

        EXAMPLES::

            sage: o = PerfectMatching([(1, 7), (2, 4), (3, 8), (5, 6)])
            sage: p = PerfectMatching([(1, 6), (2, 7), (3, 4), (5, 8)])
            sage: it = o.loops_iterator(p)
            sage: it.next()
            [1, 7, 2, 4, 3, 8, 5, 6]
            sage: it.next()
            Traceback (most recent call last):
            ...
            StopIteration
        """
        if other is None:
            other = self.parent().an_element()
        elif self.parent() != other.parent():
            s = "%s is not a matching of the ground set of %s" % (other, self)
            raise ValueError(s)
        remain = flatten(self.value)
        while len(remain) > 0:
            a = remain.pop(0)
            b = self.partner(a)
            remain.remove(b)
            loop = [a, b]
            c = other.partner(b)
            while c != a:
                b = self.partner(c)
                remain.remove(c)
                loop.append(c)
                remain.remove(b)
                loop.append(b)
                c = other.partner(b)
            yield loop

    def module_generator(self):
        """
        Return the module generator (or highest weight element) of ``self``.

        The module generator is the unique tableau of shape `(n-1, \ldots, 2,
        1)` with weight `0`.

        EXAMPLES::

            sage: T = crystals.infinity.Tableaux(['D',4])
            sage: T.module_generator()
            [[1, 1, 1], [2, 2], [3]]
        """
        n = self._cartan_type.rank()
        p = Partition([x for x in reversed(range(1, n))])
        # The column canonical tableau, read by columns
        module_generator = flatten([[p[j]-i for i in range(p[j])] for j in range(n-1)])
        return self(list=[self.letters(x) for x in module_generator])

    def _build_module_generators(self):
        r"""
        Build the module generators.

        There is only one module generator which corresponds to a single
        `r \times s` rectangle.

        EXAMPLES::

            sage: KRT = KirillovReshetikhinTableaux(['A', 4, 1], 2, 3)
            sage: KRT._build_module_generators()
            ([[1, 1, 1], [2, 2, 2]],)
        """
        tableau = []
        for i in range(self._s):
            tableau.append( [self._r - j for j in range(self._r)] )

        return (self([self.letters(x) for x in flatten(tableau)]),)

def generateCode (f, pars):
	'''
	Creates list1 and list 2 from the right side function f of an ODE:
	input: 
		f -> right side function for ODE system
		pars -> list with the parameters on f
	output:
		C code for Automatic Differentiation


	Example with Lorenz Equation
	sage: var('t, x, y, z')		# variables for lorenz equations
	sage: var('s, r, b')		# parameters for lorenz equations
	sage: f(t,x,y,z) = [s*(y-x), x*(r-z) - y, x*y - b*z]	# Right side function for Lorenz equation
	sage: generateCode (f, [s, r, b])

	'''
	list1, list2 = createLists (f, pars)
	removeRepeated (list1, list2)
	constList1, constList2 = removeConstants (list1, list2, pars)
	list3 = createCodeList(list1, list2, constList1, f, pars)
	constList3 = createConstCodeList (constList1, constList2, pars)

	parsedConstList = createParsedConstList (constList3)
	parsedList = createParsedList (list1, list3, f)
	vars = f[0].arguments ()
	auxSet = set(flatten([i.variables() for i in f]))	# set of variables in f
	if set ([vars[0]]).issubset (auxSet):				# non autonomous system
		print '\tdouble T[order];'
		print '\tfor (i=2; i<order; i++) T[i] = 0.0;'
		print '\tT[0] = t;'
		print '\tT[1] = 1.0;'

	if len(list1) > 0:	# checks if there are links
		print '\tdouble l[{}][order];'.format (len (list1))
	if len(constList1) > 0: # checks if there are constant expresions
		print '\tdouble c[{}];'.format (len (constList1))
	for s in parsedConstList:
        	print s
	print '\tfor (i=0; i<order; i++) {'
	for s in parsedList:
		print s
	print '\t}'

    def plot3d(self,depth=None):
        # FIXME: refactor this before publishing
        from sage.plot.graphics import Graphics
        from sage.plot.point import point
        from sage.misc.flatten import flatten
        from sage.plot.plot3d.shapes2 import sphere
        if self._n !=3:
            raise ValueError("Can only 3d plot fans.")
        if depth == None:
            depth = self._depth
        if not self.is_finite() and depth==infinity:
            raise ValueError("For infinite algebras you must specify the depth.")

        colors = dict([(0,'red'),(1,'green'),(2,'blue'),(3,'cyan')])
        G = Graphics()

        roots = self.d_vectors(depth=depth)
        compatible = []
        while roots:
            x = roots.pop()
            for y in roots:
                if self.compatibility_degree(x,y) == 0:
                    compatible.append((x,y))
        for (u,v) in compatible:
            G += _arc3d((_normalize(vector(u)),_normalize(vector(v))),thickness=0.5,color='black')

        for i in range(3):
            orbit = self.ith_orbit(i,depth=depth)
            for j in orbit:
                G += point(_normalize(vector(orbit[j])),color=colors[i],size=10,zorder=len(G.all))

        if self.is_affine():
            tube_vectors=map(vector,flatten(self.affine_tubes()))
            tube_vectors=map(_normalize,tube_vectors)
            for v in tube_vectors:
                G += point(v,color=colors[3],size=10,zorder=len(G.all))
            G += _arc3d((tube_vectors[0],tube_vectors[1]),thickness=5,color='gray',zorder=0)
        
        G += sphere((0,0,0),opacity=0.1,zorder=0)
        G._extra_kwds['frame']=False
        G._extra_kwds['aspect_ratio']=1 
        return G

 def _tube_support(self, alpha):
     gck = self.gamma().associated_coroot()
     if gck.scalar(alpha) != 0:
         raise ValueError("Root not in U_c")
     tubes = list(self.affine_tubes())
     ack = alpha.associated_coroot()
     while tubes:
         tube = tubes.pop()
         if any([ack.scalar(x) != 0 for x in tube[0]]):
             sup_a = []
             roots = flatten(tube)+[0]
             basis = tube[0]
             a = copy(alpha)
             while a != 0:
                 edges = [ x for x in basis if a-x in roots ]
                 sup_a += edges
                 while edges:
                     a = a-edges.pop()
             return tuple(sup_a)
     raise ValueError("Unable to compute support of root")

    def nestings_iterator(self):
        r"""
        INPUT:

            A perfect matching on a *totally ordered* ground set.

        OUTPUT:

            We place the element of a ground set and draw the perfect matching
            by linking the elements of the same pair in the upper
            half-plane. This function returns an iterator over the pairs of
            nesting lines (as a line correspond to a pair, the iterator
            produces pairs of pairs).

        EXAMPLES::

            sage: n = PerfectMatching([(1, 6), (2, 7), (3, 5), (4, 8)])
            sage: it = n.nestings_iterator();
            sage: it.next()
            ((1, 6), (3, 5))
            sage: it.next()
            ((2, 7), (3, 5))
            sage: it.next()
            Traceback (most recent call last):
            ...
            StopIteration
        """
        x = self.value[:]
        if len(x) == 0:
            return
        (i, j) = x.pop(0)
        for (a, b) in x:
            # if (i<a<j<b) or (i<b<j<a) or (j<a<i<b) or (j<b<i<a) or (
            #        a<i<b<j) or (a<j<b<i) or (b<i<a<j) or (b<j<a<i):
            labij = sorted([a, b, i, j])
            posij = sorted([labij.index(i), labij.index(j)])
            if posij == [0, 3] or posij == [1, 2]:
                yield ((i, j), (a, b))
        for nest in PerfectMatchings(flatten(x))(x).nestings_iterator():
            yield nest

    def face_poset(self):
        r"""
        The face poset of this cell complex, the poset of
        nonempty cells, ordered by inclusion.

        This uses the :meth:`cells` method, and also assumes that for
        each cell ``f``, all of ``f.faces()``, ``tuple(f)``, and
        ``f.dimension()`` make sense.  (If this is not the case in
        some derived class, as happens with `\Delta`-complexes, then
        override this method.)

        EXAMPLES::

            sage: P = SimplicialComplex([[0, 1], [1,2], [2,3]]).face_poset(); P
            Finite poset containing 7 elements
            sage: P.list()
            [(3,), (2,), (2, 3), (1,), (0,), (0, 1), (1, 2)]

            sage: S2 = cubical_complexes.Sphere(2)
            sage: S2.face_poset()
            Finite poset containing 26 elements
        """
        from sage.combinat.posets.posets import Poset
        from sage.misc.flatten import flatten

        covers = {}
        # The code for posets seems to work better if each cell is
        # converted to a tuple.
        all_cells = flatten([list(f) for f in self.cells().values()])

        for C in all_cells:
            if C.dimension() >= 0:
                covers[tuple(C)] = []
        for C in all_cells:
            for face in C.faces():
                if face.dimension() >= 0:
                    covers[tuple(face)].append(tuple(C))
        return Poset(covers)

    def __classcall_private__(cls, p):
        r"""
        This function tries to recognize the input (it can be either a list or
        a tuple of pairs, or a fix-point free involution given as a list or as
        a permutation), constructs the parent (enumerated set of
        PerfectMatchings of the ground set) and calls the __init__ function to
        construct our object.

        EXAMPLES::

            sage: m = PerfectMatching([('a','e'),('b','c'),('d','f')]);m
            [('a', 'e'), ('b', 'c'), ('d', 'f')]
            sage: isinstance(m,PerfectMatching)
            True
            sage: n = PerfectMatching([3, 8, 1, 7, 6, 5, 4, 2]);n
            [(1, 3), (2, 8), (4, 7), (5, 6)]
            sage: n.parent()
            Set of perfect matchings of {1, 2, 3, 4, 5, 6, 7, 8}
            sage: PerfectMatching([(1, 4), (2, 3), (5, 6)]).is_non_crossing()
            True

        The function checks that the given list or permutation is a valid perfect
        matching (i.e. a list of pairs with pairwise disjoint elements  or a
        fixpoint-free involution) and raises a ValueError otherwise:

            sage: PerfectMatching([(1, 2, 3), (4, 5)])
            Traceback (most recent call last):
            ...
            ValueError: [(1, 2, 3), (4, 5)] is not a valid perfect matching: all elements of the list must be pairs

        If you know your datas are in a good format, use directly
        ``PerfectMatchings(objects)(data)``.

        TESTS::

             sage: m = PerfectMatching([('a','e'),('b','c'),('d','f')])
             sage: TestSuite(m).run()
             sage: m = PerfectMatching([])
             sage: TestSuite(m).run()
             sage: PerfectMatching(6)
             Traceback (most recent call last):
             ...
             ValueError: cannot convert p (= 6) to a PerfectMatching
             sage: PerfectMatching([(1,2,3)])
             Traceback (most recent call last):
             ...
             ValueError: [(1, 2, 3)] is not a valid perfect matching:
             all elements of the list must be pairs
             sage: PerfectMatching([(1,1)])
             Traceback (most recent call last):
             ...
             ValueError: [(1, 1)] is not a valid perfect matching:
             there are some repetitions
             sage: PerfectMatching(Permutation([4,2,1,3]))
             Traceback (most recent call last):
             ...
             ValueError: The permutation p (= [4, 2, 1, 3]) is not a fixed point free involution
        """
        # we have to extract from the argument p the set of objects of the
        # matching and the list of pairs.
        # First case: p is a list (resp tuple) of lists (resp tuple).
        if (isinstance(p, list) or isinstance(p, tuple)) and (
                all([isinstance(x, list) or isinstance(x, tuple) for x in p])):
            objects = Set(flatten(p))
            data = (map(tuple, p))
            #check if the data are correct
            if not all([len(t) == 2 for t in data]):
                raise ValueError("%s is not a valid perfect matching:\n"
                                 "all elements of the list must be pairs" % p)
            if len(objects) < 2*len(data):
                raise ValueError("%s is not a valid perfect matching:\n"
                                 "there are some repetitions" % p)
        # Second case: p is a permutation or a list of integers, we have to
        # check if it is a fix-point-free involution.
        elif ((isinstance(p, list) and
               all(map(lambda x: (isinstance(x, Integer) or isinstance(x, int)), p)))
              or isinstance(p, Permutation)):
            p = Permutation(p)
            n = len(p)
            if not(p.cycle_type() == [2 for i in range(n//2)]):
                raise ValueError("The permutation p (= %s) is not a "
                                 "fixed point free involution" % p)
            objects = Set(range(1, n+1))
            data = p.to_cycles()
        # Third case: p is already a perfect matching, we return p directly
        elif isinstance(p, PerfectMatching):
            return p
        else:
            raise ValueError("cannot convert p (= %s) to a PerfectMatching" % p)
        # Finally, we create the parent and the element using the element
        # class of the parent. Note: as this function is private, when we
        # create an object via parent.element_class(...), __init__ is directly
        # executed and we do not have an infinite loop.
        return PerfectMatchings(objects)(data)