diff --git a/Sources/AngouriMath/Functions/Evaluation/Evaluation.Definition.cs b/Sources/AngouriMath/Functions/Evaluation/Evaluation.Definition.cs index 11bad1196..724c1bfec 100644 --- a/Sources/AngouriMath/Functions/Evaluation/Evaluation.Definition.cs +++ b/Sources/AngouriMath/Functions/Evaluation/Evaluation.Definition.cs @@ -333,8 +333,8 @@ private static Entity CollectLikeTerms(Entity expanded) public Entity Factorize(int level = 2) => level <= 1 // InnerSimplified, so that the factors come back finished: the rules leave // x ^ 1 where they mean x, and sqrt(4) where they mean 2. - ? this.Replace(Patterns.FactorizeRules).InnerSimplified - : this.Replace(Patterns.FactorizeRules).InnerSimplified.Factorize(level - 1); + ? this.Replace(Patterns.PerfectSquareRules).Replace(Patterns.FactorizeRules).InnerSimplified + : this.Replace(Patterns.PerfectSquareRules).Replace(Patterns.FactorizeRules).InnerSimplified.Factorize(level - 1); /// /// Simplifies an equation ( e.g. (x - y) * (x + y) -> x^2 - y^2, but 3 * x + y * x = (3 + y) * x ) diff --git a/Sources/AngouriMath/Functions/Simplification/Patterns/Patterns.ExpandFactorize.cs b/Sources/AngouriMath/Functions/Simplification/Patterns/Patterns.ExpandFactorize.cs index 3637ff616..34b418076 100644 --- a/Sources/AngouriMath/Functions/Simplification/Patterns/Patterns.ExpandFactorize.cs +++ b/Sources/AngouriMath/Functions/Simplification/Patterns/Patterns.ExpandFactorize.cs @@ -5,6 +5,7 @@ // Website: https://am.angouri.org. // +using PeterO.Numbers; using static AngouriMath.Entity; namespace AngouriMath.Functions @@ -64,6 +65,116 @@ Sumf or Minusf when CollectCommonFactors(x) is { } collected => collected, _ => x }; + + /// + /// A three-term sum that is a perfect square with a radical in it. + /// + /// + /// Its own pass, run before rather than as an arm of + /// it. Replace walks bottom-up, so the common-factor rule reaches the inner + /// 4 + 4*sqrt(x) of 4 + 4*sqrt(x) + x first and rewrites it to + /// 4 * (1 + sqrt(x)) -- by the time the outer sum is looked at, the three + /// terms this needs are no longer there. Ordering the arms within the switch does + /// not help, because the two rules are looking at different nodes. + /// https://github.com/asc-community/AngouriMath/issues/176 + /// + internal static Entity PerfectSquareRules(Entity x) + => x is Sumf or Minusf && CollapseToPerfectSquare(x) is { } square ? square : x; + + /// + /// u + 2*sqrt(u)*sqrt(v) + v collapses to (sqrt(u) + sqrt(v))^2, which + /// is what 1 + sqrt(2x) + x/2 is. + /// https://github.com/asc-community/AngouriMath/issues/176 + /// + /// + /// + /// The identity itself is unconditional: (sqrt(u) + sqrt(v))^2 expands to + /// u + 2*sqrt(u)*sqrt(v) + v because sqrt(u)^2 is u for every + /// complex u. It is the *square of a principal root*, not the root of a + /// square -- sqrt(u^2) = u is the false one, and holds only for a + /// non-negative u. https://github.com/asc-community/AngouriMath/issues/752 + /// + /// + /// What cannot be trusted is the test for whether the cross term matches. Deciding + /// it needs , and the simplifier equates + /// sqrt(x)*sqrt(y) with sqrt(x*y), which is false on the branch cuts: + /// at x = y = -1, x + 2*sqrt(x*y) + y is 0 while + /// (sqrt(x) + sqrt(y))^2 is -4. Asked symbolically, this rule fired on that + /// sum and produced a wrong answer. + /// + /// + /// So the symbolic match only proposes, and a numeric check at sample points + /// disposes. The points include negative values, which is where a branch-cut error + /// shows and nowhere else, and every free variable is given a different one so that + /// x and y cannot coincide into a case that happens to hold. The rule + /// withdraws unless every sampled point agrees, so a variable it cannot evaluate at + /// simply means no collapse. + /// + /// + /// Restricted to sums that contain a radical, both because that is where the gap is + /// -- a polynomial trinomial is collapsed by the rules above -- and to keep the cost + /// of that Simplify off every three-term sum in every tree. + /// + /// + /// when the sum is not a square, so the rule does not fire. + private static Entity? CollapseToPerfectSquare(Entity expr) + { + if (!expr.Nodes.Any(node => node is Powf(_, Rational and not Integer))) + return null; + var terms = Sumf.LinearChildren(expr).ToList(); + if (terms.Count != 3) + return null; + + for (var cross = 0; cross < 3; cross++) + { + var w = terms[cross]; + var p = new Powf(terms[(cross + 1) % 3], Rational.Create(1, 2)); + var q = new Powf(terms[(cross + 2) % 3], Rational.Create(1, 2)); + var product = (2 * p * q).Simplify(); + var candidate = + (product - w).Simplify() == 0 ? new Powf((p + q).Simplify(), 2) : + (product + w).Simplify() == 0 ? new Powf((p - q).Simplify(), 2) : + null; + if (candidate is { } square && AgreesNumerically(expr, square)) + return square; + } + return null; + } + + /// Sample points, negative first: a branch-cut error shows nowhere else. + [ConstantField] private static readonly EDecimal[] PerfectSquareSamplePoints = + { EDecimal.FromString("-1.7"), EDecimal.FromString("-0.6"), + EDecimal.FromString("0.8"), EDecimal.FromString("2.3") }; + + /// + /// Whether a proposed collapse is the same number as what it replaces, checked at + /// with each free variable offset from the + /// last so that two of them never take the same value. + /// + private static bool AgreesNumerically(Entity original, Entity candidate) + { + var variables = original.Vars.ToList(); + if (variables.Count == 0) + return true; + for (var i = 0; i < PerfectSquareSamplePoints.Length; i++) + { + Entity before = original, after = candidate; + for (var v = 0; v < variables.Count; v++) + { + var point = PerfectSquareSamplePoints[(i + v) % PerfectSquareSamplePoints.Length]; + before = before.Substitute(variables[v], Real.Create(point)); + after = after.Substitute(variables[v], Real.Create(point)); + } + if (before.Evaled is not Complex left || after.Evaled is not Complex right) + return false; + var difference = (left - right).Abs(); + if (difference is not Real real || real.EDecimal.Abs() + .CompareTo(EDecimal.FromString("1e-12")) > 0) + return false; + } + return true; + } + /// /// Pulls a shared factor out of the terms of a sum: a*c + a*d + b*c + b*d /// becomes a*(c + d) + b*(c + d), which the pairwise rules above close to diff --git a/Sources/Tests/UnitTests/PatternsTest/CollapsePerfectSquareTest.cs b/Sources/Tests/UnitTests/PatternsTest/CollapsePerfectSquareTest.cs new file mode 100644 index 000000000..524b473dd --- /dev/null +++ b/Sources/Tests/UnitTests/PatternsTest/CollapsePerfectSquareTest.cs @@ -0,0 +1,116 @@ +// +// 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 System.Linq; +using AngouriMath; +using AngouriMath.Extensions; +using Xunit; + +namespace AngouriMath.Tests.PatternsTest +{ + /// + /// A three-term sum that is a perfect square was left as written when one of its terms + /// is a radical: 1 + sqrt(2x) + x/2 is (1 + sqrt(x/2))^2. + /// https://github.com/asc-community/AngouriMath/issues/176 + /// https://github.com/asc-community/AngouriMath/issues/203 + /// + [Trait("Area", "PatternsTest")] + public sealed class CollapsePerfectSquareTest + { + private static readonly double[] Points = { 0.3, 1.4, 2.7, 5.1, 8.6 }; + + /// + /// Sampled rather than compared as text: the square can be spelled several equal + /// ways -- sqrt(x/2) and sqrt(2x)/2 are one number -- and what matters is that the + /// collapsed form is the same function, not how it prints. + /// + private static void AssertSameValue(string original, Entity collapsed) + { + var compared = 0; + foreach (var point in Points) + { + double before, after; + try + { + before = original.ToEntity().Substitute("x", point).EvalNumerical() + .RealPart.EDecimal.ToDouble(); + after = collapsed.Substitute("x", point).EvalNumerical() + .RealPart.EDecimal.ToDouble(); + } + catch { continue; } + if (double.IsNaN(before) || double.IsNaN(after)) continue; + compared++; + var scale = Math.Max(Math.Max(Math.Abs(before), Math.Abs(after)), 1e-12); + Assert.True(Math.Abs(before - after) <= 1e-9 * scale, + $"{original} collapsed to {collapsed.Stringize()}, which at x = {point} " + + $"is {after} rather than {before}"); + } + Assert.True(compared > 0, $"{original}: no point was comparable, so nothing was checked"); + } + + [Theory] + [InlineData("1 + sqrt(2 * x) + x / 2")] + [InlineData("1 + 2 * sqrt(x) + x")] + [InlineData("4 + 4 * sqrt(x) + x")] + public void ASumWithARadicalCollapsesToASquare(string expr) + { + var collapsed = expr.ToEntity().Factorize(); + Assert.Contains(collapsed.Nodes, + node => node is Entity.Powf(_, Entity.Number.Integer(2))); + AssertSameValue(expr, collapsed); + } + + /// + /// A sum that is not a square must be left alone -- the cross term has to be + /// exactly twice the product of the two roots, and a rule that rounded that off + /// would be inventing an identity. + /// + [Theory] + [InlineData("1 + 3 * sqrt(x) + x")] + [InlineData("1 + sqrt(x) + x")] + [InlineData("1 + sqrt(2 * x) + x / 3")] + public void ASumThatIsNotASquareIsLeftAlone(string expr) + => AssertSameValue(expr, expr.ToEntity().Factorize()); + + /// + /// sqrt(u)^2 is u for every complex u, but sqrt(u^2) is + /// not u -- it is u only where u is non-negative. So a + /// polynomial trinomial must not be collapsed through this route, which would be + /// #752 all over again: at x = -3, x^2 + 2x + 1 is 4 and + /// (sqrt(x^2) + 1)^2 is 16. + /// + /// + /// Not a square over the complex plane, though it reads like one: sqrt(x)*sqrt(y) + /// is sqrt(x*y) only for suitable branches. At x = y = -1 the sum is 0 while + /// (sqrt(x) + sqrt(y))^2 is -4, so collapsing it would be #752 in a new place. + /// + [Fact] + public void ATwoVariableSumThatOnlyLooksLikeASquareIsLeftAlone() + { + var collapsed = "x + 2 * sqrt(x * y) + y".ToEntity().Factorize(); + var at = collapsed.Substitute("x", -1).Substitute("y", -1).EvalNumerical(); + Assert.True(Math.Abs(at.RealPart.EDecimal.ToDouble()) < 1e-9, + $"x + 2*sqrt(x*y) + y collapsed to {collapsed.Stringize()}, which at x = y = -1 " + + $"is {at.Stringize()} rather than 0"); + } + + [Theory] + [InlineData("x ^ 2 + 2 * x + 1", -3.0)] + [InlineData("x ^ 2 - 2 * x + 1", -3.0)] + public void APolynomialTrinomialKeepsItsValueAtANegativePoint(string expr, double point) + { + var collapsed = expr.ToEntity().Factorize(); + var before = expr.ToEntity().Substitute("x", point).EvalNumerical() + .RealPart.EDecimal.ToDouble(); + var after = collapsed.Substitute("x", point).EvalNumerical() + .RealPart.EDecimal.ToDouble(); + Assert.True(Math.Abs(before - after) < 1e-9, + $"{expr} collapsed to {collapsed.Stringize()}, which at x = {point} is {after} rather than {before}"); + } + } +}