def cardinality(self):
        """
        Returns the number of Lyndon words with the evaluation e.

        EXAMPLES::

            sage: LyndonWords([]).cardinality()
            0
            sage: LyndonWords([2,2]).cardinality()
            1
            sage: LyndonWords([2,3,2]).cardinality()
            30

        Check to make sure that the count matches up with the number of
        Lyndon words generated.

        ::

            sage: comps = [[],[2,2],[3,2,7],[4,2]]+Compositions(4).list()
            sage: lws = [ LyndonWords(comp) for comp in comps]
            sage: all( [ lw.cardinality() == len(lw.list()) for lw in lws] )
            True
        """
        evaluation = self.e
        le = __builtin__.list(evaluation)
        if len(evaluation) == 0:
            return 0

        n = sum(evaluation)

        return sum([moebius(j)*factorial(n/j) / prod([factorial(ni/j) for ni in evaluation]) for j in divisors(gcd(le))])/n

    def cardinality(self):
        """
        Returns the number of integer necklaces with the evaluation e.
        
        EXAMPLES::
        
            sage: Necklaces([]).cardinality()
            0
            sage: Necklaces([2,2]).cardinality()
            2
            sage: Necklaces([2,3,2]).cardinality()
            30
        
        Check to make sure that the count matches up with the number of
        Lyndon words generated.
        
        ::
        
            sage: comps = [[],[2,2],[3,2,7],[4,2]]+Compositions(4).list()
            sage: ns = [ Necklaces(comp) for comp in comps]
            sage: all( [ n.cardinality() == len(n.list()) for n in ns] )
            True
        """
        evaluation = self.e
        le = list(evaluation)
        if len(le) == 0:
            return 0

        n = sum(le)

        return sum([euler_phi(j)*factorial(n/j) / prod([factorial(ni/j) for ni in evaluation]) for j in divisors(gcd(le))])/n

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

        The number of ordered set partitions of a set of length `k` with
        composition shape `\mu` is equal to

        .. MATH::

            \frac{k!}{\prod_{\mu_i \neq 0} \mu_i!}.

        EXAMPLES::

            sage: OrderedSetPartitions(5,[2,3]).cardinality()
            10
            sage: OrderedSetPartitions(0, []).cardinality()
            1
            sage: OrderedSetPartitions(0, [0]).cardinality()
            1
            sage: OrderedSetPartitions(0, [0,0]).cardinality()
            1
            sage: OrderedSetPartitions(5, [2,0,3]).cardinality()
            10
        """
        return factorial(len(self._set))/prod([factorial(i) for i in self.c])

    def zero_eigenvectors(self, tg, flags):

        if self._use_symmetry:
            return self.symm_zero_eigenvectors(tg, flags)

        cn = self._graph.n
        s = tg.n
        k = flags[0].n  # assume all flags the same order

        rows = []

        for tv in Tuples(range(1, cn + 1), s):

            it = self._graph.degenerate_induced_subgraph(tv)

            using_phantom_edge = False
            phantom_edge = None

            if hasattr(self, "_phantom_edge") and it.ne == tg.ne - 1:
                extra_edges = [e for e in tg if not e in it]
                if len(extra_edges) == 1:
                    phantom_edge = extra_edges[0]
                    if all(tv[phantom_edge[i] - 1] == self._phantom_edge[i] for i in range(tg.r)):
                        it.add_edge(phantom_edge)
                        using_phantom_edge = True

            if not (using_phantom_edge or it.is_labelled_isomorphic(tg)):
                continue

            total = Integer(0)
            row = [0] * len(flags)

            for ov in UnorderedTuples(range(1, cn + 1), k - s):

                factor = factorial(k - s)
                for i in range(1, cn + 1):
                    factor /= factorial(ov.count(i))

                if self._weights:
                    for v in ov:
                        factor *= self._weights[v - 1]

                ig = self._graph.degenerate_induced_subgraph(tv + ov)
                if using_phantom_edge:
                    ig.add_edge(phantom_edge)
                ig.t = s
                ig.make_minimal_isomorph()

                for j in range(len(flags)):
                    if ig.is_labelled_isomorphic(flags[j]):
                        row[j] += factor
                        total += factor
                        break

            for j in range(len(flags)):
                row[j] /= total
            rows.append(row)

        return matrix_of_independent_rows(self._field, rows, len(flags))

    def cardinality(self):
        r"""
        Returns the number of subwords of w of length k.

        EXAMPLES::

            sage: Subwords([1,2,3], 2).cardinality()
            3
        """
        w = self.w
        k = self.k
        return factorial(len(w))/(factorial(k)*factorial(len(w)-k))

def volume_hamming(n,q,r):
    r"""
    Returns number of elements in a Hamming ball of radius r in `\GF{q}^n`.
    Agrees with Guava's SphereContent(n,r,GF(q)).

    EXAMPLES::
    
        sage: volume_hamming(10,2,3)
        176
    """
    ans=sum([factorial(n)/(factorial(i)*factorial(n-i))*(q-1)**i for i in range(r+1)])
    return ans

    def by_taylor_expansion(self, fs, k, is_integral=False) :
        r"""
        We combine the theta decomposition and the heat operator as in [Sko].
        This yields a bijections of Jacobi forms of weight `k` and
        `M_k \times S_{k+2} \times .. \times S_{k+2m}`.
        """
        ## we introduce an abbreviations
        if is_integral :
            PS = self.integral_power_series_ring()
        else :
            PS = self.power_series_ring()
            
        if not len(fs) == self.__precision.jacobi_index() + 1 :
            raise ValueError("fs must be a list of m + 1 elliptic modular forms or their fourier expansion")
        
        qexp_prec = self._qexp_precision()
        if qexp_prec is None : # there are no forms below the precision
            return dict()
        
        f_divs = dict()
        for (i, f) in enumerate(fs) :
            f_divs[(i, 0)] = PS(f(qexp_prec), qexp_prec)
        
        if self.__precision.jacobi_index() == 1 :
            return self._by_taylor_expansion_m1(f_divs, k, is_integral)
        
        for i in xrange(self.__precision.jacobi_index() + 1) :
            for j in xrange(1, self.__precision.jacobi_index() - i + 1) :
                f_divs[(i,j)] = f_divs[(i, j - 1)].derivative().shift(1)
            
        phi_divs = list()
        for i in xrange(self.__precision.jacobi_index() + 1) :
            ## This is the formula in Skoruppas thesis. He uses d/ d tau instead of d / dz which yields
            ## a factor 4 m
            phi_divs.append( sum( f_divs[(j, i - j)] * (4 * self.__precision.jacobi_index())**i
                                  * binomial(i,j) / 2**i#2**(self.__precision.jacobi_index() - i + 1)
                                  * prod(2*(i - l) + 1 for l in xrange(1, i))
                                  / factorial(i + k + j - 1)
                                  * factorial(2*self.__precision.jacobi_index() + k - 1) 
                                  for j in xrange(i + 1) ) )
            
        phi_coeffs = dict()
        for r in xrange(self.__precision.jacobi_index() + 1) :
            series = sum( map(operator.mul, self._theta_factors()[r], phi_divs) )
            series = self._eta_factor() * series

            for n in xrange(qexp_prec) :
                phi_coeffs[(n, r)] = series[n]

        return phi_coeffs

    def cardinality(self):
        """
        Returns the number of permutations of k things from a list of n
        things.

        EXAMPLES::

            sage: from sage.combinat.permutation_nk import PermutationsNK
            sage: PermutationsNK(3,2).cardinality()
            6
            sage: PermutationsNK(5,4).cardinality()
            120
        """
        n, k = self._n, self._k
        return factorial(n)/factorial(n-k)

    def cardinality(self):
        """
        EXAMPLES::

            sage: OrderedSetPartitions(5,[2,3]).cardinality()
            10
            sage: OrderedSetPartitions(0, []).cardinality()
            1
            sage: OrderedSetPartitions(0, [0]).cardinality()
            1
            sage: OrderedSetPartitions(0, [0,0]).cardinality()
            1
            sage: OrderedSetPartitions(5, [2,0,3]).cardinality()
            10
        """
        return factorial(len(self.s))/prod([factorial(i) for i in self.c])

    def subgraph_densities(self, n):

        if n < 0:
            raise ValueError

        if self._use_symmetry:
            return self.symm_subgraph_densities(n)

        cn = self._graph.n
        total = Integer(0)
        sharp_graph_counts = {}
        sharp_graphs = []

        for P in UnorderedTuples(range(1, cn + 1), n):

            factor = factorial(n)
            for i in range(1, cn + 1):
                factor /= factorial(P.count(i))

            if self._weights:
                for v in P:
                    factor *= self._weights[v - 1]

            ig = self._graph.degenerate_induced_subgraph(P)
            igc = copy(ig)  # copy for phantom edge
            ig.make_minimal_isomorph()

            ghash = hash(ig)
            if ghash in sharp_graph_counts:
                sharp_graph_counts[ghash] += factor
            else:
                sharp_graphs.append(ig)
                sharp_graph_counts[ghash] = factor

            total += factor

            if hasattr(self, "_phantom_edge") and all(x in P for x in self._phantom_edge):
                phantom_edge = [P.index(x) + 1 for x in self._phantom_edge]
                igc.add_edge(phantom_edge)
                igc.make_minimal_isomorph()

                ghash = hash(igc)
                if not ghash in sharp_graph_counts:
                    sharp_graphs.append(igc)
                    sharp_graph_counts[ghash] = Integer(0)

        return [(g, sharp_graph_counts[hash(g)] / total) for g in sharp_graphs]

 def _test__jacobi_corrected_taylor_expansions(nu, phi, weight) :
     r"""
     Return the ``2 nu``-th corrected Taylor coefficient.
     INPUT:
     
     - ``nu`` -- An integer.  
     
     - ``phi`` -- A Fourier expansion of a Jacobi form.
     
     - ``weight`` -- An integer.
     
     OUTPUT:
     
     - A power series in `q`.
     
     ..TODO:
     
     Implement this for all Taylor coefficients.
     
     TESTS::
     
         sage: from psage.modform.jacobiforms.jacobiformd1nn_fegenerators import *
         sage: from psage.modform.jacobiforms.jacobiformd1nn_types import *
         sage: nu_bound = 10
         sage: precision = 100
         sage: weight = 10; index = 7
         sage: phis = [jacobi_form_by_taylor_expansion(i, JacobiFormD1NNFilter(precision, index), weight) for i in range(JacobiFormD1NN_Gamma(weight, index)._rank(QQ))]
         sage: fss = [ [JacobiFormD1NNFactory_class._test__jacobi_corrected_taylor_expansions(nu, phi, weight) for phi in phis] for nu in range(nu_bound) ]
         sage: fss_vec = [ [vector(f.padded_list(precision)) for f in fs] for fs in fss ]
         sage: mf_spans = [ span([vector(b.qexp(precision).padded_list(precision)) for b in ModularForms(1, weight + 2 * nu).basis()]) for nu in range(nu_bound) ] 
         sage: all(f_vec in mf_span for (fs_vec, mf_span) in zip(fss_vec, mf_spans) for f_vec in fs_vec)
         True
     """
     ## We use EZ85, p.29 (3), the factorial in one of the factors is missing
     factors = [ (-1)**mu * factorial(2 * nu) * factorial(weight + 2 * nu - mu - 2) / ZZ(factorial(mu) * factorial(2 * nu - 2 * mu) * factorial(weight + nu - 2))
                 for mu in range(nu + 1) ]
     gegenbauer = lambda n, r: sum( f * r**(2 * nu - 2 * mu) * n**mu 
                                    for (mu,f) in enumerate(factors) )
     ch = JacobiFormD1WeightCharacter(weight)
     jacobi_index = phi.precision().jacobi_index()
     
     coeffs = dict( (n, QQ(0)) for n in range(phi.precision().index()) )
     for (n, r) in phi.precision().monoid_filter() :
         coeffs[n] += gegenbauer(jacobi_index * n, r) * phi[(ch, (n,r))]
     
     return PowerSeriesRing(QQ, 'q')(coeffs)

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

        The number of fully packed loops on  `n \times n` grid

        .. MATH::

            \prod_{k=0}^{n-1} \frac{(3k+1)!}{(n+k)!} = \frac{1! 4! 7! 10!
            \cdots (3n-2)!}{n! (n+1)! (n+2)! (n+3)! \cdots (2n-1)!}.

        EXAMPLES::

            sage: [AlternatingSignMatrices(n).cardinality() for n in range(0, 11)]
            [1, 1, 2, 7, 42, 429, 7436, 218348, 10850216, 911835460, 129534272700]
        """
        return Integer(prod( [ factorial(3*k+1)/factorial(self._n+k)
                       for k in range(self._n)] ))

