def _add_splines(c, b1, d, b2):
    """Construct c*b1 + d*b2."""
    if b1 == S.Zero or c == S.Zero:
        return expand(piecewise_fold(d * b2))
    if b2 == S.Zero or d == S.Zero:
        return expand(piecewise_fold(c * b1))
    new_args = []
    n_intervals = len(b1.args)
    assert n_intervals == len(b2.args)
    new_args.append((expand(c * b1.args[0].expr), b1.args[0].cond))
    for i in range(1, n_intervals - 1):
        new_args.append((expand(c * b1.args[i].expr + d * b2.args[i - 1].expr), b1.args[i].cond))
    new_args.append((expand(d * b2.args[-2].expr), b2.args[-2].cond))
    new_args.append(b2.args[-1])
    return Piecewise(*new_args)

def _add_splines(c, b1, d, b2):
    """Construct c*b1 + d*b2."""
    if b1 == S.Zero or c == S.Zero:
        rv = piecewise_fold(d*b2)
    elif b2 == S.Zero or d == S.Zero:
        rv = piecewise_fold(c*b1)
    else:
        new_args = []
        n_intervals = len(b1.args)
        if n_intervals != len(b2.args):
            raise ValueError("Args of b1 and b2 are not equal")
        new_args.append((c*b1.args[0].expr, b1.args[0].cond))
        for i in range(1, n_intervals - 1):
            new_args.append((
                c*b1.args[i].expr + d*b2.args[i - 1].expr,
                b1.args[i].cond
            ))
        new_args.append((d*b2.args[-2].expr, b2.args[-2].cond))
        new_args.append(b2.args[-1])
        rv = Piecewise(*new_args)

    return rv.expand()

def test_deltasummation_mul_add_x_kd_add_y_kd():
    assert ds((x + KD(i, k))*(y + KD(i, j)), (j, 1, 3)) == piecewise_fold(
        Piecewise((KD(i, k) + x, And(S(1) <= i, i <= 3)), (0, True)) +
        3*(KD(i, k) + x)*y)
    assert ds((x + KD(i, k))*(y + KD(i, j)), (j, 1, 1)) == piecewise_fold(
        Piecewise((KD(i, k) + x, Eq(i, 1)), (0, True)) +
        (KD(i, k) + x)*y)
    assert ds((x + KD(i, k))*(y + KD(i, j)), (j, 2, 2)) == piecewise_fold(
        Piecewise((KD(i, k) + x, Eq(i, 2)), (0, True)) +
        (KD(i, k) + x)*y)
    assert ds((x + KD(i, k))*(y + KD(i, j)), (j, 3, 3)) == piecewise_fold(
        Piecewise((KD(i, k) + x, Eq(i, 3)), (0, True)) +
        (KD(i, k) + x)*y)
    assert ds((x + KD(i, k))*(y + KD(i, j)), (j, 1, k)) == piecewise_fold(
        Piecewise((KD(i, k) + x, And(S(1) <= i, i <= k)), (0, True)) +
        k*(KD(i, k) + x)*y)
    assert ds((x + KD(i, k))*(y + KD(i, j)), (j, k, 3)) == piecewise_fold(
        Piecewise((KD(i, k) + x, And(k <= i, i <= 3)), (0, True)) +
        (4 - k)*(KD(i, k) + x)*y)
    assert ds((x + KD(i, k))*(y + KD(i, j)), (j, k, l)) == piecewise_fold(
        Piecewise((KD(i, k) + x, And(k <= i, i <= l)), (0, True)) +
        (l - k + 1)*(KD(i, k) + x)*y)

