def test_arithmetic_sums():
    assert summation(1, (n, a, b)) == b-a+1
    assert summation(1, (n, 1, 10)) == 10
    assert summation(2*n, (n, 0, 10**10)) == 100000000010000000000
    assert summation(4*n*m, (n, a, 1), (m, 1, d)).expand() == \
        2*d + 2*d**2 + a*d + a*d**2 - d*a**2 - a**2*d**2
    assert summation(cos(n), (n, -2, 1)) == cos(-2)+cos(-1)+cos(0)+cos(1)

def test_polynomial_sums():
    assert summation(n**2, (n, 3, 8)) == 199
    assert summation(n, (n, a, b)) == \
        ((a + b)*(b - a + 1)/2).expand()
    assert summation(n**2, (n, 1, b)) == \
        ((2*b**3 + 3*b**2 + b)/6).expand()
    assert summation(n**3, (n, 1, b)) == \
        ((b**4 + 2*b**3 + b**2)/4).expand()
    assert summation(n**6, (n, 1, b)) == \
        ((6*b**7 + 21*b**6 + 21*b**5 - 7*b**3 + b)/42).expand()

def test_arithmetic_sums():
    assert summation(1, (n, a, b)) == b - a + 1
    assert Sum(S.NaN, (n, a, b)) is S.NaN
    assert Sum(x, (n, a, a)).doit() == x
    assert Sum(x, (x, a, a)).doit() == a
    assert Sum(x, (n, 1, a)).doit() == a*x
    assert Sum(x, (x, Range(1, 11))).doit() == 55
    assert Sum(x, (x, Range(1, 11, 2))).doit() == 25
    assert Sum(x, (x, Range(1, 10, 2))) == Sum(x, (x, Range(9, 0, -2)))
    lo, hi = 1, 2
    s1 = Sum(n, (n, lo, hi))
    s2 = Sum(n, (n, hi, lo))
    assert s1 != s2
    assert s1.doit() == 3 and s2.doit() == 0
    lo, hi = x, x + 1
    s1 = Sum(n, (n, lo, hi))
    s2 = Sum(n, (n, hi, lo))
    assert s1 != s2
    assert s1.doit() == 2*x + 1 and s2.doit() == 0
    assert Sum(Integral(x, (x, 1, y)) + x, (x, 1, 2)).doit() == \
        y**2 + 2
    assert summation(1, (n, 1, 10)) == 10
    assert summation(2*n, (n, 0, 10**10)) == 100000000010000000000
    assert summation(4*n*m, (n, a, 1), (m, 1, d)).expand() == \
        2*d + 2*d**2 + a*d + a*d**2 - d*a**2 - a**2*d**2
    assert summation(cos(n), (n, -2, 1)) == cos(-2) + cos(-1) + cos(0) + cos(1)
    assert summation(cos(n), (n, x, x + 2)) == cos(x) + cos(x + 1) + cos(x + 2)
    assert isinstance(summation(cos(n), (n, x, x + S.Half)), Sum)
    assert summation(k, (k, 0, oo)) == oo
    assert summation(k, (k, Range(1, 11))) == 55