def gamma__exact(n):
    """
    Evaluates the exact value of the gamma function at an integer or
    half-integer argument.

    EXAMPLES::

        sage: gamma__exact(4)
        6
        sage: gamma__exact(3)
        2
        sage: gamma__exact(2)
        1
        sage: gamma__exact(1)
        1

        sage: gamma__exact(1/2)
        sqrt(pi)
        sage: gamma__exact(3/2)
        1/2*sqrt(pi)
        sage: gamma__exact(5/2)
        3/4*sqrt(pi)
        sage: gamma__exact(7/2)
        15/8*sqrt(pi)

        sage: gamma__exact(-1/2)
        -2*sqrt(pi)
        sage: gamma__exact(-3/2)
        4/3*sqrt(pi)
        sage: gamma__exact(-5/2)
        -8/15*sqrt(pi)
        sage: gamma__exact(-7/2)
        16/105*sqrt(pi)

    """
    from sage.all import sqrt
    ## SANITY CHECK
    if (not n in QQ) or (denominator(n) > 2):
        raise TypeError, "Oops!  You much give an integer or half-integer argument."

    if (denominator(n) == 1):
        if n <= 0:
            return infinity
        if n > 0:
            return factorial(n-1) 
    else:
        ans = QQ(1)
        while (n != QQ(1)/2):
            if (n<0): 
                ans *= QQ(1)/n
                n = n + 1
            elif (n>0):
                n = n - 1
                ans *= n 
            
        ans *= sqrt(pi)
        return ans

    def cardinality(self):
        r"""
        Return the number of integer necklaces with the evaluation ``content``.

        The formula for the number of necklaces of content `\alpha`
        a composition of `n` is:

        .. MATH::

            \sum_{d|gcd(\alpha)} \phi(d)
            \binom{n/d}{\alpha_1/d, \ldots, \alpha_\ell/d},

        where `\phi(d)` is the Euler `\phi` function.

        EXAMPLES::

            sage: Necklaces([]).cardinality()
            0
            sage: Necklaces([2,2]).cardinality()
            2
            sage: Necklaces([2,3,2]).cardinality()
            30
            sage: Necklaces([0,3,2]).cardinality()
            2

        Check to make sure that the count matches up with the number of
        necklace words generated.

        ::

            sage: comps = [[],[2,2],[3,2,7],[4,2],[0,4,2],[2,0,4]]+Compositions(4).list()
            sage: ns = [ Necklaces(comp) for comp in comps]
            sage: all( [ n.cardinality() == len(n.list()) for n in ns] )
            True
        """
        evaluation = self._content
        le = list(evaluation)
        if not le:
            return 0

        n = sum(le)

        return sum(euler_phi(j)*factorial(n/j) / prod(factorial(ni/j)
                    for ni in evaluation) for j in divisors(gcd(le))) / n

