Merge pull request sympy#23608 from eagleoflqj/manualintegrate_poly_d…

…iv_sqrt_quadratic

manualintegrate poly/sqrt(a+b*x+c*x**2)

sympy/core/basic.py

@@ -116,6 +116,8 @@ class Basic(Printable, metaclass=ManagedProperties):
     is_Point = False
     is_MatAdd = False
     is_MatMul = False
+    is_real: bool | None
+    is_zero: bool | None
     is_negative: bool | None
     is_commutative: bool | None
 

sympy/integrals/manualintegrate.py

@@ -88,6 +88,7 @@ def __eq__(self, other):
 ArcsinRule = Rule("ArcsinRule")
 ArcsinhRule = Rule("ArcsinhRule")
 ReciprocalSqrtQuadraticRule = Rule("ReciprocalSqrtQuadraticRule", "a b c")
+SqrtQuadraticDenomRule = Rule("SqrtQuadraticDenomRule", "a b c coeffs")
 SqrtQuadraticRule = Rule("SqrtQuadraticRule", "a b c")
 AlternativeRule = Rule("AlternativeRule", "alternatives")
 DontKnowRule = Rule("DontKnowRule")
@@ -136,19 +137,22 @@ def _evaluates(func):
         return func
     return _evaluates
 
+
 def contains_dont_know(rule):
     if isinstance(rule, DontKnowRule):
         return True
-    else:
-        for val in rule:
-            if isinstance(val, tuple):
-                if contains_dont_know(val):
-                    return True
-            elif isinstance(val, list):
-                if any(contains_dont_know(i) for i in val):
-                    return True
+    if not isinstance(rule, tuple):
+        return False
+    for val in rule:
+        if isinstance(val, tuple):
+            if contains_dont_know(val):
+                return True
+        elif isinstance(val, list):
+            if any(contains_dont_know(i) for i in val):
+                return True
     return False
 