def test_deltasummation_mul_add_x_y_add_kd_kd():
    assert ds((x + y)*(KD(i, k) + KD(j, k)), (k, 1, 3)) == piecewise_fold(
        Piecewise((x + y, And(S(1) <= i, i <= 3)), (0, True)) +
        Piecewise((x + y, And(S(1) <= j, j <= 3)), (0, True)))
    assert ds((x + y)*(KD(i, k) + KD(j, k)), (k, 1, 1)) == piecewise_fold(
        Piecewise((x + y, Eq(i, 1)), (0, True)) +
        Piecewise((x + y, Eq(j, 1)), (0, True)))
    assert ds((x + y)*(KD(i, k) + KD(j, k)), (k, 2, 2)) == piecewise_fold(
        Piecewise((x + y, Eq(i, 2)), (0, True)) +
        Piecewise((x + y, Eq(j, 2)), (0, True)))
    assert ds((x + y)*(KD(i, k) + KD(j, k)), (k, 3, 3)) == piecewise_fold(
        Piecewise((x + y, Eq(i, 3)), (0, True)) +
        Piecewise((x + y, Eq(j, 3)), (0, True)))
    assert ds((x + y)*(KD(i, k) + KD(j, k)), (k, 1, l)) == piecewise_fold(
        Piecewise((x + y, And(S(1) <= i, i <= l)), (0, True)) +
        Piecewise((x + y, And(S(1) <= j, j <= l)), (0, True)))
    assert ds((x + y)*(KD(i, k) + KD(j, k)), (k, l, 3)) == piecewise_fold(
        Piecewise((x + y, And(l <= i, i <= 3)), (0, True)) +
        Piecewise((x + y, And(l <= j, j <= 3)), (0, True)))
    assert ds((x + y)*(KD(i, k) + KD(j, k)), (k, l, m)) == piecewise_fold(
        Piecewise((x + y, And(l <= i, i <= m)), (0, True)) +
        Piecewise((x + y, And(l <= j, j <= m)), (0, True)))

def deltasummation(f, limit, no_piecewise=False):
    """
    Handle summations containing a KroneckerDelta.

    The idea for summation is the following:

    - If we are dealing with a KroneckerDelta expression, i.e. KroneckerDelta(g(x), j),
      we try to simplify it.

      If we could simplify it, then we sum the resulting expression.
      We already know we can sum a simplified expression, because only
      simple KroneckerDelta expressions are involved.

      If we couldn't simplify it, there are two cases:

      1) The expression is a simple expression: we return the summation,
         taking care if we are dealing with a Derivative or with a proper
         KroneckerDelta.

      2) The expression is not simple (i.e. KroneckerDelta(cos(x))): we can do
         nothing at all.

    - If the expr is a multiplication expr having a KroneckerDelta term:

      First we expand it.

      If the expansion did work, then we try to sum the expansion.

      If not, we try to extract a simple KroneckerDelta term, then we have two
      cases:

      1) We have a simple KroneckerDelta term, so we return the summation.

      2) We didn't have a simple term, but we do have an expression with
         simplified KroneckerDelta terms, so we sum this expression.

    Examples
    ========

    >>> from sympy import oo
    >>> from sympy.abc import i, j, k
    >>> from sympy.concrete.delta import deltasummation
    >>> from sympy import KroneckerDelta, Piecewise
    >>> deltasummation(KroneckerDelta(i, k), (k, -oo, oo))
    1
    >>> deltasummation(KroneckerDelta(i, k), (k, 0, oo))
    Piecewise((1, i >= 0), (0, True))
    >>> deltasummation(KroneckerDelta(i, k), (k, 1, 3))
    Piecewise((1, And(1 <= i, i <= 3)), (0, True))
    >>> deltasummation(k*KroneckerDelta(i, j)*KroneckerDelta(j, k), (k, -oo, oo))
    j*KroneckerDelta(i, j)
    >>> deltasummation(j*KroneckerDelta(i, j), (j, -oo, oo))
    i
    >>> deltasummation(i*KroneckerDelta(i, j), (i, -oo, oo))
    j

    See Also
    ========

    deltaproduct
    sympy.functions.special.tensor_functions.KroneckerDelta
    sympy.concrete.sums.summation
    """
    from sympy.concrete.summations import summation
    from sympy.solvers import solve

    if ((limit[2] - limit[1]) < 0) is True:
        return S.Zero

    if not f.has(KroneckerDelta):
        return summation(f, limit)

    x = limit[0]

    g = _expand_delta(f, x)
    if g.is_Add:
        return piecewise_fold(
            g.func(*[deltasummation(h, limit, no_piecewise) for h in g.args]))

    # try to extract a simple KroneckerDelta term
    delta, expr = _extract_delta(g, x)

    if not delta:
        return summation(f, limit)

    solns = solve(delta.args[0] - delta.args[1], x)
    if len(solns) == 0:
        return S.Zero
    elif len(solns) != 1:
        return Sum(f, limit)
    value = solns[0]
    if no_piecewise:
        return expr.subs(x, value)
    return Piecewise(
        (expr.subs(x, value), Interval(*limit[1:3]).as_relational(value)),
        (S.Zero, True)
    )

