From 5ec850145ee8df9824d4955c3dc0743da850cca9 Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Sun, 27 Sep 2026 15:04:49 +0000 Subject: [PATCH 1/3] The Gaussian, its even moments and the error functions are integrated int F^(a x^2 + b x + c) dx is F^(c - b^2/(4a)) times the Gaussian of A = a ln F in u = x + b/(2a), and int e^(A u^2) du is sqrt(pi)/(2 sqrt(-A)) erf(sqrt(-A) u), written with erfi and the real root where A is decidably positive. Beside an even power of x, on either side of 0, parts reduce it two powers at a time to the Gaussian; an odd power ends at an elementary integral or at the exponential integral and is not taken. erf, erfc and erfi are integrated by parts against 1. Part of #1501. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura --- BREAKING-CHANGES.md | 18 +++ Sources/.editorconfig | 3 + .../Integration/IntegralPatterns.cs | 136 ++++++++++++++++++ .../Calculus/GaussianIntegralTest.cs | 132 +++++++++++++++++ 4 files changed, 289 insertions(+) create mode 100644 Sources/Tests/UnitTests/Calculus/GaussianIntegralTest.cs diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 244c1e78d..6a1801003 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -1618,6 +1618,24 @@ antiderivative with no limit at the bound still gives `NaN`, as `-cos(x)` does a | `integral(1/x^2, x, 0, 1)` | `NaN` | `+oo` | | `integral(sin(x), x, 0, +oo)` | `NaN` | `NaN` | +### The Gaussian and the error functions are integrated + +`e^(-x^2)` was left unintegrated, and its antiderivative is `sqrt(pi)/2 erf(x)`. An exponential +of a quadratic is now integrated to an error function of the square it completes, written with +`erfi` and a real root where the square's coefficient is decidably positive. Beside an even power +of `x`, on either side of 0, it is integrated by parts down to that, and `erf`, `erfc` and `erfi` +themselves by parts against 1. A symbolic exponent is answered for the generic case, as +`F^(a x)/(a ln F)` already was ([#1501](https://github.com/asc-community/AngouriMath/issues/1501)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"e^(-x^2)".Integrate("x")` | `integral(e ^ (-x ^ 2), x)` | `sqrt(pi) / 2 * erf(x) + C` | +| `"e^(x^2)".Integrate("x")` | `integral(e ^ x ^ 2, x)` | `sqrt(pi) / 2 * erfi(x) + C` | +| `"x^2*e^(-x^2)".Integrate("x")` | `integral(x ^ 2 * e ^ (-x ^ 2), x)` | `x * e ^ (-x ^ 2) / (-2) + 1/2 * sqrt(pi) / 2 * erf(x) + C` | +| `"f^(a+b*x^2)".Integrate("x")` | `integral(f ^ (a + b * x ^ 2), x)` | `f ^ a * sqrt(pi) / (2 * sqrt(-b * ln(f))) * erf(sqrt(-b * ln(f)) * x) + C` | +| `"erf(x)".Integrate("x")` | `UnrecognizedFunctionParseException` | `x * erf(x) + e ^ (-x ^ 2) / sqrt(pi) + C` | +| `integral(e^(-x^2), x, -oo, +oo)` | `integral(e ^ (-x ^ 2), x, -oo, +oo)` | `sqrt(pi)` | + ### `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 ea6433622..74f47fa3a 100644 --- a/Sources/.editorconfig +++ b/Sources/.editorconfig @@ -624,3 +624,6 @@ file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed [Tests/UnitTests/Calculus/ImproperIntegralTest.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/GaussianIntegralTest.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/IntegralPatterns.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IntegralPatterns.cs index 82052f481..fcd3f7752 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IntegralPatterns.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IntegralPatterns.cs @@ -58,6 +58,9 @@ internal static Entity AntiderivativeLog(Entity arg) Entity.Cotanf(var arg) => arg, Entity.Absf(var arg) => arg, Entity.Signumf(var arg) => arg, + Entity.Erff(var arg) => arg, + Entity.Erfcf(var arg) => arg, + Entity.Erfif(var arg) => arg, Entity.Arcsinf(var arg) => arg, Entity.Arccosf(var arg) => arg, Entity.Arctanf(var arg) => arg, @@ -67,6 +70,111 @@ internal static Entity AntiderivativeLog(Entity arg) _ => null }; + /// + /// int F^(a x^2 + b x + c) dx for a non-zero a: the square completed is + /// a (x + b/(2a))^2 + c - b^2/(4a), so the integral is F^(c - b^2/(4a)) times + /// the Gaussian of A = a ln F in u = x + b/(2a). + /// + private static Entity Gaussian(Entity @base, Entity a, Entity b, Entity c, Entity.Variable x) + { + // Built without the terms a zero linear part would contribute: b/(2a) is 0 only + // provided a is not zero, and that condition, stated twice over, was all a symbolic + // F^(a x^2) got for it. The answer is for the generic case, as F^(a x)/(a ln F) is. + if (TreeAnalyzer.IsZero(b)) + return MathS.Pow(@base, c) * UnitGaussian(a * MathS.Ln(@base), x); + return MathS.Pow(@base, c - b * b / (4 * a)) * UnitGaussian(a * MathS.Ln(@base), x + b / (2 * a)); + } + + /// + /// 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. + /// + /// + /// The sign is read from Evaled, as the integrator's other sign tests are; they move + /// together to the interval evaluation of https://github.com/asc-community/AngouriMath/pull/1497. + /// Finiteness is asked as well, because holds + /// for NaN and +oo. + /// + private static Entity UnitGaussian(Entity A, Entity u) + => A.Evaled is Entity.Number.Real { EDecimal.IsFinite: true, IsPositive: true } + ? MathS.Sqrt(MathS.pi) / (2 * MathS.Sqrt(A)) * MathS.Erfi(MathS.Sqrt(A) * u) + : MathS.Sqrt(MathS.pi) / (2 * MathS.Sqrt(-A)) * MathS.Erf(MathS.Sqrt(-A) * u); + + /// + /// int k x^m F^(a x^2 + c) dx for an even whole m, not zero, and a non-zero + /// a: the Gaussian moments. With A = a ln F and I_m the integral of + /// x^m e^(A x^2), parts give I_m = x^(m - 1) e^(A x^2)/(2A) - (m - 1)/(2A) I_(m - 2), + /// and read the other way I_m = x^(m + 1) e^(A x^2)/(m + 1) - 2A/(m + 1) I_(m + 2) for a + /// negative m, each two powers nearer the Gaussian I_0. An odd m ends at + /// I_1 = e^(A x^2)/(2A), which is elementary and answered elsewhere, or at I_(-1), + /// which is the exponential integral, so it is not taken here. + /// https://github.com/asc-community/AngouriMath/issues/1501 + /// + private static Entity? GaussianMoment(Entity expr, Entity.Variable x) + { + // Asked of every integrand that reaches the table, so a type test comes first, and + // nothing here captures a local: a closure is allocated on entry, before any return. + if (expr is not (Entity.Mulf or Entity.Divf)) + return null; + var (numerator, denominator) = expr is Entity.Divf(var n, var d) ? (n, d) : (expr, (Entity)Entity.Number.Integer.One); + Entity coefficient = Entity.Number.Integer.One; + int? power = null; + Entity.Powf? gaussian = null; + foreach (var factor in Entity.Mulf.LinearChildren(numerator)) + if (!TakeMomentFactor(factor, inverted: false, x, ref coefficient, ref power, ref gaussian)) + return null; + foreach (var factor in Entity.Mulf.LinearChildren(denominator)) + if (!TakeMomentFactor(factor, inverted: true, x, ref coefficient, ref power, ref gaussian)) + return null; + if (power is not { } m || m == 0 || m % 2 != 0 || System.Math.Abs(m) > 32 || gaussian is null + || !TreeAnalyzer.TryGetPolyQuadratic(gaussian.Exponent, x, out var a, out var b, out var c) + || TreeAnalyzer.IsZero(a) || !TreeAnalyzer.IsZero(b)) + return null; + var bell = MathS.Pow(gaussian.Base, a * MathS.Sqr(x)); + return coefficient * MathS.Pow(gaussian.Base, c) * Moment(m, a * MathS.Ln(gaussian.Base), bell, x); + } + + /// + /// Reads one factor of a Gaussian moment into the constant coefficient, the power of + /// or the exponential, and says whether it was one of the three. + /// + private static bool TakeMomentFactor(Entity factor, bool inverted, Entity.Variable x, + ref Entity coefficient, ref int? power, ref Entity.Powf? gaussian) + { + if (!factor.ContainsNode(x)) + { + coefficient = inverted ? coefficient / factor : coefficient * factor; + return true; + } + if (power is null && factor == x) + { + power = inverted ? -1 : 1; + return true; + } + if (power is null && factor is Entity.Powf(var b, Entity.Number.Integer n) && b == x && n.EInteger.CanFitInInt32()) + { + power = inverted ? -n.EInteger.ToInt32Checked() : n.EInteger.ToInt32Checked(); + return true; + } + if (!inverted && gaussian is null && factor is Entity.Powf(var @base, _) exponential && !@base.ContainsNode(x)) + { + gaussian = exponential; + return true; + } + return false; + } + + /// + /// I_k, the integral of x^k e^(A x^2) for an even k, where + /// is that exponential as the integrand writes it. + /// + private static Entity Moment(int k, Entity A, Entity bell, Entity.Variable x) => k switch + { + 0 => UnitGaussian(A, x), + > 0 => MathS.Pow(x, k - 1) * bell / (2 * A) + ((1 - k) / (2 * A)).InnerSimplified * Moment(k - 2, A, bell, x), + _ => MathS.Pow(x, k + 1) * bell / (k + 1) + (-2 * A / (k + 1)).InnerSimplified * Moment(k + 2, A, bell, x), + }; + 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 @@ -126,6 +234,34 @@ _ when TryReadSineCosinePowers(expr, x, out var trigArg, out var sinePower, out !@base.ContainsNode(x) && TreeAnalyzer.TryGetPolyLinear(power, x, out var a, out _) => MathS.Pow(@base, power) / (a * MathS.Ln(@base)), + // The Gaussian: an exponential of a quadratic is an error function of the square it + // completes. With A u^2 the exponent's square part in u = x + b/(2a), the integral + // of e^(A u^2) is sqrt(pi)/(2 sqrt(-A)) erf(sqrt(-A) u), which differentiates back for + // every A that is not zero, whichever root is taken. Where A is decidably positive it + // is written with erfi and the real root instead, so that e^(x^2) is sqrt(pi)/2 erfi(x) + // rather than an error function of i x. https://github.com/asc-community/AngouriMath/issues/1501 + Entity.Powf(var @base, var power) when + !@base.ContainsNode(x) && TreeAnalyzer.TryGetPolyQuadratic(power, x, out var a, out var b, out var c) + && !TreeAnalyzer.IsZero(a) => + Gaussian(@base, a, b, c, x), + + // And the Gaussian beside an even power of x, to either side of it. + _ when GaussianMoment(expr, x) is { } moment => moment, + + // By parts against 1, the way the inverse trigonometric functions are: + // int erf(u) = u erf(u) + e^(-u^2)/sqrt(pi), and the same for the other two. + Entity.Erff(var arg) when + TreeAnalyzer.TryGetPolyLinear(arg, x, out var a, out _) => + (arg * MathS.Erf(arg) + MathS.Pow(MathS.e, -MathS.Sqr(arg)) / MathS.Sqrt(MathS.pi)) / a, + + Entity.Erfcf(var arg) when + TreeAnalyzer.TryGetPolyLinear(arg, x, out var a, out _) => + (arg * MathS.Erfc(arg) - MathS.Pow(MathS.e, -MathS.Sqr(arg)) / MathS.Sqrt(MathS.pi)) / a, + + Entity.Erfif(var arg) when + TreeAnalyzer.TryGetPolyLinear(arg, x, out var a, out _) => + (arg * MathS.Erfi(arg) - MathS.Pow(MathS.e, MathS.Sqr(arg)) / MathS.Sqrt(MathS.pi)) / a, + Entity.Absf(var arg) when TreeAnalyzer.TryGetPolyLinear(arg, x, out var a, out _) => // ∫ |ax + b| dx = sgn(ax + b) * (ax + b)^2 / (2a) MathS.Signum(arg) * MathS.Pow(arg, 2) / (2 * a), diff --git a/Sources/Tests/UnitTests/Calculus/GaussianIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/GaussianIntegralTest.cs new file mode 100644 index 000000000..4fb730a28 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/GaussianIntegralTest.cs @@ -0,0 +1,132 @@ +// +// 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 +{ + /// + /// The integrals that answer in the error functions: the Gaussian F^(a x^2 + b x + c) + /// by the square it completes, the Gaussian beside an even power of x by parts two + /// powers at a time, and the error functions themselves by parts against 1. + /// https://github.com/asc-community/AngouriMath/issues/1501 + /// + [Trait("Area", "Calculus")] + public sealed class GaussianIntegralTest + { + /// Off 0, where the negative powers are undefined. + private static readonly double[] Points = { -1.7, -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; + } + + [Theory] + [InlineData("e^(-x^2)")] + [InlineData("e^(-x^2 + 2*x + 1)")] + [InlineData("e^(-3*x^2 + x)")] + [InlineData("(1/2)^(x^2)")] + public void TheGaussianIsAnErrorFunction(string integrand) + => Assert.Contains(DifferentiatesBack(integrand).Nodes, node => node is Entity.Erff); + + /// + /// Where the square's coefficient is decidably positive the answer is written with + /// erfi and a real root, not as an error function of an imaginary argument. + /// + [Theory] + [InlineData("e^(x^2)")] + [InlineData("2^(x^2)")] + [InlineData("e^(2*x^2 - x)")] + public void AGrowingGaussianIsAnImaginaryErrorFunction(string integrand) + { + var integral = DifferentiatesBack(integrand); + Assert.Contains(integral.Nodes, node => node is Entity.Erfif); + Assert.DoesNotContain(integral.Nodes, node => node == MathS.i); + } + + /// + /// A symbolic exponent is answered for the generic case, as F^(a x)/(a ln F) is, + /// and the answer differentiates back whichever sign the square's coefficient turns out + /// to have: the root of -b ln f is imaginary for a positive b ln f, and the + /// answer is still right. + /// + [Theory] + [InlineData("f^(a + b*x^2)", "-3/2")] + [InlineData("f^(a + b*x^2)", "3/2")] + [InlineData("x^2*f^(a + b*x^2)", "-3/2")] + [InlineData("x^8*f^(a + b*x^2)", "-3/2")] + [InlineData("f^(a + b*x^2)/x^4", "3/2")] + public void ASymbolicGaussian(string integrand, string b) + => DifferentiatesBack(integrand, ("f", "2"), ("a", "1/3"), ("b", b)); + + [Fact] + public void ASymbolicGaussianWithALinearTerm() + => DifferentiatesBack("f^(a + b*x + c*x^2)", ("f", "3"), ("a", "1/2"), ("b", "2/3"), ("c", "-5/4")); + + /// A symbolic exponent carries no condition that its linear part is zero. + [Fact] + public void ASymbolicGaussianIsUnconditional() + => Assert.DoesNotContain(DifferentiatesBack("f^(a + b*x^2)", ("f", "2"), ("a", "1/3"), ("b", "-1")).Nodes, + node => node is Entity.Providedf); + + [Theory] + [InlineData("x^2*e^(-x^2)")] + [InlineData("x^4*e^(x^2)")] + [InlineData("x^8*e^(-x^2)")] + [InlineData("e^(-x^2)/x^2")] + [InlineData("e^(x^2)/x^4")] + [InlineData("3*x^6*e^(-2*x^2)/5")] + public void TheGaussianBesideAnEvenPower(string integrand) => DifferentiatesBack(integrand); + + /// + /// An odd negative power ends at int e^(A x^2)/x, which is the exponential + /// integral, so it is not taken. + /// + [Fact] + public void TheGaussianOverAnOddPowerIsNotAnErrorFunction() + => Assert.Contains("integral(", "e^(-x^2)/x".ToEntity().Integrate("x").Stringize()); + + [Theory] + [InlineData("erf(x)")] + [InlineData("erf(2*x + 1)")] + [InlineData("erfc(3*x)")] + [InlineData("erfi(x/2)")] + [InlineData("5*erfc(1 - x)")] + public void TheErrorFunctionsByParts(string integrand) => DifferentiatesBack(integrand); + + /// Over the whole line, through the error functions' values at +-oo. + [Theory] + [InlineData("e^(-x^2)", "sqrt(pi)")] + [InlineData("e^(-x^2/2)", "sqrt(2 * pi)")] + [InlineData("e^(-(x - 1)^2)", "sqrt(pi)")] + public void TheGaussianOverTheWholeLine(string integrand, string value) + => Assert.Equal(value.ToEntity().Simplify(), MathS.Integral(integrand.ToEntity(), "x", "-oo".ToEntity(), "+oo".ToEntity()).Simplify()); + } +} From 1427915f0dc1b32aa97e2b4fdb8226f052a4afa9 Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Sun, 27 Sep 2026 17:16:10 +0000 Subject: [PATCH 2/3] e^(x^2) is answered, so the examples of an unanswered integral are x^x and e^(x^3) The documented example of Integrate(string), the transformation tests and the exponential ansatz's declined row used e^(x^2) as the integral with no answer; it is sqrt(pi)/2 erfi(x) now. They use x^x, which has no antiderivative in any function the library will have, and e^(x^3) where the row is about the ansatz's shape. y' + y = e^(x^2) moves from what the ODE solver declines to what it solves. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura --- .../AngouriMath/Convenience/AngouriMathExtensions.cs | 10 ++++++---- .../Tests/UnitTests/Calculus/ExponentialAnsatzTest.cs | 8 +++++--- .../Calculus/ExponentialSubstitutionIntegralTest.cs | 3 +++ .../Calculus/OrdinaryDifferentialEquationTest.cs | 8 +++++--- .../UnitTests/Convenience/ExtensionsExampleTest.cs | 3 ++- .../Core/Transformations/TransformationTest.cs | 10 +++++----- 6 files changed, 26 insertions(+), 16 deletions(-) diff --git a/Sources/AngouriMath/Convenience/AngouriMathExtensions.cs b/Sources/AngouriMath/Convenience/AngouriMathExtensions.cs index 9c2dc1a9d..fdfe1d497 100644 --- a/Sources/AngouriMath/Convenience/AngouriMathExtensions.cs +++ b/Sources/AngouriMath/Convenience/AngouriMathExtensions.cs @@ -784,17 +784,19 @@ public static Entity Differentiate(this string str, Variable x) /// An integrated expression. It might remain the same or be transformed into nodes with no integrals. /// /// - /// The constant of integration is carried as the variable C; an integral that - /// has no elementary antiderivative comes back as the node itself, which is how this - /// says it could not settle the question: + /// The constant of integration is carried as the variable C; an integral with no + /// antiderivative written in the library's functions comes back as the node itself, + /// which is how this says it could not settle the question: /// /// Console.WriteLine("1 / x".Integrate("x")); /// Console.WriteLine("e ^ (x ^ 2)".Integrate("x")); + /// Console.WriteLine("x ^ x".Integrate("x")); /// /// Prints /// /// ln(x) + C - /// integral(e ^ x ^ 2, x) + /// sqrt(pi) / 2 * erfi(x) + C + /// integral(x ^ x, x) /// /// public static Entity Integrate(this string str, Variable x) diff --git a/Sources/Tests/UnitTests/Calculus/ExponentialAnsatzTest.cs b/Sources/Tests/UnitTests/Calculus/ExponentialAnsatzTest.cs index afa98cdbc..d0bab5823 100644 --- a/Sources/Tests/UnitTests/Calculus/ExponentialAnsatzTest.cs +++ b/Sources/Tests/UnitTests/Calculus/ExponentialAnsatzTest.cs @@ -28,8 +28,10 @@ namespace AngouriMath.Tests.Calculus /// /// /// Every answer is differentiated back and compared to the integrand numerically; the - /// non-elementary neighbours — e^x/x, e^(x^2) — are pinned as declined, since - /// answering them would be a wrong answer and not a missing one. + /// non-elementary neighbours — e^x/x, e^(x^3) — are pinned as declined, since + /// an elementary answer to them would be a wrong answer and not a missing one. + /// (e^(x^2) is answered, with erfi, since + /// https://github.com/asc-community/AngouriMath/issues/1501.) /// /// [Trait("Area", "Calculus")] @@ -181,7 +183,7 @@ public void AnExponentialTimesARootOfAQuadraticIsDeclinedWhereThereIsNone() /// [Theory] [InlineData("e^x/x")] - [InlineData("e^(x^2)")] + [InlineData("e^(x^3)")] [InlineData("e^x/(1 + x)")] [InlineData("e^(1/x)")] public void ANonElementaryOneIsDeclined(string integrand) diff --git a/Sources/Tests/UnitTests/Calculus/ExponentialSubstitutionIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/ExponentialSubstitutionIntegralTest.cs index 398685a62..65d69950d 100644 --- a/Sources/Tests/UnitTests/Calculus/ExponentialSubstitutionIntegralTest.cs +++ b/Sources/Tests/UnitTests/Calculus/ExponentialSubstitutionIntegralTest.cs @@ -333,6 +333,9 @@ public void ASymbolicBase() /// elementary antiderivative at all, and its exponent is not linear, which is the check /// that stops it — a rewrite would leave an x standing beside the new variable. /// Should it later be answered by something else, this moves rather than being deleted. + /// It now is, with erfi, by the Gaussian's own rule + /// (https://github.com/asc-community/AngouriMath/issues/1501), and what this row pins, that + /// the substitution did not write it, still holds. /// [Theory] [InlineData("e^(x^2)")] diff --git a/Sources/Tests/UnitTests/Calculus/OrdinaryDifferentialEquationTest.cs b/Sources/Tests/UnitTests/Calculus/OrdinaryDifferentialEquationTest.cs index 38b0fc8d6..610842029 100644 --- a/Sources/Tests/UnitTests/Calculus/OrdinaryDifferentialEquationTest.cs +++ b/Sources/Tests/UnitTests/Calculus/OrdinaryDifferentialEquationTest.cs @@ -88,6 +88,7 @@ private static void Satisfies(string equation) [InlineData("derivative(apply(y, x), x) + 2*x*apply(y, x) - x")] [InlineData("derivative(apply(y, x), x) - apply(y, x) - e^x")] [InlineData("derivative(apply(y, x), x) + apply(y, x) - sin(x)")] + [InlineData("derivative(apply(y, x), x) + apply(y, x) - e^(x^2)")] public void WithARightHandSide(string equation) => Satisfies(equation); /// @@ -125,8 +126,9 @@ public void TheSameEquationWrittenTwoWaysHasOneSolution() /// /// What it declines. The first four are not linear in the unknown and its derivative, so /// the method does not apply at all; the last is linear and declined because - /// int e^(x^2) dx has no closed form, which is a fact about the integral rather - /// than about this solver. + /// int e^(x^3 + x) dx has no closed form, which is a fact about the integral rather + /// than about this solver. With e^(x^2) in its place it is answered, through + /// erfi. /// /// /// Declining is the right answer to each. The alternative — reading a nonlinear equation @@ -139,7 +141,7 @@ public void TheSameEquationWrittenTwoWaysHasOneSolution() [InlineData("derivative(apply(y, x), x) + apply(y, x)^2")] [InlineData("sin(derivative(apply(y, x), x)) - 1")] [InlineData("derivative(apply(y, x), x)^2 - apply(y, x)")] - [InlineData("derivative(apply(y, x), x) + apply(y, x) - e^(x^2)")] + [InlineData("derivative(apply(y, x), x) + apply(y, x) - e^(x^3)")] public void WhatItDeclines(string equation) => Assert.Null(MathS.SolveOde(equation.ToEntity(), "y", "x")); diff --git a/Sources/Tests/UnitTests/Convenience/ExtensionsExampleTest.cs b/Sources/Tests/UnitTests/Convenience/ExtensionsExampleTest.cs index b1b1576de..4e3800e61 100644 --- a/Sources/Tests/UnitTests/Convenience/ExtensionsExampleTest.cs +++ b/Sources/Tests/UnitTests/Convenience/ExtensionsExampleTest.cs @@ -259,7 +259,8 @@ [Fact] public void DifferentiateString() [Fact] public void IntegrateIndefiniteString() { Printed("ln(x) + C", "1 / x".Integrate("x")); - Printed("integral(e ^ x ^ 2, x)", "e ^ (x ^ 2)".Integrate("x")); + Printed("sqrt(pi) / 2 * erfi(x) + C", "e ^ (x ^ 2)".Integrate("x")); + Printed("integral(x ^ x, x)", "x ^ x".Integrate("x")); } [Fact] public void IntegrateDefiniteString() diff --git a/Sources/Tests/UnitTests/Core/Transformations/TransformationTest.cs b/Sources/Tests/UnitTests/Core/Transformations/TransformationTest.cs index 93c8497aa..e3a2e08f2 100644 --- a/Sources/Tests/UnitTests/Core/Transformations/TransformationTest.cs +++ b/Sources/Tests/UnitTests/Core/Transformations/TransformationTest.cs @@ -85,9 +85,9 @@ public void AFixedPointSucceedsWithoutChangingAnything() [Fact] public void ApplyOrKeepHandsBackTheInputWhereThereIsNoAnswer() { - var unanswerable = Transformation.Integration("x").Apply(Parse("e ^ (x ^ 2)")); + var unanswerable = Transformation.Integration("x").Apply(Parse("x ^ x")); Assert.False(unanswerable.Succeeded); - Assert.Equal(Parse("e ^ (x ^ 2)"), unanswerable.OutputOrInput); + Assert.Equal(Parse("x ^ x"), unanswerable.OutputOrInput); } [Fact] @@ -140,7 +140,7 @@ public void ANoAnswerAnywhereInAChainIsANoAnswerForTheChain() // The integral has no closed form, so nothing downstream of it can have run. var chain = Transformation.Integration("x").Then(Transformation.Simplification); - Assert.False(chain.Apply(Parse("e ^ (x ^ 2)")).Succeeded); + Assert.False(chain.Apply(Parse("x ^ x")).Succeeded); } [Fact] @@ -342,14 +342,14 @@ public void SubstituteIsItsTransformation(string raw) [Fact] public void AnIntegralWithNoClosedFormHasNoAnswerRatherThanAWrongOne() { - var result = Transformation.Integration("x").Apply(Parse("e ^ (x ^ 2)")); + var result = Transformation.Integration("x").Apply(Parse("x ^ x")); Assert.False(result.Succeeded); Assert.Null(result.Output); // The 1.x method makes the same claim in the shape its callers expect: an // unevaluated node, which is "I could not settle this" and not NaN. - var legacy = Parse("e ^ (x ^ 2)").Integrate("x"); + var legacy = Parse("x ^ x").Integrate("x"); Assert.IsType(legacy); Assert.False(legacy.IsNaN); } From 30662d3873049e5c73d8cbe660a9c2c28d83d849 Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Sun, 27 Sep 2026 17:17:31 +0000 Subject: [PATCH 3/3] The Gaussian and its moments over the half line With #1509, a definite integral takes the limit where its antiderivative is undefined at a bound, so the moments' x^k e^(-x^2) at +oo is 0 and the half line gives sqrt(pi)/2, sqrt(pi)/4 and 3 sqrt(pi)/8. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_012sonx8iAspMiwRwokT1Ura --- .../Tests/UnitTests/Calculus/GaussianIntegralTest.cs | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/Sources/Tests/UnitTests/Calculus/GaussianIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/GaussianIntegralTest.cs index 4fb730a28..1a4f2847a 100644 --- a/Sources/Tests/UnitTests/Calculus/GaussianIntegralTest.cs +++ b/Sources/Tests/UnitTests/Calculus/GaussianIntegralTest.cs @@ -128,5 +128,17 @@ public void TheGaussianOverAnOddPowerIsNotAnErrorFunction() [InlineData("e^(-(x - 1)^2)", "sqrt(pi)")] public void TheGaussianOverTheWholeLine(string integrand, string value) => Assert.Equal(value.ToEntity().Simplify(), MathS.Integral(integrand.ToEntity(), "x", "-oo".ToEntity(), "+oo".ToEntity()).Simplify()); + + /// + /// And over the half line, where the moments' elementary terms are x^k e^(-x^2) at + /// +oo, which is a limit rather than a value to substitute + /// (https://github.com/asc-community/AngouriMath/issues/1507). + /// + [Theory] + [InlineData("e^(-x^2)", "sqrt(pi) / 2")] + [InlineData("x^2*e^(-x^2)", "sqrt(pi) / 4")] + [InlineData("x^4*e^(-x^2)", "3 * sqrt(pi) / 8")] + public void TheGaussianOverTheHalfLine(string integrand, string value) + => Assert.Equal(value.ToEntity().Simplify(), MathS.Integral(integrand.ToEntity(), "x", "0".ToEntity(), "+oo".ToEntity()).Simplify()); } }