diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index bcf774609..957d4c65e 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -476,6 +476,20 @@ Pythagorean identity with a coefficient is now a rule of its own, found wherever | `"2 + (-2)*cos(x)^2 + (-2)*sin(x)^2".Simplify()` | `2 + (-2) * cos(x) ^ 2 + (-2) * sin(x) ^ 2` | `0` | | `"2 + (-2) * cos(2 * x) + (-4) * sin(x) ^ 2".Simplify()` | unchanged | `0` | +### A compiled odd root of a negative real is real + +**Silent: compiled, a negative real to a rational power with an odd denominator below ten is the real root to the +numerator, as evaluation has it, where both compilers took the principal power, and `Compile` gave +`NaN`.** `(-8)^(1/3)` evaluates to `-2`, `(-8)^(2/3)` to `4`. Newton's method, which `Solve` falls back to, runs on the +compiled function, so a root needing the real root was not found either. A base off the real line, a positive one, an +even denominator or one of ten or more is the principal power, as before. + +| Input | Was (2.6.0) | Now | +|---|---|---| +| `"x^(1/3)".ToEntity().Compile("x").Call(-8)` | `<1.0000000000000002; 1.7320508075688772>` | `<-2; 0>` | +| `"x^(1/3)".ToEntity().Compile("x")(-8)` | `NaN` | `-2` | +| `"(3 - 4x)^(1/3) + sqrt(4 + 4x) = 3".ToEntity().Solve("x")` | `{ (-2.466…) / 4, 2.952… / 4 }` | `{ (-2.466…) / 4, 2.952… / 4, 33.514… / 4 }` | + ### A compiled division by zero over `double` has no value **Silent: compiled over `double` or `float`, a division by zero and 0 to a power that is not positive are `NaN`, as in diff --git a/Sources/.editorconfig b/Sources/.editorconfig index ece1ed53a..01f7dd40c 100644 --- a/Sources/.editorconfig +++ b/Sources/.editorconfig @@ -1034,6 +1034,9 @@ file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed [Tests/UnitTests/Common/HyperbolicTangentPastTheRangeTest.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/Common/CompiledOddRootTest.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/Common/RefuseNonRealCompilationTest.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/Compilation/IntoFE/Compiler.cs b/Sources/AngouriMath/Functions/Compilation/IntoFE/Compiler.cs index c33fe5e8e..efde9d334 100644 --- a/Sources/AngouriMath/Functions/Compilation/IntoFE/Compiler.cs +++ b/Sources/AngouriMath/Functions/Compilation/IntoFE/Compiler.cs @@ -166,11 +166,16 @@ public partial record Powf { private protected override bool CompilesItself => true; + // A rational exponent with an odd denominator below ten takes the real root of a negative real, + // as evaluation does: (-8)^(1/3) is -2, which the principal power made 1 + 1.73i. private protected override void CompileNode(Compiler compiler) { Exponent.InnerCompile(compiler); Base.InnerCompile(compiler); - compiler.Instructions.Add(new(InstructionType.CALL_POW)); + compiler.Instructions.Add( + Exponent is Rational { ERational: var ratio } && Core.Compilation.RealProducts.TakesTheRealRootOfANegative(ratio) + ? new(InstructionType.CALL_REAL_ROOT_POW, ratio.Numerator.IsEven ? 0 : 1) + : new(InstructionType.CALL_POW)); } } diff --git a/Sources/AngouriMath/Functions/Compilation/IntoFE/FastExpression.cs b/Sources/AngouriMath/Functions/Compilation/IntoFE/FastExpression.cs index 2f5fb9f93..37ce14fd9 100644 --- a/Sources/AngouriMath/Functions/Compilation/IntoFE/FastExpression.cs +++ b/Sources/AngouriMath/Functions/Compilation/IntoFE/FastExpression.cs @@ -72,6 +72,8 @@ internal enum InstructionType CALL_DIV, CALL_MOD, CALL_POW, + // A power to a rational with an odd denominator below ten; Reference is 1 where its numerator is odd. + CALL_REAL_ROOT_POW, CALL_LOG, CALL_MAX, CALL_MIN, @@ -231,6 +233,9 @@ public System.Numerics.Complex Substitute(params System.Numerics.Complex[] value case InstructionType.CALL_POW: stack.Push(AngouriMath.Core.Compilation.RealProducts.Power(stack.Pop(), stack.Pop())); break; + case InstructionType.CALL_REAL_ROOT_POW: + stack.Push(AngouriMath.Core.Compilation.RealProducts.RealRootPower(stack.Pop(), stack.Pop(), instruction.Reference == 1)); + break; case InstructionType.CALL_SIN: stack.Push(System.Numerics.Complex.Sin(stack.Pop())); break; diff --git a/Sources/AngouriMath/Functions/Compilation/IntoLinq/CompilationProtocol.cs b/Sources/AngouriMath/Functions/Compilation/IntoLinq/CompilationProtocol.cs index 7643d85c8..62cd9670e 100644 --- a/Sources/AngouriMath/Functions/Compilation/IntoLinq/CompilationProtocol.cs +++ b/Sources/AngouriMath/Functions/Compilation/IntoLinq/CompilationProtocol.cs @@ -262,6 +262,14 @@ Divf when ShouldBeAtLeastDouble(left) is var newLeft && ShouldBeAtLeastDouble(ri Modf => Operator(ExpressionType.Modulo, Operator(ExpressionType.Add, Operator(ExpressionType.Modulo, left, right), right), right), + // A rational exponent with an odd denominator below ten takes the real root of a negative + // real, as evaluation does: (-8)^(1/3) is -2, which was 1 + 1.73i, and NaN over double. + Powf { Exponent: Rational { ERational: var ratio } } when Core.Compilation.RealProducts.TakesTheRealRootOfANegative(ratio) + && ShouldBeAtLeastDouble(left) is var newLeft + && ShouldBeAtLeastDouble(right) is var newRight + && (newLeft.Type == typeof(System.Numerics.Complex) || newLeft.Type == typeof(double)) + => Expression.Call(newLeft.Type == typeof(double) ? realRootPowerOfReals : realRootPower, + newLeft, newRight, Expression.Constant(!ratio.Numerator.IsEven)), Powf when ShouldBeAtLeastDouble(left) is var newLeft && ShouldBeAtLeastDouble(right) is var newRight @@ -577,6 +585,10 @@ internal Expression ConvertMembership(Expression element, AngouriMath.Numerics.N : Expression.Call(realMembership, ConvertType(widened, typeof(double)), Expression.Constant((int)set)); } + [ConstantField] private static readonly MethodInfo realRootPower = + ((MethodCallExpression)((Expression>)((a, b, odd) => Core.Compilation.RealProducts.RealRootPower(a, b, odd))).Body).Method; + [ConstantField] private static readonly MethodInfo realRootPowerOfReals = + ((MethodCallExpression)((Expression>)((a, b, odd) => Core.Compilation.RealProducts.RealRootPower(a, b, odd))).Body).Method; [ConstantField] private static readonly MethodInfo complexRemainder = ((MethodCallExpression)((Expression>)((a, b) => RealOnly.Mod(a, b))).Body).Method; diff --git a/Sources/AngouriMath/Functions/Compilation/RealProducts.cs b/Sources/AngouriMath/Functions/Compilation/RealProducts.cs index e7ae386e1..9fa4c7004 100644 --- a/Sources/AngouriMath/Functions/Compilation/RealProducts.cs +++ b/Sources/AngouriMath/Functions/Compilation/RealProducts.cs @@ -79,6 +79,33 @@ public static NumericsComplex Power(NumericsComplex a, NumericsComplex b) private static bool IsWhole(double y) => !double.IsInfinity(y) && Math.Floor(y) == y; + /// + /// A power to a rational exponent whose denominator is odd and below ten, as evaluation takes it: + /// of a negative real, the real root to the numerator, (-8)^(1/3) = -2 and (-8)^(2/3) = 4, + /// as Number.Pow's search for a good root finds it; of any other number, . + /// The principal power gave 1 + 1.73i and -2 + 3.46i, and Math.Pow NaN over + /// double, where evaluation has a real value. Whether the exponent is such a rational is only + /// known where it is written, so the compilers choose this where they compile the power. + /// + public static NumericsComplex RealRootPower(NumericsComplex a, NumericsComplex b, bool oddNumerator) + => a.Imaginary == 0 && a.Real < 0 && b.Imaginary == 0 && !double.IsNaN(b.Real) + ? new NumericsComplex((oddNumerator ? -1 : 1) * Math.Pow(-a.Real, b.Real), 0) + : Power(a, b); + + /// over . + public static double RealRootPower(double a, double b, bool oddNumerator) + => a < 0 && !double.IsNaN(b) + ? (oddNumerator ? -1 : 1) * Math.Pow(-a, b) + : IntoLinq.MathAllMethods.Pow(a, b); + + /// + /// Whether a power to takes the real root of a negative real, in evaluation: + /// a rational that is not whole, with an odd denominator below ten. + /// + internal static bool TakesTheRealRootOfANegative(PeterO.Numbers.ERational exponent) + => !exponent.Denominator.Equals(PeterO.Numbers.EInteger.One) && !exponent.Denominator.IsEven + && exponent.Denominator.CompareTo(PeterO.Numbers.EInteger.FromInt32(10)) < 0; + /// /// A whole power of a number that is not real, as a product of its factors, which is how evaluation takes it. /// takes e^(n ln z), diff --git a/Sources/Tests/UnitTests/Common/CompiledOddRootTest.cs b/Sources/Tests/UnitTests/Common/CompiledOddRootTest.cs new file mode 100644 index 000000000..c1589d580 --- /dev/null +++ b/Sources/Tests/UnitTests/Common/CompiledOddRootTest.cs @@ -0,0 +1,74 @@ +// +// 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 AngouriMath.Extensions; +using Xunit; +using static AngouriMath.Entity; +using static AngouriMath.Entity.Set; +using NumericsComplex = System.Numerics.Complex; + +namespace AngouriMath.Tests.Common +{ + /// + /// A negative real to a rational power with an odd denominator below ten is the real root to the + /// numerator, in evaluation: (-8)^(1/3) is -2. Both compilers took the principal + /// power, 1 + 1.73i, and NaN over double; and since Newton's method runs on + /// the compiled function, the roots that need the real root were not found: + /// (3 - 4x)^(1/3) + sqrt(4 + 4x) = 3 was answered without 8.379. Found by a + /// cross-check of equations in roots against a Newton search; no issue was filed. + /// + [Trait("Area", "Compilation")] + public sealed class CompiledOddRootTest + { + [Theory] + [InlineData("x^(1/3)", -8)] + [InlineData("x^(2/3)", -8)] + [InlineData("x^(-1/3)", -8)] + [InlineData("x^(4/3)", -8)] + [InlineData("x^(3/5)", -32)] + [InlineData("(1 + x)^(1/3) + 2", -9)] + public void TheCompiledFormsTakeTheRealRoot(string expression, double at) + { + var entity = expression.ToEntity(); + var value = entity.Substitute("x", (Entity)(int)at).EvalNumerical(); + Assert.True(value.ImaginaryPart.IsZero, $"{expression} at {at} evaluates to {value}"); + var expected = value.RealPart.EDecimal.ToDouble(); + Assert.Equal(expected, entity.Compile("x").Call(new NumericsComplex(at, 0)).Real, 12); + Assert.Equal(0, entity.Compile("x").Call(new NumericsComplex(at, 0)).Imaginary, 12); + Assert.Equal(expected, entity.Compile("x")(new NumericsComplex(at, 0)).Real, 12); + Assert.Equal(0, entity.Compile("x")(new NumericsComplex(at, 0)).Imaginary, 12); + Assert.Equal(expected, entity.Compile("x")(at), 12); + } + + /// + /// The rest is the principal power, as in evaluation: a base off the real line, a positive one, + /// an even denominator and one of ten or more. + /// + [Theory] + [InlineData("x^(1/3)", -8, 1)] + [InlineData("x^(1/3)", 8, 0)] + [InlineData("x^(1/2)", -4, 0)] + [InlineData("x^(1/11)", -1, 0)] + public void TheRestIsThePrincipalPower(string expression, double re, double im) + { + var entity = expression.ToEntity(); + var value = entity.Substitute("x", Number.Complex.Create(re, im)).EvalNumerical().ToNumerics(); + var compiled = entity.Compile("x").Call(new NumericsComplex(re, im)); + Assert.True((compiled - value).Magnitude < 1e-12, $"{expression} at {re} + {im}i: evaluated {value}, compiled {compiled}"); + var linq = entity.Compile("x")(new NumericsComplex(re, im)); + Assert.True((linq - value).Magnitude < 1e-12, $"{expression} at {re} + {im}i: evaluated {value}, LINQ {linq}"); + } + + /// The root Newton's method could not see. + [Fact] + public void TheSearchFindsTheRootThatNeedsTheRealRoot() + { + var answer = Assert.IsType("(3 - 4x)^(1/3) + sqrt(4 + 4x) = 3".ToEntity().Solve("x")); + Assert.Contains(answer, root => System.Math.Abs(root.EvalNumerical().RealPart.EDecimal.ToDouble() - 8.3785553513) < 1e-8); + } + } +}