+
 def manual_diff(f, symbol):
     """Derivative of f in form expected by find_substitutions
 
@@ -444,7 +448,26 @@ def special_function_rule(integral):
                     return p[3](*args)
 
 
-def inverse_trig_rule(integral):
+def _add_degenerate_step(generic_cond, generic_step, degenerate_step):
+    if degenerate_step is None:
+        return generic_step
+    if isinstance(generic_step, PiecewiseRule):
+        subfunctions = [(substep, (cond & generic_cond).simplify())
+                        for substep, cond in generic_step.subfunctions]
+    else:
+        subfunctions = [(generic_step, generic_cond)]
+    if isinstance(degenerate_step, PiecewiseRule):
+        subfunctions += degenerate_step.subfunctions
+    else:
+        subfunctions.append((degenerate_step, S.true))
+    return PiecewiseRule(subfunctions, generic_step.context, generic_step.symbol)
+
+
+def inverse_trig_rule(integral: IntegralInfo, degenerate=True):
+    """
+    Set degenerate=False on recursive call where coefficient of quadratic term
+    is assumed non-zero.
+    """
     integrand, symbol = integral
     base, exp = integrand.as_base_exp()
     a = Wild('a', exclude=[symbol])
@@ -475,10 +498,17 @@ def make_inverse_trig(RuleClass, a, sign_a, c, sign_c, h):
         return substep
 
     a, b, c = [match.get(i, S.Zero) for i in (a, b, c)]
+    generic_cond = Ne(c, 0)
+    if not degenerate or generic_cond is S.true:
+        degenerate_step = None
+    elif b.is_zero:
+        degenerate_step = ConstantRule(a ** exp, a ** exp, symbol)
+    else:
+        degenerate_step = sqrt_linear_rule(IntegralInfo((a + b * symbol) ** exp, symbol))
 
     if simplify(2*exp + 1) == 0:
         h, k = -b/(2*c), a - b**2/(4*c)  # rewrite base to k + c*(symbol-h)**2
-        general_rule = ReciprocalSqrtQuadraticRule(a, b, c, integrand, symbol)
+        step = general_rule = ReciprocalSqrtQuadraticRule(a, b, c, integrand, symbol)
         if k.is_real and c.is_real:
             # list of ((rule, base_exp, a, sign_a, b, sign_b), condition)
             possibilities = []
@@ -494,10 +524,13 @@ def make_inverse_trig(RuleClass, a, sign_a, c, sign_c, h):
                 rules = [(make_inverse_trig(*args), cond) for args, cond in possibilities]
                 if not k.is_positive:  # conditions are not thorough, need fall back rule
                     rules.append((general_rule, S.true))
-                return PiecewiseRule(rules, integrand, symbol)
-        return general_rule
+                step = PiecewiseRule(rules, integrand, symbol)
+            else:
+                step = general_rule
+        return _add_degenerate_step(generic_cond, step, degenerate_step)
     if exp == S.Half:
-        return SqrtQuadraticRule(a, b, c, integrand, symbol)
+        step = SqrtQuadraticRule(a, b, c, integrand, symbol)
+        return _add_degenerate_step(generic_cond, step, degenerate_step)
 
 
 def add_rule(integral):
@@ -865,26 +898,36 @@ def sqrt_linear_rule(integral: IntegralInfo):
     substep = integral_steps(substituted*u**(q0-1)*q0/b0, u)
     if not contains_dont_know(substep):
         step = URule(u, u_x, None, substep, integrand, x)
-        generate_cond = Ne(b0, 0)
-        if generate_cond is not S.true:  # possible degenerate case
+        generic_cond = Ne(b0, 0)
+        if generic_cond is not S.true:  # possible degenerate case
             simplified = integrand.subs({b: 0 for b in bs})
             degenerate_step = integral_steps(simplified, x)
-            step = PiecewiseRule([(step, generate_cond), (degenerate_step, S.true)], integrand, x)
+            step = PiecewiseRule([(step, generic_cond), (degenerate_step, S.true)], integrand, x)
         return step
 
 
-def sqrt_quadratic_denom_rule(integral: IntegralInfo):
+def sqrt_quadratic_denom_rule(integral: IntegralInfo, degenerate=True):
     integrand, x = integral
-    numer, denom = integrand.as_numer_denom()
     a = Wild('a', exclude=[x])
     b = Wild('b', exclude=[x])
     c = Wild('c', exclude=[x, 0])
-    match = denom.match(sqrt(a+b*x+c*x**2))
+    numer = Wild('numer')
+    match = integrand.match(numer/sqrt(a+b*x+c*x**2))
     if match:
+        a, b, c, numer = match[a], match[b], match[c], match[numer]
         numer_poly = numer.as_poly(x)
         if numer_poly is None:
             return
-        a, b, c = match[a], match[b], match[c]
+
+        generic_cond = Ne(c, 0)
+        if not degenerate or generic_cond is S.true:
+            degenerate_step = None
+        elif b.is_zero:
+            degenerate_step = integral_steps(numer/sqrt(a), x)
+        else:
+            degenerate_step = sqrt_linear_rule(IntegralInfo(numer/sqrt(a+b*x), x))
+
+        denom = sqrt(a+b*x+c*x**2)
         deg = numer_poly.degree()
         if deg <= 1:
             # integrand == (d+e*x)/sqrt(a+b*x+c*x**2)
@@ -901,13 +944,18 @@ def sqrt_quadratic_denom_rule(integral: IntegralInfo):
                 if A != 1:
                     linear_step = ConstantTimesRule(A, pre_substitute, linear_step, A*pre_substitute, x)
             if B != 0:
-                constant_step = inverse_trig_rule(IntegralInfo(1/denom, x))
+                constant_step = inverse_trig_rule(IntegralInfo(1/denom, x), degenerate=False)
                 if B != 1:
                     constant_step = ConstantTimesRule(B, 1/denom, constant_step, B/denom, x)
             if linear_step and constant_step:
                 add = Add(A*pre_substitute, B/denom, evaluate=False)
-                return RewriteRule(add, AddRule([linear_step, constant_step], add, x), integrand, x)
-            return linear_step or constant_step
+                step = RewriteRule(add, AddRule([linear_step, constant_step], add, x), integrand, x)
+            else:
+                step = linear_step or constant_step
+        else:
+            coeffs = numer_poly.all_coeffs()
+            step = SqrtQuadraticDenomRule(a, b, c, coeffs, integrand, x)
+        return _add_degenerate_step(generic_cond, step, degenerate_step)
 
 
 def hyperbolic_rule(integral: tuple[Expr, Symbol]):
@@ -1539,9 +1587,38 @@ def eval_reciprocal_sqrt_quadratic(a, b, c, integrand, x):
     return log(2*sqrt(c)*sqrt(a+b*x+c*x**2)+b+2*c*x)/sqrt(c)
 
 
+@evaluates(SqrtQuadraticDenomRule)
+def eval_sqrt_quadratic_denom(a, b, c, coeffs: list[Expr], integrand, x):
+    # Integrate poly/sqrt(a+b*x+c*x**2) using recursion.
+    # coeffs are coefficients of the polynomial.
+    # Let I_n = x**n/sqrt(a+b*x+c*x**2), then
+    # I_n = A * x**(n-1)*sqrt(a+b*x+c*x**2) - B * I_{n-1} - C * I_{n-2}
+    # where A = 1/(n*c), B = (2*n-1)*b/(2*n*c), C = (n-1)*a/(n*c)
+    # See https://github.com/sympy/sympy/pull/23608 for proof.
+    result_coeffs = []
+    coeffs = coeffs.copy()
+    for i in range(len(coeffs)-2):
+        n = len(coeffs)-1-i
+        coeff = coeffs[i]/(c*n)
+        result_coeffs.append(coeff)
+        coeffs[i+1] -= (2*n-1)*b/2*coeff
+        coeffs[i+2] -= (n-1)*a*coeff
+    d, e = coeffs[-1], coeffs[-2]
+    s = sqrt(a+b*x+c*x**2)
+    constant = d-b*e/(2*c)
+    if constant == 0:
+        I0 = 0
+    else:
+        step = inverse_trig_rule(IntegralInfo(1/s, x), degenerate=False)
+        I0 = constant*_manualintegrate(step)
+    return Add(*(result_coeffs[i]*x**(len(coeffs)-2-i)
+                 for i in range(len(result_coeffs))), e/c)*s + I0
+
+
 @evaluates(SqrtQuadraticRule)
 def eval_sqrt_quadratic(a, b, c, integrand, x):
-    return x*integrand/2 + _manualintegrate(sqrt_quadratic_denom_rule(IntegralInfo((2*a+b*x)/integrand, x)))/4
+    step = sqrt_quadratic_denom_rule(IntegralInfo((2*a+b*x)/integrand, x), degenerate=False)
+    return x*integrand/2 + _manualintegrate(step)/4
 
 
 @evaluates(AlternativeRule)

sympy/integrals/tests/test_integrals.py

@@ -1676,7 +1676,6 @@ def test_issue_13112():
     assert integrate(sin(t)**2 / (5 - 4*cos(t)), [t, 0, 2*pi]) == pi / 4
 
 
-@slow
 def test_issue_14709b():
     h = Symbol('h', positive=True)
     i = integrate(x*acos(1 - 2*x/h), (x, 0, h))
@@ -1948,6 +1947,12 @@ def test_issue_11254d():
     assert integrate((sech(x)**2).rewrite(sinh), x) == 2*tanh(x/2)/(tanh(x/2)**2 + 1)
 
 
+def test_issue_22863():
+    i = integrate((3*x**3-x**2+2*x-4)/sqrt(x**2-3*x+2), (x, 0, 1))
+    assert i == -101*sqrt(2)/8 - 135*log(3 - 2*sqrt(2))/16
+    assert math.isclose(i.n(), -2.98126694400554)
+
+
 def test_hyperbolic():
     assert integrate(coth(x)) == x - log(tanh(x) + 1) + log(tanh(x))
     assert integrate(sech(x)) == 2*atan(tanh(x/2))
@@ -1958,7 +1963,9 @@ def test_sqrt_quadratic():
     assert integrate(1/sqrt(3*x**2+4*x+5)) == sqrt(3)*asinh(3*sqrt(11)*(x + S(2)/3)/11)/3
     assert integrate(1/sqrt(-3*x**2+4*x+5)) == sqrt(3)*asin(3*sqrt(19)*(x - S(2)/3)/19)/3
     assert integrate(1/sqrt(3*x**2+4*x-5)) == sqrt(3)*log(6*x + 2*sqrt(3)*sqrt(3*x**2 + 4*x - 5) + 4)/3
-    assert integrate(1/sqrt(a+b*x+c*x**2), x) == log(2*sqrt(c)*sqrt(a+b*x+c*x**2)+b+2*c*x)/sqrt(c)
+    assert integrate(1/sqrt(a+b*x+c*x**2), x) == \
+        Piecewise((log(b + 2*sqrt(c)*sqrt(a + b*x + c*x**2) + 2*c*x)/sqrt(c), Ne(c, 0)),
+                  (2*sqrt(a + b*x)/b, Ne(b, 0)), (x/sqrt(a), True))
 
     assert integrate((7*x+6)/sqrt(3*x**2+4*x+5)) == \
            7*sqrt(3*x**2 + 4*x + 5)/3 + 4*sqrt(3)*asinh(3*sqrt(11)*(x + S(2)/3)/11)/9
@@ -1967,10 +1974,19 @@ def test_sqrt_quadratic():
     assert integrate((7*x+6)/sqrt(3*x**2+4*x-5)) == \
            7*sqrt(3*x**2 + 4*x - 5)/3 + 4*sqrt(3)*log(6*x + 2*sqrt(3)*sqrt(3*x**2 + 4*x - 5) + 4)/9
     assert integrate((d+e*x)/sqrt(a+b*x+c*x**2), x) == \
-           e*sqrt(a + b*x + c*x**2)/c + (-b*e/(2*c) + d)*log(b + 2*sqrt(c)*sqrt(a + b*x + c*x**2) + 2*c*x)/sqrt(c)
+        Piecewise((e*sqrt(a + b*x + c*x**2)/c +
+                   (-b*e/(2*c) + d)*log(b + 2*sqrt(c)*sqrt(a + b*x + c*x**2) + 2*c*x)/sqrt(c), Ne(c, 0)),
+                  ((2*d*sqrt(a + b*x) + 2*e*(-a*sqrt(a + b*x) + (a + b*x)**(S(3)/2)/3)/b)/b, Ne(b, 0)),
+                  ((d*x + e*x**2/2)/sqrt(a), True))
+
+    assert integrate((3*x**3-x**2+2*x-4)/sqrt(x**2-3*x+2)) == \
+           sqrt(x**2 - 3*x + 2)*(x**2 + 13*x/4 + S(101)/8) + 135*log(2*x + 2*sqrt(x**2 - 3*x + 2) - 3)/16
 
     assert integrate(sqrt(53225*x**2-66732*x+23013)) == \
            x*sqrt(53225*x**2 - 66732*x + 23013)/2 - 16683*sqrt(53225*x**2 - 66732*x + 23013)/53225 + \
            111576969*sqrt(2129)*asinh(53225*x/10563 - S(11122)/3521)/1133160250
-    assert integrate(sqrt(a+b*x+c*x**2), x) == b*sqrt(a + b*x + c*x**2)/(4*c) + x*sqrt(a + b*x + c*x**2)/2 + \
-           (2*a - b**2/(2*c))*log(b + 2*sqrt(c)*sqrt(a + b*x + c*x**2) + 2*c*x)/(4*sqrt(c))
+    assert integrate(sqrt(a+b*x+c*x**2), x) == \
+        Piecewise((b*sqrt(a + b*x + c*x**2)/(4*c) + x*sqrt(a + b*x + c*x**2)/2 +
+                   (2*a - b**2/(2*c))*log(b + 2*sqrt(c)*sqrt(a + b*x + c*x**2) + 2*c*x)/(4*sqrt(c)), Ne(c, 0)),
+                  (2*(a + b*x)**(S(3)/2)/(3*b), Ne(b, 0)),
+                  (sqrt(a)*x, True))

sympy/integrals/tests/test_manual.py

@@ -172,7 +172,7 @@ def test_manualintegrate_inversetrig():
     assert manualintegrate(1/sqrt(4*x**2 + 1), x) == \
         asinh(2*x)/2
     assert manualintegrate(1/sqrt(ra*x**2 + 1), x) == \
-        Piecewise((asin(x*sqrt(-ra))/sqrt(-ra), ra < 0), (asinh(sqrt(ra)*x)/sqrt(ra), ra > 0))
+        Piecewise((asin(x*sqrt(-ra))/sqrt(-ra), ra < 0), (asinh(sqrt(ra)*x)/sqrt(ra), ra > 0), (x, True))
     assert manualintegrate(1/sqrt(ra + x**2), x) == \
         Piecewise((asinh(x*sqrt(1/ra)), ra > 0), (log(2*x + 2*sqrt(ra + x**2)), True))
 
@@ -181,19 +181,26 @@ def test_manualintegrate_inversetrig():
     assert manualintegrate(1/sqrt(x**2 - 4), x) == log(2*x + 2*sqrt(x**2 - 4))
     assert manualintegrate(1/sqrt(4*x**2 - 4), x) == log(8*x + 4*sqrt(4*x**2 - 4))/2
     assert manualintegrate(1/sqrt(9*x**2 - 1), x) == log(18*x + 6*sqrt(9*x**2 - 1))/3
-    assert manualintegrate(1/sqrt(ra*x**2 - 4), x) == log(2*sqrt(ra)*sqrt(ra*x**2 - 4) + 2*ra*x)/sqrt(ra)
+    assert manualintegrate(1/sqrt(ra*x**2 - 4), x) == \
+           Piecewise((log(2*sqrt(ra)*sqrt(ra*x**2 - 4) + 2*ra*x)/sqrt(ra), Ne(ra, 0)), (-I*x/2, True))
     assert manualintegrate(1/sqrt(-ra + 4*x**2), x) == \
         Piecewise((asinh(2*x*sqrt(-1/ra))/2, ra < 0), (log(8*x + 4*sqrt(-ra + 4*x**2))/2, True))
 
     # From https://www.wikiwand.com/en/List_of_integrals_of_inverse_trigonometric_functions
     # asin
     assert manualintegrate(asin(x), x) == x*asin(x) + sqrt(1 - x**2)
     assert manualintegrate(asin(a*x), x) == Piecewise(((a*x*asin(a*x) + sqrt(-a**2*x**2 + 1))/a, Ne(a, 0)), (0, True))
-    assert manualintegrate(x*asin(a*x), x) == -a*Integral(x**2/sqrt(-a**2*x**2 + 1), x)/2 + x**2*asin(a*x)/2
+    assert manualintegrate(x*asin(a*x), x) == \
+           -a*Piecewise((-x*sqrt(-a**2*x**2 + 1)/(2*a**2) +
+                         log(-2*a**2*x + 2*sqrt(-a**2)*sqrt(-a**2*x**2 + 1))/(2*a**2*sqrt(-a**2)), Ne(a**2, 0)),
+                        (x**3/3, True))/2 + x**2*asin(a*x)/2
     # acos
     assert manualintegrate(acos(x), x) == x*acos(x) - sqrt(1 - x**2)
     assert manualintegrate(acos(a*x), x) == Piecewise(((a*x*acos(a*x) - sqrt(-a**2*x**2 + 1))/a, Ne(a, 0)), (pi*x/2, True))
-    assert manualintegrate(x*acos(a*x), x) == a*Integral(x**2/sqrt(-a**2*x**2 + 1), x)/2 + x**2*acos(a*x)/2
+    assert manualintegrate(x*acos(a*x), x) == \
+           a*Piecewise((-x*sqrt(-a**2*x**2 + 1)/(2*a**2) +
+                        log(-2*a**2*x + 2*sqrt(-a**2)*sqrt(-a**2*x**2 + 1))/(2*a**2*sqrt(-a**2)), Ne(a**2, 0)),
+                       (x**3/3, True))/2 + x**2*acos(a*x)/2
     # atan
     assert manualintegrate(atan(x), x) == x*atan(x) - log(x**2 + 1)/2
     assert manualintegrate(atan(a*x), x) == Piecewise(((a*x*atan(a*x) - log(a**2*x**2 + 1)/2)/a, Ne(a, 0)), (0, True))
@@ -215,11 +222,13 @@ def test_manualintegrate_inversetrig():
     assert manualintegrate(1/sqrt(ra-rb*x**2), x) == \
         Piecewise((asin(x*sqrt(rb/ra))/sqrt(rb), And(-rb < 0, ra > 0)),
                   (asinh(x*sqrt(-rb/ra))/sqrt(-rb), And(-rb > 0, ra > 0)),
-                  (log(-2*rb*x + 2*sqrt(-rb)*sqrt(ra - rb*x**2))/sqrt(-rb), True))
+                  (log(-2*rb*x + 2*sqrt(-rb)*sqrt(ra - rb*x**2))/sqrt(-rb), Ne(rb, 0)),
+                  (x/sqrt(ra), True))
     assert manualintegrate(1/sqrt(ra + rb*x**2), x) == \
         Piecewise((asin(x*sqrt(-rb/ra))/sqrt(-rb), And(ra > 0, rb < 0)),
                   (asinh(x*sqrt(rb/ra))/sqrt(rb), And(ra > 0, rb > 0)),
-                  (log(2*sqrt(rb)*sqrt(ra + rb*x**2) + 2*rb*x)/sqrt(rb), True))
+                  (log(2*sqrt(rb)*sqrt(ra + rb*x**2) + 2*rb*x)/sqrt(rb), Ne(rb, 0)),
+                  (x/sqrt(ra), True))
 
 
 def test_manualintegrate_trig_substitution():
@@ -604,7 +613,9 @@ def test_manualintegrate_sqrt_quadratic():
     assert_is_integral_of(1/sqrt(3*x**2+4*x+5), sqrt(3)*asinh(3*sqrt(11)*(x + S(2)/3)/11)/3)
     assert_is_integral_of(1/sqrt(-3*x**2+4*x+5), sqrt(3)*asin(3*sqrt(19)*(x - S(2)/3)/19)/3)
     assert_is_integral_of(1/sqrt(3*x**2+4*x-5), sqrt(3)*log(6*x + 2*sqrt(3)*sqrt(3*x**2 + 4*x - 5) + 4)/3)
-    assert manualintegrate(1/sqrt(a+b*x+c*x**2), x) == log(2*sqrt(c)*sqrt(a+b*x+c*x**2)+b+2*c*x)/sqrt(c)
+    assert manualintegrate(1/sqrt(a+b*x+c*x**2), x) == \
+        Piecewise((log(b + 2*sqrt(c)*sqrt(a + b*x + c*x**2) + 2*c*x)/sqrt(c), Ne(c, 0)),
+                  (2*sqrt(a + b*x)/b, Ne(b, 0)), (x/sqrt(a), True))
 
     assert_is_integral_of((7*x+6)/sqrt(3*x**2+4*x+5),
                           7*sqrt(3*x**2 + 4*x + 5)/3 + 4*sqrt(3)*asinh(3*sqrt(11)*(x + S(2)/3)/11)/9)
@@ -613,12 +624,22 @@ def test_manualintegrate_sqrt_quadratic():
     assert_is_integral_of((7*x+6)/sqrt(3*x**2+4*x-5),
                           7*sqrt(3*x**2 + 4*x - 5)/3 + 4*sqrt(3)*log(6*x + 2*sqrt(3)*sqrt(3*x**2 + 4*x - 5) + 4)/9)
     assert manualintegrate((d+e*x)/sqrt(a+b*x+c*x**2), x) == \
-           e*sqrt(a + b*x + c*x**2)/c + (-b*e/(2*c) + d)*log(b + 2*sqrt(c)*sqrt(a + b*x + c*x**2) + 2*c*x)/sqrt(c)
+        Piecewise((e*sqrt(a + b*x + c*x**2)/c +
+                   (-b*e/(2*c) + d)*log(b + 2*sqrt(c)*sqrt(a + b*x + c*x**2) + 2*c*x)/sqrt(c), Ne(c, 0)),
+                  ((2*d*sqrt(a + b*x) + 2*e*(-a*sqrt(a + b*x) + (a + b*x)**(S(3)/2)/3)/b)/b, Ne(b, 0)),
+                  ((d*x + e*x**2/2)/sqrt(a), True))
+
+    assert manualintegrate((3*x**3-x**2+2*x-4)/sqrt(x**2-3*x+2), x) == \
+        sqrt(x**2 - 3*x + 2)*(x**2 + 13*x/4 + S(101)/8) + 135*log(2*x + 2*sqrt(x**2 - 3*x + 2) - 3)/16
 
     assert_is_integral_of(sqrt(53225*x**2-66732*x+23013),
                           x*sqrt(53225*x**2 - 66732*x + 23013)/2 - 16683*sqrt(53225*x**2 - 66732*x + 23013)/53225 +
                           111576969*sqrt(2129)*asinh(53225*x/10563 - S(11122)/3521)/1133160250)
     assert manualintegrate(sqrt(a+c*x**2), x) == \
-           x*sqrt(a + c*x**2)/2 + 2*a*log(2*sqrt(c)*sqrt(a + c*x**2) + 2*c*x)/(4*sqrt(c))
-    assert manualintegrate(sqrt(a+b*x+c*x**2), x) == b*sqrt(a + b*x + c*x**2)/(4*c) + x*sqrt(a + b*x + c*x**2)/2 + \
-           (2*a - b**2/(2*c))*log(b + 2*sqrt(c)*sqrt(a + b*x + c*x**2) + 2*c*x)/(4*sqrt(c))
+        Piecewise((a*log(2*sqrt(c)*sqrt(a + c*x**2) + 2*c*x)/(2*sqrt(c)) + x*sqrt(a + c*x**2)/2, Ne(c, 0)),
+                  (sqrt(a)*x, True))
+    assert manualintegrate(sqrt(a+b*x+c*x**2), x) == \
+        Piecewise((b*sqrt(a + b*x + c*x**2)/(4*c) + x*sqrt(a + b*x + c*x**2)/2 +
+                   (2*a - b**2/(2*c))*log(b + 2*sqrt(c)*sqrt(a + b*x + c*x**2) + 2*c*x)/(4*sqrt(c)), Ne(c, 0)),
+                  (2*(a + b*x)**(S(3)/2)/(3*b), Ne(b, 0)),
+                  (sqrt(a)*x, True))