def test_geometric_sums():
    assert summation(pi**n, (n, 0, b)) == (1-pi**(b+1)) / (1-pi)
    assert summation(2 * 3**n, (n, 0, b)) == 3**(b+1) - 1
    assert summation(Rational(1,2)**n, (n, 1, oo)) == 1
    assert summation(2**n, (n, 0, b)) == 2**(b+1) - 1
    assert summation(2**n, (n, 1, oo)) == oo
    assert summation(2**(-n), (n, 1, oo)) == 1
    assert summation(3**(-n), (n, 4, oo)) == Rational(1,54)
    assert summation(2**(-4*n+3), (n, 1, oo)) == Rational(8,15)
    assert summation(2**(n+1), (n, 1, b)).expand() == 4*(2**b-1)

    def expectation(self, expr, var, evaluate=True, **kwargs):
        """ Expectation of expression over distribution """
        # TODO: support discrete sets with non integer stepsizes

        if evaluate:
            try:
                p = poly(expr, var)

                t = Dummy('t', real=True)

                mgf = self.moment_generating_function(t)
                deg = p.degree()
                taylor = poly(series(mgf, t, 0, deg + 1).removeO(), t)
                result = 0
                for k in range(deg+1):
                    result += p.coeff_monomial(var ** k) * taylor.coeff_monomial(t ** k) * factorial(k)

                return result

            except PolynomialError:
                return summation(expr * self.pdf(var),
                                 (var, self.set.inf, self.set.sup), **kwargs)

        else:
            return Sum(expr * self.pdf(var),
                         (var, self.set.inf, self.set.sup), **kwargs)

 def to_sympy(self, expr, **kwargs):
     if expr.has_form('Sum', 2) and expr.leaves[1].has_form('List', 3):
         index = expr.leaves[1]
         result = sympy.summation(expr.leaves[0].to_sympy(), (
             index.leaves[0].to_sympy(), index.leaves[1].to_sympy(),
             index.leaves[2].to_sympy()))            
         return result

 def to_sympy(self, expr, **kwargs):
     if expr.has_form('Sum', 2) and expr.leaves[1].has_form('List', 3):
         index = expr.leaves[1]
         arg = expr.leaves[0].to_sympy()
         bounds = (index.leaves[0].to_sympy(), index.leaves[1].to_sympy(), index.leaves[2].to_sympy())
         if arg is not None and None not in bounds:
             return sympy.summation(arg, bounds)

def test_hypersum():
    from sympy import simplify, sin, hyper
    assert simplify(summation(x**n/fac(n), (n, 1, oo))) == -1 + exp(x)
    assert summation((-1)**n * x**(2*n) / fac(2*n), (n, 0, oo)) == cos(x)
    assert simplify(summation((-1)**n*x**(2*n+1)/factorial(2*n+1),
                              (n, 3, oo))) \
           == -x + sin(x) + x**3/6 - x**5/120

    # TODO to get this without hyper need to improve hyperexpand
    assert summation(1/(n+2)**3, (n, 1, oo)) == \
           hyper([3, 3, 3, 1], [4, 4, 4], 1)/27

    s = summation(x**n*n, (n, -oo, 0))
    assert s.is_Piecewise
    assert s.args[0].args[0] == -1/(x*(1-1/x)**2)
    assert s.args[0].args[1] == (abs(1/x) < 1)

def test_Sum_doit():
    assert Sum(n*Integral(a**2), (n, 0, 2)).doit() == a**3
    assert Sum(n*Integral(a**2), (n, 0, 2)).doit(deep=False) == \
        3*Integral(a**2)
    assert summation(n*Integral(a**2), (n, 0, 2)) == 3*Integral(a**2)

    # test nested sum evaluation
    s = Sum( Sum( Sum(2,(z,1,n+1)), (y,x+1,n)), (x,1,n))
    assert 0 == (s.doit() - n*(n+1)*(n-1)).factor()

    assert Sum(Sum(KroneckerDelta(m, n), (m, 1, 3)), (n, 1, 3)).doit() == 3
    assert Sum(Sum(KroneckerDelta(k, m), (m, 1, 3)), (n, 1, 3)).doit() == \
        3*Piecewise((1, And(S(1) <= k, k <= 3)), (0, True))
    assert Sum(f(n)*Sum(KroneckerDelta(m, n), (m, 0, oo)), (n, 1, 3)).doit() == \
        f(1) + f(2) + f(3)
    assert Sum(f(n)*Sum(KroneckerDelta(m, n), (m, 0, oo)), (n, 1, oo)).doit() == \
        Sum(Piecewise((f(n), n >= 0), (0, True)), (n, 1, oo))
    l = Symbol('l', integer=True, positive=True)
    assert Sum(f(l)*Sum(KroneckerDelta(m, l), (m, 0, oo)), (l, 1, oo)).doit() == \
        Sum(f(l), (l, 1, oo))

    # issue 2597
    nmax = symbols('N', integer=True, positive=True)
    pw = Piecewise((1, And(S(1) <= n, n <= nmax)), (0, True))
    assert Sum(pw, (n, 1, nmax)).doit() == Sum(pw, (n, 1, nmax))

    def __init__(self, states, interval, differential_order):
        """
        :param states: tuple of states in beginning and end of interval
        :param interval: time interval (tuple)
        :param differential_order: grade of differential flatness :math:`\\gamma`
        """
        self.yd = states
        self.t0 = interval[0]
        self.t1 = interval[1]
        self.dt = interval[1] - interval[0]
        gamma = differential_order  # + 1 # TODO check this against notes

        # setup symbolic expressions
        tau, k = sp.symbols('tau, k')

        alpha = sp.factorial(2 * gamma + 1)

        f = sp.binomial(gamma, k) * (-1) ** k * tau ** (gamma + k + 1) / (gamma + k + 1)
        phi = alpha / sp.factorial(gamma) ** 2 * sp.summation(f, (k, 0, gamma))

        # differentiate phi(tau), index in list corresponds to order
        dphi_sym = [phi]  # init with phi(tau)
        for order in range(differential_order):
            dphi_sym.append(dphi_sym[-1].diff(tau))

        # lambdify
        self.dphi_num = []
        for der in dphi_sym:
            self.dphi_num.append(sp.lambdify(tau, der, 'numpy'))