def quadratic_L_function__exact(n, d):
    r"""
    Returns the exact value of a quadratic twist of the Riemann Zeta function
    by `\chi_d(x) = \left(\frac{d}{x}\right)`.

    The input `n` must be a critical value.

    EXAMPLES::

        sage: quadratic_L_function__exact(1, -4)
        1/4*pi
        sage: quadratic_L_function__exact(-4, -4)
        5/2

    TESTS::

        sage: quadratic_L_function__exact(2, -4)
        Traceback (most recent call last):
        ...
        TypeError: n must be a critical value (i.e. odd > 0 or even <= 0)

    REFERENCES:

    - [Iwasawa]_, pp 16-17, Special values of `L(1-n, \chi)` and `L(n, \chi)`
    - [IreRos]_
    - [WashCyc]_
    """
    from sage.all import SR, sqrt

    if n <= 0:
        return QuadraticBernoulliNumber(1 - n, d) / (n - 1)
    elif n >= 1:
        # Compute the kind of critical values (p10)
        if kronecker_symbol(fundamental_discriminant(d), -1) == 1:
            delta = 0
        else:
            delta = 1

        # Compute the positive special values (p17)
        if (n - delta) % 2 == 0:
            f = abs(fundamental_discriminant(d))
            if delta == 0:
                GS = sqrt(f)
            else:
                GS = I * sqrt(f)
            ans = SR(ZZ(-1) ** (1 + (n - delta) / 2))
            ans *= (2 * pi / f) ** n
            ans *= GS  # Evaluate the Gauss sum here! =0
            ans *= 1 / (2 * I ** delta)
            ans *= QuadraticBernoulliNumber(n, d) / factorial(n)
            return ans
        else:
            if delta == 0:
                raise TypeError("n must be a critical value (i.e. even > 0 or odd < 0)")
            if delta == 1:
                raise TypeError("n must be a critical value (i.e. odd > 0 or even <= 0)")