def _solveset(f, symbol, domain, _check=False):
    """Helper for solveset to return a result from an expression
    that has already been sympify'ed and is known to contain the
    given symbol."""
    # _check controls whether the answer is checked or not

    from sympy.simplify.simplify import signsimp
    orig_f = f
    f = together(f)
    if f.is_Mul:
        _, f = f.as_independent(symbol, as_Add=False)
    if f.is_Add:
        a, h = f.as_independent(symbol)
        m, h = h.as_independent(symbol, as_Add=False)
        f = a/m + h  # XXX condition `m != 0` should be added to soln
    f = piecewise_fold(f)

    # assign the solvers to use
    solver = lambda f, x, domain=domain: _solveset(f, x, domain)
    if domain.is_subset(S.Reals):
        inverter_func = invert_real
    else:
        inverter_func = invert_complex
    inverter = lambda f, rhs, symbol: inverter_func(f, rhs, symbol, domain)

    result = EmptySet()

    if f.expand().is_zero:
        return domain
    elif not f.has(symbol):
        return EmptySet()
    elif f.is_Mul and all(_is_finite_with_finite_vars(m, domain)
            for m in f.args):
        # if f(x) and g(x) are both finite we can say that the solution of
        # f(x)*g(x) == 0 is same as Union(f(x) == 0, g(x) == 0) is not true in
        # general. g(x) can grow to infinitely large for the values where
        # f(x) == 0. To be sure that we are not silently allowing any
        # wrong solutions we are using this technique only if both f and g are
        # finite for a finite input.
        result = Union(*[solver(m, symbol) for m in f.args])
    elif _is_function_class_equation(TrigonometricFunction, f, symbol) or \
            _is_function_class_equation(HyperbolicFunction, f, symbol):
        result = _solve_trig(f, symbol, domain)
    elif f.is_Piecewise:
        dom = domain
        result = EmptySet()
        expr_set_pairs = f.as_expr_set_pairs()
        for (expr, in_set) in expr_set_pairs:
            if in_set.is_Relational:
                in_set = in_set.as_set()
            if in_set.is_Interval:
                dom -= in_set
            solns = solver(expr, symbol, in_set)
            result += solns
    else:
        lhs, rhs_s = inverter(f, 0, symbol)
        if lhs == symbol:
            # do some very minimal simplification since
            # repeated inversion may have left the result
            # in a state that other solvers (e.g. poly)
            # would have simplified; this is done here
            # rather than in the inverter since here it
            # is only done once whereas there it would
            # be repeated for each step of the inversion
            if isinstance(rhs_s, FiniteSet):
                rhs_s = FiniteSet(*[Mul(*
                    signsimp(i).as_content_primitive())
                    for i in rhs_s])
            result = rhs_s
        elif isinstance(rhs_s, FiniteSet):
            for equation in [lhs - rhs for rhs in rhs_s]:
                if equation == f:
                    if any(_has_rational_power(g, symbol)[0]
                           for g in equation.args) or _has_rational_power(
                           equation, symbol)[0]:
                        result += _solve_radical(equation,
                                                 symbol,
                                                 solver)
                    elif equation.has(Abs):
                        result += _solve_abs(f, symbol, domain)
                    else:
                        result += _solve_as_rational(equation, symbol, domain)
                else:
                    result += solver(equation, symbol)
        else:
            result = ConditionSet(symbol, Eq(f, 0), domain)

    if _check:
        if isinstance(result, ConditionSet):
            # it wasn't solved or has enumerated all conditions
            # -- leave it alone
            return result

        # whittle away all but the symbol-containing core
        # to use this for testing
        fx = orig_f.as_independent(symbol, as_Add=True)[1]
        fx = fx.as_independent(symbol, as_Add=False)[1]

        if isinstance(result, FiniteSet):
            # check the result for invalid solutions
