diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 733433eba..33458eb6e 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -1725,6 +1725,26 @@ the exponential integral, as `x e^(-1/x^2)` does, is still left unintegrated | `"f^(a + b/(c + d*x)^2)/(c + d*x)^4".Integrate("x")` | `integral(f ^ (a + b / (c + d * x) ^ 2) / (c + d * x) ^ 4, x)` | the antiderivative | | `"sinh(a + b/x^2)".Integrate("x")` | `integral((e ^ (a + b / x ^ 2) - e ^ (-(a + b / x ^ 2))) / 2, x)` | the antiderivative | +### Exponentials of quadratics beside sines, cosines and hyperbolic functions of quadratics are integrated + +`e^(x^2) sin(x)` was left unintegrated. A product of exponentials of polynomials of degree at most +two, whole powers of sines, cosines and hyperbolic functions of such polynomials, and a polynomial +is now integrated wherever some exponent is a quadratic. Each sine or cosine is written as +exponentials of `i` times its argument, and every term of the product multiplied out is the +Gaussian's. So the answer for a real integrand with a sine or cosine in it is written with the +imaginary unit, as Rubi's is: `e^(x^2) sin(x)` integrates to +`-i/2 e^(1/4) sqrt(pi)/2 erfi(x + i/2) + i/2 e^(1/4) sqrt(pi)/2 erfi(x - i/2)`. A sine or cosine of a +quadratic alone, or beside a polynomial only, is left to the rules that answered it before +([#1501](https://github.com/asc-community/AngouriMath/issues/1501)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"e^(x^2)*sin(x)".Integrate("x")` | `integral(e ^ x ^ 2 * sin(x), x)` | the antiderivative, with `erfi` | +| `"f^(a + b*x + c*x^2)*cos(d + p*x + q*x^2)".Integrate("x")` | `integral(f ^ (a + b * x + c * x ^ 2) * cos(d + p * x + q * x ^ 2), x)` | the antiderivative, with `erf` | +| `"f^(a + c*x^2)*cosh(d + p*x)^2".Integrate("x")` | `integral(f ^ (a + c * x ^ 2) * ((e ^ (d + p * x) + e ^ (-(d + p * x))) / 2) ^ 2, x)` | the antiderivative | +| `"x^2*sinh(a + b*x + c*x^2)^2".Integrate("x")` | `integral(x ^ 2 * ((e ^ (a + b * x + c * x ^ 2) - e ^ (-(a + b * x + c * x ^ 2))) / 2) ^ 2, x)` | the antiderivative | +| `"2^(x^2)*3^(x + 1)".Integrate("x")` | `integral(2 ^ x ^ 2 * 3 ^ (x + 1), x)` | the antiderivative, with `erfi` | + ### `binomial(n, k)` is a function **Addition, not silent.** The binomial coefficient is a node, `Entity.Binomialf`, spelled diff --git a/Sources/.editorconfig b/Sources/.editorconfig index 3614f5364..47bf3f25e 100644 --- a/Sources/.editorconfig +++ b/Sources/.editorconfig @@ -639,3 +639,6 @@ file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed [Tests/UnitTests/Calculus/HalfOddPowerOfALogarithmTest.cs] file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n + +[Tests/UnitTests/Calculus/ExponentialsOfQuadraticsTest.cs] +file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index e995820b9..3ef20c341 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -3225,6 +3225,374 @@ private static bool TryReadAMonomial(Entity term, Entity.Variable x, out int deg } } + /// + /// A product of exponentials of polynomials of degree at most two, whole powers of sines + /// and cosines of such polynomials, whole powers of sums of those exponentials -- which is + /// how a hyperbolic function is written -- and a polynomial, where some exponent is left a + /// quadratic. A sine or cosine is written as exponentials, sin(q) = (e^(i q) - e^(-i q))/(2i), + /// and the product multiplied out is a sum of terms k P(x) e^(Q(x)), with Q of + /// degree at most two: the Gaussian's moments where Q is a quadratic, and elementary + /// where it is not. Rubi's 4.7.6, 6.1.5 and 6.2.5, f^(a + b x + c x^2) sin(d + e x + f x^2)^n + /// and f^(a + b x + c x^2) cosh(d + e x + f x^2)^n, 6.1.4's x^2 sinh(a + b x + c x^2)^2, + /// and 2.3's f^(a + b x + c x^2) g^(d + e x + f x^2). + /// https://github.com/asc-community/AngouriMath/issues/1501 + /// + /// + /// + /// A sine or cosine is taken only beside an exponential. Alone, or beside a polynomial, + /// the sine of a quadratic is a Fresnel integral, and x sin(x^2) is -cos(x^2)/2, + /// which the substitution answers with no imaginary unit in it. A sum of exponentials is + /// taken beside another exponential, or in a power: x^2 sinh(q) alone is a sum the + /// integrator splits, and each of its terms is the Gaussian's. + /// + /// + /// A term's exponent is kept as whole multiples of the exponents the integrand writes, not + /// as their sum: e^(i q) times e^(-i q) is then the multiple 0, decidably, + /// where the sum i q + (-i) q is not read as zero, and was taken for a Gaussian + /// whose square's coefficient is zero. + /// + /// + internal static Entity? SolveAProductOfExponentialsOfQuadratics(Entity expr, Entity.Variable x) + { + if (expr is not (Mulf or Divf or Powf)) + return null; + // By the written shape first, which reads and allocates nothing: most products asked + // here have one exponential in them, or none, and no sine. + int exponentials = 0, trigonometrics = 0; + var aPowerOfASum = false; + CountFactorsByShape(expr, x, ref exponentials, ref trigonometrics, ref aPowerOfASum); + if (trigonometrics > 0 ? exponentials == 0 : exponentials < 2 && !aPowerOfASum) + return null; + (exponentials, trigonometrics, aPowerOfASum) = (0, 0, false); + Entity constant = Number.Integer.One; + Entity? polynomial = null; + var exponents = new List(); + var terms = new List<(Entity Coefficient, int[] Multiples)> { (Number.Integer.One, new int[MostExponents]) }; + foreach (var (factor, underneath) in FactorsOfTheIntegrand(expr)) + { + if (!factor.ContainsNode(x)) + { + constant = underneath ? constant / factor : constant * factor; + continue; + } + if (AsExponentialsOfQuadratics(factor, x, exponents) is { } read) + { + // Below the bar, a single exponential only, whose reciprocal is one too. + if (underneath) + { + if (read.Count != 1) + return null; + var (coefficient, multiples) = read[0]; + for (var i = 0; i < multiples.Length; i++) + multiples[i] = -multiples[i]; + read[0] = (1 / coefficient, multiples); + } + if (factor.Nodes.Any(node => node is Sinf or Cosf)) + trigonometrics++; + else + exponentials++; + aPowerOfASum |= read.Count > 1 && factor is Powf(_, Number.Integer { EInteger: var power }) && power.CompareTo(EInteger.One) > 0; + if (Multiplied(terms, read) is not { } product) + return null; + terms = product; + continue; + } + if (!underneath && TreeAnalyzer.TryGetPolynomial(factor, x, out var inX) + && inX.Keys.All(degree => degree.Sign >= 0 && degree.CompareTo(EInteger.FromInt32(MostPowerOfAPolynomialBesideExponentials)) <= 0)) + { + polynomial = polynomial is null ? factor : polynomial * factor; + continue; + } + return null; + } + if (trigonometrics > 0 ? exponentials == 0 : exponentials < 2 && !aPowerOfASum) + return null; + var quadratics = new List<(Entity Coefficient, Entity A, Entity B, Entity C)>(terms.Count); + foreach (var (coefficient, multiples) in terms) + { + var (a, b, c) = TheExponent(multiples, exponents); + quadratics.Add((coefficient, a, b, c)); + } + // Some exponent left a quadratic: with none, the product is the elementary + // exponential's, and the rules for that answer it without an imaginary unit. + if (!quadratics.Any(term => !TreeAnalyzer.IsZero(term.A))) + return null; + Dictionary? monomials = null; + if (polynomial is not null && !TreeAnalyzer.TryGetPolynomial(polynomial, x, out monomials)) + return null; + Entity total = Number.Integer.Zero; + foreach (var (coefficient, a, b, c) in quadratics) + { + Entity? integral; + if (!TreeAnalyzer.IsZero(a)) + integral = IntegralPatterns.APolynomialTimesTheGaussian(polynomial, MathS.e, a, b, c, x); + else if (!TreeAnalyzer.IsZero(b)) + integral = MathS.Pow(MathS.e, b * x + c) * APolynomialAgainstAnExponential(monomials, b, x); + else + integral = MathS.Pow(MathS.e, c) * APolynomialIntegrated(monomials, x); + if (integral is null) + return null; + total += coefficient * integral; + } + return (constant * total).InnerSimplified; + } + + /// + /// An exponent the integrand writes, A x^2 + B x + C once read, of the exponential + /// Base^Written, or of e^(i Written) for a sine or cosine of it. + /// + private sealed record WrittenExponent(Entity Base, Entity Written, bool Imaginary, Entity A, Entity B, Entity C); + + private const int MostExponents = 6; + private const int MostPowerOfAPolynomialBesideExponentials = 32; + private const int MostPowerOfASumOfExponentials = 8; + private const int MostExponentialTerms = 64; + + /// + /// Counts the factors of by their written shape alone: a sine or + /// cosine, or a power of one; an exponential of ; and a sum, which may + /// be one of exponentials, or a whole power of a sum. An over-count, which the reading + /// that follows corrects: (1 + x) is counted as a sum. + /// + private static void CountFactorsByShape(Entity expr, Entity.Variable x, ref int exponentials, ref int trigonometrics, ref bool aPowerOfASum) + { + switch (expr) + { + case Mulf(var left, var right): + CountFactorsByShape(left, x, ref exponentials, ref trigonometrics, ref aPowerOfASum); + CountFactorsByShape(right, x, ref exponentials, ref trigonometrics, ref aPowerOfASum); + break; + case Divf(var numerator, var denominator): + CountFactorsByShape(numerator, x, ref exponentials, ref trigonometrics, ref aPowerOfASum); + CountFactorsByShape(denominator, x, ref exponentials, ref trigonometrics, ref aPowerOfASum); + break; + case Sinf or Cosf or Powf(Sinf or Cosf, Number.Integer): + trigonometrics++; + break; + case Powf(var @base, var exponent) when !@base.ContainsNode(x) && exponent.ContainsNode(x): + case Sumf or Minusf: + exponentials++; + break; + case Powf(Sumf or Minusf or Mulf or Divf, Number.Integer { EInteger: var power }) when power.CompareTo(EInteger.One) > 0: + exponentials++; + aPowerOfASum = true; + break; + } + } + + /// + /// as a sum of terms k e^(m_1 Q_1 + m_2 Q_2 + ...), where the + /// Q_j are and the m_j whole, or + /// : an exponential of a polynomial of degree at most two, a sine or + /// cosine of one, and sums, products, quotients by a constant and small whole powers of those. + /// + private static List<(Entity Coefficient, int[] Multiples)>? AsExponentialsOfQuadratics(Entity expr, Entity.Variable x, List exponents) + { + switch (expr) + { + case var constant when !constant.ContainsNode(x): + return new() { (constant, new int[MostExponents]) }; + case Powf(var @base, var exponent) when !@base.ContainsNode(x): + { + // e^(-q), as sinh(q) writes it, is -1 times the exponent of its e^q. + var (written, multiple) = exponent switch + { + Mulf(Number.Integer k, var inner) when k.EInteger.Abs().CompareTo(EInteger.FromInt32(MostPowerOfASumOfExponentials)) <= 0 => (inner, k.EInteger.ToInt32Checked()), + Mulf(var inner, Number.Integer k) when k.EInteger.Abs().CompareTo(EInteger.FromInt32(MostPowerOfASumOfExponentials)) <= 0 => (inner, k.EInteger.ToInt32Checked()), + _ => (exponent, 1) + }; + if (TheWrittenExponent(@base, written, imaginary: false, x, exponents) is not { } index) + return null; + var multiples = new int[MostExponents]; + multiples[index] = multiple; + return new() { (Number.Integer.One, multiples) }; + } + // sin(q) = (e^(i q) - e^(-i q))/(2i) and cos(q) = (e^(i q) + e^(-i q))/2. + case Sinf(var argument): + return AnImaginaryPair(-MathS.i / 2, MathS.i / 2, argument, x, exponents); + case Cosf(var argument): + return AnImaginaryPair(Number.Rational.Create(1, 2), Number.Rational.Create(1, 2), argument, x, exponents); + case Sumf(var left, var right): + { + if (AsExponentialsOfQuadratics(left, x, exponents) is not { } l || AsExponentialsOfQuadratics(right, x, exponents) is not { } r) + return null; + l.AddRange(r); + return l.Count > MostExponentialTerms ? null : l; + } + case Minusf(var left, var right): + { + if (AsExponentialsOfQuadratics(left, x, exponents) is not { } l || AsExponentialsOfQuadratics(right, x, exponents) is not { } r) + return null; + foreach (var (coefficient, multiples) in r) + l.Add((-coefficient, multiples)); + return l.Count > MostExponentialTerms ? null : l; + } + case Mulf(var left, var right): + { + if (AsExponentialsOfQuadratics(left, x, exponents) is not { } l || AsExponentialsOfQuadratics(right, x, exponents) is not { } r) + return null; + return Multiplied(l, r); + } + case Divf(var numerator, var denominator) when !denominator.ContainsNode(x): + { + if (AsExponentialsOfQuadratics(numerator, x, exponents) is not { } n) + return null; + for (var i = 0; i < n.Count; i++) + n[i] = (n[i].Coefficient / denominator, n[i].Multiples); + return n; + } + case Powf(var inner, Number.Integer power) when power.EInteger.Sign > 0 + && power.EInteger.CompareTo(EInteger.FromInt32(MostPowerOfASumOfExponentials)) <= 0: + { + if (AsExponentialsOfQuadratics(inner, x, exponents) is not { } once) + return null; + var result = once; + for (var k = power.EInteger.ToInt32Checked(); k > 1 && result is not null; k--) + result = Multiplied(result, once); + return result; + } + default: + return null; + } + } + + /// + /// The index of Base^Written, or of e^(i Written), among + /// , read and added where it is not there yet, or + /// where it is not a polynomial of degree at most two. + /// + private static int? TheWrittenExponent(Entity @base, Entity written, bool imaginary, Entity.Variable x, List exponents) + { + for (var i = 0; i < exponents.Count; i++) + if (exponents[i].Imaginary == imaginary && exponents[i].Base == @base && exponents[i].Written == written) + return i; + if (exponents.Count == MostExponents || !written.ContainsNode(x) + || !TreeAnalyzer.TryGetPolyQuadratic(written, x, out var a, out var b, out var c)) + return null; + // A zero coefficient is left the integer zero: 0 ln(f), simplified, is 0 provided f is + // not 0, which is not read as zero, and a Gaussian was taken of 0 x^2. + Entity scale = imaginary ? MathS.i : @base == MathS.e ? Number.Integer.One : MathS.Ln(@base); + exponents.Add(new(@base, written, imaginary, Scaled(a, scale), Scaled(b, scale), Scaled(c, scale))); + return exponents.Count - 1; + + static Entity Scaled(Entity coefficient, Entity scale) + => TreeAnalyzer.IsZero(coefficient) ? Number.Integer.Zero + : scale == Number.Integer.One ? coefficient : (scale * coefficient).InnerSimplified; + } + + /// + /// k_+ e^(i q) + k_- e^(-i q) for the q, or + /// where it is not a polynomial of degree at most two. + /// + private static List<(Entity Coefficient, int[] Multiples)>? AnImaginaryPair(Entity plus, Entity minus, Entity argument, Entity.Variable x, List exponents) + { + if (TheWrittenExponent(MathS.e, argument, imaginary: true, x, exponents) is not { } index) + return null; + var (up, down) = (new int[MostExponents], new int[MostExponents]); + (up[index], down[index]) = (1, -1); + return new() { (plus, up), (minus, down) }; + } + + /// + /// The product of two sums of exponentials, with the terms of one exponent gathered, and + /// past . + /// + private static List<(Entity Coefficient, int[] Multiples)>? Multiplied(List<(Entity Coefficient, int[] Multiples)> left, List<(Entity Coefficient, int[] Multiples)> right) + { + if (left.Count * right.Count > MostExponentialTerms) + return null; + var gathered = new List<(Entity Coefficient, int[] Multiples)>(); + foreach (var l in left) + foreach (var r in right) + { + var multiples = new int[MostExponents]; + for (var i = 0; i < multiples.Length; i++) + multiples[i] = l.Multiples[i] + r.Multiples[i]; + var coefficient = l.Coefficient * r.Coefficient; + var at = gathered.FindIndex(term => term.Multiples.SequenceEqual(multiples)); + if (at < 0) + gathered.Add((coefficient, multiples)); + else + gathered[at] = (gathered[at].Coefficient + coefficient, multiples); + } + var product = new List<(Entity Coefficient, int[] Multiples)>(gathered.Count); + foreach (var (coefficient, multiples) in gathered) + { + var simplified = coefficient.InnerSimplified; + if (!TreeAnalyzer.IsZero(simplified)) + product.Add((simplified, multiples)); + } + return product; + } + + /// + /// (A, B, C) of the exponent sum_j m_j Q_j, with the zero multiples left out, + /// so that a term with none is decidably e^0. + /// + private static (Entity A, Entity B, Entity C) TheExponent(int[] multiples, List exponents) + { + Entity a = Number.Integer.Zero, b = Number.Integer.Zero, c = Number.Integer.Zero; + for (var i = 0; i < exponents.Count; i++) + { + if (multiples[i] == 0) + continue; + var (m, exponent) = (multiples[i], exponents[i]); + a = Plus(a, m, exponent.A); + b = Plus(b, m, exponent.B); + c = Plus(c, m, exponent.C); + } + return (a.InnerSimplified, b.InnerSimplified, c.InnerSimplified); + + static Entity Plus(Entity sum, int multiple, Entity coefficient) + { + if (coefficient == Number.Integer.Zero) + return sum; + var term = multiple == 1 ? coefficient : multiple * coefficient; + return sum == Number.Integer.Zero ? term : sum + term; + } + } + + /// + /// int P(x) e^(b x) dx over e^(b x): sum_j (-1)^j P^(j)(x)/b^(j + 1), by + /// parts, for P read into , or 1 where they are absent. + /// + private static Entity APolynomialAgainstAnExponential(Dictionary? monomials, Entity rate, Entity.Variable x) + { + if (monomials is null) + return 1 / rate; + Entity sum = Number.Integer.Zero; + foreach (var monomial in monomials) + { + var power = monomial.Key.ToInt32Checked(); + // The j-th derivative of x^k is k!/(k - j)! x^(k - j). + Entity falling = Number.Integer.One; + for (var j = 0; j <= power; j++) + { + var term = monomial.Value * falling * MathS.Pow(x, power - j) / MathS.Pow(rate, j + 1); + sum = j % 2 == 0 ? sum + term : sum - term; + falling *= power - j; + } + } + return sum; + } + + /// + /// int P(x) dx for P read into , or x where + /// they are absent. + /// + private static Entity APolynomialIntegrated(Dictionary? monomials, Entity.Variable x) + { + if (monomials is null) + return x; + Entity sum = Number.Integer.Zero; + foreach (var monomial in monomials) + { + var power = monomial.Key.ToInt32Checked() + 1; + sum += monomial.Value * MathS.Pow(x, power) / power; + } + return sum; + } + /// /// A polynomial times an exponential times a sine or a cosine — /// P(x) e^(a x) cos(b x) and its kin — integrated by the repeated by-parts that diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IntegralPatterns.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IntegralPatterns.cs index 19acdee79..4c9e9d394 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IntegralPatterns.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IntegralPatterns.cs @@ -85,6 +85,15 @@ private static Entity Gaussian(Entity @base, Entity a, Entity b, Entity c, Entit return MathS.Pow(@base, c - b * b / (4 * a)) * UnitGaussian(a * MathS.Ln(@base), x + b / (2 * a)); } + /// + /// int P(x) F^(a x^2 + b x + c) dx for a non-zero a and a polynomial P, + /// which may be absent: the Gaussian, or the sum of its moments at the square's centre. + /// + internal static Entity? APolynomialTimesTheGaussian(Entity? polynomial, Entity @base, Entity a, Entity b, Entity c, Entity.Variable x) + => polynomial is null + ? Gaussian(@base, a, b, c, x) + : ShiftedMoments(Entity.Number.Integer.One, 0, polynomial, @base, a, b, c, x); + /// /// int e^(A u^2) du for a non-zero A: sqrt(pi)/(2 sqrt(-A)) erf(sqrt(-A) u), /// or with erfi and the real root where A is decidably positive. diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs index 5e5bef213..a41f92922 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs @@ -582,6 +582,9 @@ private static Entity Normalized(Entity expr, Entity.Variable x) => if ((answer = IndefiniteIntegralSolver.SolveByWritingAnImaginaryTangentAsAnExponential(expr, x, integrateByParts)) is { }) return answer; // And `A cos(z) + i A sin(z)`, which is `A e^(i z)`, where no rotation is real. if ((answer = IndefiniteIntegralSolver.SolveByWritingAnImaginarySumOfACosineAndASineAsAnExponential(expr, x, integrateByParts)) is { }) return answer; + // Exponentials of quadratics, with sines, cosines and sums of them multiplied out: + // each term is the Gaussian's, where the search would take the product by parts. + if ((answer = IndefiniteIntegralSolver.SolveAProductOfExponentialsOfQuadratics(expr, x)) is { }) return answer; // A fractional power of a perfect square is the power of the modulus, sgn(P) P^(2r). if ((answer = IndefiniteIntegralSolver.SolveByTakingARootOfAPerfectSquare(expr, x, integrateByParts)) is { }) return answer; // x^(n - 1) g(x^n) with a symbolic n is g(u)/n under u = x^n. diff --git a/Sources/Tests/UnitTests/Calculus/ExponentialsOfQuadraticsTest.cs b/Sources/Tests/UnitTests/Calculus/ExponentialsOfQuadraticsTest.cs new file mode 100644 index 000000000..4f37a3e28 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/ExponentialsOfQuadraticsTest.cs @@ -0,0 +1,111 @@ +// +// Copyright (c) 2019-2026 Angouri. +// AngouriMath is licensed under MIT. +// Details: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md. +// Website: https://am.angouri.org. +// + +using System; +using System.Linq; +using AngouriMath.Extensions; +using Xunit; + +namespace AngouriMath.Tests.Calculus +{ + /// + /// A product of exponentials of quadratics, with sines, cosines and hyperbolic functions of + /// quadratics beside them, multiplied out into exponentials: each term is the Gaussian's. + /// https://github.com/asc-community/AngouriMath/issues/1501 + /// + [Trait("Area", "Calculus")] + public sealed class ExponentialsOfQuadraticsTest + { + private static readonly double[] Points = { -1.3, -0.6, 0.35, 0.9, 1.45 }; + + /// + /// Integrates, pins the parameters, and compares the derivative of the answer with the + /// integrand at . The parameters are pinned after integrating, so the + /// rule is asked the symbolic question. + /// + private static Entity DifferentiatesBack(string integrand, params (string Name, string Value)[] pins) + { + var integral = integrand.ToEntity().Integrate("x"); + Assert.DoesNotContain("integral(", integral.Stringize()); + Entity Pinned(Entity e) => pins.Aggregate(e, (current, pin) => current.Substitute(pin.Name, pin.Value.ToEntity())); + var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x"); + var original = Pinned(integrand.ToEntity()); + foreach (var at in Points) + { + var got = derivative.Substitute("x", at).EvalNumerical(); + var want = original.Substitute("x", at).EvalNumerical(); + var difference = Math.Abs((double)(got - want).RealPart) + Math.Abs((double)(got - want).ImaginaryPart); + var scale = Math.Max(1.0, Math.Abs((double)want.RealPart)); + Assert.True(difference / scale < 1e-9, + $"d/dx of the antiderivative of {integrand} is {got} at x = {at}, where the integrand is {want}"); + } + return integral; + } + + private static readonly (string, string)[] Pins = + { ("a", "1/3"), ("b", "1/2"), ("c", "-1"), ("d", "2/3"), ("f", "2"), ("g", "3"), ("p", "3/4"), ("q", "1/5") }; + + /// + /// Rubi's 4.7.6: a sine or cosine written as exponentials of the imaginary unit times its + /// argument, so the answer holds an error function of a complex argument. + /// + [Theory] + [InlineData("e^(x^2)*sin(x)")] + [InlineData("e^(x^2)*cos(2*x + 1)")] + [InlineData("e^x*sin(x^2)")] + [InlineData("e^x*cos(1 + x + x^2)")] + [InlineData("2^(x^2 + x)*cos(x^2)")] + [InlineData("e^(-x^2)*sin(x)^2")] + [InlineData("e^x*cos(x^2 + x)^3")] + [InlineData("x*e^(x^2)*sin(x)")] + [InlineData("f^(a + b*x + c*x^2)*sin(d + p*x + q*x^2)")] + [InlineData("f^(a + b*x)*cos(d + q*x^2)^2")] + [InlineData("f^(a + c*x^2)*sin(d + p*x)^3")] + public void ASineOrCosineBesideAnExponential(string integrand) + => DifferentiatesBack(integrand, Pins); + + /// Rubi's 6.1.5 and 6.2.5, where every exponent is real. + [Theory] + [InlineData("e^x*sinh(x^2)")] + [InlineData("2^(x^2)*cosh(x + 1)")] + [InlineData("e^(x^2)*sinh(x)^2")] + [InlineData("f^(a + b*x + c*x^2)*sinh(d + p*x + q*x^2)")] + [InlineData("f^(a + c*x^2)*cosh(d + p*x)^2")] + [InlineData("f^(a + b*x)*sinh(d + q*x^2)^3")] + public void AHyperbolicBesideAnExponential(string integrand) + => DifferentiatesBack(integrand, Pins); + + /// + /// Rubi's 6.1.4 and 6.2.4: a power of a hyperbolic function of a quadratic beside a + /// polynomial, which was taken by parts past the budget. + /// + [Theory] + [InlineData("x^2*sinh(1 + x + x^2)^2")] + [InlineData("(1 + 2*x)^2*cosh(x - x^2)^2")] + [InlineData("x^2*cosh(a + b*x + c*x^2)^2")] + public void APowerOfAHyperbolicBesideAPolynomial(string integrand) + => DifferentiatesBack(integrand, Pins); + + /// Rubi's 2.3: two exponentials of quadratics, of different bases. + [Theory] + [InlineData("2^(x^2)*3^(x + 1)")] + [InlineData("f^(a + b*x + c*x^2)*g^(d + p*x + q*x^2)")] + public void TwoExponentials(string integrand) + => DifferentiatesBack(integrand, Pins); + + /// + /// A sine or cosine beside a polynomial only, or of a linear argument beside an + /// exponential of one, is left to the rules that answer it without an imaginary unit. + /// + [Theory] + [InlineData("x*sin(x^2)")] + [InlineData("e^x*sin(2*x)")] + [InlineData("e^(2*x)*cos(x)^2")] + public void ARealAnswerStaysReal(string integrand) + => Assert.DoesNotContain(DifferentiatesBack(integrand).Nodes, node => node == MathS.i); + } +}