def zeta__exact(n):
    r"""
    Returns the exact value of the Riemann Zeta function

    The argument must be a critical value, namely either positive even
    or negative odd.

    See for example [Iwasawa]_, p13, Special value of `\zeta(2k)`

    EXAMPLES:

    Let us test the accuracy for negative special values::

        sage: RR = RealField(100)
        sage: for i in range(1,10):
        ...       print "zeta(" + str(1-2*i) + "): ", RR(zeta__exact(1-2*i)) - zeta(RR(1-2*i))
        zeta(-1):  0.00000000000000000000000000000
        zeta(-3):  0.00000000000000000000000000000
        zeta(-5):  0.00000000000000000000000000000
        zeta(-7):  0.00000000000000000000000000000
        zeta(-9):  0.00000000000000000000000000000
        zeta(-11):  0.00000000000000000000000000000
        zeta(-13):  0.00000000000000000000000000000
        zeta(-15):  0.00000000000000000000000000000
        zeta(-17):  0.00000000000000000000000000000

    Let us test the accuracy for positive special values::

        sage: all(abs(RR(zeta__exact(2*i))-zeta(RR(2*i))) < 10**(-28) for i in range(1,10))
        True

    TESTS::

        sage: zeta__exact(5)
        Traceback (most recent call last):
        ...
        TypeError: n must be a critical value (i.e. even > 0 or odd < 0)

    REFERENCES:

    .. [Iwasawa] Iwasawa, *Lectures on p-adic L-functions*
    .. [IreRos] Ireland and Rosen, *A Classical Introduction to Modern Number Theory*
    .. [WashCyc] Washington, *Cyclotomic Fields*
    """
    if n < 0:
        return bernoulli(1-n)/(n-1)
    elif n > 1:
        if (n % 2 == 0):
            return ZZ(-1)**(n/2 + 1) * ZZ(2)**(n-1) * pi**n * bernoulli(n) / factorial(n)
        else:
            raise TypeError("n must be a critical value (i.e. even > 0 or odd < 0)")
    elif n==1:
        return infinity
    elif n==0:
        return -1/2