#.........这里部分代码省略.........

#.........这里部分代码省略.........

    original_expr = expr = signsimp(expr)

    from sympy.simplify.hyperexpand import hyperexpand
    from sympy.functions.special.bessel import BesselBase
    from sympy import Sum, Product

    if not isinstance(expr, Basic) or not expr.args:  # XXX: temporary hack
        return expr

    if not isinstance(expr, (Add, Mul, Pow, ExpBase)):
        if isinstance(expr, Function) and hasattr(expr, "inverse"):
            if len(expr.args) == 1 and len(expr.args[0].args) == 1 and \
               isinstance(expr.args[0], expr.inverse(argindex=1)):
                return simplify(expr.args[0].args[0], ratio=ratio,
                                measure=measure, fu=fu)
        return expr.func(*[simplify(x, ratio=ratio, measure=measure, fu=fu)
                         for x in expr.args])

    # TODO: Apply different strategies, considering expression pattern:
    # is it a purely rational function? Is there any trigonometric function?...
    # See also https://github.com/sympy/sympy/pull/185.

    def shorter(*choices):
        '''Return the choice that has the fewest ops. In case of a tie,
        the expression listed first is selected.'''
        if not has_variety(choices):
            return choices[0]
        return min(choices, key=measure)

    expr = bottom_up(expr, lambda w: w.normal())
    expr = Mul(*powsimp(expr).as_content_primitive())
    _e = cancel(expr)
    expr1 = shorter(_e, _mexpand(_e).cancel())  # issue 6829
    expr2 = shorter(together(expr, deep=True), together(expr1, deep=True))

    if ratio is S.Infinity:
        expr = expr2
    else:
        expr = shorter(expr2, expr1, expr)
    if not isinstance(expr, Basic):  # XXX: temporary hack
        return expr

    expr = factor_terms(expr, sign=False)

    # hyperexpand automatically only works on hypergeometric terms
    expr = hyperexpand(expr)

    expr = piecewise_fold(expr)

    if expr.has(BesselBase):
        expr = besselsimp(expr)

    if expr.has(TrigonometricFunction) and not fu or expr.has(
            HyperbolicFunction):
        expr = trigsimp(expr, deep=True)

    if expr.has(log):
        expr = shorter(expand_log(expr, deep=True), logcombine(expr))

    if expr.has(CombinatorialFunction, gamma):
        expr = combsimp(expr)

    if expr.has(Sum):
        expr = sum_simplify(expr)

    if expr.has(Product):
        expr = product_simplify(expr)

    short = shorter(powsimp(expr, combine='exp', deep=True), powsimp(expr), expr)
    short = shorter(short, factor_terms(short), expand_power_exp(expand_mul(short)))
    if short.has(TrigonometricFunction, HyperbolicFunction, ExpBase):
        short = exptrigsimp(short, simplify=False)

    # get rid of hollow 2-arg Mul factorization
    hollow_mul = Transform(
        lambda x: Mul(*x.args),
        lambda x:
        x.is_Mul and
        len(x.args) == 2 and
        x.args[0].is_Number and
        x.args[1].is_Add and
        x.is_commutative)
    expr = short.xreplace(hollow_mul)

    numer, denom = expr.as_numer_denom()
    if denom.is_Add:
        n, d = fraction(radsimp(1/denom, symbolic=False, max_terms=1))
        if n is not S.One:
            expr = (numer*n).expand()/d

    if expr.could_extract_minus_sign():
        n, d = fraction(expr)
        if d != 0:
            expr = signsimp(-n/(-d))

    if measure(expr) > ratio*measure(original_expr):
        expr = original_expr

    return expr