def repeated(n, i, e, a, b):
    # let n_a = n; n_{a+k+1} = e(n=n_{a+k}, i=a+k+1)
    # return n_b
    if e.has(i):
        if e.has(n):
            if not (e - n).has(n):
                term = e - n
                return n + sympy.summation(term, (i, a, b))
            raise NotImplementedError("has i and n")
        else:
            return e.subs(i, b)
    else:
        if e.has(n):
            if not (e - n).has(n):
                term = e - n
                return n + term * (b - a + 1)
            c, args = e.as_coeff_add(n)
            arg, = args
            if not (arg / n).simplify().has(n):
                coeff = arg / n
                # print("Coefficient is %s, iterations is %s" %
                #       (coeff, (b-a+1)))
                return n * coeff ** (b - a + 1)
            raise NotImplementedError
        else:
            return e

def test_totient():
    assert [totient(k) for k in range(1, 12)] == \
        [1, 1, 2, 2, 4, 2, 6, 4, 6, 4, 10]
    assert totient(5005) == 2880
    assert totient(5006) == 2502
    assert totient(5009) == 5008
    assert totient(2**100) == 2**99

    raises(ValueError, lambda: totient(30.1))
    raises(ValueError, lambda: totient(20.001))

    m = Symbol("m", integer=True)
    assert totient(m)
    assert totient(m).subs(m, 3**10) == 3**10 - 3**9
    assert summation(totient(m), (m, 1, 11)) == 42

    n = Symbol("n", integer=True, positive=True)
    assert totient(n).is_integer

    x=Symbol("x", integer=False)
    raises(ValueError, lambda: totient(x))

    y=Symbol("y", positive=False)
    raises(ValueError, lambda: totient(y))

    z=Symbol("z", positive=True, integer=True)
    raises(ValueError, lambda: totient(2**(-z)))

def test_Sum_doit():
    assert Sum(n * Integral(a ** 2), (n, 0, 2)).doit() == a ** 3
    assert Sum(n * Integral(a ** 2), (n, 0, 2)).doit(deep=False) == 3 * Integral(a ** 2)
    assert summation(n * Integral(a ** 2), (n, 0, 2)) == 3 * Integral(a ** 2)

    # test nested sum evaluation
    S = Sum(Sum(Sum(2, (z, 1, n + 1)), (y, x + 1, n)), (x, 1, n))
    assert 0 == (S.doit() - n * (n + 1) * (n - 1)).factor()