def zeta__exact(n):
    r"""
    Returns the exact value of the Riemann Zeta function

    References:

    - Iwasawa's "Lectures on p-adic L-functions", p13, "Special value of `\zeta(2k)`"
    - Ireland and Rosen's "A Classical Introduction to Modern Number Theory"
    - Washington's "Cyclotomic Fields"

    EXAMPLES::

        sage: ## Testing the accuracy of the negative special values
        sage: RR = RealField(100)
        sage: for i in range(1,10):
        ...       print "zeta(" + str(1-2*i) + "): ", RR(zeta__exact(1-2*i)) - zeta(RR(1-2*i))
        zeta(-1):  0.00000000000000000000000000000
        zeta(-3):  0.00000000000000000000000000000
        zeta(-5):  0.00000000000000000000000000000
        zeta(-7):  0.00000000000000000000000000000
        zeta(-9):  0.00000000000000000000000000000
        zeta(-11):  0.00000000000000000000000000000
        zeta(-13):  0.00000000000000000000000000000
        zeta(-15):  0.00000000000000000000000000000
        zeta(-17):  0.00000000000000000000000000000

        sage: RR = RealField(100)
        sage: for i in range(1,10):
        ...       print "zeta(" + str(1-2*i) + "): ", RR(zeta__exact(1-2*i)) - zeta(RR(1-2*i))
        zeta(-1):  0.00000000000000000000000000000
        zeta(-3):  0.00000000000000000000000000000
        zeta(-5):  0.00000000000000000000000000000
        zeta(-7):  0.00000000000000000000000000000
        zeta(-9):  0.00000000000000000000000000000
        zeta(-11):  0.00000000000000000000000000000
        zeta(-13):  0.00000000000000000000000000000
        zeta(-15):  0.00000000000000000000000000000
        zeta(-17):  0.00000000000000000000000000000

    """
    if n<0:
        k = 1-n
        return -bernoulli(k)/k
    elif n>1:
        if (n % 2 == 0):
            return ZZ(-1)**(n/2 + 1) * ZZ(2)**(n-1) * pi**n * bernoulli(n) / factorial(n)
        else:
            raise TypeError, "n must be a critical value! (I.e. even > 0 or odd < 0.)"
    elif n==1:
        return infinity
    elif n==0:
        return -1/2