def solveset_real(f, symbol):
    """ Solves a real valued equation.

    Parameters
    ==========

    f : Expr
        The target equation
    symbol : Symbol
        The variable for which the equation is solved

    Returns
    =======

    Set
        A set of values for `symbol` for which `f` is equal to
        zero. An `EmptySet` is returned if no solution is found.

    `solveset_real` claims to be complete in the set of the solution it
    returns.

    Raises
    ======

    NotImplementedError
        The algorithms for to find the solution of the given equation are
        not yet implemented.
    ValueError
        The input is not valid.
    RuntimeError
        It is a bug, please report to the github issue tracker.


    See Also
    =======

    solveset_complex : solver for complex domain

    Examples
    ========

    >>> from sympy import Symbol, exp, sin, sqrt, I
    >>> from sympy.solvers.solveset import solveset_real
    >>> x = Symbol('x', real=True)
    >>> a = Symbol('a', real=True, finite=True, positive=True)
    >>> solveset_real(x**2 - 1, x)
    {-1, 1}
    >>> solveset_real(sqrt(5*x + 6) - 2 - x, x)
    {-1, 2}
    >>> solveset_real(x - I, x)
    EmptySet()
    >>> solveset_real(x - a, x)
    {a}
    >>> solveset_real(exp(x) - a, x)
    {log(a)}

    In case the equation has infinitely many solutions an infinitely indexed
    `ImageSet` is returned.

    >>> solveset_real(sin(x) - 1, x)
    ImageSet(Lambda(_n, 2*_n*pi + pi/2), Integers())

    If the equation is true for any arbitrary value of the symbol a `S.Reals`
    set is returned.

    >>> solveset_real(x - x, x)
    (-oo, oo)

    """
    if not symbol.is_Symbol:
        raise ValueError(" %s is not a symbol" % (symbol))

    f = sympify(f)
    if not isinstance(f, (Expr, Number)):
        raise ValueError(" %s is not a valid sympy expression" % (f))

    original_eq = f
    f = together(f)

    if f.has(Piecewise):
        f = piecewise_fold(f)
    result = EmptySet()

    if f.expand().is_zero:
        return S.Reals
    elif not f.has(symbol):
        return EmptySet()
    elif f.is_Mul and all([_is_finite_with_finite_vars(m) for m in f.args]):
        # if f(x) and g(x) are both finite we can say that the solution of
        # f(x)*g(x) == 0 is same as Union(f(x) == 0, g(x) == 0) is not true in
        # general. g(x) can grow to infinitely large for the values where
        # f(x) == 0. To be sure that we not are silently allowing any
        # wrong solutions we are using this technique only if both f and g and
        # finite for a finite input.
        result = Union(*[solveset_real(m, symbol) for m in f.args])
    elif _is_function_class_equation(C.TrigonometricFunction, f, symbol) or \
            _is_function_class_equation(C.HyperbolicFunction, f, symbol):
        result = _solve_real_trig(f, symbol)
    elif f.is_Piecewise:
        result = EmptySet()
#.........这里部分代码省略.........

