Repository: sympy/sympy

Issue:
Slow inverse_laplace_transform
Running an inverse Laplace Transform can take a very long time. Take the following code:

    a, t, s = sympy.symbols('a, t, s', real=True)
    
    f = -15-sympy.exp(6*t)-sympy.sin(3*t)
![image](https://user-images.githubusercontent.com/17097491/173278322-f2e6a39d-6d33-4b84-9ef3-99474a4afa93.png)

    F =  sympy.laplace_transform(f, t, s, noconds=True)
![image](https://user-images.githubusercontent.com/17097491/173278442-16815fad-6d4a-404a-be71-e94becafea39.png)

The above code completes in a few seconds but when executing the following to "undo" the previous transform:

    invF = sympy.inverse_laplace_transform(F, s, t)

The process takes many minutes (at least15 on my machine) to complete.

Discussion:

Patch:
diff --git a/sympy/integrals/laplace.py b/sympy/integrals/laplace.py
index 69d5e1d358..a88bafe263 100644
--- a/sympy/integrals/laplace.py
+++ b/sympy/integrals/laplace.py
@@ -1,6 +1,7 @@
 """Laplace Transforms"""
 from sympy.core import S, pi, I
 from sympy.core.add import Add
+from sympy.core.cache import cacheit
 from sympy.core.function import (
     AppliedUndef, Derivative, expand, expand_complex, expand_mul, expand_trig,
     Lambda, WildFunction, diff)
@@ -21,10 +22,11 @@
 from sympy.functions.special.gamma_functions import digamma, gamma, lowergamma
 from sympy.integrals import integrate, Integral
 from sympy.integrals.transforms import (
-    _simplify, _noconds_, IntegralTransform, IntegralTransformError)
+    _simplify, IntegralTransform, IntegralTransformError)
 from sympy.logic.boolalg import to_cnf, conjuncts, disjuncts, Or, And
 from sympy.matrices.matrices import MatrixBase
 from sympy.polys.matrices.linsolve import _lin_eq2dict
+from sympy.polys.polyerrors import PolynomialError
 from sympy.polys.polyroots import roots
 from sympy.polys.polytools import Poly
 from sympy.polys.rationaltools import together
@@ -266,7 +268,8 @@ def _laplace_deep_collect(f, t):
             return func(*args)
 
 
-def _laplace_build_rules(t, s):
+@cacheit
+def _laplace_build_rules():
     """
     This is an internal helper function that returns the table of Laplace
     transform rules in terms of the time variable `t` and the frequency
@@ -283,12 +286,16 @@ def _laplace_build_rules(t, s):
     to the expression before matching. For most rules it should be
     ``_laplace_deep_collect``.
     """
+    t = Dummy('t')
+    s = Dummy('s')
     a = Wild('a', exclude=[t])
     b = Wild('b', exclude=[t])
     n = Wild('n', exclude=[t])
     tau = Wild('tau', exclude=[t])
     omega = Wild('omega', exclude=[t])
     dco = lambda f: _laplace_deep_collect(f, t)
+    debug('_laplace_build_rules is building rules')
+
     laplace_transform_rules = [
     (a, a/s,
      S.true, S.Zero, dco), # 4.2.1
@@ -463,7 +470,7 @@ def _laplace_build_rules(t, s):
     (besselk(0, a*t), log((s + sqrt(s**2-a**2))/a)/(sqrt(s**2-a**2)),
      S.true, -re(a), dco) # 4.16.23
     ]
-    return laplace_transform_rules
+    return laplace_transform_rules, t, s
 
 
 def _laplace_rule_timescale(f, t, s):
@@ -755,12 +762,12 @@ def _laplace_apply_simple_rules(f, t, s):
     This function applies all simple rules and returns the result if one
     of them gives a result.
     """
-    simple_rules = _laplace_build_rules(t, s)
+    simple_rules, t_, s_ = _laplace_build_rules()
     prep_old = ''
     prep_f = ''
     for t_dom, s_dom, check, plane, prep in simple_rules:
-        if not prep_old==prep:
-            prep_f = prep(f)
+        if not prep_old == prep:
+            prep_f = prep(f.subs({t: t_}))
             prep_old = prep
         ma = prep_f.match(t_dom)
         if ma:
@@ -775,7 +782,8 @@ def _laplace_apply_simple_rules(f, t, s):
                 debug('      f:     %s'%(f,))
                 debug('      rule:  %s o---o %s'%(t_dom, s_dom))
                 debug('      match: %s'%(ma, ))
-                return (s_dom.xreplace(ma), plane.xreplace(ma), c)
+                return (s_dom.xreplace(ma).subs({s_: s}),
+                        plane.xreplace(ma), S.true)
     return None
 
 
@@ -868,7 +876,7 @@ def doit(self, **hints):
 
         Standard hints are the following:
         - ``noconds``:  if True, do not return convergence conditions. The
-        default setting is `False`.
+        default setting is `True`.
         - ``simplify``: if True, it simplifies the final result. This is the
         default behaviour
         """
@@ -1006,8 +1014,7 @@ def laplace_transform(f, t, s, legacy_matrix=True, **hints):
         return LT[0]
 
 
-@_noconds_(True)
-def _inverse_laplace_transform(F, s, t_, plane, simplify=True):
+def _inverse_laplace_transform_integration(F, s, t_, plane, simplify=True):
     """ The backend function for inverse Laplace transforms. """
     from sympy.integrals.meijerint import meijerint_inversion, _get_coeff_exp
     from sympy.integrals.transforms import inverse_mellin_transform
@@ -1034,7 +1041,7 @@ def pw_simp(*args):
         F = F.apart(s)
 
     if F.is_Add:
-        f = Add(*[_inverse_laplace_transform(X, s, t, plane, simplify)\
+        f = Add(*[_inverse_laplace_transform_integration(X, s, t, plane, simplify)\
                      for X in F.args])
         return _simplify(f.subs(t, t_), simplify), True
 
@@ -1089,6 +1096,246 @@ def simp_exp(arg):
     return _simplify(f.subs(t, t_), simplify), cond
 
 
+def _complete_the_square_in_denom(f, s):
+    from sympy.simplify.radsimp import fraction
+    [n, d] = fraction(f)
+    if d.is_polynomial(s):
+        cf = d.as_poly(s).all_coeffs()
+        if len(cf)==3:
+            a, b, c = cf
+            d = a*((s+b/(2*a))**2+c/a-(b/(2*a))**2)
+    return n/d
+
+
+@cacheit
+def _inverse_laplace_build_rules():
+    """
+    This is an internal helper function that returns the table of inverse
+    Laplace transform rules in terms of the time variable `t` and the
+    frequency variable `s`.  It is used by `_inverse_laplace_apply_rules`.
+    """
+    s = Dummy('s')
+    t = Dummy('t')
+    a = Wild('a', exclude=[s])
+    b = Wild('b', exclude=[s])
+    c = Wild('c', exclude=[s])
+
+    debug('_inverse_laplace_build_rules is building rules')
+    def _frac(f, s):
+        try:
+            return f.factor(s)
+        except PolynomialError:
+            return f
+    same = lambda f: f
+    # This list is sorted according to the prep function needed.
+    _ILT_rules = [
+        (a/s, a, S.true, same, 1),
+        (b*(s+a)**(-c), t**(c-1)*exp(-a*t)/gamma(c), c>0, same, 1),
+    ]
+    return _ILT_rules, s, t
+
+
+def _inverse_laplace_apply_simple_rules(f, s, t):
+    """
+    Helper function for the class InverseLaplaceTransform.
+    """
+    if f==1:
+        debug('_inverse_laplace_apply_simple_rules match:')
+        debug('      f:    %s'%(1,))
+        debug('      rule: 1 o---o DiracDelta(%s)'%(t,))
+        return DiracDelta(t), S.true
+
+    _ILT_rules, s_, t_ = _inverse_laplace_build_rules()
+    _prep = ''
+    fsubs = f.subs({s: s_})
+
+    for s_dom, t_dom, check, prep, fac in _ILT_rules:
+        if not _prep == (prep, fac):
+            _F = prep(fsubs*fac)
+            _prep = (prep, fac)
+        ma = _F.match(s_dom)
+        if ma:
+            try:
+                c = check.xreplace(ma)
+            except TypeError:
+                continue
+            if c:
+                debug('_inverse_laplace_apply_simple_rules match:')
+                debug('      f:    %s'%(f,))
+                debug('      rule: %s o---o %s'%(s_dom, t_dom))
+                debug('      ma:   %s'%(ma,))
+                return Heaviside(t)*t_dom.xreplace(ma).subs({t_: t}), S.true
+    return None
+
+
+def _inverse_laplace_time_shift(F, s, t, plane):
+    """
+    Helper function for the class InverseLaplaceTransform.
+    """
+    a = Wild('a', exclude=[s])
+    g = Wild('g')
+
+    if not F.has(s):
+        return F*DiracDelta(t), S.true
+    ma1 = F.match(exp(a*s))
+    if ma1:
+        if ma1[a].is_negative:
+            debug('_inverse_laplace_time_shift match:')
+            debug('      f:    %s'%(F,))
+            debug('      rule: exp(-a*s) o---o DiracDelta(t-a)')
+            debug('      ma:   %s'%(ma1,))
+            return DiracDelta(t+ma1[a]), S.true
+        else:
+            debug('_inverse_laplace_time_shift match: negative time shift')
+            return InverseLaplaceTransform(F, s, t, plane), S.true
+
+    ma1 = F.match(exp(a*s)*g)
+    if ma1:
+        if ma1[a].is_negative:
+            debug('_inverse_laplace_time_shift match:')
+            debug('      f:    %s'%(F,))
+            debug('      rule: exp(-a*s)*F(s) o---o Heaviside(t-a)*f(t-a)')
+            debug('      ma:   %s'%(ma1,))
+            return _inverse_laplace_transform(ma1[g], s, t+ma1[a], plane)
+        else:
+            debug('_inverse_laplace_time_shift match: negative time shift')
+            return InverseLaplaceTransform(F, s, t, plane), S.true
+    return None
+
+
+def _inverse_laplace_apply_prog_rules(F, s, t, plane):
+    """
+    Helper function for the class InverseLaplaceTransform.
+    """
+    prog_rules = [_inverse_laplace_time_shift]
+
+    for p_rule in prog_rules:
+        if (r := p_rule(F, s, t, plane)) is not None:
+            return r
+    return None
+
+
+def _inverse_laplace_expand(fn, s, t, plane):
+    """
+    Helper function for the class InverseLaplaceTransform.
+    """
+    if fn.is_Add:
+        return None
+    r = expand(fn, deep=False)
+    if r.is_Add:
+        return _inverse_laplace_transform(r, s, t, plane)
+    r = expand_mul(fn)
+    if r.is_Add:
+        return _inverse_laplace_transform(r, s, t, plane)
+    r = expand(fn)
+    if r.is_Add:
+        return _inverse_laplace_transform(r, s, t, plane)
+    if fn.is_rational_function(s):
+        r = fn.apart(s).doit()
+    if r.is_Add:
+        return _inverse_laplace_transform(r, s, t, plane)
+    return None
+
+
+def _inverse_laplace_rational(fn, s, t, plane, simplify):
+    """
+    Helper function for the class InverseLaplaceTransform.
+    """
+    debug('[ILT _i_l_r] (%s, %s, %s)'%(fn, s, t))
+    x_ = symbols('x_')
+    f = fn.apart(s)
+    terms = Add.make_args(f)
+    terms_t = []
+    conditions = [S.true]
+    for term in terms:
+        [n, d] = term.as_numer_denom()
+        dc = d.as_poly(s).all_coeffs()
+        dc_lead = dc[0]
+        dc = [ x/dc_lead for x in dc ]
+        nc = [ x/dc_lead for x in n.as_poly(s).all_coeffs() ]
+        if len(dc)==1:
+            r = nc[0]*DiracDelta(t)
+            terms_t.append(r)
+        elif len(dc)==2:
+            r = nc[0]*exp(-dc[1]*t)
+            terms_t.append(Heaviside(t)*r)
+        elif len(dc)==3:
+            a = dc[1]/2
+            b = (dc[2]-a**2).factor()
+            if len(nc)==1:
+                nc = [S.Zero] + nc
+            l, m = tuple(nc)
+            if b==0:
+                r = (m*t+l*(1-a*t))*exp(-a*t)
+            else:
+                hyp = False
+                if b.is_negative:
+                    b=-b
+                    hyp = True
+                b2 = list(roots(x_**2-b, x_).keys())[0]
+                bs = sqrt(b).simplify()
+                if hyp:
+                    r = l*exp(-a*t)*cosh(b2*t) + (m-a*l)/bs*exp(-a*t)*sinh(bs*t)
+                else:
+                    r = l*exp(-a*t)*cos(b2*t) + (m-a*l)/bs*exp(-a*t)*sin(bs*t)
+            terms_t.append(Heaviside(t)*r)
+        else:
+            ft, cond = _inverse_laplace_transform(fn, s, t, plane,
+                                               simplify=True, dorational=False)
+            terms_t.append(ft)
+            conditions.append(cond)
+
+    result = Add(*terms_t)
+    if simplify:
+        result = result.simplify(doit=False)
+    debug('[ILT _i_l_r]   returns %s'%(result,))
+    return result, And(*conditions)
+
+
+def _inverse_laplace_transform(fn, s_, t_, plane, simplify=True, dorational=True):
+    """
+    Front-end function of the inverse Laplace transform. It tries to apply all
+    known rules recursively.  If everything else fails, it tries to integrate.
+    """
+    terms = Add.make_args(fn)
+    terms_t = []
+    conditions = []
+
+    debug('[ILT _i_l_t] (%s, %s, %s)'%(fn, s_, t_))
+
+    for term in terms:
+        k, f = term.as_independent(s_, as_Add=False)
+        if dorational and term.is_rational_function(s_) and \
+            (r := _inverse_laplace_rational(f, s_, t_, plane, simplify)) is not None:
+            pass
+        elif (r := _inverse_laplace_apply_simple_rules(f, s_, t_)) is not None:
+            pass
+        elif (r := _inverse_laplace_apply_prog_rules(f, s_, t_, plane)) is not None:
+            pass
+        elif (r := _inverse_laplace_expand(f, s_, t_, plane)) is not None:
+            pass
+        elif any(undef.has(s_) for undef in f.atoms(AppliedUndef)):
+            # If there are undefined functions f(t) then integration is
+            # unlikely to do anything useful so we skip it and given an
+            # unevaluated LaplaceTransform.
+            r = (InverseLaplaceTransform(f, s_, t_, plane), S.true)
+        elif (r := _inverse_laplace_transform_integration(f, s_, t_, plane,
+                                      simplify=simplify)) is not None:
+            pass
+        else:
+            r = (InverseLaplaceTransform(f, s_, t_, plane), S.true)
+        (ri_, ci_) = r
+        terms_t.append(k*ri_)
+        conditions.append(ci_)
+
+    result = Add(*terms_t)
+    if simplify:
+        result = result.simplify(doit=False)
+    condition = And(*conditions)
+
+    return result, condition
+
+
 class InverseLaplaceTransform(IntegralTransform):
     """
     Class representing unevaluated inverse Laplace transforms.
@@ -1116,12 +1363,45 @@ def fundamental_plane(self):
         return plane
 
     def _compute_transform(self, F, s, t, **hints):
-        return _inverse_laplace_transform(F, s, t, self.fundamental_plane, **hints)
+        return _inverse_laplace_transform_integration(F, s, t,
+                                            self.fundamental_plane, **hints)
 
     def _as_integral(self, F, s, t):
         c = self.__class__._c
         return Integral(exp(s*t)*F, (s, c - S.ImaginaryUnit*S.Infinity,
-                                     c + S.ImaginaryUnit*S.Infinity))/(2*S.Pi*S.ImaginaryUnit)
+                    c + S.ImaginaryUnit*S.Infinity))/(2*S.Pi*S.ImaginaryUnit)
+
+    def doit(self, **hints):
+        """
+        Try to evaluate the transform in closed form.
+
+        Explanation
+        ===========
+
+        Standard hints are the following:
+        - ``noconds``:  if True, do not return convergence conditions. The
+        default setting is `True`.
+        - ``simplify``: if True, it simplifies the final result. This is the
+        default behaviour
+        """
+        _noconds = hints.get('noconds', True)
+        _simplify = hints.get('simplify', True)
+
+        debug('[ILT doit] (%s, %s, %s)'%(self.function,
+                                         self.function_variable,
+                                         self.transform_variable))
+
+        s_ = self.function_variable
+        t_ = self.transform_variable
+        fn = self.function
+        plane = self.fundamental_plane
+
+        r = _inverse_laplace_transform(fn, s_, t_, plane, simplify=_simplify)
+
+        if _noconds:
+            return r[0]
+        else:
+            return r
 
 
 def inverse_laplace_transform(F, s, t, plane=None, **hints):
diff --git a/sympy/polys/matrices/_typing.py b/sympy/polys/matrices/_typing.py
deleted file mode 100644
index 45dfe80647..0000000000
--- a/sympy/polys/matrices/_typing.py
+++ /dev/null
@@ -1,12 +0,0 @@
-from typing import TypeVar, Protocol
-
-
-T = TypeVar('T')
-
-
-class RingElement(Protocol):
-    def __add__(self: T, other: T, /) -> T: ...
-    def __sub__(self: T, other: T, /) -> T: ...
-    def __mul__(self: T, other: T, /) -> T: ...
-    def __pow__(self: T, other: int, /) -> T: ...
-    def __neg__(self: T, /) -> T: ...
diff --git a/sympy/polys/matrices/dense.py b/sympy/polys/matrices/dense.py
index 6c56c77f83..26411af968 100644
--- a/sympy/polys/matrices/dense.py
+++ b/sympy/polys/matrices/dense.py
@@ -34,62 +34,55 @@
     [[1, 0, -1], [0, 1, 2]]
 
 """
-from __future__ import annotations
 from operator import mul
+
 from .exceptions import (
     DMShapeError,
     DMNonInvertibleMatrixError,
     DMNonSquareMatrixError,
-)
-from typing import Sequence, TypeVar
-from sympy.polys.matrices._typing import RingElement
-
-
-T = TypeVar('T')
-R = TypeVar('R', bound=RingElement)
+    )
 
 
-def ddm_transpose(matrix: Sequence[Sequence[T]]) -> list[list[T]]:
+def ddm_transpose(a):
     """matrix transpose"""
-    return list(map(list, zip(*matrix)))
+    aT = list(map(list, zip(*a)))
+    return aT
 
 
-def ddm_iadd(a: list[list[R]], b: Sequence[Sequence[R]]) -> None:
+def ddm_iadd(a, b):
     """a += b"""
     for ai, bi in zip(a, b):
         for j, bij in enumerate(bi):
             ai[j] += bij
 
 
-def ddm_isub(a: list[list[R]], b: Sequence[Sequence[R]]) -> None:
+def ddm_isub(a, b):
     """a -= b"""
     for ai, bi in zip(a, b):
         for j, bij in enumerate(bi):
             ai[j] -= bij
 
 
-def ddm_ineg(a: list[list[R]]) -> None:
+def ddm_ineg(a):
     """a  <--  -a"""
     for ai in a:
         for j, aij in enumerate(ai):
             ai[j] = -aij
 
 
-def ddm_imul(a: list[list[R]], b: R) -> None:
+def ddm_imul(a, b):
     for ai in a:
         for j, aij in enumerate(ai):
             ai[j] = aij * b
 
 
-def ddm_irmul(a: list[list[R]], b: R) -> None:
+def ddm_irmul(a, b):
     for ai in a:
         for j, aij in enumerate(ai):
             ai[j] = b * aij
 
 
-def ddm_imatmul(
-    a: list[list[R]], b: Sequence[Sequence[R]], c: Sequence[Sequence[R]]
-) -> None:
+def ddm_imatmul(a, b, c):
     """a += b @ c"""
     cT = list(zip(*c))
 


Test patch:
diff --git a/sympy/integrals/tests/test_laplace.py b/sympy/integrals/tests/test_laplace.py
index 703743c673..1dfaf0d6fe 100644
--- a/sympy/integrals/tests/test_laplace.py
+++ b/sympy/integrals/tests/test_laplace.py
@@ -5,7 +5,6 @@
 from sympy.core import EulerGamma, Subs, Derivative, diff
 from sympy.core.exprtools import factor_terms
 from sympy.core.numbers import I, oo, pi
-from sympy.core.relational import Eq
 from sympy.core.singleton import S
 from sympy.core.symbol import Symbol, symbols
 from sympy.simplify.simplify import simplify
@@ -13,7 +12,7 @@
 from sympy.functions.elementary.exponential import exp, log, exp_polar
 from sympy.functions.elementary.hyperbolic import cosh, sinh, coth, asinh
 from sympy.functions.elementary.miscellaneous import sqrt
-from sympy.functions.elementary.trigonometric import atan, atan2, cos, sin
+from sympy.functions.elementary.trigonometric import atan, cos, sin
 from sympy.functions.special.gamma_functions import lowergamma, gamma
 from sympy.functions.special.delta_functions import DiracDelta, Heaviside
 from sympy.functions.special.zeta_functions import lerchphi
@@ -364,9 +363,11 @@ def test_laplace_transform():
 
 @slow
 def test_inverse_laplace_transform():
+    from sympy.functions.special.delta_functions import DiracDelta
     ILT = inverse_laplace_transform
-    a, b, c, = symbols('a b c', positive=True)
-    t = symbols('t')
+    a, b, c, d = symbols('a b c d', positive=True)
+    r = Symbol('r', real=True)
+    t, z = symbols('t z')
 
     def simp_hyp(expr):
         return factor_terms(expand_mul(expr)).rewrite(sin)
@@ -375,13 +376,14 @@ def simp_hyp(expr):
     assert ILT(1/s, s, t) == Heaviside(t)
     assert ILT(a/(a + s), s, t) == a*exp(-a*t)*Heaviside(t)
     assert ILT(s/(a + s), s, t) == -a*exp(-a*t)*Heaviside(t) + DiracDelta(t)
+    assert ILT((s-a)**(-b), s, t) == t**(b - 1)*exp(a*t)*Heaviside(t)/gamma(b)
     assert ILT((a + s)**(-2), s, t) == t*exp(-a*t)*Heaviside(t)
     assert ILT((a + s)**(-5), s, t) == t**4*exp(-a*t)*Heaviside(t)/24
     assert ILT(a/(a**2 + s**2), s, t) == sin(a*t)*Heaviside(t)
     assert ILT(s/(s**2 + a**2), s, t) == cos(a*t)*Heaviside(t)
     assert ILT(b/(b**2 + (a + s)**2), s, t) == exp(-a*t)*sin(b*t)*Heaviside(t)
-    assert ILT(b*s/(b**2 + (a + s)**2), s, t) +\
-        (a*sin(b*t) - b*cos(b*t))*exp(-a*t)*Heaviside(t) == 0
+    assert ILT(b*s/(b**2 + (a + s)**2), s, t) == \
+        (-a*sin(b*t) + b*cos(b*t))*exp(-a*t)*Heaviside(t)
     assert ILT(exp(-a*s)/s, s, t) == Heaviside(-a + t)
     assert ILT(exp(-a*s)/(b + s), s, t) == exp(b*(a - t))*Heaviside(-a + t)
     assert ILT((b + s)/(a**2 + (b + s)**2), s, t) == \
@@ -396,6 +398,9 @@ def simp_hyp(expr):
     assert ILT(1 - 1/(s**2 + 1), s, t) == -sin(t)*Heaviside(t) + DiracDelta(t)
     assert ILT(1/s**2, s, t) == t*Heaviside(t)
     assert ILT(1/s**5, s, t) == t**4*Heaviside(t)/24
+    # Issue #24424
+    assert ILT((s + 8)/((s + 2)*(s**2 + 2*s + 10)), s, t) ==\
+        ((8*sin(3*t) - 9*cos(3*t))*exp(t) + 9)*exp(-2*t)*Heaviside(t)/15
     assert simp_hyp(ILT(a/(s**2 - a**2), s, t)) == sinh(a*t)*Heaviside(t)
     assert simp_hyp(ILT(s/(s**2 - a**2), s, t)) == cosh(a*t)*Heaviside(t)
     # TODO sinh/cosh shifted come out a mess. also delayed trig is a mess
@@ -419,42 +424,40 @@ def simp_hyp(expr):
 
     assert ILT( (s * eye(2) - Matrix([[1, 0], [0, 2]])).inv(), s, t) ==\
         Matrix([[exp(t)*Heaviside(t), 0], [0, exp(2*t)*Heaviside(t)]])
-
-
-def test_inverse_laplace_transform_delta():
-    from sympy.functions.special.delta_functions import DiracDelta
-    ILT = inverse_laplace_transform
-    t = symbols('t')
+    # New tests for rules
+    assert ILT(b/(s**2-a**2), s, t) == b*sinh(a*t)*Heaviside(t)/a
+    assert ILT(b*s/(s**2-a**2), s, t) == b*cosh(a*t)*Heaviside(t)
+    assert ILT(b/(s*(s+a)), s, t) == b*(exp(a*t) - 1)*exp(-a*t)*Heaviside(t)/a
+    assert ILT(b*s/(s+a)**2, s, t) == b*(-a*t + 1)*exp(-a*t)*Heaviside(t)
+    assert ILT(c/((s+a)*(s+b)), s, t) ==\
+        c*(exp(a*t) - exp(b*t))*exp(-t*(a + b))*Heaviside(t)/(a - b)
+    assert ILT(c*s/((s+a)*(s+b)), s, t) ==\
+        c*(a*exp(b*t) - b*exp(a*t))*exp(-t*(a + b))*Heaviside(t)/(a - b)
+    assert ILT(c*s/(d**2*(s+a)**2+b**2), s, t) ==\
+        c*(-a*d*sin(b*t/d) + b*cos(b*t/d))*exp(-a*t)*Heaviside(t)/(b*d**2)
+    # Rules for testing different DiracDelta cases
     assert ILT(2, s, t) == 2*DiracDelta(t)
     assert ILT(2*exp(3*s) - 5*exp(-7*s), s, t) == \
-        2*DiracDelta(t + 3) - 5*DiracDelta(t - 7)
+        2*InverseLaplaceTransform(exp(3*s), s, t, None) - 5*DiracDelta(t - 7)
     a = cos(sin(7)/2)
     assert ILT(a*exp(-3*s), s, t) == a*DiracDelta(t - 3)
-    assert ILT(exp(2*s), s, t) == DiracDelta(t + 2)
+    assert ILT(exp(2*s), s, t) == InverseLaplaceTransform(exp(2*s), s, t, None)
     r = Symbol('r', real=True)
-    assert ILT(exp(r*s), s, t) == DiracDelta(t + r)
-
-
-def test_inverse_laplace_transform_delta_cond():
-    from sympy.functions.elementary.complexes import im
-    from sympy.functions.special.delta_functions import DiracDelta
-    ILT = inverse_laplace_transform
-    t = symbols('t')
-    r = Symbol('r', real=True)
-    assert ILT(exp(r*s), s, t, noconds=False) == (DiracDelta(t + r), True)
-    z = Symbol('z')
-    assert ILT(exp(z*s), s, t, noconds=False) == \
-        (DiracDelta(t + z), Eq(im(z), 0))
+    assert ILT(exp(r*s), s, t) == InverseLaplaceTransform(exp(r*s), s, t, None)
+    # Rules from the previous test_inverse_laplace_transform_delta_cond():
+    assert ILT(exp(r*s), s, t, noconds=False) ==\
+        (InverseLaplaceTransform(exp(r*s), s, t, None), True)
     # inversion does not exist: verify it doesn't evaluate to DiracDelta
     for z in (Symbol('z', extended_real=False),
               Symbol('z', imaginary=True, zero=False)):
         f = ILT(exp(z*s), s, t, noconds=False)
         f = f[0] if isinstance(f, tuple) else f
         assert f.func != DiracDelta
-    # issue 15043
-    assert ILT(1/s + exp(r*s)/s, s, t, noconds=False) == (
-        Heaviside(t) + Heaviside(r + t), True)
-
+    # old test for Issue 8514, is not important anymore since this function
+    # is not solved by integration anymore
+    assert ILT(1/(a*s**2+b*s+c),s, t) ==\
+        2*exp(-b*t/(2*a))*sin(t*sqrt(4*a*c - b**2)/(2*a))*\
+            Heaviside(t)/sqrt(4*a*c - b**2)
 
 @slow
 def test_expint():
@@ -478,19 +481,3 @@ def test_expint():
     assert inverse_laplace_transform((s - log(s + 1))/s**2, s,
                 x).rewrite(expint).expand() == \
         (expint(2, x)*Heaviside(x)).rewrite(Ei).rewrite(expint).expand()
-
-
-@slow
-def test_issue_8514():
-    a, b, c, = symbols('a b c', positive=True)
-    t = symbols('t', positive=True)
-    ft = simplify(inverse_laplace_transform(1/(a*s**2+b*s+c),s, t))
-    assert ft == (I*exp(t*cos(atan2(0, -4*a*c + b**2)/2)*sqrt(Abs(4*a*c -
-                  b**2))/a)*sin(t*sin(atan2(0, -4*a*c + b**2)/2)*sqrt(Abs(
-                  4*a*c - b**2))/(2*a)) + exp(t*cos(atan2(0, -4*a*c + b**2)
-                  /2)*sqrt(Abs(4*a*c - b**2))/a)*cos(t*sin(atan2(0, -4*a*c
-                  + b**2)/2)*sqrt(Abs(4*a*c - b**2))/(2*a)) + I*sin(t*sin(
-                  atan2(0, -4*a*c + b**2)/2)*sqrt(Abs(4*a*c - b**2))/(2*a))
-                  - cos(t*sin(atan2(0, -4*a*c + b**2)/2)*sqrt(Abs(4*a*c -
-                  b**2))/(2*a)))*exp(-t*(b + cos(atan2(0, -4*a*c + b**2)/2)
-                  *sqrt(Abs(4*a*c - b**2)))/(2*a))/sqrt(-4*a*c + b**2)