def quadratic_L_function__exact(n, d):
    r"""
    Returns the exact value of a quadratic twist of the Riemann Zeta function
    by `\chi_d(x) = \left(\frac{d}{x}\right)`.  

    References:

    - Iwasawa's "Lectures on p-adic L-functions", p16-17, "Special values of
      `L(1-n, \chi)` and `L(n, \chi)`
    - Ireland and Rosen's "A Classical Introduction to Modern Number Theory"
    - Washington's "Cyclotomic Fields"

    EXAMPLES::

        sage: bool(quadratic_L_function__exact(1, -4) == pi/4)
        True

    """
    from sage.all import SR, sqrt
    if n<=0:
        k = 1-n
        return -QuadraticBernoulliNumber(k,d)/k
    elif n>=1:
        ## Compute the kind of critical values (p10)
        if kronecker_symbol(fundamental_discriminant(d), -1) == 1:
            delta = 0
        else:
            delta = 1        

        ## Compute the positive special values (p17)
        if ((n - delta) % 2 == 0):
            f = abs(fundamental_discriminant(d))
            if delta == 0:                            
                GS = sqrt(f)
            else:
                GS = I * sqrt(f)
            ans = SR(ZZ(-1)**(1+(n-delta)/2))
            ans *= (2*pi/f)**n
            ans *= GS     ## Evaluate the Gauss sum here! =0
            ans *= 1/(2 * I**delta)
            ans *= QuadraticBernoulliNumber(n,d)/factorial(n)
            return ans
        else:
            if delta == 0:
                raise TypeError, "n must be a critical value!\n" + "(I.e. even > 0 or odd < 0.)"
            if delta == 1:
                raise TypeError, "n must be a critical value!\n" + "(I.e. odd > 0 or even <= 0.)"

    def __init__(self,b,N,iprec,x0=0):
        self.bsym = b
        self.N = N
        self.iprec = iprec
        self.x0sym = x0
        self.prec = None


        bname = repr(b).strip('0').replace('.',',')
        if b == sqrt(2):
           bname = "sqrt2"
        if b == e**(1/e):
           bname = "eta"
           
        x0name = repr(x0)
        if x0name.find('.') > -1:
            if x0.is_real():
                x0name = repr(float(real(x0))).strip('0').replace('.',',')
            else:
                x0name = repr(complex(x0)).strip('0').replace('.',',')
        # by some reason save does not work with additional . inside the path

        self.path = "savings/msexp_%s"%bname + "_N%04d"%N + "_iprec%05d"%iprec + "_a%s"%x0name

        b = b.n(iprec)
        self.b = b
        x0 = x0.n(iprec)
        self.x0 = x0
        if isinstance(x0,RealNumber):
            R = RealField(iprec)
	else:
            R = ComplexField(iprec)    

        #symbolic carleman matrix
        if x0 == 0:
            #C = Matrix([[ m**n*log(b)**n/factorial(n) for n in range(N)] for m in range(N)])
            coeffs = [log(b)**n/factorial(n) for n in xrange(N)]
        else:
            #too slow
            #c = b**x0
            #C = Matrix([ [log(b)**n/factorial(n)*sum([binomial(m,k)*k**n*c**k*(-x0)**(m-k) for k in range(m+1)]) for n in range(N)] for m in range(N)])
            coeffs = [b**x0-x0]+[b**x0*log(b)**n/factorial(n) for n in xrange(1,N)]

        def psmul(A,B):
            N = len(B)
            return [sum([A[k]*B[n-k] for k in xrange(n+1)]) for n in xrange(N)]
        
        C = Matrix(R,N)
        row = vector(R,[1]+(N-1)*[0])
        C[0] = row
        for m in xrange(1,N):
            row = psmul(row,coeffs)
            C[m] = row
  
        print "Carleman matrix created."

        #numeric matrix and eigenvalues
        #self.CM = C.n(iprec) #n seems to reduce to a RealField
        self.CM = C

        self.eigenvalues = self.CM.eigenvalues()
        print "Eigenvalues computed."

        self.IM = None
        self.calc_IM()
        print "Iteration matrix computed."
        self.coeffs_1 = self.IM * vector([1]+[(1-x0)**n for n in range(1,N)])
        if x0 == 0:
            self.coeffs_0 = self.IM.column(0)
        else:
            self.coeffs_0 = self.IM * vector([1]+[(-x0)**n for n in range(1,N)])
        #there would also be a method to only coeffs_0 with x0,
        #this would require to make the offset -1-slog(0)
       
        self.L = None

 def f(n):
     return diff(expr,var,n).substitute({var:at})/factorial(n)