def test_deltasummation():
    ds = deltasummation
    assert ds(x, (j, 1, 0)) == 0
    assert ds(x, (j, 1, 3)) == 3*x
    assert ds(x + y, (j, 1, 3)) == 3*(x + y)
    assert ds(x*y, (j, 1, 3)) == 3*x*y
    assert ds(KD(i, j), (k, 1, 3)) == 3*KD(i, j)
    assert ds(x*KD(i, j), (k, 1, 3)) == 3*x*KD(i, j)
    assert ds(x*y*KD(i, j), (k, 1, 3)) == 3*x*y*KD(i, j)

    # return unevaluated, until it gets implemented
    assert ds(KD(i**2, j**2), (j, -oo, oo)) == \
        Sum(KD(i**2, j**2), (j, -oo, oo))

    assert Piecewise((KD(i, k), And(S(1) <= i, i <= 3)), (0, True)) == \
        ds(KD(i, j)*KD(j, k), (j, 1, 3)) == \
        ds(KD(j, k)*KD(i, j), (j, 1, 3))

    assert ds(KD(i, k), (k, -oo, oo)) == 1
    assert ds(KD(i, k), (k, 0, oo)) == Piecewise((1, i >= 0), (0, True))
    assert ds(KD(i, k), (k, 1, 3)) == \
        Piecewise((1, And(S(1) <= i, i <= 3)), (0, True))
    assert ds(k*KD(i, j)*KD(j, k), (k, -oo, oo)) == j*KD(i, j)
    assert ds(j*KD(i, j), (j, -oo, oo)) == i
    assert ds(i*KD(i, j), (i, -oo, oo)) == j
    assert ds(x, (i, 1, 3)) == 3*x
    assert ds((i + j)*KD(i, j), (j, -oo, oo)) == 2*i

    assert ds(KD(i, j), (j, 1, 3)) == \
        Piecewise((1, And(S(1) <= i, i <= 3)), (0, True))
    assert ds(KD(i, j), (j, 1, 1)) == Piecewise((1, Eq(i, 1)), (0, True))
    assert ds(KD(i, j), (j, 2, 2)) == Piecewise((1, Eq(i, 2)), (0, True))
    assert ds(KD(i, j), (j, 3, 3)) == Piecewise((1, Eq(i, 3)), (0, True))
    assert ds(KD(i, j), (j, 1, k)) == \
        Piecewise((1, And(S(1) <= i, i <= k)), (0, True))
    assert ds(KD(i, j), (j, k, 3)) == \
        Piecewise((1, And(k <= i, i <= 3)), (0, True))
    assert ds(KD(i, j), (j, k, l)) == \
        Piecewise((1, And(k <= i, i <= l)), (0, True))

    assert ds(x*KD(i, j), (j, 1, 3)) == \
        Piecewise((x, And(S(1) <= i, i <= 3)), (0, True))
    assert ds(x*KD(i, j), (j, 1, 1)) == Piecewise((x, Eq(i, 1)), (0, True))
    assert ds(x*KD(i, j), (j, 2, 2)) == Piecewise((x, Eq(i, 2)), (0, True))
    assert ds(x*KD(i, j), (j, 3, 3)) == Piecewise((x, Eq(i, 3)), (0, True))
    assert ds(x*KD(i, j), (j, 1, k)) == \
        Piecewise((x, And(S(1) <= i, i <= k)), (0, True))
    assert ds(x*KD(i, j), (j, k, 3)) == \
        Piecewise((x, And(k <= i, i <= 3)), (0, True))
    assert ds(x*KD(i, j), (j, k, l)) == \
        Piecewise((x, And(k <= i, i <= l)), (0, True))

    assert ds((x + y)*KD(i, j), (j, 1, 3)) == \
        Piecewise((x + y, And(S(1) <= i, i <= 3)), (0, True))
    assert ds((x + y)*KD(i, j), (j, 1, 1)) == \
        Piecewise((x + y, Eq(i, 1)), (0, True))
    assert ds((x + y)*KD(i, j), (j, 2, 2)) == \
        Piecewise((x + y, Eq(i, 2)), (0, True))
    assert ds((x + y)*KD(i, j), (j, 3, 3)) == \
        Piecewise((x + y, Eq(i, 3)), (0, True))
    assert ds((x + y)*KD(i, j), (j, 1, k)) == \
        Piecewise((x + y, And(S(1) <= i, i <= k)), (0, True))
    assert ds((x + y)*KD(i, j), (j, k, 3)) == \
        Piecewise((x + y, And(k <= i, i <= 3)), (0, True))
    assert ds((x + y)*KD(i, j), (j, k, l)) == \
        Piecewise((x + y, And(k <= i, i <= l)), (0, True))

    assert ds(KD(i, k) + KD(j, k), (k, 1, 3)) == piecewise_fold(
        Piecewise((1, And(S(1) <= i, i <= 3)), (0, True)) +
        Piecewise((1, And(S(1) <= j, j <= 3)), (0, True)))
    assert ds(KD(i, k) + KD(j, k), (k, 1, 1)) == piecewise_fold(
        Piecewise((1, Eq(i, 1)), (0, True)) +
        Piecewise((1, Eq(j, 1)), (0, True)))
    assert ds(KD(i, k) + KD(j, k), (k, 2, 2)) == piecewise_fold(
        Piecewise((1, Eq(i, 2)), (0, True)) +
        Piecewise((1, Eq(j, 2)), (0, True)))
    assert ds(KD(i, k) + KD(j, k), (k, 3, 3)) == piecewise_fold(
        Piecewise((1, Eq(i, 3)), (0, True)) +
        Piecewise((1, Eq(j, 3)), (0, True)))
    assert ds(KD(i, k) + KD(j, k), (k, 1, l)) == piecewise_fold(
        Piecewise((1, And(S(1) <= i, i <= l)), (0, True)) +
        Piecewise((1, And(S(1) <= j, j <= l)), (0, True)))
    assert ds(KD(i, k) + KD(j, k), (k, l, 3)) == piecewise_fold(
        Piecewise((1, And(l <= i, i <= 3)), (0, True)) +
        Piecewise((1, And(l <= j, j <= 3)), (0, True)))
    assert ds(KD(i, k) + KD(j, k), (k, l, m)) == piecewise_fold(
        Piecewise((1, And(l <= i, i <= m)), (0, True)) +
        Piecewise((1, And(l <= j, j <= m)), (0, True)))

    assert ds(x*KD(i, k) + KD(j, k), (k, 1, 3)) == piecewise_fold(
        Piecewise((x, And(S(1) <= i, i <= 3)), (0, True)) +
        Piecewise((1, And(S(1) <= j, j <= 3)), (0, True)))
    assert ds(x*KD(i, k) + KD(j, k), (k, 1, 1)) == piecewise_fold(
        Piecewise((x, Eq(i, 1)), (0, True)) +
        Piecewise((1, Eq(j, 1)), (0, True)))
    assert ds(x*KD(i, k) + KD(j, k), (k, 2, 2)) == piecewise_fold(
        Piecewise((x, Eq(i, 2)), (0, True)) +
        Piecewise((1, Eq(j, 2)), (0, True)))
    assert ds(x*KD(i, k) + KD(j, k), (k, 3, 3)) == piecewise_fold(
        Piecewise((x, Eq(i, 3)), (0, True)) +
#.........这里部分代码省略.........

def _add_splines(c, b1, d, b2):
    """Construct c*b1 + d*b2."""
    if b1 == S.Zero or c == S.Zero:
        rv = piecewise_fold(d*b2)
    elif b2 == S.Zero or d == S.Zero:
        rv = piecewise_fold(c*b1)
    else:
        new_args = []
        n_intervals = len(b1.args)
        if n_intervals != len(b2.args):
            # Args of b1 and b2 are not equal. Just combining the Piecewise without any fancy optimisation
            p1 = piecewise_fold(c*b1)
            p2 = piecewise_fold(d*b2)

            # Search all Piecewise arguments except (0, True)
            p2args = list(p2.args[:-1])

            # This merging algorithm assume the conditions in p1 and p2 are sorted
            for arg in p1.args[:-1]:
                # Conditional of Piecewise are And objects
                # the args of the And object is a tuple of two Relational objects
                # the numerical value is in the .rhs of the Relational object
                expr = arg.expr
                cond = arg.cond

                lower = cond.args[0].rhs

                # Check p2 for matching conditions that can be merged
                for i, arg2 in enumerate(p2args):
                    expr2 = arg2.expr
                    cond2 = arg2.cond

                    lower_2 = cond2.args[0].rhs
                    upper_2 = cond2.args[1].rhs

                    if cond2 == cond:
                        # Conditions match, join expressions
                        expr += expr2
                        # Remove matching element
                        del p2args[i]
                        # No need to check the rest
                        break
                    elif lower_2 < lower and upper_2 <= lower:
                        # Check if arg2 condition smaller than arg1, add to new_args by itself (no match expected in p1)
                        new_args.append(arg2)
                        del p2args[i]
                        break

                # Checked all, add expr and cond
                new_args.append((expr, cond))

            # Add remaining items from p2args
            new_args.extend(p2args)

            # Add final (0, True)
            new_args.append((0, True))
        else:
            new_args.append((c*b1.args[0].expr, b1.args[0].cond))
            for i in range(1, n_intervals - 1):
                new_args.append((
                    c*b1.args[i].expr + d*b2.args[i - 1].expr,
                    b1.args[i].cond
                ))
            new_args.append((d*b2.args[-2].expr, b2.args[-2].cond))
            new_args.append(b2.args[-1])

        rv = Piecewise(*new_args)

    return rv.expand()