def test_composite_sums():
    f = Rational(1, 2)*(7 - 6*n + Rational(1, 7)*n**3)
    s = summation(f, (n, a, b))
    assert not isinstance(s, Sum)
    A = 0
    for i in range(-3, 5):
        A += f.subs(n, i)
    B = s.subs(a, -3).subs(b, 4)
    assert A == B

 def expectation(self, expr, var, evaluate=True, **kwargs):
     """ Expectation of expression over distribution """
     # TODO: support discrete sets with non integer stepsizes
     if evaluate:
         return summation(expr * self.pdf(var),
                      (var, self.set.inf, self.set.sup), **kwargs)
     else:
         return Sum(expr * self.pdf(var),
                      (var, self.set.inf, self.set.sup), **kwargs)

def test_Sum_interface():
    assert isinstance(Sum(0, (n, 0, 2)), Sum)
    assert Sum(nan, (n, 0, 2)) is nan
    assert Sum(nan, (n, 0, oo)) is nan
    assert Sum(0, (n, 0, 2)).doit() == 0
    assert isinstance(Sum(0, (n, 0, oo)), Sum)
    assert Sum(0, (n, 0, oo)).doit() == 0
    raises(ValueError, lambda: Sum(1))
    raises(ValueError, lambda: summation(1))

    def compute_characteristic_function(self, **kwargs):
        """ Compute the characteristic function from the PDF

        Returns a Lambda
        """
        x, t = symbols('x, t', real=True, finite=True, cls=Dummy)
        pdf = self.pdf(x)
        cf = summation(exp(I*t*x)*pdf, (x, self.set.inf, self.set.sup))
        return Lambda(t, cf)

    def _eval_product(self, a, n, term):
        from sympy import summation, Sum
        k = self.index

        if not term.has(k):
            return term**(n-a+1)
        elif term.is_polynomial(k):
            poly = term.as_poly(k)

            A = B = Q = S.One
            C_= poly.LC()

            all_roots = roots(poly, multiple=True)

            for r in all_roots:
                A *= C.RisingFactorial(a-r, n-a+1)
                Q *= n - r

            if len(all_roots) < poly.degree():
                B = Product(quo(poly, Q.as_poly(k)), (k, a, n))

            return poly.LC()**(n-a+1) * A * B
        elif term.is_Add:
            p, q = term.as_numer_denom()

            p = self._eval_product(a, n, p)
            q = self._eval_product(a, n, q)

            return p / q
        elif term.is_Mul:
            exclude, include = [], []

            for t in term.args:
                p = self._eval_product(a, n, t)

                if p is not None:
                    exclude.append(p)
                else:
                    include.append(t)

            if not exclude:
                return None
            else:
                A, B = Mul(*exclude), term._new_rawargs(*include)
                return A * Product(B, (k, a, n))
        elif term.is_Pow:
            if not term.base.has(k):
                s = summation(term.exp, (k, a, n))

                if not isinstance(s, Sum):
                    return term.base**s
            elif not term.exp.has(k):
                p = self._eval_product(a, n, term.base)

                if p is not None:
                    return p**term.exp

def test_totient():
    assert [totient(k) for k in range(1, 12)] == \
        [1, 1, 2, 2, 4, 2, 6, 4, 6, 4, 10]
    assert totient(5005) == 2880
    assert totient(5006) == 2502
    assert totient(5009) == 5008
    assert totient(2**100) == 2**99

    m = Symbol("m", integer=True)
    assert totient(m)
    assert totient(m).subs(m, 3**10) == 3**10 - 3**9
    assert summation(totient(m), (m, 1, 11)) == 42

def test_primeomega():
    assert primeomega(2) == 1
    assert primeomega(2 * 2) == 2
    assert primeomega(2 * 2 * 3) == 3
    assert primeomega(3 * 25) == primeomega(3) + primeomega(25)
    assert [primeomega(p) for p in primerange(1, 10)] == [1, 1, 1, 1]
    assert primeomega(fac(50)) == 108
    assert primeomega(2 ** 9941 - 1) == 1
    n = Symbol('n', integer=True)
    assert primeomega(n)
    assert primeomega(n).subs(n, 2 ** 31 - 1) == 1
    assert summation(primeomega(n), (n, 2, 30)) == 59