diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IntegralPatterns.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IntegralPatterns.cs index fd5ce97ec..41fe74903 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IntegralPatterns.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IntegralPatterns.cs @@ -11,8 +11,53 @@ namespace AngouriMath.Functions.Algebra { internal static class IntegralPatterns { + /// + /// The argument through which each rule below reads its integrand's dependence on + /// , or for a shape none of them matches. + /// + /// + /// The conditions on the two two-argument nodes are the ones their own rules carry: + /// a base that mentions x is a different integral, not this one. + /// + private static Entity? RateBearingArgument(Entity expr, Entity.Variable x) => expr switch + { + Entity.Sinf(var arg) => arg, + Entity.Cosf(var arg) => arg, + Entity.Secantf(var arg) => arg, + Entity.Cosecantf(var arg) => arg, + Entity.Tanf(var arg) => arg, + Entity.Cotanf(var arg) => arg, + Entity.Absf(var arg) => arg, + Entity.Signumf(var arg) => arg, + Entity.Arcsinf(var arg) => arg, + Entity.Arccosf(var arg) => arg, + Entity.Arctanf(var arg) => arg, + Entity.Arccotanf(var arg) => arg, + Entity.Logf(var @base, var arg) when !@base.ContainsNode(x) => arg, + Entity.Powf(var @base, var power) when !@base.ContainsNode(x) => power, + _ => null + }; + internal static Entity? TryStandardIntegrals(Entity expr, Entity.Variable x) => expr switch { + // Every rule below divides by the linear rate of its integrand's argument, and + // an argument that mentions x without depending on it has a rate of zero. That + // put a literal division by zero into the answer, so `e ^ (x + -x)` integrated + // to `e ^ (x + -x) / (0 * ln(e))`, which evaluates to NaN. Each of those + // integrands depends on x only through that argument, so a zero rate means the + // integrand is a constant and integrates to itself times x. + // + // The rate has to be decidably zero: a symbolic one, as in sin(a * x), is not, + // and answering sin(a * x) * x there would be wrong for every non-zero a. + // + // Reachable from ordinary input once two powers of one base are gathered: + // e^x * e^(-x) becomes e^(x + -x), whose rate is zero while `x + -x` is still + // written out. https://github.com/asc-community/AngouriMath/issues/785 + _ when RateBearingArgument(expr, x) is { } constant + && TreeAnalyzer.TryGetPolyLinear(constant, x, out var rate, out _) + && rate.Evaled is Entity.Number.Complex { IsZero: true } + => expr * x, + Entity.Sinf(var arg) when TreeAnalyzer.TryGetPolyLinear(arg, x, out var a, out _) => -MathS.Cos(arg) / a, diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs index e3cd35982..0799f9ff2 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs @@ -46,9 +46,33 @@ namespace AngouriMath.Functions.Algebra { internal static partial class Integration { + /// + /// Brings the two powers of one base in a product together, so that an integrand + /// which is a power in disguise is recognised as one. + /// + /// + /// + /// This has to happen on every recursive call rather than once on the way in. The + /// shape is produced by the integrator itself -- distributing a product over a sum + /// turns sin(x)^4 * (5 - 6*sin(x)^2) into a term sin(x)^4 * (-6) * sin(x)^2 + /// -- so normalising only the caller's input would miss every case the integrator + /// generates for itself. https://github.com/asc-community/AngouriMath/issues/781 + /// + /// + /// Deliberately not followed by InnerSimplified: that rewrites x^(-2) + /// back into 1/x^2, and + /// rewrites a 1/x^n it is handed into Pow(x, -n), so the pair recurse + /// into each other until the stack runs out. The gathering folds the exponents it + /// builds and leaves the rest of the tree alone. + /// + /// + private static Entity Normalized(Entity expr) => + expr.Replace(Patterns.GatherPowersOfOneBase); + /// Does not add the constant of integration because this is called recursively. internal static Entity? ComputeIndefiniteIntegral(Entity expr, Entity.Variable x, bool integrateByParts = true) { + expr = Normalized(expr); if (!expr.ContainsNode(x)) return expr * x; // base case, handle here if ((IntegralPatterns.TryStandardIntegrals(expr, x)) is { } answer) return answer; // The flag has to be handed on. Every one of these recurses, and a solver that diff --git a/Sources/AngouriMath/Functions/Simplification/Patterns/Patterns.Power.cs b/Sources/AngouriMath/Functions/Simplification/Patterns/Patterns.Power.cs index c2d9c3825..27f154766 100644 --- a/Sources/AngouriMath/Functions/Simplification/Patterns/Patterns.Power.cs +++ b/Sources/AngouriMath/Functions/Simplification/Patterns/Patterns.Power.cs @@ -143,6 +143,110 @@ when ReduceRadical(radicand, power) is { } reduced => reduced, _ => x }; + /// + /// Gathers the factors of a product that are powers of one base, wherever they sit + /// in it: a^n * c * a^m becomes a^(n+m) * c. + /// + /// + /// + /// already has {}^n * {}^m = {}^(n+m), but it pairs + /// two sibling nodes, and a product is a tree rather than a list -- in + /// sin(x)^4 * (-6) * sin(x)^2 the constant sits between the two powers, so + /// they are never siblings and the rule never fires. Full + /// gets these because it reassociates and sorts the factors first; nothing that + /// normalises without sorting does. + /// + /// + /// a^n * a^m = a^(n+m) needs no condition: with a^n read as + /// e^(n Log a) on the principal branch, the two sides are + /// e^(n Log a) * e^(m Log a) and e^((n+m) Log a), which are equal for + /// every complex n and m. This is what makes it unlike + /// (a^b)^c = a^(b*c), which moves the branch and is guarded above. + /// The one point it does not cover is a = 0 with a negative exponent, where + /// 0^2 * 0^(-1) is undefined and 0^1 is 0 -- so this is applied to + /// integrands, where an antiderivative may differ on a measure-zero set, and not in + /// . https://github.com/asc-community/AngouriMath/issues/781 + /// + /// + internal static Entity GatherPowersOfOneBase(Entity x) + { + if (x is not (Mulf or Divf)) + return x; + + // Exponent null marks a factor that is carried through untouched. + var factors = new List<(Entity Base, Entity? Exponent)>(); + var merged = false; + foreach (var factor in Mulf.LinearChildren(x)) + { + var (@base, exponent) = Decompose(factor); + // Numeric factors are left to evaluation, which folds 2 * 2 into 4. Gathering + // them here would write it as 2^2 and call that progress. + if (@base is Number) + { + factors.Add((factor, null)); + continue; + } + var found = false; + for (var i = 0; i < factors.Count; i++) + if (factors[i].Exponent is { } already && factors[i].Base == @base) + { + // Only the exponent is folded, never the product around it. + // InnerSimplified on the whole expression rewrites x^(-2) back into + // 1/x^2, and SolveAsPolynomialTerm turns a 1/x^n it is handed into + // Pow(x, -n) again -- the two normalisations chase each other until + // the stack runs out. + factors[i] = (@base, (already + exponent).InnerSimplified); + merged = found = true; + break; + } + if (!found) + factors.Add((@base, exponent)); + } + + // Rebuilding unconditionally would rewrite every quotient in the tree as a + // negative power, since LinearChildren flattens a / b into a * b^(-1). + if (!merged) + return x; + + Entity? result = null; + foreach (var (@base, exponent) in factors) + { + var factor = exponent switch + { + null => @base, + Integer(1) => @base, + _ => new Powf(@base, exponent) + }; + result = result is null ? factor : result * factor; + } + return result ?? x; + } + + /// + /// Reads a factor of a product as a base and an exponent. + /// + /// + /// The nested case is not cosmetic: writes a + /// divisor as (...)^(-1), so x^3 / x^2 arrives as + /// x^3 * (x^2)^(-1), whose second factor has base x^2 rather than + /// x and would be read as an unrelated base. Unwrapping it is + /// (a^b)^n = a^(b*n), which holds for a whole n whatever the sign of + /// a -- the same guard the ({}^{})^{} rule above carries, and for the + /// same reason. https://github.com/asc-community/AngouriMath/issues/752 + /// + private static (Entity Base, Entity Exponent) Decompose(Entity factor) + { + if (factor is not Powf(var @base, var exponent)) + return (factor, 1); + while (@base is Powf(var inner, var innerExponent) + && (exponent is Integer || inner.Evaled is Real { IsPositive: true })) + { + exponent = innerExponent * exponent; + @base = inner; + } + return (@base, exponent); + } + /// /// Largest divisor tried when reducing a radical. Reducing sqrt(n) exactly would /// mean factoring n, which is not something a simplification pass can afford to diff --git a/Sources/Tests/UnitTests/Calculus/ConstantRateIntegrandTest.cs b/Sources/Tests/UnitTests/Calculus/ConstantRateIntegrandTest.cs new file mode 100644 index 000000000..b35c7f614 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/ConstantRateIntegrandTest.cs @@ -0,0 +1,97 @@ +// +// 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 +{ + /// + /// Every standard-integral rule reads a linear rate out of its integrand's argument and + /// divides by it. An argument that mentions x without depending on it has a rate + /// of zero, and the division was written into the answer unguarded: + /// + /// e ^ (x + -x) -> e ^ (x + -x) / (0 * ln(e)) -- NaN, where it is x + /// + /// https://github.com/asc-community/AngouriMath/issues/785 + /// Found while fixing https://github.com/asc-community/AngouriMath/issues/781, whose + /// gathering of two powers of one base turns e^x * e^(-x) into exactly this + /// shape -- but the defect predates it and is reachable by writing the exponent out. + /// + public sealed class ConstantRateIntegrandTest + { + /// + /// A zero rate means the integrand is a constant, so the antiderivative is that + /// constant times x. Checked by differentiating back rather than by its printed + /// form, which still carries the unreduced argument. + /// + [Theory] + [InlineData("e ^ (x + -x)", "1")] + [InlineData("sin(x + -x)", "0")] + [InlineData("cos(x - x)", "1")] + [InlineData("2 ^ (x - x)", "1")] + public void AZeroRateIntegrandIsAConstant(string integrand, string value) + { + var antiderivative = integrand.Integrate("x"); + Assert.False(antiderivative is Entity.Integralf, + $"{integrand} was declined: {antiderivative.Stringize()}"); + var difference = (antiderivative.Differentiate("x") - value.ToEntity()).Simplify(); + while (difference is Entity.Providedf(var inner, _)) difference = inner; + Assert.Equal(Entity.Number.Integer.Create(0), difference); + } + + /// + /// The guard must fire on a zero rate and on nothing else -- a rate of one or of + /// any other constant is what the rules below it are for, and reading those as + /// constants would replace every one of these answers with the integrand times x. + /// + [Theory] + [InlineData("e ^ (2 * x)")] + [InlineData("sin(3 * x + 1)")] + [InlineData("cos(x)")] + [InlineData("ln(2 * x)")] + public void ANonZeroRateStillIntegratesByItsRule(string integrand) + { + var antiderivative = integrand.Integrate("x"); + Assert.False(antiderivative is Entity.Integralf, + $"{integrand} was declined: {antiderivative.Stringize()}"); + var difference = (antiderivative.Differentiate("x") - integrand.ToEntity()).Simplify(); + while (difference is Entity.Providedf(var inner, _)) difference = inner; + Assert.Equal(Entity.Number.Integer.Create(0), difference); + } + + /// + /// abs takes the guard too, and is checked on the antiderivative rather than + /// by differentiating it back: the library differentiates abs(u) as + /// sgn(u) * u' and answers NaN at u = 0, so the round trip + /// says nothing about this rule either way. + /// + [Fact] + public void AZeroRateAbsIntegrandIsZero() + { + var antiderivative = "abs(x + -x)".Integrate("x"); + Assert.False(antiderivative is Entity.Integralf); + Assert.Equal(Entity.Number.Integer.Create(0), + (antiderivative - "C".ToEntity()).Simplify()); + } + + /// + /// A symbolic rate is not decidably zero, so the guard must leave it to the rules -- + /// answering sin(a*x) * x here would be wrong for every non-zero a. + /// + [Fact] + public void ASymbolicRateIsLeftToTheRules() + { + var antiderivative = "sin(a * x)".Integrate("x"); + Assert.False(antiderivative is Entity.Integralf); + var difference = (antiderivative.Differentiate("x") - "sin(a * x)".ToEntity()).Simplify(); + while (difference is Entity.Providedf(var inner, _)) difference = inner; + Assert.Equal(Entity.Number.Integer.Create(0), difference); + } + } +} diff --git a/Sources/Tests/UnitTests/Calculus/PowersOfOneBaseMergedTest.cs b/Sources/Tests/UnitTests/Calculus/PowersOfOneBaseMergedTest.cs new file mode 100644 index 000000000..c03abb9d6 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/PowersOfOneBaseMergedTest.cs @@ -0,0 +1,119 @@ +// +// 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 System; +using AngouriMath; +using AngouriMath.Extensions; +using Xunit; + +namespace AngouriMath.Tests.Calculus +{ + /// + /// Two powers of one base sitting in the same product were never merged, so an + /// integrand that is a power in disguise was declined: + /// + /// x ^ 2 / x -> integral(x ^ 2 / x, x) -- it is x + /// sin(x) ^ 4 * (-6) * sin(x) ^ 2 -> declined -- it is -6 sin(x)^6 + /// + /// The rule that merges them exists in Patterns.PowerRules, but it pairs two + /// sibling nodes, and a constant factor sitting between the two powers makes them + /// non-siblings. Integrate normalises with InnerSimplified, which does + /// not run PowerRules at all, so even the sibling case was missed. + /// https://github.com/asc-community/AngouriMath/issues/781 + /// + public sealed class PowersOfOneBaseMergedTest + { + /// Sample points, chosen to include negative x and to avoid the zeros of sin. + private static readonly double[] Points = { -1.3, -0.4, 0.25, 0.7, 1.9, 3.3 }; + + /// Relative, because these integrands range over several orders of magnitude. + private const double RelativeTolerance = 1e-9; + + private static Entity AssertAnswered(string integrand) + { + var antiderivative = integrand.Integrate("x"); + Assert.False(antiderivative is Entity.Integralf, + $"{integrand} was declined: {antiderivative.Stringize()}"); + return antiderivative; + } + + /// + /// The antiderivative is checked by differentiating it back, rather than against a + /// printed form -- the point is that an answer exists and is right, not how it is + /// spelled. + /// + private static void AssertIntegrates(string integrand, string expectedDerivative) + { + var difference = (AssertAnswered(integrand).Differentiate("x") - expectedDerivative.ToEntity()).Simplify(); + while (difference is Entity.Providedf(var inner, _)) difference = inner; + Assert.Equal(Entity.Number.Integer.Create(0), difference); + } + + /// + /// The same check for integrands whose differentiated-back form is a trigonometric + /// identity that cannot close -- for those, asserting + /// a symbolic zero would be testing the simplifier's reach rather than the + /// integrator's answer. The antiderivative is sampled instead. + /// + private static void AssertIntegratesNumerically(string integrand, string expectedDerivative) + { + var derivative = AssertAnswered(integrand).Differentiate("x"); + foreach (var point in Points) + { + var actual = derivative.Substitute("x", point).Substitute("C", 0) + .EvalNumerical().RealPart.EDecimal.ToDouble(); + var expected = expectedDerivative.ToEntity().Substitute("x", point) + .EvalNumerical().RealPart.EDecimal.ToDouble(); + var scale = Math.Max(Math.Max(Math.Abs(expected), Math.Abs(actual)), 1e-12); + Assert.True(Math.Abs(expected - actual) <= RelativeTolerance * scale, + $"{integrand} at x = {point}: differentiated back to {actual}, expected {expected}"); + } + } + + [Theory] + [InlineData("x ^ 2 / x", "x")] + [InlineData("x ^ 2 * (1 / x)", "x")] + [InlineData("x ^ 3 / x ^ 2", "x")] + public void PowersOfOneBaseAreMerged(string integrand, string expectedDerivative) + => AssertIntegrates(integrand, expectedDerivative); + + [Theory] + [InlineData("sin(x) ^ 4 * (-6) * sin(x) ^ 2", "-6 * sin(x) ^ 6")] + [InlineData("sin(x) ^ 4 * 3 * sin(x) ^ 2", "3 * sin(x) ^ 6")] + public void PowersOfOneTrigBaseAreMerged(string integrand, string expectedDerivative) + => AssertIntegratesNumerically(integrand, expectedDerivative); + + /// + /// Every argument other than a bare x already answered before the fix, and + /// must keep answering. These are the control cases from the issue. + /// + [Theory] + [InlineData("x ^ 2 * x", "x ^ 3")] + public void AlreadyAnsweredCasesStillAnswer(string integrand, string expectedDerivative) + => AssertIntegrates(integrand, expectedDerivative); + + [Theory] + [InlineData("sin(2 * x) ^ 4 * (-6) * sin(2 * x) ^ 2", "-6 * sin(2 * x) ^ 6")] + [InlineData("sin(x + 1) ^ 4 * (-6) * sin(x + 1) ^ 2", "-6 * sin(x + 1) ^ 6")] + public void AlreadyAnsweredTrigCasesStillAnswer(string integrand, string expectedDerivative) + => AssertIntegratesNumerically(integrand, expectedDerivative); + + /// + /// Merging powers of one base must not disturb a product whose factors have + /// different bases -- the gathering is keyed on the base, and a wrong key here + /// would silently rewrite unrelated factors together. + /// + [Theory] + [InlineData("x ^ 2 * a ^ 3")] + [InlineData("sin(x) ^ 2 * cos(x) ^ 3")] + public void DifferentBasesAreLeftAlone(string expr) + { + var gathered = expr.ToEntity().InnerSimplified; + Assert.Equal(MathS.Boolean.True, gathered.EqualTo(expr.ToEntity()).Simplify()); + } + } +}