diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IntegralPatterns.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IntegralPatterns.cs index ec0f921da..6211b028c 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IntegralPatterns.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IntegralPatterns.cs @@ -5,6 +5,8 @@ // Website: https://am.angouri.org. // +using PeterO.Numbers; + namespace AngouriMath.Functions.Algebra { internal static class IntegralPatterns @@ -19,17 +21,12 @@ internal static class IntegralPatterns TreeAnalyzer.TryGetPolyLinear(arg, x, out var a, out _) => MathS.Sin(arg) / a, - // By power reduction: sin(u)^2 = (1 - cos(2u)) / 2, so the integral is - // x/2 - sin(2u)/(4a). Without this, integrating sin(x)^2 fell through to - // integration by parts and cycled there. - Entity.Powf(Entity.Sinf(var arg), Entity.Number.Integer(2)) when - TreeAnalyzer.TryGetPolyLinear(arg, x, out var a, out _) => - x / 2 - MathS.Sin(2 * arg) / (4 * a), - - // cos(u)^2 = (1 + cos(2u)) / 2 - Entity.Powf(Entity.Cosf(var arg), Entity.Number.Integer(2)) when - TreeAnalyzer.TryGetPolyLinear(arg, x, out var a, out _) => - x / 2 + MathS.Sin(2 * arg) / (4 * a), + // sin(u)^n * cos(u)^m for whole n and m, which covers sin(x)^2 and cos(x)^2 + // as much as sin(x)^3 or sin(x)^2 * cos(x)^2. + _ when TryReadSineCosinePowers(expr, x, out var trigArg, out var sinePower, out var cosinePower) + && sinePower + cosinePower >= 2 + && TreeAnalyzer.TryGetPolyLinear(trigArg, x, out var trigRate, out _) => + IntegrateSineCosinePowers(trigArg, sinePower, cosinePower, trigRate, x), Entity.Secantf(var arg) when TreeAnalyzer.TryGetPolyLinear(arg, x, out var a, out _) => @@ -73,9 +70,247 @@ internal static class IntegralPatterns && TreeAnalyzer.TryGetPolyQuadratic(denominator, x, out var a, out var b, out var c) // ∫ k/(ax^2 + bx + c) dx => IntegrateRationalQuadratic(numerator, a, b, c, x), + // The inverse trigonometric functions, each of which is integration by parts + // against 1 -- a shape the by-parts solver does not look for, since there is no + // product to split. + Entity.Arcsinf(var arg) when + TreeAnalyzer.TryGetPolyLinear(arg, x, out var a, out _) => + (arg * MathS.Arcsin(arg) + MathS.Sqrt(1 - arg * arg)) / a, + + Entity.Arccosf(var arg) when + TreeAnalyzer.TryGetPolyLinear(arg, x, out var a, out _) => + (arg * MathS.Arccos(arg) - MathS.Sqrt(1 - arg * arg)) / a, + + Entity.Arctanf(var arg) when + TreeAnalyzer.TryGetPolyLinear(arg, x, out var a, out _) => + (arg * MathS.Arctan(arg) - MathS.Ln(MathS.Abs(1 + arg * arg)) / 2) / a, + + Entity.Arccotanf(var arg) when + TreeAnalyzer.TryGetPolyLinear(arg, x, out var a, out _) => + (arg * MathS.Arccotan(arg) + MathS.Ln(MathS.Abs(1 + arg * arg)) / 2) / a, + + // ∫ B^(px + q) * sin(mx + n) dx and its cosine twin. Integrating by parts + // twice returns the integral it started from, so the usual machinery cycles + // rather than terminating; solving that equation for the integral once gives + // the closed form below, which is what goes in the table. + Entity.Mulf(var exponential, Entity.Sinf(var wave)) when + IsExponentialRate(exponential, x, out var rate) + && TreeAnalyzer.TryGetPolyLinear(wave, x, out var frequency, out _) => + exponential * (rate * MathS.Sin(wave) - frequency * MathS.Cos(wave)) + / (rate * rate + frequency * frequency), + + Entity.Mulf(Entity.Sinf(var wave), var exponential) when + IsExponentialRate(exponential, x, out var rate) + && TreeAnalyzer.TryGetPolyLinear(wave, x, out var frequency, out _) => + exponential * (rate * MathS.Sin(wave) - frequency * MathS.Cos(wave)) + / (rate * rate + frequency * frequency), + + Entity.Mulf(var exponential, Entity.Cosf(var wave)) when + IsExponentialRate(exponential, x, out var rate) + && TreeAnalyzer.TryGetPolyLinear(wave, x, out var frequency, out _) => + exponential * (rate * MathS.Cos(wave) + frequency * MathS.Sin(wave)) + / (rate * rate + frequency * frequency), + + Entity.Mulf(Entity.Cosf(var wave), var exponential) when + IsExponentialRate(exponential, x, out var rate) + && TreeAnalyzer.TryGetPolyLinear(wave, x, out var frequency, out _) => + exponential * (rate * MathS.Cos(wave) + frequency * MathS.Sin(wave)) + / (rate * rate + frequency * frequency), + + // ∫ sqrt(ax^2 + bx + c) dx, which is one integration by parts away from the + // reciprocal form below and is written in terms of it. + Entity.Powf(var radicand, Entity.Number.Rational(Entity.Number.Integer(1), Entity.Number.Integer(2))) when + TreeAnalyzer.TryGetPolyQuadratic(radicand, x, out var qa, out var qb, out var qc) + && qa.Evaled is Entity.Number.Complex { IsZero: false } + => IntegrateRootOfQuadratic(qa, qb, qc, radicand, x), + + // ∫ k / sqrt(ax^2 + bx + c) dx -- the arcsine and logarithm forms. Without + // these, 1/sqrt(1 - x^2) had no antiderivative at all. + Entity.Divf(var numerator, + Entity.Powf(var radicand, Entity.Number.Rational(Entity.Number.Integer(1), Entity.Number.Integer(2)))) when + !numerator.ContainsNode(x) + && TreeAnalyzer.TryGetPolyQuadratic(radicand, x, out var ra, out var rb, out var rc) + => IntegrateOverRootOfQuadratic(numerator, ra, rb, rc, radicand, x), + _ => null }; + + /// + /// Reads an expression as sin(arg)^n * cos(arg)^m, where either power may be + /// zero, or reports that it is not of that shape. Both factors have to be functions + /// of the same argument. + /// + private static bool TryReadSineCosinePowers( + Entity expr, Entity.Variable x, out Entity arg, out int sinePower, out int cosinePower) + { + arg = 0; + sinePower = cosinePower = 0; + switch (expr) + { + case Entity.Sinf(var a): + arg = a; sinePower = 1; return a.ContainsNode(x); + case Entity.Cosf(var a): + arg = a; cosinePower = 1; return a.ContainsNode(x); + case Entity.Powf(Entity.Sinf(var a), Entity.Number.Integer power) + when power.EInteger.Sign > 0 && power.EInteger.CanFitInInt32(): + arg = a; sinePower = power.EInteger.ToInt32Checked(); return a.ContainsNode(x); + case Entity.Powf(Entity.Cosf(var a), Entity.Number.Integer power) + when power.EInteger.Sign > 0 && power.EInteger.CanFitInInt32(): + arg = a; cosinePower = power.EInteger.ToInt32Checked(); return a.ContainsNode(x); + case Entity.Mulf(var left, var right) + when TryReadSineCosinePowers(left, x, out var leftArg, out var leftSine, out var leftCosine) + && TryReadSineCosinePowers(right, x, out var rightArg, out var rightSine, out var rightCosine) + && leftArg == rightArg: + arg = leftArg; + sinePower = leftSine + rightSine; + cosinePower = leftCosine + rightCosine; + return true; + default: + return false; + } + } + + /// + /// The antiderivative of sin(u)^n * cos(u)^m where u is linear in x + /// with slope . + /// + /// + /// With an odd power there is a substitution: peel one factor off to be the + /// differential and write what is left in the other function, which turns the + /// integral into a polynomial. With both powers even there is no such factor to + /// peel, so the halved-angle identities go in instead and the result is integrated + /// again -- the total power halves each time, so this ends. + /// + private static Entity IntegrateSineCosinePowers( + Entity arg, int sinePower, int cosinePower, Entity rate, Entity.Variable x) + { + if (sinePower == 0 && cosinePower == 0) + return x; + + // cos^(2k+1) * sin^n: let s = sin(u), ds = cos(u) du + // = (1/rate) * sum_j C(k, j) (-1)^j s^(n + 2j + 1) / (n + 2j + 1) + // + // Tried before the sine substitution, so that where both powers are odd the + // answer comes back in sin rather than cos: sin(x)cos(x) integrates to + // sin(x)^2/2, which is the form everyone writes, rather than the equally + // correct -cos(x)^2/2 that differs from it by a constant. + if (cosinePower % 2 == 1) + { + var k = (cosinePower - 1) / 2; + Entity sum = 0; + for (var j = 0; j <= k; j++) + { + var power = sinePower + 2 * j + 1; + sum += Binomial(k, j) * (j % 2 == 0 ? 1 : -1) * MathS.Sin(arg).Pow(power) / power; + } + return sum / rate; + } + + // sin^(2k+1) * cos^m: let c = cos(u), dc = -sin(u) du + if (sinePower % 2 == 1) + { + var k = (sinePower - 1) / 2; + Entity sum = 0; + for (var j = 0; j <= k; j++) + { + var power = cosinePower + 2 * j + 1; + sum += Binomial(k, j) * (j % 2 == 0 ? 1 : -1) * MathS.Cos(arg).Pow(power) / power; + } + return -sum / rate; + } + + // Both even: sin^2 = (1 - cos(2u))/2 and cos^2 = (1 + cos(2u))/2, expanded into + // powers of cos(2u), each of which is integrated by the same rules. + var p = sinePower / 2; + var q = cosinePower / 2; + var doubled = 2 * arg; + Entity result = 0; + // (1 - t)^p (1 + t)^q with t = cos(2u), over 2^(p+q) + for (var i = 0; i <= p; i++) + for (var j = 0; j <= q; j++) + { + var coefficient = Binomial(p, i) * Binomial(q, j) * (i % 2 == 0 ? 1 : -1); + var term = IntegrateSineCosinePowers(doubled, 0, i + j, 2 * rate, x); + result += coefficient * term; + } + return result / Entity.Number.Integer.Create(EInteger.One.ShiftLeft(p + q)); + } + + private static Entity.Number.Integer Binomial(int n, int k) + { + var result = EInteger.One; + for (var i = 0; i < k; i++) + result = result * EInteger.FromInt32(n - i) / EInteger.FromInt32(i + 1); + return Entity.Number.Integer.Create(result); + } + + /// + /// Whether is an exponential in , and + /// at what rate: B^(px + q) grows as e^(rate * x) with + /// rate = p * ln(B). The constant factor B^q needs no separating out, + /// because the antiderivative is written in terms of the original expression. + /// + private static bool IsExponentialRate(Entity expr, Entity.Variable x, out Entity rate) + { + rate = 0; + if (expr is not Entity.Powf(var @base, var exponent) + || @base.ContainsNode(x) + || !TreeAnalyzer.TryGetPolyLinear(exponent, x, out var perX, out _)) + return false; + rate = perX * MathS.Ln(@base); + return true; + } + + /// + /// The antiderivative of sqrt(a x^2 + b x + c): + /// (2ax + b) sqrt(Q) / (4a) + ((4ac - b^2) / (8a)) times the integral of + /// 1/sqrt(Q) -- integration by parts once, leaving the reciprocal form that + /// already knows. + /// + /// + /// Only where the leading coefficient is a number other than zero. With a = 0 this + /// is the square root of something linear, which the ordinary power rule already + /// integrates, and dividing by a would not be allowed anyway. + /// + private static Entity IntegrateRootOfQuadratic( + Entity a, Entity b, Entity c, Entity radicand, Entity.Variable x) + => (2 * a * x + b) * MathS.Sqrt(radicand) / (4 * a) + + (4 * a * c - b * b) / (8 * a) * IntegrateOverRootOfQuadratic(1, a, b, c, radicand, x); + + /// + /// The antiderivative of k / sqrt(a x^2 + b x + c), which takes one of two + /// forms depending on the sign of the leading coefficient: + /// + /// a < 0, an arc of a circle: -k/sqrt(-a) * arcsin((2ax + b) / sqrt(b^2 - 4ac)) + /// a > 0, a hyperbolic arc: k/sqrt(a) * ln|2ax + b + 2 sqrt(a) sqrt(a x^2 + b x + c)| + /// + /// Returned as a piecewise on that sign, the way the rational quadratic below is, + /// since which one applies is not known until a and the coefficients are. + /// + private static Entity IntegrateOverRootOfQuadratic( + Entity numerator, Entity a, Entity b, Entity c, Entity radicand, Entity.Variable x) + { + var twoAxPlusB = 2 * a * x + b; + + // a < 0: the radicand is a downward parabola, positive between its roots + var arcsinCase = + -numerator * MathS.Arcsin(twoAxPlusB / MathS.Sqrt(b * b - 4 * a * c)) / MathS.Sqrt(-a); + + // a > 0 + var logarithmCase = + numerator * MathS.Ln(MathS.Abs(twoAxPlusB + 2 * MathS.Sqrt(a) * MathS.Sqrt(radicand))) / MathS.Sqrt(a); + + // a = 0: sqrt(bx + c), which integrates as an ordinary power + var linearCase = 2 * numerator * MathS.Sqrt(b * x + c) / b; + + return MathS.Piecewise([ + new Entity.Providedf(linearCase, a.EqualTo(0)), + new Entity.Providedf(arcsinCase, a < 0), + new Entity.Providedf(logarithmCase, a > 0) + ]); + } + private static Entity IntegrateRationalQuadratic(Entity numerator, Entity a, Entity b, Entity c, Entity.Variable x) { // The formula depends on whether it's linear (a = 0) or quadratic (a ≠ 0) diff --git a/Sources/Tests/UnitTests/Calculus/IntegrationTest.cs b/Sources/Tests/UnitTests/Calculus/IntegrationTest.cs index 99e5b7aa1..55865161c 100644 --- a/Sources/Tests/UnitTests/Calculus/IntegrationTest.cs +++ b/Sources/Tests/UnitTests/Calculus/IntegrationTest.cs @@ -268,13 +268,16 @@ public void TestLnAbsSquared() } } - [Theory(Skip = "TODO: integration by parts multiple times")] - [InlineData("ln(abs(x)) ^ 3", "C + x * (ln(abs(x)) ^ 3 - ln(abs(x)) ^ 2 - ln(abs(x)) ^ 2 - ln(abs(x)) ^ 2) + 6 * (x * (ln(abs(x)) - 1) + -x)")] // Triple integration by parts - [InlineData("e^x * sin(x)", "-1/2 * cos(x) * e ^ x + 1/2 * sin(x) * e ^ x + C")] // Classic integration by parts - [InlineData("e^x * cos(x)", "1/2 * cos(x) * e ^ x + 1/2 * sin(x) * e ^ x + C")] // Classic integration by parts - [InlineData("arctan(x)", "x * arctan(x) - 1/2 * ln(abs(x ^ 2 + 1)) + C")] // Integration by parts with 1 * arctan(x) - [InlineData("arcsin(x)", "x * arcsin(x) + sqrt(1 - x ^ 2) + C")] // Integration by parts with 1 * arcsin(x) - [InlineData("arccos(x)", "x * arccos(x) - sqrt(1 - x ^ 2) + C")] // Integration by parts with 1 * arccos(x) + // These five no longer need integration by parts at all. e^x*sin(x) and its + // cosine twin cycle under by parts, so they are solved as the closed form that + // cycle resolves to; the inverse trigonometric ones are by parts against 1, which + // has no product for the by-parts solver to split, so they are table entries. + [Theory] + [InlineData("e^x * sin(x)", "-1/2 * cos(x) * e ^ x + 1/2 * sin(x) * e ^ x + C")] + [InlineData("e^x * cos(x)", "1/2 * cos(x) * e ^ x + 1/2 * sin(x) * e ^ x + C")] + [InlineData("arctan(x)", "x * arctan(x) - 1/2 * ln(abs(x ^ 2 + 1)) + C")] + [InlineData("arcsin(x)", "x * arcsin(x) + sqrt(1 - x ^ 2) + C")] + [InlineData("arccos(x)", "x * arccos(x) - sqrt(1 - x ^ 2) + C")] public void TestIntegrationByPartsNonPolynomial(string initial, string expected) { var result = initial.Integrate("x").InnerSimplified; @@ -282,6 +285,16 @@ public void TestIntegrationByPartsNonPolynomial(string initial, string expected) Assert.Equal(MathS.Boolean.True, result.EqualTo(expectedResult).Simplify()); } + // Still open: this one wants by parts applied three times over. + [Theory(Skip = "TODO: integration by parts multiple times")] + [InlineData("ln(abs(x)) ^ 3", "C + x * (ln(abs(x)) ^ 3 - ln(abs(x)) ^ 2 - ln(abs(x)) ^ 2 - ln(abs(x)) ^ 2) + 6 * (x * (ln(abs(x)) - 1) + -x)")] + public void TestTripleIntegrationByParts(string initial, string expected) + { + var result = initial.Integrate("x").InnerSimplified; + var expectedResult = expected.ToEntity().InnerSimplified; + Assert.Equal(MathS.Boolean.True, result.EqualTo(expectedResult).Simplify()); + } + [Theory(Skip = "TODO: integration by parts multiple times")] [InlineData("sin(ln(abs(x)))", "x / 2 * (sin(ln(abs(x))) - cos(ln(abs(x)))) + C")] // Integration by parts twice [InlineData("cos(ln(abs(x)))", "x / 2 * (sin(ln(abs(x))) + cos(ln(abs(x)))) + C")] // Integration by parts twice diff --git a/Sources/Tests/UnitTests/Calculus/StandardIntegralsTest.cs b/Sources/Tests/UnitTests/Calculus/StandardIntegralsTest.cs new file mode 100644 index 000000000..97c250c68 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/StandardIntegralsTest.cs @@ -0,0 +1,134 @@ +// +// Copyright (c) 2019-2022 Angouri. +// AngouriMath is licensed under MIT. +// Details: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md. +// Website: https://am.angouri.org. +// + +using AngouriMath; +using AngouriMath.Extensions; +using Xunit; + +namespace AngouriMath.Tests.Calculus +{ + /// + /// Integrals added to the table rather than reached by the general solvers. Each is + /// checked by differentiating the answer back and comparing at points, since what + /// matters is that it is an antiderivative, not what form it is written in. + /// + public sealed class StandardIntegralsTest + { + /// + /// Chosen inside the integrand's own domain. 1/sqrt(x^2 - 1) is imaginary on + /// (-1, 1), so testing it there says nothing about whether the answer is right. + /// + private static void AssertIsAntiderivative(string integrand, params double[] points) + { + var f = integrand.ToEntity(); + var antiderivative = f.Integrate("x"); + Assert.DoesNotContain("integral(", antiderivative.Stringize()); + var derivative = antiderivative.Substitute("C", 0).Differentiate("x"); + foreach (var point in points) + { + var expected = f.Substitute("x", point).EvalNumerical().RealPart.EDecimal.ToDouble(); + var actual = derivative.Substitute("x", point).EvalNumerical().RealPart.EDecimal.ToDouble(); + Assert.Equal(expected, actual, 8); + } + } + + // k / sqrt(a x^2 + b x + c), which is an arcsine where a < 0 and a logarithm where + // a > 0. Nothing integrated 1/sqrt(1 - x^2) at all before. + [Theory] + [InlineData("1 / sqrt(1 - x ^ 2)", new[] { 0.31, 0.72, -0.4 })] + [InlineData("1 / sqrt(4 - x ^ 2)", new[] { 0.31, 1.7, -1.2 })] + [InlineData("2 / sqrt(9 - x ^ 2)", new[] { 0.5, 2.2, -2.5 })] + [InlineData("1 / sqrt(x ^ 2 + 1)", new[] { 0.31, 2.4, -1.9 })] + [InlineData("1 / sqrt(x ^ 2 - 1)", new[] { 1.4, 2.7, 5.1 })] + [InlineData("1 / sqrt(2 * x + 3)", new[] { 0.5, 2.2 })] + public void RootOfAQuadraticInTheDenominator(string integrand, double[] points) => + AssertIsAntiderivative(integrand, points); + + // An exponential times a sine returns the integral it started from after being + // integrated by parts twice, so the general solver cycles. Solving that equation + // once gives a closed form, which is what the table holds. + [Theory] + [InlineData("e ^ x * sin(x)", new[] { 0.31, 1.4, -0.8 })] + [InlineData("e ^ x * cos(x)", new[] { 0.31, 1.4, -0.8 })] + [InlineData("sin(x) * e ^ x", new[] { 0.31, 1.4, -0.8 })] + [InlineData("cos(x) * e ^ x", new[] { 0.31, 1.4, -0.8 })] + [InlineData("e ^ (2 * x) * sin(3 * x)", new[] { 0.31, 1.4, -0.8 })] + [InlineData("e ^ (-x) * sin(x)", new[] { 0.31, 1.4, -0.8 })] + [InlineData("2 ^ x * cos(x)", new[] { 0.31, 1.4, -0.8 })] + public void ExponentialTimesAWave(string integrand, double[] points) => + AssertIsAntiderivative(integrand, points); + + // The inverse trigonometric functions are integration by parts against 1, and the + // by-parts solver looks for a product to split, so it never sees them. + [Theory] + [InlineData("arcsin(x)", new[] { 0.21, 0.44, -0.33 })] + [InlineData("arccos(x)", new[] { 0.21, 0.44, -0.33 })] + [InlineData("arctan(x)", new[] { 0.21, 0.44, -0.33 })] + [InlineData("arccotan(x)", new[] { 0.21, 0.44, -0.33 })] + [InlineData("arcsin(2 * x)", new[] { 0.21, 0.44, -0.33 })] + [InlineData("arctan(3 * x + 1)", new[] { 0.21, 0.44, -0.33 })] + public void InverseTrigonometricFunctions(string integrand, double[] points) => + AssertIsAntiderivative(integrand, points); + + // The shapes these sit next to in the table have to keep working. + [Theory] + [InlineData("1 / (x ^ 2 + 1)", new[] { 0.31, 1.4, -0.8 })] + [InlineData("e ^ x", new[] { 0.31, 1.4, -0.8 })] + [InlineData("sin(x)", new[] { 0.31, 1.4, -0.8 })] + [InlineData("x * e ^ x", new[] { 0.31, 1.4, -0.8 })] + [InlineData("sin(x) ^ 2", new[] { 0.31, 1.4, -0.8 })] + [InlineData("1 / x", new[] { 0.31, 1.4 })] + public void NeighbouringFormsAreUnaffected(string integrand, double[] points) => + AssertIsAntiderivative(integrand, points); + + // sin(u)^n * cos(u)^m for whole n and m. An odd power gives a factor to peel off + // as the differential, which turns the integral into a polynomial; with both even + // the halved-angle identities go in and the result is integrated again. + [Theory] + [InlineData("sin(x) ^ 3", new[] { 0.31, 1.4, -0.8 })] + [InlineData("cos(x) ^ 3", new[] { 0.31, 1.4, -0.8 })] + [InlineData("sin(x) ^ 4", new[] { 0.31, 1.4, -0.8 })] + [InlineData("cos(x) ^ 5", new[] { 0.31, 1.4, -0.8 })] + [InlineData("sin(x) ^ 6", new[] { 0.31, 1.4, -0.8 })] + [InlineData("sin(x) ^ 2 * cos(x) ^ 2", new[] { 0.31, 1.4, -0.8 })] + [InlineData("sin(x) ^ 3 * cos(x) ^ 2", new[] { 0.31, 1.4, -0.8 })] + [InlineData("sin(x) ^ 2 * cos(x) ^ 3", new[] { 0.31, 1.4, -0.8 })] + [InlineData("sin(2 * x) ^ 2", new[] { 0.31, 1.4, -0.8 })] + [InlineData("sin(3 * x + 1) ^ 3", new[] { 0.31, 1.4, -0.8 })] + public void PowersOfSineAndCosine(string integrand, double[] points) => + AssertIsAntiderivative(integrand, points); + + // Where both powers are odd either substitution works, and they differ by a + // constant. The sine one is the form everyone writes, and two existing tests + // assert it, so it is the one tried first. + [Fact] + public void BothPowersOddGivesTheSineForm() => + Assert.Equal(MathS.Boolean.True, + "sin(x) * cos(x)".Integrate("x").InnerSimplified + .EqualTo("sin(x) ^ 2 / 2 + C".ToEntity().InnerSimplified).Simplify()); + + // sqrt(a x^2 + b x + c), which is one integration by parts away from the + // reciprocal form above and is written in terms of it. + [Theory] + [InlineData("sqrt(1 - x ^ 2)", new[] { 0.31, 0.72, -0.4 })] + [InlineData("sqrt(4 - x ^ 2)", new[] { 0.31, 1.7, -1.2 })] + [InlineData("sqrt(x ^ 2 + 1)", new[] { 0.31, 2.4, -1.9 })] + [InlineData("sqrt(x ^ 2 - 1)", new[] { 1.4, 2.7, 5.1 })] + [InlineData("sqrt(2 * x ^ 2 + 3 * x + 5)", new[] { 0.31, 1.4 })] + public void RootOfAQuadratic(string integrand, double[] points) => + AssertIsAntiderivative(integrand, points); + + // With no quadratic term this is the square root of something linear, which the + // ordinary power rule already integrates -- and dividing by the leading + // coefficient would not be allowed. It must not be taken over. + [Theory] + [InlineData("sqrt(x)", new[] { 0.31, 1.4 })] + [InlineData("sqrt(2 * x + 3)", new[] { 0.5, 2.2 })] + public void RootOfSomethingLinearIsLeftToThePowerRule(string integrand, double[] points) => + AssertIsAntiderivative(integrand, points); + } +}